# NumericalSimulat_省略_ntaminantsinSoil_.docx

Pedosphere 112 131 136, 2001 ISSN 1002-0160/CN 32-1315/P 2001 SCIENCE PRESS, BEIJING 131 Numerical Simulation of Preferential Flow of Contaminants in Soil*1 XU SHAOHUI, DU ENHAO and ZHANG JIABAO Institute of Soil Science the Chinese Academy of Sciencesf Nanjing 210008 China Received December 2000; revised February 20, 2001 ABSTRACT A simple modeling approach was suggested to simulate preferential transport of water and contaminants in soil. After saturated hydraulic conductivity was interpolated by means of Krige interpolation or scaling , and then zoned, the locations where saturated hydraulic conductivity Weis larger represented regions where preferential flow occurred, because heterogeneity of soil, one of the mechanisms resulting in preferential flow, could be reflected through the difference in saturated hydraulic conductivity. The modeling approach was validated through numerical simulation of contaminant transport in a two-dimensional hypothetical soil profile. The results of the numerical simulation showed that the approach suggested in this study was feasible. Key Words contaminant, numerical simulation, preferential flow, soil INTRODUCTION Preferential flow refers to the rapid transport of water and solutes through some small portions of a soil volume that receives over its entire inlet boundary. It is the general phenomenon rather than the exception in soils. As a result of preferential flow, the utilization of water and nutrients by plant is reduced and the groundwater recharge is speeded up. However, because of the short time and small areas they contact with soil matrix, some of contaminants will remain undegraded and move rapidly down, increasing the possibility of contaminating groundwater. Therefore, in recent twenty years, preferential flow has become one of the central issues of research in the related fields sucli as water science, environmental science， agronomy， soil physics, etc” in the world. At the same time, owing to the great variability in time and space, preferential flow is also a difficult problem of research. Study on preferential flow in China has just begun. A number of experiments in the laboratory and field have been conducted to extensively study influencing factors, flow mechanisms, observation s， etc” of preferential flow. But，study on how to describe quantitatively preferential flow is in its infancy. Many models have ever been proposed to simulate preferential flow, for example, mobile-immobile model Project supported by the National Natural Science Foundation of China No. 49971041, the National Key Basic Research and Development Program of China G1999011803 and the Director Foundation of Institute of Soil Science, the Chinese Academy of Sciences ISSDF0004. 132 S. H. XU et al van Genuchten and Wierenga, 1976, two-region model Skopp et al. y 1981, double-pore model Gerke and van Genuchten, 1993, kinematic wave model Germann and Beven, 1985, numerical model of piecewise linear approximation of hydraulic conductivity Steenhuis et 1990, two-phase model Hosang, 1993, etc. There were some defects or deficiencies in all of these models, and they could not be applied in practice. Therefore, the objective of this study was to find a more simple approach to simulate preferenticd flow and validate the approach through numerical simulation. MODELING APPROACH SUGGESTED It is known that one of the mechanisms resulting in preferential flow is heterogeneity of soil, whereas the heterogeneity can be characterized by saturated hydraulic conductivity. Provided that saturated hydraulic conductivity values of a number of points in the study area were measured, the values of many other points would also be obtained by means of Krige interpolation or scaling . Furthermore, the soil heterogeneity might be zoned; namely, the points where the values of saturated hydraulic conductivity were close were divided into one region and an identical vahie of hydraulic conductivity was prescribed for the region. The zone with larger parameter value, where water and contaminants transported rapidly, was that of preferential flow. Many researchers e.g.y Roth et a/., 1991 suggested that the Richards equation and convection-dispersion equation based on the Darcy law and Fick law, respectively, are not suitable for describing the preferential flow of water and solute in soil. This was because an average or identical parameter value was utilized in the whole study area. If different values of saturated hydraulic conductivity were given in discrete points of space according to the practical situation, or the values were given according to zoning of saturated hydraulic conductivity, modeling preferential transport of water and solute in soil would be completed by solving convection-dispersion equation after calculating the pressure head and water content through solving Richards equation and obtaining the seepage velocity by using the of Yeh 1981. NUMERICAL SIMULATION A 90 cm x 200 cm two-dimensional hypothetical vertical soil profile X-Z plane was discretized into 400 triangular elements and 231 nodes, and the discrete map is shown in Fig.L The governing equation and initial and boundary conditions to describe water flow in the soil profile may be written as la SIMULATION OF PREFERENTIAL FLOW OF CONTAMINANTS 133 where hxz,t is the pressure head L; Kh the unsaturated hydraulic conductivity L T 一工 ；C*/i the specific capacity L-1; and z the coordinates in X and Z directions Z direction downward positive, respectively L and t the time T. Xm 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0,8 Fig. 1 Discrete triangular elements of the soil profile in the modeling region and zonation of soil saturated hydraulic conductivity. I and II represent zones of smaller and larger soil hydraulic conductivity, respectively. Fig. 2 Contaminant concentration distribution in soil profile at t 1 h a and t 12 h b. The governing equation to describe contaminant transport in the soil profile and the corresponding initial and boundary conditions may be written as where cx， 之， t is the concentration of contaminant M L 一 3; 0 the volumetric water content L3 L-3; Vx and Vz the velocities in X direction and Z direction, respectively L T_1, V yJV 乃 xx and Zzz the dispersion coefficients in X direction and Z direction, 134 S. H. XU et al. respectively L2 T--1; and x and 2 the coordinates in X direction and Z direction, respectively L. The modeling region was divided into two subregions based on the soil saturated hydraulic conductivity one with the saturated hydraulic conductivity value K8I of 0.26 cm h*-1 Zone I and the other with the saturated hydraulic conductivity value 2 29.7 cm h-1 Zone II . The zoning map of the soil saturated hydraulic conductivity is presented in Fig. 1. The functions to characterize soil hydraulic property axe the van Genuchten-Mualem model t1 3 Kh - 1 - 5e1/mm2 where Se is the degree of saturation, dimensionless; 0 the volumetric water content L3 L-3; 6 the saturated water content L3 L-3, 0r the residual water content L3 L-3; a the parameter to represent the pore size distribution L-1; n and m the parameters to describe the shape of water retention curve, m 1 1/n， dimensionless; Z the parameter to characterize soil pore counection generally taken as 05 dimensionless. a,n and m are all larger than zero. In this study, the values of parameters in the van Genuchten model are as follows Leij et a/., 1996 9ri 0.095 cm3 cm3, dr2 0.045 cm3 cm-3; 06 0.41 cm3 cm-3, 0e2 -43 cm3 cm-3; ai 0.019 cm-1, a2 0.145 cm-1 and ni 1.31, n2 2.68 for Zone I and Zone II, respectively. The values of contaminant transport parameters in the soil were longitudinal dispersivity 0.5 cm, transverse dispersivity ax l cm and time step At 0.01 h. The solution was obtained in four steps 1 Step I. In order to keep mass conservation, Equation la was rewritten as the mixed of water content 9x,z,t and pressure head hx,z,t Celia et a/., 1990 By using the finite element , pressure head hx,zt and water content 6xz,t may be get by solvingEquation 5 and giving the corresponding initial and boundary conditions. 2 Step II. If the seepage velocity Vx and Vz were calculated based on the Darcy law， Vx and Vz 火 / - 1， and the pressure head gradient 石 and were approximated by the difference , then, velocity Vx and Vz would be dis- continuous at the finite element boundaries, and the bigger error would be caused when the convection-dispersion equation was solved. To overcome this defect, the of Yeh 1981 was used to calculate velocity Vx and Vz, that is, by means of solving an equation group GVx Qx and GV Z Qz. In the equation group, G is N order sym- M rr metric positively definite matrix and its element is gij ifiiipjdxdz, and Qx and 1 e SIMULATION OF PREFERENTIAL FLOW OF CONTAMINANTS 135 Qz 3X6 the N dimensional vectors, their elements are qxi M -E 1 e Kh dh 6S dx fiidxdz and 如 _E|甲 造 -1 respectively, where N is the total number of nodes, M is the total number of element, and fii and are the base functions of triangular element and weight function, respectively. 3 Step III. According to the relationship between disper- rn 1 J1, , ahV aTV aTV aLV - sion coefficient D and velocity Vy namely, Dzz ------- L - ---- and Dx --------- -------- , where Dzz is the longitudinal dispersion coefficient, Dxx is the transverse dispersion coeffi- cient and the molecular diffusion coefficients are ignored, the concentration of contaminant cx, z, t could be calculated by solving Equation 2a with the finite element Simunek et aln 1994. The distribution of contaminant concentrations in the soil profile at time t 1 h and t 12 h shown in Fig. 2 indicated that contaminant transport in Subregion II where the saturated hydraulic conductivity was larger, or preferential flow occurred, was much more rapid than that in Subregion I where the saturated hydraulic conductivity was smaller, or no preferential flow occurred, which reflects the nature of preferential flow. CONCLUSIONS Transport of contaminant in the subregion where the saturated hydraulic conductivity was larger, or preferential flow occurred, was much more rapid than that in the subregion where the saturated hydraulic conductivity was smaller, or no preferential flow occurred, which reflects the nature of preferential flow, and indicates that the issue of modeling preferential flow may be resolved more efficiently by using the suggested in this study. Since it is greatly difficult to study preferential flow in soil, a number of comprehensive experiments in laboratory and field, and novel observation s will be needed to obtain a vast number of data, and to accurately describe the occurrence and development of preferential flow as well as to enhance the modeling techniques, so that the theory of preferential flow could be perfected to solve some practical problems such as utilization of water and nutrient by plants, groundwater recharge, contamination of groundwater, etc. REFERENCES Celia, M. A., Bououtas, E. T. and R. L. Zarba, 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res, 20 14831496, Gerke, H. H. and van Genuchten, M. T. 1993. A dual-porosity model for simulating preferential movement of watr and solutes in structured porous media. Water Resour. Res. 29 3053L9. Germann, P. F. and Beven, K. 1985. Kinematic wave approximation to infiltration into soils with sorbing macropores. Water Resour. Res. 21 990996. Hosang, J. 1993. Modelling preferential flow of water in soils a two-phase approach for field conditions. Geoderma. 58 149163. Leij, F. J., Alves, W. J., van Genuchten, M. T. and Williams, J. R. 1996. The UNSODA Unsaturated Soil Hydraulic Database. Users Manual Version 1.0. EPA, USA. 113pp. 136 S. H. XU et al Roth, K., Jury, W. A., Flhler, H. and Attinger, W. 1991. Transport of chloride through an unsaturated field soil.iyater/wotir1Res.272 5332 541. Simiinek, J., Vogel, T. and van Genuchten, M. T. 1994. The SWMS.2D Code for Simulating Water Flow and Solute Transport in Two-Dimensional Variably Saturated Media. Version 1.1. Research Report No. 132. U. S. Salinity Laboratory, USDA, ARS, Riverside, CA. 216pp. Skopp, J., Gardner, W. R. and Tyler, E. J. 1981. Solute movement in structured soils Two-region model with small intersiction. Soil Sci. Soc. Am. J, 45 837842. Steenhuis, T. S., Parlange, J. Y. and Andreini, M. S. 1990. A numerical model for preferential solute movement in structured soils. 7eoderma. 46 193 20 van Genuchten, M. T. and Wierenga, P. J. 1976. MJISS transfer in sorbing transfer porous media I. Analytical solutions. Soil Sci. Soc. Am. J, 40 473480. Yeh, G. T. 1981. On the computation of dsurcia velocity and mass balance the finite element modeling of groundwater flow. Heater 17 1529 1534.