Investigation of solute transport in nonstationary unsaturated flow fields

Most theoretical work on analyzing plume spreading at the field scale in a partially saturated heterogeneous formation assumes weak stationarity of velocity field. While this assumption is not applicable to the case of a bounded flow domain, the nonstationarity in fluctuations of unsaturated velocity fields is induced through the presence of the boundary conditions in the flow domain. In this work, we attempt to quantify the large-time macrodispersion in nonstationary unsaturated velocity fields caused by the presence of a fixed head boundary condition. Application of the perturbation-based nonstationary spectral approach leads to an analytical expression for the macrodispersion describing the field-scale dispersive solute flux in terms of the statistical moments of two formation parameters: the Gardener’s parameter ( α) and the saturated hydraulic conductivity (Ks). The results predicted from the expression indicate that the enhanced unsaturated plume spreading can arise from a larger correlation scale of ln Ks or ln α fields. In addition, theα-parameter takes the role in reducing the field-scale plume spreading.


Introduction
The spatial variability in the specific discharge fields, arising from the heterogeneities of the geologic formations, enhances the spreading of nonreactive field-scale plumes in heterogeneous formations.The common theme in analyzing the transport processes of field-scale plumes for various stochastic methods is therefore to relate the statistics of the solute displacement or the solute dispersive flux to the statistics of the specific discharge fields.As such, the fieldscale spreading of solute plumes may be described through the macrodispersion tensor (e.g.Dagan, 1989;Gelhar, 1993;Rubin, 2003).
Most of stochastic analyses of field-scale solute transport processes have been performed under the assumption of weak stationarity of velocity field, which is the consequence of the assumptions of an unbounded flow domain and uniformity of mean flow.However, some real-world problems of solute transport in heterogeneous subsurface formations demand predictions over a finite flow domain.It has been recognized that the conditions of finite boundary could cause flow nonuniformity in the mean and in turn the nonstationarity in specific discharge field.The flow nonstationarity significantly impacts the spreading of the solute plume (e.g.Sun and Zhang, 2000;Wu et al., 2003;Dai et al., 2007;Lu et al., 2010).The practical implication of that is clear: the use of the approach appropriate for the unbounded domain may result in a significant error in quantifying the field-scale dispersion process in the nonstationary flow field.The field-scale transport processes are obviously influenced in more complicated ways by the nonstationary velocity fields.That is why the problem of plume transport by nonstationary groundwater flow at the field scale has so far attracted only limited attention in the groundwater hydrology literature.
The Eulerian analysis of field-scale solute transport in heterogeneous media is built based on the solution of the stochastic convection-dispersion perturbation equation.As such, the macrodispersive flux (Gelhar and Axness, 1983), an outcome of the correlation between the velocity field and concentration fluctuations, is introduced to quantify the spreading of the solute plume at the field scale.Comprehensive overviews of the construction of the Eulerian approach and its application to the analysis of the solute transport in heterogeneous media were given by Gelhar (1993) and Published by Copernicus Publications on behalf of the European Geosciences Union.Rubin (2003).It is well known that the correlation between the fluctuations of the velocity and concentration enhances the degree of spreading of solute plumes (field-scale dispersion or macrodispersion) in heterogeneous aquifers that is much greater than what would occur by local dispersion (e.g.Gelhar and Axness, 1983).The quantification of the field-scale dispersion created by the heterogeneities of the geologic formations is therefore rather important in the analysis of unsaturated solute transport.
In this paper, we attempt to quantify the field-scale dispersion at large times in nonstationary unsaturated velocity fields analytically in terms of the statistical moments of two formation parameters: the Gardner's parameter (Gardner, 1958) and the saturated hydraulic conductivity.Eulerian approaches are well suited for modeling large-time transport of solutes in heterogeneous media.Therefore, the quantification of the field-scale dispersion in this study will be carried out within the Eulerian framework.Note that in this study the nonstationarity in unsaturated velocity fields is introduced due to the presence of a fixed head boundary condition at an arbitrary depth in the flow domain and a constant flux at the land surface.The present investigation of field-scale unsaturated transport processes at large times can be considered as an extended work to the previous study by Chang and Yeh (2009) for bounded unsaturated flow processes in heterogeneous aquifers.The large-time macrodispersive flux is quantified by applying Fourier-Stieltjes integral representations of the randomly nonstationary fluctuations (Li and McLaughlin, 1991).The process of solute spreading in field porous formations is dominated by the spatial heterogeneity of the formation properties.We therefore conclude by analyzing the influence of the correlation scales of these two parameters on the spreading of field-scale unsaturated plume, which is the main focus of this work.
Existing stochastic studies (e.g.Rubin and Bellin, 1994;Rubin and Seong, 1994;Indelman and Rubin, 1996;Zhang, 1999;Sun and Zhang, 2000;Foussereau et al., 2000;Destouni et al., 2001;Wu et al., 2003;Hu, 2006;Dai et al., 2007;Russo and Fiori, 2009;Lu et al., 2010) on the issue of transport process of field-scale plumes in nonstationary groundwater flow fields were mostly carried out using the Lagrangian perspective.The analytical analysis of the fieldscale plumes in nonstationary unsaturated groundwater flow fields using the Eulerian concept has not been presented so far.This study attempts to analyze the spreading of plumes in field-scale unsaturated formation based on the Eulerian framework and nonstationary spectral approach.The findings presented in this paper may be of interest to researchers in seeking further research.

Formulation of stochastic unsaturated solute transport equation
We consider the case of the steady-state transport of nonreactive conservative solute plumes in variably saturated, heterogeneous formations.The flow system we investigate is of infinite extent along the horizontal direction, where a constant flux is introduced on the top and a prescribed capillary pressure head is specified at the lower boundary.The steadystate transport of plumes at a local scale modeled through the mass conservation equation is given by (e.g.Bear, 1979) where C is the solute concentration, q i is the i-th component of the specific discharge q, D ij are the components of the local dispersion coefficient tensor, and θ is the soil moisture content.
The effects of the pore-scale dispersion and the spatial variations in moisture content are negligible when compared with those of varying hydraulic conductivity.This simplifies Eq. ( 1) to Note that it has been mentioned in the works of Russo (1998) and Harter and Zhang (1999) that the field-scale dispersion at the large time is insensitive to variability in water content compared to the variability of ln K s .Harter and Zhang (1999) also mentioned that the preasymptotic macrodispersion is significantly larger in soils with variable water content than in those assuming constant content.The impact of spatial variability in water content is negligible for moist soils, but it could be significant for very dry soils (Harter and Zhang, 1999).The development of equations for the stochastic mean concentration and perturbation starts from the local transport equation (Eq.2).Consider that C and q i are realizations of the random variables represented, respectively, by a small perturbation expansion (Gelhar and Axness, 1983): (3) where the expected value operator is denoted by the angle bracket.The substitution of the perturbation expansions of Eqs.
(3) and (4) into Eq.( 2) produces the stochastic transport equation.The following form for the mean equation of solute transport is obtained after expanding terms in the stochastic transport equation and taking the ensemble average of them: where H denotes the average of capillary tension head.
The validity of the small perturbation approximation (i.e. the convergence of Eqs. 3 and 4) is preserved by σ 2 f (the variance of log saturated hydraulic conductivity) that should be small compared to unity (Gutjahr and Gelhar, 1981).However, the study of Monte Carlo simulations of two-dimensional flow through heterogeneous formations by Zhang and Winter (1999) confirmed the accuracy of the head moment solutions obtained from the application of the small perturbation approximation in σ 2 f at the variance up to 4.38.Similar comparison with Monte Carlo simulations reported in Guadagnini and Neuman (1999) yields accurate results (namely the statistics of hydraulic head) even for σ 2 f as high as 4 to 5 (strongly heterogeneous media).
Equations ( 5) and ( 6) are simplified by taking the mean specific discharge in the X 1 direction (i.e.q = ( q 1 , q 2 , q 3 ) = (q, 0, 0)) and approximating the pore-scale dispersion tensor in the form (Bear, 1979;Gelhar and Axness, 1983) where α L and α T denote the longitudinal and transverse porescale dispersivities, respectively, and |q| denotes the magnitude of q.As such, where the conservation for the fluid mass has been applied in the development of Eqs. ( 8) and ( 9).The cross-correlation term in Eq. ( 8), referred to as macrodispersive flux by Gelhar and Axness (1983), has been used to quantify the field-scale dispersion effect in saturated heterogeneous media.It creates additional mass transport arising from the correlation between the variation in specific discharge and concentration.Note that the term in Eq. ( 9), involving the product of specific discharge perturbation and mean concentration gradient, serves as the source in contributing to the output perturbations in concentration field.In other words, Eq. ( 9) links the output variation in concentration fields to the input variation in specific discharge fields.
In the following section, we attempt to quantify the macrodispersive solute flux in Eq. ( 8), describing the fieldscale spreading behavior, in terms of statistical properties of the specific discharge fields.This may be obtained from solving Eq. ( 9).

Concentration perturbation
It is recognized that there is a disparity in scale between the mean and the perturbation in the concentration field.Generally, the mean concentration field is slowly varying in space.The scale on which the perturbations in the concentration field fluctuate is much smaller than that related to the variation in the mean concentration field (e.g.Gelhar and Axness, 1983;Vomvoris and Gelhar, 1990).It is then possible to simplify the perturbation equation (Eq.9), by approximating the mean concentration gradient, the coefficient in Eq. ( 9), as a constant when solving Eq. ( 9).And within this framework, the solution of Eq. ( 9) allows for quantifying the macrodispersive flux term in Eq. ( 8).It must be noticed that the assumption of spatially uniform mean concentration cannot be made in the prediction of the spreading process of field-scale plumes near the source of plume where a sharp concentration gradient exists.In other words, the asymptotic transport relationship is applicable only after a substantial displacement distance.In addition, it is according to previous studies (Russo, 1996(Russo, , 1998) ) that the spreading of the fieldscale plume would therefore reach its large-time behavior as long as the lateral length scale, which is used to characterize the size of the solute body, is much larger than the scale of heterogeneity.
To solve Eq. ( 9), the coefficients q and q i in Eq. ( 9) must be specified first.Those can be determined simply from the work of Chang and Yeh (2009) in which the general expressions for the mean and perturbation of the specific discharge were developed.They considered the case where a constant specific flux q 0 is introduced on the soil surface X 1 = X L , and the prescribed pressure head φ 0 is specified at X 1 = X 0 .One can obtain the mean specific flux, q, from Chang and Yeh (2009) (by substituting their Eqs.40 and 41 into Eqs.37 and 38, respectively) as where q is negative for infiltration and In Eq. ( 11), α is the soil pore-size distribution parameter, α g = exp (< ln α >), R i are the components of the wave number vector R, dZ f (R) and dZ β (R) are the zero-mean random Fourier-Stieltjes increments of ln K s and ln α fields, respectively, and where F =< ln K s >.Note that the Gardner exponential constitutive model (Gardner, 1958) and the statistical independence of ln K s and ln α random processes have been assumed in the development of these expressions.In addition, X in Eq. ( 11) denotes the absolute position, not the separation lag.
An efficient approach for developing the analytical solution to the perturbation transport equation is through the Fourier-Stieltjes integral representation of a stationary or stationary-increment random field (e.g.Gelhar and Axness, 1983;Vomvoris and Gelhar, 1990;Rehfelt and Gelhar, 1992).As indicated in Eq. ( 11), the specific discharge perturbation is space-dependent, which is the cause of the effect of a finite flow domain (Chang and Yeh, 2009).This spacevariant coefficient will lead the output perturbation in concentration field in Eq. ( 9) to be nonstationary.Therefore, the infinite-domain spectral theory (stationary spectral theory) may not be directly applicable to the case of solute transport in nonuniform flow fields (such as Eq.9).The solution of Eq. ( 9) may be determined from the use of the nonstationary spectral representation (Li and McLaughlin, 1991) based on the Fourier-Stieltjes integral representation of a nonstationary random field in terms of a transfer function.
Note that the spatial variability of specific discharge, which is the outcomes of the spatial variations in ln K s and ln α, contributes to the spatial variation in concentration field.Motivated by that, the response to the linear Eq. ( 9), the concentration fluctuation, is separated into two parts: one representing the outcome of the variation in ln K s field and the other the variation in the ln α field.We then replace C in Eq. ( 9) by where Cf (X, R) and Cβ (X, R) are unknown transfer functions.The components in Eq. ( 13) represent the consequences of the spatial variations in ln K s , and ln α fields, respectively.Substituting Eqs. ( 10), ( 11) and ( 13) into Eq.( 9), separating the two components and making use of the uniqueness of the spectral representation leads Eq. ( 9) to the following two differential equations: The solutions of Eqs. ( 14) and ( 15) give the following two transfer functions, respectively, as where µ = α g α L .The fluctuations in concentration field can then be determined by applying Eqs. ( 16) and ( 17) into Eq.( 13).

Field-scale macrodispersion coefficient
The cross-correlation term, macrodispersive flux, in Eq. ( 8) can be determined from multiplying Eq. ( 13) by the complex conjugate of Eq. ( 11), taking the expected value, and applying the theorem of the spectral representation.The result is where and where S ββ (R) and S ff (R) are the spectral density functions of ln α and ln K s fluctuations, respectively.From Eqs. ( 8) and ( 18), the resulting ensemble average of concentration field is thus governed by the following convection-dispersion equation: with A f 33 = A f 22 and A β33 = A β22 .
To simplify the initial analysis within the Eulerian framework, we assume that the local process is isotropic (α L = α T ) and the medium is statistically isotropic.An exponential correlation function (Zhang and Winter, 1998;Zhang, 1999) is used to represent the correlation structure of random ln K s or ln α fields, which has the spectrum for the ln K s fluctuations or for the ln α fluctuations.In Eqs. ( 22) and ( 23), λ f and σ 2 f denote the correction scale and variance of the ln K s fluctuations, respectively, while λ β and σ 2 β denote the correction scale and variance of the ln α fluctuations, respectively.
With Eq. ( 22), the computation of the macrodispersivity integral Eq. ( 19), produced by the influence of the variation in ln K s fields, over the wave number domain yields where ξ = α g λ f .In general, the evaluation of Eq. ( 20) with Eq. ( 23) cannot be performed analytically.However, to take the advantage of the analytical solution, offering a clear insight for the role of the statistics of two formation parameters in influencing the large-time behavior of macrodispersion, we focus only on the case where α g α L 1.Note that the typical values for α L and α −1 would probably be from 10 −2 to 1 m (Gelhar et al., 1979;Matheron and de Marsily, 1980) and 0.2 to 2 m (Yeh et al., 1985), respectively.Under that conduction, the longitudinal macrodispersivity, caused by the effect of the variation of α, can be estimated as where ε = α g λ β and (α g ) is defined by Eq. ( 12).The plot of Eq. ( 24) presented as the function of the ln K s correlation scale for selected values of α g is shown in Fig. 1.The larger the ln K s correlation scale, the more the plume spreads, as indicated in the figure.An increase in the correlation scale of ln K s introduces a larger spatial consistency of fluctuations in the specific discharge above or below the mean specific discharge, and consequently produces a greater specific discharge variance.The fluctuations in specific discharge contribute to the irregular character of the concentration distribution in the field-scale unsaturated plumes and enhance the field-scale plume spreading.Also, from Fig. 1, if the correlation scales of ln K s and α L remain constant, the macrodispersivity will decrease with α g .A larger α g results in the smaller unsaturated hydraulic conductivity variability and, in turn, the smaller specific discharge variability (Yeh et al., 1985).Consequently, less spreading of the solute will take place.
Figure 2 illustrates how the longitudinal macrodispersivity component in Eq. ( 24) varies with the correction scale of the log α-parameter.The log α-parameter correlation scale is of importance in increasing the variability of unsaturated hydraulic conductivity and thereby the variability of specific discharge (Chang and Yeh, 2009), which enhances the  field-scale plume spreading.As expected, the spreading of the solute plume is correlated inversely with the α-parameter for fixed values of λ β and α L .The plot of dependence of longitudinal macrodispersivity upon the position is illustrated in Fig. 3 based on Eq. ( 25).The longitudinal macrodispersivity component increases rapidly with the position at the beginning, then has practically attained its asymptotic value when the plume is close to the downstream boundary of the constant head.This behavior is the outcome of the spatial variations in specific discharge field.

Concluding remarks
The work uses the nonstationary spectral perturbation techniques to develop a closed-form expression quantifying the field-scale plume spreading in a partially saturated heterogeneous aquifer.This expression, related to the statistics of two formation parameters, i.e. ln K s and ln α, has allowed to investigate how these statistical properties influence the spreading process of the field-scale unsaturated plume.
Our results indicate that the field-scale dispersive solute flux increases with the variabilities of these two parameters.The correlation scales of these two parameters influence the spreading of the field-scale unsaturated plumes positively.In addition, the α-parameter is of great importance in reducing the field-scale plume spreading.

Fig. 1 .
Fig. 1.Dimensionless component of the longitudinal macrodispersivity as a function of dimensionless ln K s correlation scale.

Fig. 2 .
Fig. 2. Dimensionless component of the longitudinal macrodispersivity as a function of dimensionless ln α correlation scale.

Fig. 3 .
Fig. 3. Dimensionless component of the longitudinal macrodispersivity as a function of dimensionless position, where the dimensionless flow domain we investigate corresponds to 0 < (X 1 − X 0 )/α L < (X L − X 0 )/α L , and (X 1 − X 0 )/α L = 0 corresponds to the location of the bottom of the flow domain.