STOCHASTIC ANALYSIS OF WATER FLOW IN HETEROGENEOUS MEDIA

313Journal of Hydrodynamics ,Ser. B,2005,17(3) :313- 322China Ocean Press，Beijing 一Printed in ChinaSTOCHASTIC ANALYSIS OF WATER FLOW IN HETEROGENEOUSMEDIAYANG Jin-zhong，WANG Wei-ping, CAI Shu-ying，LI Shao-longNational Key Laboratory of Water Resources and Hydropower Engineering Science, W uhan University ,W uhan 430072，China, e-mail:jzyang@ whu. edu. cn( Received Jan. 5, 2004)ABSTRACT: A stochastic model for saturated- unsaturated flow and solute transport processes becomes neces-flow is developed based on the combination of the Karhunen-sary (Bresler and Dagan 1981; Dagan and BreslerLoeve expansion of the input random soil properties with a1979). Field observations show that the hydraulicperturbation method. The saturated hydraulic conductivityk,properties of soils vary significantly with special(x) is assumed to be log normal random functions， ex-pressed by f(x). f(x ) is decomposed as infinite series in alocation even for a same soil type. It has been rec-set of orthogonal normal random variables by the Karhunen- ognized that the theory of stochastic processes pro-Loeve (KL) expansion and the pressure head is expand as vides a useful method for evaluating flow andpolynomial chaos with the same set of orthogonal randomtransport uncertainties. Many stochastic theoriesvariables. With these expansions, the stochastic saturated-have been developed to study the effects of spatialunsaturated flow equation and the corresponding initial andvariability on flow and transport in both saturatedboundary conditions are represented by a series of determin-istic partial differential equations which can be solved subse-(Gelhar 1993; Dagan 1989) and unsaturated zonesquently by a suitable numerical method. Some examples are (e. g.，Yeh et al. 1985a， b; Mantoglo and Gelhargiven to show the rliability and eficieneyg of the proposed 1987 ; Mantoglo 1992; Yang et al.1996 a,b; Yangmethod.et al. 1997; Zhang and Lu 2002). In the unsatu-KEY WORDS: saturated- unsaturated flow，Karhunen-Lo-rated zone the problem is further complicated byeve (KL) expansion， perturbation method，stochastic nu-the fact that the flow equations are nonlinear bemerical modelingcause unsaturated hydraulic conductivity dependson pressure head. Liedl (1994) proposed a pertur-1.INTRODUCTIONbation model for transient unsaturated flow. Li andThe unsaturated zone is a most active zone and Yeh (1998) studied transient unsaturated flow ina buffer which connect all the processes aboveheterogeneous porous media using a vector state-ground and the groundwater. Movement of waterspace approach and investigated the behavior ofand pollutants in the vadose zone affects thehead variances for transient unsaturated flow ingrowth of vegetation and wildlife, the amount oftwo dimensions. Zhang ( 1999) studied transientrecharge and evapotranspiration, and overall waterunsaturated flow for nonstationary situation andquality. Spatial variability of hydraulic properties .derived partial differential equations governing thein the vadose zone is one important factor that con-statistical moments by perturbation expansions andtrols the migration rate and path of water and pol-lutants (Ye and Yang 1995; Y'e and Yang 1996).then implement these equations by the method ofBecause of our incomplete knowledge about thefinite differences. Lu and Zhang ( 2002 ) studiedspatial distribution of hydraulic properties, predic-the中國煤化工base of the van Genu-tion of flow and transport processes in the vadose chteYHC N M H Gclation and found thezone always involves some degree of uncertainty. impacts of the different constitutive models. Be-To address the uncertainty, stochastic modeling of cause the complicity of the soil distribution and the* Project supported by the National Natural Science Foundation of China (Grant Nos: 50279039，50379042).Biography :人瓦黎括zhong (1953-)， Male, Ph. D.， Professor314boundary conditions in the real world, the analyti- ductivity, Q,(x,t)， hr, (x,t) and hami(x) are thecal model can only be used for the analysis of the flux out of the Neumann boundary F2，the pres-simpler problem. The state-space approach (Li and sure head at the first type boundary F，and the in-Yeh 1998) and the perturbation moment equation itial pressure head distribution in the domain, reapproach (Mantoglo 1992; Zhang 1999; Zhang and spectively. For simplicity, g(x,t),Q.(x,t),hn (x,Lu 2002) can be used in more general complicated 1) and hit (x) are treated as deterministic func-situaticns. However, the challenge is the compu- tions. For the simplicity of the analysis and thetational effort even for the first order moment ise- comparison of the results with the literature revaluated.sults, K (x,t) and θ(x,t) are expressed by theIn the spectral stochastic finite element meth-Gardner-Russo model asod, the spatial random media property in elementsis expressed by the KL expansion and the random- K(x,t) = K,(x)exp[a(x)h(x,t)], h<0(3a)ness of the system is treated as an additional di-mension in which a set of basic functions is de- K(x,t)= K,(x), h≥0(3b)fined，which is referred to as Polynomial Chaos .(Ghanem and Spanos 1991; Wiener 1938). The re-0.(x,l) = (0. - 0,){exp[號a(x)h(x,l)]●sulting deterministic system of algebraic equationcan be used for the solution of the deterministic co-efficient of the polynomial chaos expansion.[1- -a(x)h(x,t)]}/(m+2), h< 0(3c)In this work, the stochastic saturated unsatu-rated flow equation with the spatial random hy-draulic conductivity and moisture capacity is solved 0.(x,t) = (0.- 0,),h≥0(3d)based on the combination of KL expansion of therandom input of the media properties and the per- where0, 0,，and 0, are the effective, residual, andturbation method. The steady and unsateady satu- saturated water content, respectively, K,(x) is therated-unsaturated flow problem are solved with saturated hydraulic conductivity, a(x) is the soilmore than one spatial random media function in- parameter related to the pore size distribution, andputs.m is a experience parameter related to tortuosity ，STOCHASTIC SATURATED-UNSATURATED which can be determined from the field experi-FLOW EQUATIONment, and the specific soil water capacity C(x,t)=The equation of transient flow in saturated- d0/ dh.unsaturated media can be expressed asThe variabilities of 0, and 0, are likely smallcompared to that of the effective water content0. .C(x,t)ah(x.t). = v[ K(x,t)Vh(x,t)+z]+In this study, 0, and 0, are assumed to be determin-)tistic constants. The soil pore size distribution pa-rameter a (x) and the log-transformed saturatedg(x,t)(1) hydraulic conductivity f(x)= lnK,(x) are treatedas spatial random functions. They are assumed towith the following initial and boundary conditionsbe normal random functions with known covari-ances.一{K(x)V[h(x,l)+z]. n}|r, = Q,(x,l) (2a) :KARHUNEN-LOEVE EXPANSION OF RAN-DOM FUNCTIONSh(x,l) Ir = hr, (x,l)(2b)The Karhunen-Loeve expansion of a spatialand中國煤化Ibased on the spectral .h(x,t) |-o = han:(x)(2c) expalMHCN MHGfunction Cw (x，y)，where, x anay are usea to denote spatial coordi-where h(x,t)十z is the total hydraulic head, h(x，nates at two different spatial locations, w indicatest) is the pressure head, g(x,t) is the source or the random nature of the corresponding quantity.sink term，C(x,t) is the specific moisture capaci- It can be shown that the covariance function hasty, and K(萬(wàn)雙病s the unsaturated hydraulic con- the following form (Ghanem and Spanos 1991):315φ(x,w) = q(x)+2√λgφ;° (x)E,(w) =h(x,l) = H'°(x,l)+ 2H,(x.l)F(5,)+q(x)+ 2q(x)&(w)(4)ZH,(x.0)r(5.5,)+. =1 j=1where φ(x) is the mean of the spatial random func-H (x,l)F(5,5,.5) + .....tion p(x,c), and {&(w)} form a set of orthogonal2221random variables. Furthermore, {p;° (x)} are the(6)eigenfunctions and {入} are the corresponding ei- whereI,(ζ, ,5 ,ζ ..5) is the nth order poly -genvalues of the covariance kernel, which can be enomial chaos( Ghanem and Spanos 1991).valuated as the solution to the following integral e-All elements in {I,(ζr，5 ,5。,.ζ )〉arequation:mutually orthogonal and form a basis in second or-der random function space. Any element of the|C(x,y)p;" (y)dy= A;pi" (x)(5) {I,(5，51e，5 ,... )} isa polynomial of(5q ,，therefore，we can rearrange the ex-where 2 denotes the spatial domain over which thepansion of h(x,t) as the summation of all possiblerandom function φ(x,w) is defined. The most im-polynomials of(ζn，5.，5portant aspect of this spectral representation is thath(x,t)= h°(x,t) + Zh;(x.t)5,+the spatial random fluctuations have been decom-posed into a set of deterministic functions multipl-ying random variables. If the random function φZ(yǔ)hr(x.055,+習ha(x.0555 +..+(x,w)is Gaussian, then the random variables,jik= 1{ξ (w)} form an orthonormal Gaussian vector.We carry out the analysis by assuming the2 (.. (x.t)[Iζ,]}+= .perfectly correlated and uncorrelated cases betweenf(x) and a(x). For the perfectly correlated case,a(x) can be expressed in terms of f (x). Thereh°(x,l)十h"(x,t)+h(2(x,l)+... (7)fore，f(x) and a(x) can be decomposed by Kar-hunen- Loeve expansion on the same set of normalA complete probability characterization of the ran-random base {ξ;(w)}. If f(x) and a(x) are uncor- dom function h(x,l) is obtained if all the determin-related, they can be decomposed by Karhunen-Lo- istic coefficients (h0) (x,t),h;(x,t),h; (x,t),hje .eve expansion based on two sets of uncorrelated (x,t) .... i, j, k... = 1, 2, .. in expression (7)random bases {6(w)} and {η:(w)} ，respectively. have been calculated.We arrange {ξ(w)} and {η:(w)} as {ζ;(w)} =PERTURBATION EQUATIONS{G(w),η(w)} and {ζ;(w)} will be used as a baseThe log- transformed unsaturated soil conduc-set series for the expansion of pressure head.tivity Y(x,t) ln K(x,t) can be written as:4.POLYNOMIAL CHAOS EXPANSION OFPRESSURE HEADY(x,t) = In[K(x,t)]= f(x) + a(x)h(x,t)，The pressure head h(x,t) is a function of thesoil properties f(x) and a(x) indicated by Eqs. (1h≤0(8)and (2)， which can be formally expressed as some中國煤化工non-linear functional of the set {ζ;(w)} used to re- The.FYH.CNMHG inputs of the systempresent the soil randomness. It has been shown with known covarlance Lfr (x，y) and Ca(x,y)，that this functional dependence can be expanded in respectively. We can defineterms of polynomials in Gaussian random varia-bles，referred to as polynomial chaos ( Cameron f(x)=f(x)+f'(x )(9a)and Mar萬(wàn)方弊據316a(x)=a(x)+a' (x )(9b)g●e-Y'0)q")(11a)where f(x) and a(x) are the means of f(x) and(-r(o(--+28. +.o)].n/n--Q...a(x) respectively， f' (x) and a' (x) are the zero(11b)mean Gaussian random functions， which can beexpressed by the K arhunen-Loeve expansion.By substituting Eqs. (3)， (9) into Eq. (1) hm》ln =hr.δm.o(11c)and Eq. (2) and rearranging， we obtain(1ld)C(x,t)e -r(x,t)h(x.t2=V*h(x,t)+dtThe above equations is the perturbation equation ofthe stochastic saturated unsaturated flow problemVY(x,t)●V[h(x,t)+z]+g(x,t)e -r(x,t)for the order m≥0.(10a)6.KARHUNEN-LOEVE AND POLYNOMIAL{er(x,l)v[h(x,t)+z]. n}Ir,=- Q,(x,t) (10b)CHAOSBASEDPERTURBATIONEQUATIONS ,When the input random soil properties f(x)h(x,t)|r =hr; (x,t)(10c)and a(x) are perfectly correlated, the pressureh(x,t)|=o=hm;(x)(10d)head can be expressed by polynomial chaos in Eq.(7) asThe pressure head h(x,l) is a random function .h(x,t)=h<0)+h(1)+ h(2) +...(12a)because it depends on the randomness of the inputsoil parameters f(x) and a(x). We may expressthe soil parameter and pressure head related h")= 2h;(x,t)ζ;(12b)quantities C(x,t)， Y(x,t)，e-Y(x,t) ，er(x,t) and .h(x,t) as:h(2)=. 2 h, (x,t)55,(12c)C(x,t)=C")+C1+C+...whereh;(x,t) and h;(x,t) are deterministicY(x,t)=Y(0) +Y)+ Y(2>+...functions. By substituting Eqs. (12) into Eqs. (11)for m=0，we havee-Y=exp[- Yo°)][q°0) +q"" +q°2) +..].h9)-=v°h°0+VY(0●Vh'°)l+g'")(13a))ter=exp[- Y"°)][p°+p°>+ p2+...]，aY(0)h(x,l)=h'0)+ h'>+h2)+...g°)=ge Y(o)en'0)+(13b)dzwhere each term of C“, Y9,q”，p”and h') isQ°=-e Y(0)Qn/ep°)-vz●n(13c)proportional to the ith order of the standarddeviation of f(x) and a(x). Substituting aboveexpressions into Eqs. (10)，we have the followingh"=hr,(13d)equations for the order m≥0(0)|中國煤化工(13e)MYHCNMHGzELe- roOCQq(nh"f=]=+(13f)alSubstituting Eq. (12) into Eq. (11) for m=1 and乙vY(2●V[hm-)+zδm -.J]+multiplying the resulting equations by 5;,j=1，2,... ,taking expectation, and recalling that <ζζ;>=31 7sequentially. This property leads to the numericalefficiency especially for the solution of the higherah;l=V°h;+vY(°)●Vh;+g;, i=l, 2, 3,..order terms in the polynomial chaos expression ofJithe pressure head.(14a)For the uncorrelated f(x) and a(x)，thepolynomial chaos expression of pressure head i[vh;●n]Ir; =Q;(14b) Eq.<7) can be expanded by two sets ofuncorrelated normal random variables {ξ;} and {η:}h;lr.=0(14c) ash;l=o=0(14d) h(x,t)=h9)十h"+ h<2>+...(16a)Q,=<-[p:"v(h°>+z)]●n}Ir, /p'(14e) .h"= 2(hξ, + h?h?)(16b)g;=g●e -Y(°q"+vY!●V(ho) +z)-h?)= 22(n55E, +hi}ξη, + h?n.n,)(16c)i=1j=1e YO(C°q( +C)q0) )ah")"(14f)diBy substituting Eqs. (16) into Eqs. (11) form=0,Equations (14) is linear because all the terms in Q;1，2，multiplying與，η，ξξ，η: η;，i, j=1, 2,takingexpectationof the resultingand g; are linear function of h;.equations,and considering theorthogonal .Similarly, form=2，we haveproperties of ξand η;，we obtaine hi=vah, +vy0) .口h,+g,. i,j=1, 2, 3，。 ho)JteClr=v2h(0)+vY0)●Vh'°+g°)(17a)dt(15a)aY(0)(15b) g°)=ge-Y(0en0) +az(17b)h;lr=0(15c) Q°=-e-Y(0)Qn/ep'°)-vz● n(17c)hj|1=o=0(15d) h=hr.(17d)Q,={-[p"vh,+ p"”v(h0)+z)]. n-hm) =hmi(17e)[p"'Vh;+ pi?'V(h0+z)]●n}Ir; /p0)/2e=e-Yo)C0) q(°)(17f)(15e)for m= 0. By suing the same procedure as theg,=g●e -Y(°q+vY"●Vh;+vY{;●perfectly correlated case, we have((h0)+z)-e YO)[(C°q$°>+C"q"0) )h +ah{=V°hi+vY(°●Vhi+g{, i=1, 2, 3,...(18a)中國煤化工(C°q9 +C"q)" +C*qo)nh,(15f)(18b)dMYHCNMH G .Note that although the zero order Eqs. (13) is (h{)Ir,=0(18c)nonlinear because the coefficient e and Y[0) dependon the dependent functionho, Eqs. (14) and (15) (h})|r-o=0(18d) .arelinethose equations are solved318Qi=-[piV(h°+z)]●n]|r, /p")(18e) approximate thespatial derivatives by thecentral-difference scheme andthe temporalgf=g. e-Y(°q;+vY?●V(h°0)+z)一derivatives by the implicit method. The zero ordermean flow equation for both perfectly correlatede row(Coq:+Cq )h?(18f)nd uncorrelated case is the same. T his equation isdtnonlinear and thus needs to be solved in aniterative manner. Once the mean pressure headah?V2h?+vY(°●Vh?+g?，i=l，2, 3...h (0) is solved, the linear equations for the otheratperturbation terms can be solved sequentially and(19a)the coefficient matrix of the resulting systemequations are the same. This behavior of the[Vh?●n]|r,=Q?(19b)resulting perturbation equations renders efficiencyin the numerical method. Due to the symmetry of(h?)Ir =0(19c) the cofficientsh;(x.t) (i, j = 1, 2,.，Np) inpolynomial chaos expansion of the pressure head, .(h?)|=o=0(19d) only a half of h,(x,t) need to be determined.Because up to second order terms areevaluated inQ!=-[p?V(h°) +z)].n]|r,/p'°)(19e) this study, only h。(x.l) (i= 1, 2,..N,) areneeded for the calculation of the moment ofg?=g●e -Y(°q?+vY?. V(h0)+z)-pressure head and water content, which reduce thecomputational efforts significantly.Jh(0)e -YO)(C"q? +C%q"> )9Jt(19f)8. NUMERICAL EXAMPLES8. 1 Comparisom with a simple analyticalfor m=1，andsol utionWe compared our numerical results with thedhinV°hn +vY°) . Vh; +gy，analytical solution of Yeh et al. (1985b) for bothdperfectly correlated and uncorrelated cases. F romtheir analytical results， we can calculate the ratioi≥j=1, 2, 3,...(20a)rum of the standard deviation of the uncorrelated case[vh;●n]Ir; =Q,(20b)to the perfectly correlated case, rme =(σn)c/(σn)ar=(o}+ Suσ2 )1/2/(σj一S.σ。)，where (σn )me and(σn )ar are the standard deviation of pressure head,(h;)|r=0(20c)Siand S. are the suction head for the uncorrelatedand correlated cases in unbounded domain ,(h,)|.=o=0(20d)respectively. By substituting o}= 2.0，2= 20m?，andS=S.= 1. 09m into above formula, wefor m=2.. have rn = 11.2. The numerical result from theEquations (13)-(15) and Eqs. (17)-(20) are KL- based model in bounded domain for perfectlythe deterministic perturbation equations of thecorrelated and uncorrelated case is rm =14.2,stochastic saturated- unsaturated flow equation,which is compatible with the analytical result ofwhich can be salved by any preferred numericalYeh et al. (1985b) for unbounded domain. We alsomethod such as finite element method or finite calculated the result with the input parameterso, =different method.0.1,中國煤化工。= 1. 93m. We haveMYHCNMHGI resultof Yehet al. .7. NUMERICAL IMPLEMENTATION(1985b). In our proposed model, we obtainedrn =The numerical implementation is facilitated by1. 60，which is very close to the result from therecognizing that all perturbation equations have theanalytical model. Both oursame format except for the driving forces. We319Table 1 Comparison of the results from KL-based model with the analytical resultsσr。"176. 00.50. 001001. 000154. 00. 8890.87199. 90.6420. 64876. 80. 5480. 54023. 10. 3900.38676.80. 000010, 7160.7160.000050. 8540.8510. 000100.001002.4202. 4300. 002003. 5603. 3600. 003504.5804.4100. 005005. 3705. 2800.10.00100.0. 7960. 7940.20. 8511.0001.01. 2101.2101.51. 3802.01. 5401.540* roe and rn represent the ratio of standard deviation of the uncorrelated case to the perfectly correlated case from theanalytical model ( Yeh et al. 1985b) and the proposed KL based model, respectively. s is the suction of waternumerical results and the analytical results of Yeh eigenfunctions can be determined analytically. Weet al. (1985b) indicate that the correlation first try to show the validity of the proposedstructures of f (x) and a (x) have large affect on stochastic models and numerical implementation bythe calculated variance of pressure head. Table 1 comparing our numerical results with those fromshows the comparison of the result from the the conventional moment- equation based stochasticanalytical result of Yeh et al. and those by our models (Lu and Zhang 2002; Zhang and Lu 2002).Karhunen- Loeve expansion and perturbation based We consider a rectangle grid of 21 by 61 nodes in astochastic model (KL-based model) for different vertical cross section of 1.2m by 3. 6m. TheparametersS=S。=S，o}，and σ2，respectively. boundary conditions are specified as follows:The results of the KL- based model are very closed constant pressure head at the bottom with theto the analytical solution, which indicated that the pressureh = 0, a constant infiltration flux Q2at theour model is reliable.top，and no-flow at the left and right sides. The8.2 Com parison with the moment based modelinput parameters are given as f= 0.0, oj=0. 1, aIn this section we attempt to demonstrate the = 3m', σ = 25m-2，0,=0.3，0,=0.0，Q2=applicability of the developed stochastic model to- 1m/中國煤化工e for both a(x) andsaturated-unsaturated flow in hypothetical soils.f (x:YHCNMHGlepicts the first twoThelog-transformedsaturatedhydraulic moments of pressure nead at the steady state fromconductivity. f (x) and the pore size distribution theKL- basedmodelparameter a (x) are assumed to be second order moment- equation- basedstochasticstationary with a separable exponential covariance. ( moment- based model). The zero order meanInthis萬(wàn)方數握tion,theeigenvaluesand equations are the same from the two stochastic320models. In the KL-based model the second ordermean pressure head is calculated and it is showns0 [that the contribution of the second order terms tothe mean pressure head is small for the selectedstatistical parameters of the random functions.100+When we expand f (x) anda (x) with 40 terms inthe K-L based model, the pressure head variancefrom the two stochastic models are almost thesame.' I his comparison indicates the correctness ofMomentour proposed model because the moment- basedmodel has been validated before by Monte Carlo40， -80 -120simulations (e.g. ，Zhang and Lu 2002). Moreimportantly，theexample demonstrates the150efficiency of the proposed model in computationalefforts compared with the moment based model, at100least for the random parameters selected in theexamples. In the moment-based model, equationsforC(x，y;t, r), Cinr(x，y;l, τ), Cin(x，y; t，50 t- n10τ ) and the mean equation of pressure head have to三- n40be solved (Zhang and Lu 2002). Because thecovariance equations are functions of arbitrary two020space point x and y in the simulation domain, thecomputation of these equations is demanding. Inthe KL-based model， instead of solving theFig.1 Comparison of the mean pressure head andcovariance equations, the deterministic coefficientsstandard deviation calculated from the KL- basedh;(x,l)cndh;(x，t)，i=1,2,...Np,inthemodel with those from moment- based model. n10polynomial chaos expansion are solved with theand n40 represent the terms of the expansion ofsimilar form of equations as the covariance equationthe input parameters in the KL based modelin the moment-based model， where N, is the The non-flow boundary conditions are prescribednumber of terms in KL-expansion of the random at two lateral boundaries. The hydraulic head isinput function f (x) anda (x). For the case of n prescribed at the left and right boundaries at 10.nodes of the simulation domain， the ratio of the 5m and 10. 0m，respectively， which produces acomputational effort from the KL- based model to mean flow from the left to the right. The inputthe moment- based model is approximatelyβ=(n;十parameters are given as f= 0.0， o}=1, 0.=0.3，2XN,+N,)/(n;+3Xn )，where n;is the number 0,=0.0 ，and the correlation scale of f (x) is 1m.of iterations for the nonlinear equation of mean At the steady state the groundwater head ispressure head, n is the number of node in the uniformly distributed in the simulation domainsimulation. In the present example, n = 1281， N。with a constant gradient from the inflow boundary= 40，and n; = 6，we have β = 0.033. It is to the outflow boundary. Figure 2 shows theindicates that the proposed KI- based model is standard variation of water head for a profile alongcomputationally effective compared with the the water flow direction calculated from theMoment-based model, which is commonly used in KL- based model, moment based model and Montethe solution the stochastic problem of water flow in Carlo simulation. The maximum value of the watersoils.head中國煤化工dle of the inflow and8.3 Stochastic analysis of water flow in outfYHCNMHGhe water headat thesaturated soilstwo bounaaries are aetermnistic. In addition， theA two- dimensional domain in a saturated figure shows that the head variance derived fromheterogeneous porous media is considered. The the KL-based model is visually identical to theflow domain is a square of a size 10m X 10m，moment based model, all the of results from theuniformlyt熬據zed into 40 X 40 square elements.two models are very closed to the Monte Carlo321simulation result with 5000 realizations.based stochastic model. The simulation resultsfrom the KL-base and moment based models are i-dentical for the same perturbation order. Howev-0.05 Eer，the proposed KL-based model is much more ef-0.04ficient because only a few terms in KL ex pansion- 0.0are required for the input random porous media0.02 EKL-based modelproperties.Moment-based model0.01” Monte Carlo simulationREFERENCESFig.2 Calculated standard deviation from the KL- based[1] BRESLER E. and DAGAN G. Convective and porescale dispersive solute transport in unsaturated hetero-model, moment-based model for saturated casegeneous fields[J]. Water Resour. Res. ，1981，17(6):and compared with result from the Monte Carlo1683- 1693.simulation[2] CAMMERON R. H. and MARTIN W. T. The or-thogonal development of nonlinear functionals in series9. SUMMARY AND CONCLUSIONSof Fourier Hermite functionals[J]. Ann. Math. ,Based on the combination of Karhunen- Loeve .1947， 48: 385 -392.expansion of the input random functions, Chaos [3] DAGAN G. Flow and transport in Porous Formationsdecomposition of the water head and the[M]. New York: Springer- Verlag，1989 ,465.perturbation method, a stochastic numerical model [4] DAGAN G. and BRESLER E. Solute dispersion inforstochastictransientflowinunsaturated heterogeneous soil at field scale, I, Theo-ry[J]. Soil Sci. Soc. Am. J. ，1979, 43: 461-467.saturated-unsaturated zones is developed， whichcan also be used for the analysis of saturated or[5] GELHAR L. W. Stochastic subsurface hydrologyC M].Prentice-Hall, Old Tappan, N. J. ，1993.unsaturated flow. The input soil parameter is the .[6] GHANEM R. G. and SPANOS P. D. Stochastic finite :soil conductivity K ( x，t ) and the soilelements~A spectral approach [ M]. New York:characteristics0(x，t )，which are expressed by theSpringer Verlag, 1991, 214.Gardner- Russo model. The randomness of the soil[7] LIB. and YEH T. C. J. Sensitivity and moment ana-propertyis relatedtothe saturated soillyses of head in variably saturated regimes[J]. Adv.conductivity K,(x) and thesoilporesizeWater Resour. ，1998， 21: 477-485.distribution parametera (x). We assume thata (x) [8] LIEDL R. A conceptual perturbation model of waterand the log -transformed saturated hydraulicmovement in stochastically heterogeneous soils [J].conductivity f (x = ln K,(x) are spatial randomAdv. Water Resour.. 1994， 17: 171-179.functions with a normal distribution and known[9] LU Zi- ming and ZHANG Dong- xiao. Stochastic analy-sis of transient flow in heterogeneous, variably satu-covariances. By using the Karh unen-Loeve ex pan-rated porous media: The van Genuchten-Mualem con-sion of the input random parameters, polynomialstitutive model[J]. Vadose Zone Journal， 2002， 1:chaos expansion of pressure head，and the pertur-137-149.bation method, the governing equation of the sto- [10] MANTOGLOU A. A theoretical approach for mode-chastic salurated unsaturated flow and the corre-ing unsaturated flow in spatially variable soils: Effec-sponding initial and boundary conditions are re-tive flow models in finite domains and nonstationaritypresented by a series of deterministic partial dif-[J]. Water Resour. Res. ，1992, 28(1): 251-267.ferential equations in which the dependent varia- [11] MANTOGL.OU A. and GELLHARL. W. Stochasticbles are coef ficients of the polynomial chaos ex-modeling of large-scale transient unsaturated flow sys-pansion. The partial differential equations can be .中國煤化工.1987. 23)()3746.eous chaos[J]. Amer. J.solved subsequently and the random water headMYHCNMHGf unction can be represented by the pol ynomial cha-[13] YANG Jin-zhong， ZHANG Ren-duo and WU Jin-os ex pansion. The moments of the water head andquan. A analytical solution of macrodispersivity forother statistical properties can be evaluated by theadsorbing solute transport in unsaturated soils [J ].polynomial chaos ex pansion. The developed stoWater Resour. Res. ，1996a, 32(2): 355-362.chastic mdtF境re compared with the moment- [14] YANG Jin- zhong，ZHANG Ren-duo and WU Jin-322quan. Stochastic analysis of adsorbing solute transportneous soils, 1, Statistically isotropic media[J]. Waterin three dimensional heterogeneous unsaturated soilsResour. Res. ，1985a，21(4): 447-456.[J]. Water Resour. Res. , 1997, 33(8): 1947-1956.[19] YEH T. C.，GELHAR L.. W. and GUTJAHR A.[15] YANG Jin- zhong. ZHANG Ren-duo and YE Zi-tong.L. Stochastic analysis of unsaturated flow in heteroge-Stochastic analysis of adsorbing solute transport inneous soils， 2，Statistically anisotropic media withheterogeneous unsaturated soil[J]. Journal of Hydro-variable a[J]. Water Resour. Res. ，1985b， 21(4):dynamics, Ser B, 1996b， 8(3): 31-43.457-464.[16] YE Zi-tong， YANG Jin-zhong. Solute transport in va- [20] ZHANG Dong xiao. Nonstationary stochastic analysisdose zone: field studies[J]. Journal of Hydrodynam-of transient unsaturated flow in randomly heterogene-ics, Ser B，1995，7(4): 54-62.ous media[J]. Water Resour. Res. ，1999， 35(4):[17] YE Zi tong, YANG Jin- zhong. Experimental studies1127-1141.of solute transport through unsaturated natural soils [21] ZHANG Dong xiao and LU Zi-ming. Stochastic analy-[J]. Journal of Hydrodynamics, Ser B，1996， 8(1):sis of flow in a heterogeneous unsaturated- saturated31-38.system[J]. Water Resour. Res. ，2002， 38(2): 137-[18] YEH T. C.. GELHAR L. w. and GUTJAHR A.149.中國煤化工MYHCNMH G .