Veena S Soraganvi

Rachid Ababou

M.S Mohan Kumar

Effective flow and transport properties of heterogeneous unsaturated soils

Keywords:

The heterogeneity of field scale soils poses a challenge to predictive large scale flow and transport modeling. The theory of effective macroscale parameters holds good and is applicable in dealing with such problems. But the validity of the analytic stochastic solutions obtained for randomly heterogeneous soils is debatable, as the test cases under which they are validated are of limited scope due to linearization and perturbation approximations. In this study, samples of heterogeneous soils are generated using sets of spatially correlated random field parameters that are either geometrically isotropic, or else, geometrically anisotropic with either horizontal or vertical stratification (perfect or imperfect). Several combinations of ratios of correlation length and capillary dispersion lengths are considered. Numerical simulations of unsaturated flow are performed on each randomly heterogeneous soil sample. The principal components K ii (Ψ) of the macroscale effective unsaturated conductivity are then obtained as a function of the mean suction Ψ of the sample. They are compared to stochastic spectral perturbation theory, and to a probabilistic semi-empirical Power Average Model (PAM). They are also compared with arithmetic, geometric and harmonic mean conductivity-suction curves. The numerically upscaled principal conductivity curves match quite well the PAM, better than the classical means (Arithmetic, Geometric, Harmonic), and also somewhat better than the curves obtained from stochastic spectral perturbation theory. It is observed that the upscaled principal components K ii ( ), obtained numerically and with the PAM along directions "i " orthogonal/parallel to perfect stratification coincide with the harmonic/arithmetic mean curves at low suctions (i.e., near saturation), but deviate from it and come closer to the geometric mean at higher suctions. The PAM appears suitable for generation of approximate upscaled conductivity curves, e.g., for obtaining the mesh-scale or block-scale conductivity curves in large scale simulation codes. Transient solute transport simulations are then performed on the detailed random velocity fields obtained from the steady state simulations of unsaturated flow in the randomly heterogeneous soil samples. Snapshots of solute concentration C(x,z,t) are taken at different times. The temporal evolution of spatial moments of concentration is analyzed in order to characterize the macroscale advection and dispersion of the unsaturated concentration plume, and in particular, its macro-dispersion coefficient (D) and dispersivity length scale (A). For the synthetic soil samples considered in this study, the macro-dispersive spreading of the solute is stronger for flow parallel to vertical stratification, compared to flow perpendicular to horizontal stratification, and also, compared to flow in statistically isotropic non-stratified soil.

Introduction

Field soils are inherently heterogeneous, exhibiting significant spatial variation of hydraulic properties even within a given type of soil [START_REF] Warrick | Spatial variability of soil physical properties in the field[END_REF]. In recent years, soil scientists are being increasingly aware of the fact that small-scale heterogeneity affects the larger-scale flow and transport phenomena. The importance of considering spatial variability and assessing flow in the unsaturated zone has been examined in several theoretical works (eg. Yeh et al., 1985aYeh et al., , 1985bYeh et al., , 1985c ; ;Mantoglou and Gelhar, 1987a, 1987b, 1987c ). To deal with these issues, one way is to include the detailed heterogeneity of the system in numerical simulations of the flow. But obtaining the detailed hydraulic properties of the geologic formation at many points in the field is expensive and sometimes impossible. Despite recent progress in massively parallel Richards flow solvers (e.g. [START_REF] Orgogozo | An open source massively parallel solver for Richards equation: mechanistic modelling of water fluxes at the watershed scale[END_REF], the computational task may be impossible in practice for large scale / basin scale simulations, requiring possibly billions of numerical cells. Thus, even with future progress in computational capabilities, the direct simulation of subsurface soil processes at basin scales of hundreds of squared kilometers horizontally and tens of meters vertically seems out of reach. power averaging coefficient for the ith direction (dimensionless exponent) q or q i Darcy areal flux density (also known as "Darcy velocity ") [m/s] cross-correlation between the random log-parameters (a(x), f(x)) S S specific storage coefficient [(m 3 /m 3 )/m] or [1/m] suction head (-h) An alternative approach is to define an equivalent homogeneous medium with the upscaled (effective) flow and transport properties that represent the heterogeneous medium (e.g., [START_REF] Yeh | One dimensional steady state infiltration in heterogeneous soils[END_REF]. In particular, the stochastic approach requires only a statistical description concerning the degree of heterogeneity and the spatial structure of the soil. Only a few statistical parameters are needed for each soil property: mean, variance, and spatial structure represented by three auto-correlation length scales (one in each direction). It is expected that upscaled hydrodynamic properties can be characterized as functions of these statistical descriptors.

Nomenclature

Several authors have analyzed the spatial distribution and heterogeneity of the nonlinear unsaturated conductivity-suction curve K( , x ) for soils and shallow sedimentary formations, for example [START_REF] Hills | The second Las Cruces Trench experiment: experimental results and two dimensional flow predictions[END_REF] , Russo (1983) , [START_REF] Warrick | Spatial variability of soil physical properties in the field[END_REF] ; [START_REF] Wierenga | The Las Cruces trench site: experimental results and one dimensional flow predictions[END_REF] , and many others. It is found generally that the conductivity-suction curves, or their logarithm, are spatially variable both in slope and intercept (LnKs). Another thing is that, the slope of the conductivity-suction curve tends to be positively cross-correlated with the saturated conductivity Ks or Ln-Ks (the intercept).

Thus, [START_REF] Polmann | Application of Stochastic Methods to Transient Flow and Transport in Heterogeneous Unsaturated Soils[END_REF] determined the cross-correlation coefficient between intercept (LnKs) and slope ( ) of log-conductivity vs. suction curves for two heterogeneous soils, the Zohar soil of Russo (1983) , and the Maddock soil of [START_REF] Carvallo | Spatial variability of in situ unsaturated hydraulic conductivity of Maddock sandy loam[END_REF] . In both the cases, the cross-correlation was found to be significantly positive ( = + 0.26 for Zohar, = + 0.83 for Maddock)., See also Mantoglou and Gelhar, 1987 ), and [START_REF] Polmann | Application of Stochastic Methods to the Simulation of Large-Scale Unsaturated Flow and Transport[END_REF] To sum up, the cross-correlation ∈ [0, 1]. Two special cases can be usefully tested in unsaturated flow-transport simulations: the "neutral " case ≈ 0 and the perfectly correlated case ≈ + 1. In this study, we will assume ≈ 0 in the numerical experiments. Yeh et al., 1985aYeh et al., , 1985bYeh et al., , 1985c ) ) solved the steady state stochastic flow equation for infiltration in porous media with random field soil parameters. The stochastic flow solution was obtained for small perturbations, using linearization approximations and Fourier spectral decomposition. Then, upscaling was obtained by averaging the resulting random field solution, hydraulic gradient vector J(x) and the Darcy "velocity " (flux density) vector q(x), leading to a macroscale Darcy equation relating mean flux to mean gradient (from which effective unsaturated conductivity is deduced). Several additional assumptions were used: infinite domain, small variability, statistical homogeneity and ergodicity. Due to linearization, the solution is valid only for a limited range of suction and moisture content.

Several numerical investigations ( Ababou et al., 1988 ;Ababou & Gelhar (1988) ; [START_REF] Mccord | Application of second-type boundaries in unsaturated flow modeling[END_REF] ; [START_REF] Polmann | Stochastic modeling of large scale flow in heterogeneous unsaturated soils[END_REF] confirmed the stochastic solutions qualitatively for unsaturated flow. Though several of the numerical studies in the literature involve ensemble averaging, single realization flow simulations can also be conducted. Thus, [START_REF] Ababou | Numerical simulation of three dimensional saturated flow in randomly heterogeneous porous media[END_REF] have shown that the statistical results from single realizations of groundwater flow are meaningful if the following criteria are satisfied: (i) mesh size smaller than spatial correlation scale ( "statistical resolution "), and (ii) domain size larger than spatial correlation scales ( "statistical sampling "). Similar criteria were applied for simulating unsaturated moisture plume migration in a 3D domain ( Ababou et al. 1988 ;Ababou et al., 1992 ). In that way, a single realization can be used, instead of ensemble averaging over many replicates.

A number of field works and numerical studies have been developed concerning the applicability of the ADE (Advection-Dispersion Equation) to unsaturated solute transport at the field scale and concerning the effects of soil heterogeneity, and the scale effects on the dispersion process. In particular, the longitudinal dispersivity length scale (named "Azz " if the mean flow is vertical) appears to be scale-dependent, i.e., its effective value depends on the space-time scale of interest.

In some studies, the soil parameters were assumed to be purely random variables, constant in space (they can be viewed as uncertain parameters). For example, [START_REF] Harter | Stochastic analysis of solute transport in heterogeneous, variably saturated soils[END_REF] performed Monte Carlo simulations to analyze the spatial structure and evolution of a solute in the unsaturated zone, assuming random but spatially constant moisture content.

But more generally, unsaturated flow-transport phenomena have been studied by considering soil heterogeneity to be spatially random, i.e., with some soil parameters being random functions of space, with a given auto-correlation structure (as will be assumed in this work). Thus, [START_REF] Russo | Stochastic analysis of simulated vadose zone solute transport in avertical cross section of heterogeneous soils during nonsteady water flow[END_REF] analyzed the movement of a conservative tracer in a hypothetical heterogeneous unsaturated domain for a single realization of random field of Ks (saturated conductivity). An analytical expression was obtained for the longitudinal macro-dispersivity, as a function of the standard deviation of the random Ks. This study, among several others, showed that longitudinal macro-dispersivity is significantly larger than transverse macro-dispersivity.

In [START_REF] Russo | Numerical analysis of flow and transport in a variably saturated bimodal heterogeneous porous media[END_REF] , the random soil properties were assumed bimodal: a spatially distributed background soil contained randomly dispersed clay lenses (with same mean and variance but different two-point covariances). In [START_REF] Khaleel | Upscaled flow and transport properties for heterogeneous unsaturated media[END_REF] , effective properties of steady flow and transient solute transport (non-reactive) were obtained by Monte Carlo simulations, for the coarse-textured sand of the Hanford site with van-Genuchten constitutive relations . Interestingly, flow upscaling was performed for a large range of tensions, up to Ψ ≈ 25 m in their work.

Finally, we note that effects other than heterogeneity can play a role, such as deformation (which is not considered in this work). Recently, [START_REF] Wu | Numerical Modelling of consolidation-induced solute transport in unsaturated soil with dynamic hydraulic conductivity and degree of saturation[END_REF] introduced a time-variable porosity in Richards' equa-tion to model the deformation (consolidation) in a homogeneous soil. Their unsaturated solute transport results are compared to the incompressible case, indicating similar results except, however, near soil surface.

Upscaling the unsaturated conductivity to obtain an "effective " or a "macroscale " conductivity-suction curve has been attempted via various approaches (stochastic or not), sometimes with semi-analytical approaches. For example, [START_REF] Warrick | Effective unsaturated hydraulic conductivity for one dimensional structured heterogeneity[END_REF] upscaled the pressure-dependent hydraulic conductivity K(h,z) for steady vertical infiltration, in a soil column with a periodic pattern of N alternating layers of sand & sandy loam (an ODE solver was used). Two limiting values of K EFF (h) were obtained for small and large cell thickness. The results were, respectively: harmonic average K H (h), and arithmetic average K A (h). Interestingly, these curves are special forms of power averages (to be discussed later in this work).

Effective flow theories are potentially very useful in soil hydrology, as they allow to focus on soil structure (e.g., stratification) rather than fine details. However the test cases under which they are validated are limited. There are few studies concerning the validity of stochastic solutions for infiltration in relatively coarse sandy soils, where gravity takes over capillary dispersion at large times or quasi-steady conditions. Numerical upscaling studies have been developed to justify the theoretical results of Yeh et al., 1985aYeh et al., , 1985bYeh et al., , 1985c and others, based on Gelhar's stochastic spectral perturbation approach, but there have not been many comparisons with other upscaling theories or models, such as the Power Averaging Model (PAM) [START_REF] Ababou | Probabilistic Analysis of Unsaturated Conductivity: Anisotropy, Crossing Point, and Upper Envelope: Section 6.3, Chapter 6 (Stochastic Analysis of Flow and Transport)[END_REF] . In this work, we focus on both the semi-analytical power average expression of [START_REF] Ababou | Probabilistic Analysis of Unsaturated Conductivity: Anisotropy, Crossing Point, and Upper Envelope: Section 6.3, Chapter 6 (Stochastic Analysis of Flow and Transport)[END_REF] and on Yeh et al., 1985aYeh et al., , 1985bYeh et al., , 1985c' s stochastic spectral perturbation solution for effective unsaturated conductivity.

Capillary dispersion plays a major role in moisture migration phenomena. The recent paper by [START_REF] Soraganvi | Analysis and upscaling of unsaturated flow through randomly heterogeneous soil[END_REF] demonstrates the effect of a dimensionless capillary parameter introduced by them (the product Z ) on the behavior of unsaturated flow in randomly heterogeneous soils with exponential conductivity curve ( is Gardner's parameter). In their work, or its mean < >, was interpreted as an inverse mean capillary dispersion length scale: = 1/ CAP . The dimensionless ratio Z / CAP was introduced, where Z is the vertical correlation length of the soil (stratified or not). They used the flow code "UFT2D " [START_REF] Soraganvi | Modeling of flow and advection dominant solute transport in variably saturated porous media[END_REF] to study unsaturated flow upscaling at steady-state for randomly heterogeneous soils, for several values of ratios. The upscaled soil behaved as an equivalent homogeneous medium with suction-dependent anisotropy, depending on parameters K S and . A recent review of this and other works on nonlinear unsaturated properties in homogeneous as well as layered or heterogeneous soils, can be found in Ababou (2018, Vol.1, Section 5.3.4 & 5.4.8).

In the present paper too, we develop a set of unsaturated flow simulations, which are analyzed based on these capillary concepts. The dimensionless product Z is denoted here , and it is named "capillary ratio ". Since the mean Gardner parameter < > is the inverse capillary length scale, = 1/ CAP , the capillary ratio is also defined here as:

= Z / CAP . We will let " " vary from ≤ 1 to > 1, in order to analyze its effect on suction and moisture distribution, and on the upscaled, effective conductivity curve.

In this work, we will introduce a novel study of unsaturated flow and dispersive solute transport in perfectly and imperfectly stratified random soils, for various values of the capillary ratio . Since the mean unsaturated flow is predominantly vertical (Z) in soil hydrology, we will neglect horizontal flow (Y), to focus essentially on the 2D case (X, Z). Thus, flow and transport in imperfectly stratified soils will be studied here in 2D vertical cross-sections. The 1D case (Z) is also relevant, as it results from the special case of "perfectly stratified " heterogeneity: thus, vertical 1D flow & transport phenomena, orthogonal and parallel to layering, will also be studied. Solute transport simulations will be analyzed through spatial moments, leading to new results on unsaturated macro-dispersion.

The paper is organized as follows.

In the first part, we try to address some of the issues mentioned above, by implementing numerical experiments considering different combinations, or different capillary ratios (correlation length / capillary dispersion length). Samples of heterogeneous soils are generated using sets of spatially correlated random field parameters that are either geometrically isotropic or else, geometrically anisotropic, the latter with either perfect or imperfect stratification. Unsaturated flow simulations are performed on the 2D or 1D soil samples, and the numerically upscaled, equivalent unsaturated conductivity of the samples, are analyzed in terms of capillary ratio and statistical anisotropy. The effective conductivity curve seems well captured by the semi-empirical Power Average Model, and also to some extent, by the theoretical spectral perturbation solutions.

Transient solute transport simulations are then performed using the detailed, steady-state, unsaturated random velocity fields obtained numerically from the previous randomly heterogeneous soil samples. Snapshots of solute concentration C(x,z,t) are taken at different times. The temporal evolution of spatial moments of concentration is analyzed, to characterize the mean advection and dispersion of the unsaturated concentration plume, and in particular, its macro-dispersion coefficient and macro-dispersivity length scale.

Theory

Effective flow equations

Infiltration through a heterogeneous unsaturated medium is governed locally by Darcy's law with pressure dependent, spatially variable hydraulic conductivity K(h,x):

( ) = -( ℎ, )∇ ( ℎ ( ) + ) (1)
where q is Darcy velocity or flux density vector [m/s], K is unsaturated hydraulic conductivity [m/s], h is pressure head [m], z is the vertical elevation (height) directed upwards, and the symbol ∇ ( "nabla ") applied to a scalar field, represents the gradient operator.

The steady-state mass conservation equation for the heterogeneous soil is then, locally:

-∇ ( ) = ∇ { ( ℎ, ) ∇ ( ℎ ( ) + ) } = 0 (2)
where the symbol ∇ ( "nabla "), applied to a vector field, represents the "divergence " operator. Note that the pressure head formulation ( "h ") can be replaced by a suction head formulation ( " ") by letting = -h everywhere.

The heterogeneous unsaturated conductivity curve K(h,x) is represented by a few spatially variable parameters p 1 (x), p 2 (x), etc., that is: = ( ℎ ; 1 ( ) , 2 ( ) , …) . If each parameter p i (x) is considered to be a random function of space, Eqs. ( 1) and (2) form a system of stochastic PDE's with random field coefficients, which are nonlinear because "K " is not only variable in space but also pressure dependent at each point in space.

The random curves K ( h, x ) can be "averaged " (homogenized) as a single "effective " conductivity curve K (h) representing the combined effects of the random curves K ( h, x ) distributed throughout the soil domain. The pressure "h " in K (h) stands for mean pressure ⟨h ⟩. If the soil structure is geometrically anisotropic, K (h) turns out to be a tensor function of pressure, as indicated in several above-discussed works [START_REF] Gelhar | Stochastic Subsurface Hydrology[END_REF] ., Yeh et al., 1985aYeh et al., , 1985bYeh et al., , 1985c ) )., [START_REF] Mantoglou | A theoretical approach for modeling unsaturated flow in spatially variable soils: effective flow models in finite domains and non stationarity[END_REF] ; Mantoglou and Gelhar, 1987a, Mantoglou and Gelhar, 1987b), Polmann et al., 1991 ). In addition, Bear and Bachmat (1986) discuss the anisotropic and the tensorial character of unsaturated permeability based on averaging flow over a bundle of capillary tubes, where tubes represent "pores ". They find that the absolute conductivity K (h) (rather than the "relative conductivity " K (h)/Ks) should be considered as the effective tensorial quantity. The tensorial effective K (h) will also be denoted K ̲ ̲ e (h) below.

Accordingly, the "effective " or "macroscale " Darcy equation is written as follows, in terms of mean flux ⟨q⟩ and the mean pressure ⟨h⟩ :

⟨ ⟩ = -( ⟨ℎ ⟩) ∇ ( ⟨ℎ ⟩ + ) (3)
where K ̲ ̲ e ( ⟨h ⟩) is the second rank tensorial, anisotropic, pressuredependent effective unsaturated hydraulic conductivity (sometimes named "permeability " for short). For steady flow, ensuring mass conservation at the macroscale (or in the mean sense) yields the following mean conservation equation:

-∇ ⟨ ⟩ = ∇ ( ⟨ℎ ⟩) ∇ ( ⟨ℎ ⟩ + ) = 0 (4)
Note: the brackets ⟨…⟩ designate mean quantities (ensemble averages or spatial averages depending on context); however, we will later drop these brackets to alleviate notations.

The previous equations can all be expressed in terms of suction =h, rather than pressure h. Large suctions correspond to very dry media. Various interpretations of suction, and the physical meaning of very large suctions ( ≈ 10 m or much more), are discussed in [START_REF] Ababou | Capillary Flows in Heterogeneous and Random Porous Media[END_REF].

Eqs. ( 3) and ( 4) , together, constitute the system of PDE's to be solved for unsaturated flow at the macroscale. The spatially variable coefficients have "disappeared " in this upscaled version, but the effective "permeability " tensor function K ̲ ̲ e ( ⟨h ⟩) now needs to be characterized quantitatively in terms of the geometry and variability of local soil parameters. That is the purpose of upscaling theories.

Spectral perturbation approach

In this method the effective unsaturated conductivity is obtained in two steps. First, perturbation and linearization approximations are implemented for Eqs. ( 3) and (4) , and these are solved approximately for the zero-mean fluctuations of the stochastic suction field in Fourier space (using the Wiener-Khinchine spectral representation of random fields). Secondly, averages are calculated to relate the mean flux to the mean suction gradient, from which an approximate effective conductivity is finally obtained (using inverse Fourier transforms for the final results). The analytical expressions derived by Yeh et al., 1985, Yeh et al., 1985 ) concern steady flow in spatially correlated random media, assuming "perfectly stratified layering ", i.e., with correlation lengths X → ∞ and Y → ∞ in horizontal directions (X and Y). The final results do not depend on X and Y , but they depend on the finite vertical correlation scale Z < ∞, which can be viewed as a statistical layer thickness (although the random field description of vertical heterogeneity is not based on discrete layers). [START_REF] Gardner | Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[END_REF] relation is used for the nonlinear local unsaturated conductivity, assuming that the two random field parameters LnK S and are statistically independent (not cross-correlated). Under these assumptions, the effective conductivity along the principal directions, vertical (V) and horizontal (H), are given by [START_REF] Khaleel | Upscaled flow and transport properties for heterogeneous unsaturated media[END_REF]

) = exp - 2 2(1 + ) - - (2 + | |) 2(1 + ) 2 | | (5) = exp + 2 2(1 + ) - - (2 -| |) 2(1 + ) 2 | | (6)
where F is the mean and 2 the variance of Ln K S ( ⃗ x ) ; B is the mean and 2 is the variance of β( ⃗ x ) , and |H| is the mean "tension " or mean suction (also expressed as ⟨ ⟩ elsewhere).

A major drawback of this spectral perturbation method is that the results, like Eqs. ( 5) and (6) , are valid only for small standard deviations of the log-parameters ( f ≤ 1, ≪ 1), and then, only for a limited range of mean suctions (|H| or ⟨ ⟩ ), up to suctions not much greater than 1/B or so, where B is the mean value of the parameter. A piecewise linear approach was suggested by [START_REF] Soraganvi | Analysis and upscaling of unsaturated flow through randomly heterogeneous soil[END_REF] for upscaling high tension ranges up to 25 m (for a coarse sand) using the van Genuchten conductivity model, with (K S, v, n v ) as random parameters.

Power averaging theory

The Power Averaging theory, called "PAM " for Power Averaging Model (or sometimes called "RA model ") was developed around 1991 by R. Ababou and collaborators, based on a probabilistic theory of random conductivity curves. The theory was described in detail in Ababou, 1991b ) and Ababou et al., 1993 also contains early numerical verification tests, and[START_REF] Bagtzoglou | Effective hydraulic property calculations for unsaturated, fractured rock with semi-analytical and direct numerical techniques: review and applications[END_REF] contains a review of analytical and probabilistic approaches, as well as numerical tests on unsaturated fractured materials.

The probabilistic theory leading to the Power Average Model was framed in ensemble space rather than Euclidian space. It was assumed that the local conductivity-suction curves of the random medium are of exponential form (Gardner's model):

( ; ) = ( ) exp { -( ) } ⇔ ( , ) = ( ) -( ( ) ) (7)
where K S ( x ) and ( x ) are random positive coefficients, considered as random variables (rather than random fields), and is the local suction head in meters (the lower case letter is used for the local suction). The coefficients K S ( x ) and ( x ) may have a log-normal distribution, or any other positive distribution. In the analytical developments below, they are assumed log-normal, so that f(x) = Ln{Ks(x)} and b(x) = Ln{ (x)} are Gaussian random variables, which can be taken cross-correlated to each other. Thus, is the cross-correlation coefficient between Ln{Ks(x)} and Ln{ (x)} at any point (as discussed earlier, we expect ∈[0, 1] for most soils).

In the Power Averaging Model "PAM ", the effective conductivitysuction curve is obtained by calculating the ensemble power average ( "PA "):

̂ ( ) = ⟨ ( ; ) ⟩ 1 ∕ (8a)
The effective "PA " conductivity function on the left-hand side was obtained from probabilistic ensemble averaging calculations, exploiting the properties of Gaussian and log-normal variables. Explicit results were obtained for Gardner's exponential model ( Eq. ( 7)), and LogNormal distributions of {K S (x), (x)}. Euclidian space (x) plays the role of ensemble space in these probabilistic calculations. The resulting power average conductivity curve can be expressed as:

(Ψ) = exp 1 2 ( Ψ 2 + Ψ + ) (8b)
where: = 2 ; B i = -2(1 + p i a f ); = 2 ; and: Ψ = G ⟨ ⟩ is the dimensionless mean suction head (capital letter Ψ).

Equivalently, (Eq. ( 8b)) can be expressed as follows, in terms of the "naïve " curve K 0 ( Ψ) ≡ K G exp( -Ψ):

̂ ( Ψ) = exp [ -Ψ] exp 1 2 2 ⟨ ⟩ 2 Ψ 2 -2 ⟨ ⟩ ln Ψ + 2 ln (8c)
Note that the second exponential expression is a multiplicative correction of the "naïve " curve K 0 ( Ψ) ≡ K G exp( -Ψ), where Ψ is the dimensionless mean suction. In Eqs. (8 a,b,c), the curve K ii ( Ψ) is the principal effective conductivity component along the direction X i . The parameter p i is the power averaging exponent used for obtaining the effective conductivity along the direction X i .This power exponent, p i , may be treated at this point as a semi-empirical adjustable parameter.

In addition, the statistical parameters in Eqs. ( 8) and ( 9) are defined as follows:

= -( , ) is the slope of the log-conductivity / suction curve, G is the geometric mean of ( x ), that is, G = exp { E (ln ( x ))}. b is the standard deviation of log-coefficient b(x) = Ln{ (x)}; f is the standard deviation of saturated log-conductivity f(x) = Ln{K S (x)};

is the cross-correlation coefficient between the random logparameters ( b( ⃗ x ) , f ( ⃗ x ) ) .

The averaging powers p i can be chosen tentatively as follows (a more complete conjecture is available for some random spatial structures [START_REF] Ababou | Random porous media flow on large 3-D grids: numerics, performance, and application to homogenization[END_REF], Appendix B)):

• p i = 0 (geometric mean) in each direction X i if the medium is perfectly isotropic;

• p i = -1 (harmonic mean) in the direction orthogonal to layers (e.g. vertically);

• p i = + 1 (arithmetic mean) in the direction parallel to layers (e.g. horizontally).

However, in this work, we will calculate the p i 's from the numerical experiments by numerically fitting the power averaged curves (instead of applying the above averaging powers).

Compared to the above-described Power Average Model, the stochastic spectral perturbation results reviewed earlier Eqs. ( 5) and ( 6) can also be expressed, similarly, as a function of a direction-dependent parameter which is analogous to p i . This new parameter is denoted as p * i (p i with a "star "), by analogy with the power exponent p i of the PAM model. With this, it can be seen that the spectral perturbation result of Eqs. ( 5), ( 6) can be expressed equivalently as:

K ii ( Ψ) = K G exp p * i σ 2 b Ψ 2 -2 1 + p * i ρσ b σ f Ψ + p * i σ 2 f ∕2 (9)
where p * i is now the directional parameter of the spectral perturbation theory. Remarkably, the Power Averaging Model of Eqs. (8) yields an effective conductivity curve K ii ( Ψ) of the same form as the spectral perturbation theory (Eq. ( 9)), although the coefficients of the quadratic Ψ-polynomial are somewhat different (e.g., comparing Eq. (8b) and Eq. ( 9)).

In the general case of 3D statistically anisotropic soils, the semiempirical parameters of "PAM " are the 3 averaging powers p i ( i = 1, 2, 3) which can be taken different horizontally and vertically. For a horizontally stratified medium, taking p i = + 1 horizontally and p i = -1 vertically, yields the arithmetic and harmonic means of the nonlinear K( ,x) curves. With this choice, the PAM becomes similar (although not exactly identical) to the effective conductivity model initially proposed in [START_REF] Mualem | A modified dependent domain theory of hysteresis[END_REF] , and further investigated by [START_REF] Mccord | Application of second-type boundaries in unsaturated flow modeling[END_REF] , . However, [START_REF] Mualem | A modified dependent domain theory of hysteresis[END_REF] did not include the cross-correlation parameter , which is physically meaningful, and can play an important role on the results.

It should be recognized that the PAM theory is based on unstructured, purely random soil parameters, with a prescribed variability, but without a prescribed spatial auto-correlation (as opposed to spatially correlated random fields). In theory, the soil parameters of the PAM have zero spatial correlation scales. In some other works, on the contrary, the unsaturated flow/transport parameters are assumed to be "random constants ", which is tantamount to having spatially correlated parameters over long or infinite distances (an example previously cited in the literature review is [START_REF] Harter | Stochastic analysis of solute transport in heterogeneous, variably saturated soils[END_REF] , but there are many others in the literature). The stochastic spectral theory is the more general case, where the soil is spatially correlated over finite distances, but restrictive assumptions are required (infinite domain, statistical homogeneity, small variability, etc., as already evoked earlier). In a different approach, [START_REF] Indelman | Upscaling of conductivity of heterogeneous formations: general approach and application to isotropic media[END_REF] consider the spatial moments of block conductivities in a numerical flow model, where they relate theoretically the moments of the block-conductivities to the moments of the underlying point-wise, continuous random field conductivities. This leads them to average block conductivities which depend on the statistics of the underlying K(x,y,z) but also on the size and shapes of the blocks. Some other methods are more empirical, such as the spatial power averaging approach by [START_REF] Desbarats | Spatial averaging of hydraulic conductivity in three-dimensional heterogeneous porous media[END_REF] . For a more extensive survey, see also the upscaling review by [START_REF] Dagan | Upscaling of flow in heterogeneous porous formations: critical examination and issues of principle[END_REF] To sum up, it should be recognized that each of the existing methods is based on some assumptions and has its own limitations. Furthermore, the nonlinearity of unsaturated flow complicates upscaling methodology compared to the case of linear saturated flow. The heuristic PAM, above, depends explicitly on soil variability, but it depends also on statistical spatial structure with its semi-empirically fitted power exponents. This makes it a flexible, appealing model, easily applicable to site-specific problems. In spite of its simplicity, it yields, also, better macroscale conductivity curves than classical averages like arithmetic or harmonic mean conductivity curves.

Flow and transport simulation models

Flow model

In terms of pressure head "h ", the mixed form of Richards' equation [START_REF] Bear | Macroscopic modeling of transport phenomenon in porous media. 2: applications to mass, momentum and energy transport[END_REF] , [START_REF] Ababou | Capillary Flows in Heterogeneous and Random Porous Media[END_REF] ) is used for modeling variably saturated flow through unsaturated and/or saturated zones. In the 2D code used here (UFT2D "Unsaturated Flow and Transport 2D ": see [START_REF] Soraganvi | Modeling of flow and advection dominant solute transport in variably saturated porous media[END_REF], the pressure-based flow equation takes the form, in vertical cross-section (x, z):

( ℎ ) ℎ + ( ℎ ) ℎ + ( ℎ ) = ( ℎ ) ℎ + (10)
The right-hand side contains a compressibility term which is neglected here (let Ss = 0), and the transient moisture term / t goes to zero for steady-state flow (this is the case analyzed here, although timestepping will be used for reaching steady state). In this flow Eq. ( 10) , the spatial variability of coefficients is not shown explicitly to alleviate notations. Coefficients K xx (h) and K zz (h) [m/s] are the pressure-dependent hydraulic conductivity components in the x and z direction respectively; z) [m] are the horizontal and vertical coordinates respectively, with z upwards.

h is pressure head [m]; is volumetric moisture content [m 3 /m 3 ]; is porosity [m 3 /m 3 ]; S S is the specific storage coefficient [m -1 ]; and (x,
Note also that Eq. ( 10) can be reformulated in terms of suction head = -h [m]. At the local scale, the conductivity curves K xx (h) and K zz (h) can be assumed equal (local isotropy). On the other hand, the local conductivity curve can also depend on spatial location (as is the case here). Thus, we let K xx = K zz = K (h; x, z) for the local conductivity, when solving for detailed flow in a randomly heterogeneous soil sample.

The flux boundary condition q z = q 0 used at the top soil surface ( z = L ) is expressed as follows using the unsaturated version of Darcy's pressure-based flux law, with z upwards:

= - ( ℎ, ) ℎ + 1 ⇒ 0 = - ( ℎ, ) ℎ | | | | = + 1 = , ∀ > 0
where q 0 is the water flux density [m/s] imposed at the top, which may be specified as constant (here) or as a function of time. More generally, in the UFT2D code, several other types of Boundary Conditions can be considered, such as: Dirichlet (imposed pressure or suction), Neumann (prescribed flux, and also, gravitational flux), and Seepage Face (which combines Dirichlet and Neumann). Finally, the initial condition can be specified as either distributed moisture content or a distributed pressure head throughout the domain.

Constitutive relations (Hydraulic conductivity and moisture vs. suction): The unsaturated constitutive relations are hypothesized as follows. The exponential expression of [START_REF] Gardner | Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[END_REF] is used for the conductivity-suction curve, expressed here for a heterogeneous soil:

( , ) = exp { -( ) } ... for ≥ 0 (11)
The empirical relation of [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] is used for the moisture retention curve ( ):

Θ = - - = 1 1 + | | (12)
where: = -h is suction head ( > 0 in the unsaturated zone); is volumetric water content; r is the residual and s the saturated moisture content; Θ( ) is effective saturation vs. suction; v is the Van-Genuchten / Mualem shape factor [1/m]; n V and m V = (1 -1/n V ) are the Van-Genuchten / Mualem exponents (dimensionless).

Transport model

The Advection Dispersion Equation (ADE) for the transport of a passive solute (tracer) in the unsaturated zone is given by [START_REF] Freeze | Groundwater[END_REF] :

+ ( ) + ( ) = + ( 13
)
where c is the concentration of the solute in the liquid phase, in Kg/m 3 of water. The ADE can also be written alternatively as follows:

+ ( ) + ( ) = ( ∕ ) + ( ∕ ) (14)
The concentration C = c , is the concentration in Kg/m 3 of space, and is the volumetric water content (possibly variable in both space and time). The coefficients D xx and D zz are the diffusion-dispersion coefficients [L 2 T -1 ] in the x and z direction respectively, and they may generally depend on moisture content and (in some models) on the unsaturated flux density (q). Note that the transport Eqs. ( 13) and ( 14) can be considered at two different scales: either at the local scale, with spatially variable coefficients like moisture content (x) and Darcy velocity q(x); or at the macro-scale, where Darcy velocity and dispersion coefficients are interpreted as spatially constant macro-scale quantities (mean Darcy velocity vector, and macro-dispersion coefficients).

For the transient unsaturated transport experiments to be developed in the present work, we are using the local scale Eqs. ( 13), ( 14) with some simplifications, as explained below.

Neglecting molecular diffusion, the local dispersion coefficients D xx and D zz [m 2 /s] are assumed to be of the form:

= = (15)
where, v X and v Z [m/s] are the pore velocity components in the x and z directions respectively; "b " is a constant dimensionless exponent; and "a " is a constant coefficient similar to the dispersivity length scale (however, the units of "a " are meters only when the exponent b = 1).

In our implementation, D zz is calculated from Eq. ( 15) , using units of [cm] and [hours], with D zz = 0.018|v z | 1.15 cm 2 /h for "large " velocities v z z > 1.586 cm/h, and a constant dispersion D zz = 0.036 cm 2 /h for smaller velocities v z < 1.586 cm/h. The transverse dispersion coefficient D xx is neglected, as it was found to be small compared to D zz for both velocity ranges.

The unsaturated Darcy fluxes q x and q z in the x and z directions are given by Darcy's law (as before):

= -xx ( ) ℎ = -zz ( ) 1 + ℎ (16)
and the corresponding velocity components of the tracer in the x and z directions are given by:

v X = q X ∕θ v Z = q Z ∕θ (17)
The boundary condition employed at the soil surface is a fixed advectivedispersive solute flux:

-zz + = - ... at = , ∀ > 0 (18)
where the water flux q o is taken as positive and the minus sign (-q o C o ) forces the advective solute flux downwards The initial concentration of the solute in the soil at t = 0 is specified as a constant concentration:

( , , 0 ) = Concerning the numerical flow and transport model, the reader is referred to the code UFT2D implemented in [START_REF] Soraganvi | Modeling of flow and advection dominant solute transport in variably saturated porous media[END_REF] . The unsaturated flow model is based, numerically, on 2D implicit finite difference approximation of the mixed form Richards equation with modified Picard iterations for linearization, solved with the 2D SIP (Strongly Implicit Procedure) solver. The unsaturated transport model, on the other hand, is based on finite volume discretization of the advection-dispersion equation (ADE), with different handling of the advective and dispersive components, with a flux limiter in the advective part. This numerical scheme can handle transport problems with large Peclet numbers, i.e., high-velocity regions.

Methodology

Numerical flow experiments

Steady-state gravity drainage flow simulations are performed by a time marching approach, using the 2D numerical model presented in [START_REF] Soraganvi | Numerical and stochastic modelling of flow and transport in the vadose zone[END_REF][START_REF] Soraganvi | Numerical and stochastic modelling of flow and transport in the vadose zone[END_REF][START_REF] Soraganvi | Modeling of flow and advection dominant solute transport in variably saturated porous media[END_REF] . A uniform infiltration rate q is specified over the top boundary of the domain, and a unit hydraulic gradient is imposed at the bottom boundary, while lateral boundaries are impervious. The initial condition is a uniform pressure head.

The numerical flow simulations are performed for a series of infiltration cases, that is, for different values of infiltration flux "q " imposed at the soil surface, using gravity drainage condition at the bottom. Each transient flow simulation is performed till the steady-state is reached, i.e. by time marching, using the transient mixed form of Richard's equation presented earlier. To accelerate the time marching procedure, the initial suction 0 used for each flow experiment was specified by inverting the relation K( 0 ) = q 0 , where q 0 is the imposed flux and K( ) is the Gardner conductivity curve with mean parameters. Mass balance errors ( t ) were checked during each flow experiment.

To sum up, a series of gravitational drainage experiments were performed using the above-described method, on single realizations of randomly heterogeneous 2D and 1D soil samples. The stochastic flow fields thus obtained at steady state were then analyzed to obtain effective conductivity curves, as follows.

The spatially averaged suction, or "mean " suction, was calculated for each imposed flux value q 0 , which is also the mean vertical flux q Z through the sample. The mean hydraulic gradient through the soil sample is unity (purely gravitational), and the mean pressure gradient is zero, due to the imposed gravitational drainage condition at the bottom ( Eq. ( 3) ). Therefore, at steady state, the effective vertical conductivity K ZZ ( ⟨Ψ⟩) is equal to the imposed flux q 0 = ⟨q Z ⟩, for the given mean suction ⟨ ⟩ of the soil sample. That is how the "upscaled " conductivity curve is constructed point by point, from the different flow experiments.

Note that the numerical simulations are carried out based on the single realization approach, with a unique realization of the Gaussian random fields Ln{K S (x,z)} and { (x, z) or Ln (x, z)}. Single realization simulations are carried out for different water flux rates to obtain the effective unsaturated conductivity curve K e vs. mean suction (as explained just above ) . The macro-conductivity curves K e ( ⟨Ψ⟩) obtained numerically are then compared with the analytical results presented earlier, namely: the stochastic spectral perturbation results of ( Yeh et al., 1985a( Yeh et al., , 1985b ) ) given in the previous Eqs. ( 5) and ( 6), and the Power Average Model "PAM " of [START_REF] Ababou | Probabilistic Analysis of Unsaturated Conductivity: Anisotropy, Crossing Point, and Upper Envelope: Section 6.3, Chapter 6 (Stochastic Analysis of Flow and Transport)[END_REF] ) given in the previous Eqn 9

For the heterogeneous soil, the relationship between the hydraulic conductivity K ( x ) and suction ( x ), for each grid point x within the discretized flow domain, is represented by [START_REF] Gardner | Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[END_REF] model given by the previous Eqn 11 .

In this study, both parameters (Ks, ) are considered as random fields (i.e., random functions of space). On the other hand, the moisture retention curve ( ) of [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] is used as the Soil Moisture Characteristic (SMC), and it is assumed to be uniform throughout the domain. To sum up, the emphasis here is on the spatial variability of the conductive properties of the soil (rather than its water retention properties). Upscaled constitutive relationships (effective K( )) are then obtained by repeated numerical flow experiments, for a series of given infiltration rates, as explained earlier.

Unsaturated flow upscaling

The effective conductivity is calculated numerically by interpreting the steady-state flow simulations in terms of global quantities: mean flux and mean gradient. The resulting "upscaled " K is plotted versus mean suction ⟨ ⟩ (denoted here for simplicity), and it is then compared with two available analytical results: the stochastic spectral perturba-tion solution, and the power-averaging results (Power Averaging Model "PAM ").

In the latter case (PAM), the upscaled numerical conductivity curves are fitted to the power averaging formula in Eq. ( 9) to obtain the un-known power exponent p i using a nonlinear least-squares optimization method (Levenberg-Marquardt iterative algorithm, implemented with the M function ). The numerically upscaled K( ) curve is also compared to the Arithmetic, Harmonic and Geometric mean K( ) curves, as well as to the more general Power Average Method (PAM).

The " " parameter plays an important role in the resulting value of the averaging powers p i . It is thought to be related to pore size distribution, and its inverse represents a capillary length scale: 1/ = CAP . The length scale CAP is comparable to layer thickness, or to spatial correlation length Z . The cases Z ~ 1, Z < 1, Z > 1 are explored.

Transport simulations

The two-dimensional steady velocity field, obtained by flow simulations, is used as input to transport simulations. The simulated contaminant plume is being advected and dispersed in the steady flow field. Its space-time distribution is then analyzed in terms of macro-scale moments. For unsaturated transport simulations on the randomly heterogeneous soil samples, the local ADE (Advection Dispersion Equation) given by Eq. ( 14) is used. This ADE has two randomly heterogeneous coefficients, obtained from the steady flow simulations, namely: (i) the 2D moisture content field θ( x ⃗) , and (ii) the 2D Darcy velocity vector q ⃗( x ⃗) .

The transient transport simulations are carried out on the 2D unsaturated velocity fields that were generated by the steady flow simulations. The transient concentration field C(x,z,t) is thus obtained numerically, and it is then analyzed and "upscaled " using spatial averages and spatial moments. An equivalent homogenized ADE (or macroscale ADE) governing the mean concentration ⟨C⟩ is hypothesized to be of the form:

∇ • ∇ -( ) • ∇ = ∕ (19)
where ⟨C ⟩ is mean concentration; ⟨V ⟩ is mean velocity [m/s]; and D e stands for the "effective " or macroscale dispersion coefficient [m 2 /s], which represents the mean tensorial hydrodynamic dispersion caused by spatial variations of the velocity field. Note that the mean pore velocity is vertical,

V = [ ⟨V X ⟩, ⟨V Z ⟩] = [0, ⟨V Z ⟩],
and we are interested here essentially in vertical longitudinal macro-dispersion, . This longitudinal macro-dispersion coefficient (renamed simply "D ") is assumed proportional to the mean absolute vertical velocity |V|, that is D = A |V|. The proportionality constant "A " is the so-called macro-dispersivity length scale.

In this solute transport study, we aim at (i) analyzing the migration of the unsaturated solute plume in the randomly heterogeneous soil; and (ii) upscaling the dispersion process (macro-dispersion "D ", macrodispersivity "A ").

Moments of concentration plume, and macro-dispersivity analysis

The mean concentration and variance of the solute plume for a slug tracer are analyzed at different times after the release of the slug at the top of the soil (vertical cross-section). Average solute concentration profiles C(z,t) are obtained as a function of depth, by averaging horizontally along the x-direction. The average concentrations C(x,z,t) are then used to analyze the spatial moments of the concentration plume C(x,z,t). The plume is analyzed only with respect to the Z-direction, as it is the dominating flow direction. The average concentration C(z,t) is the average concentration of the row given by:

( , ) = max ∑ =1 ( , ) ( )Δ Δ
The spatial moments are then defined as follows:

( ) = max ∑ =1 ( , ) ( )Δ Δ (20) M Z ( t ) = 1 max ∑ =1 ( , ) ( )Δ Δ ( ) (21)
( ) = 1 max ∑ =1 ( , ) ( )Δ Δ ( ) 2 - ( ) 2 (22)
where M o ( t ) is the total mass of the solute residing in the flow domain, which consists of imax rows and jmax columns, with length and depth as Δx and Δz respectively.

Recall that "c " (lower case) is the mass concentration per volume of solvent (water), while C ( x i ,t ) = c ( x i ,t ) ( x i ) is the mass concentration per volume of space (i.e., per volume of porous medium). The zeroth moment M 0 (t) is the total mass of the solute in the domain. The first moment, M Z (t), is the vertical position of the center of mass of the plume with respect to the reference frame. The second moment M ZZ (t) can be viewed as the vertical moment of inertia of the plume: it measures its mean square spread around its center of mass in the vertical direction. The time evolution of these moments, M 0 (t), M Z (t), M ZZ (t), is then studied for each transient simulation of unsaturated solute transport. The different cases involve 2D samples of statistically isotropic and anisotropic, randomly heterogeneous soils.

In practice, the moments are calculated in that order: M 0 (t), M Z (t), M ZZ (t). Since the zero-order moment is the total mass of resident solute in the domain, it can be used as an indicator of mass balance error by comparing it to the net solute flux through all boundaries, cumulated over time. The first moment M Z (t), normalized by mass M 0 (t), provides the mean plume displacement Z(t):

Z ( t ) = M ZZ ( ) ∕ M o ( ) (23)
The time derivative of Z(t) indicates, therefore, the mean effective Lagrangian velocity Uz(t) of the center of mass of the plume:

U ( t ) = ( t ) = ∕ dt (24)
There is no reason that the calculated velocity Uz = dZ/dt be constant in time from the numerical simulations. However, a constant average velocity U 0 = ⟨U z (t) ⟩ can be obtained by averaging the calculated Lagrangian velocity U z (t) over time. Another constant average velocity V o = ⟨V z (z) ⟩ can also be obtained by averaging the Eulerian velocity V z (z) in space (vertically).

Finally, the second moment M ZZ (t) represents the spatial spreading of the plume. After normalization by the resident mass M 0 (t), it can be considered as the dispersion variance of the plume in the "z " direction ( zz 2 (t)):

σ 2 ZZ ( t ) = M ZZ ( t ) ∕ M o ( t ) (25)
This leads to a natural evaluation of the macro-dispersion coefficient "D " [m 2 /s] by analogy with Fickian diffusion, as follows:

D ( t ) = D zz ( t ) = 1 2 dσ 2 ZZ dt (26)
The calculated macro-dispersion D(t) in (Eq. ( 26) may not be constant in time (unlike the ideal case of Fickian dispersion). However, a constant average value of the macro-dispersion coefficient, D o = ⟨D(t) ⟩, can be obtained by taking the time average of the macro-dispersion coefficient D(t) defined just above. Finally, the vertical macro-dispersivity coefficient "A " [m] is defined from the above-defined macro-dispersion D, as follows:

A ( t ) = A zz ( t ) = D ( t ) U ( t ) = 1 2 dσ 2 ZZ dt dZ ∕ dt (27)
where U(t) is again the Lagrangian vertical velocity of the center of mass of the plume. Again, the calculated A(t), above, depends generally on time if macro-dispersion is non-Fickian. However, a constant mean macro-dispersivity "A 0 " may be defined using the constant mean diffusion coefficient Do (above), and the constant mean Eulerian velocity V 0 (vertically averaged flow field velocity). This yields finally the constant macro-dispersivity length scale A 0 [m], defined by:

A 0 = A zz 0 = D 0 U 0 = ∫ Tmax Tmin D ( t ) dt ∫ Tmax Tmin U ( t ) dt (28)
where:

D(t) = 1 2 dσ 2
ZZ dt ; and U(t) = dZ ∕ dt . Note: the spatial moments obtained numerically at early times and at large times are ignored in order to avoid the possible boundary effects at the top and bottom of the soil sample.

Results and discussion

Statistical representation of the soil parameters

In this subsection, we present the parameters of the soils in terms of their spatial and statistical distribution. These parameters are to be used for the flow simulations in this section, and also for the transport simulations in the next section.

In this work, we perform flow and transport experiments on single realizations of randomly heterogeneous soil samples. Therefore, we generate single realizations of the spatially distributed parameters (Ln-Ks and ) of the permeability-suction curve. The means and standard deviations of the parameters are taken from [START_REF] Khaleel | Upscaled flow and transport properties for heterogeneous unsaturated media[END_REF] , corresponding to coarse-textured sand from the Upper Hanford formation at the Hanford site, Washington, USA ( Tables 1 and2 ). The Ks parameter is taken lognormal in all the cases; the parameter is taken Gaussian or, also, lognormal in a few cases (the difference between these two choices becomes significant only for high variability of ). The two random parameters, {Ln(Ks), } or {Ln(Ks), Ln( )}, are generated independently. The 2-point spatial correlation structures and correlation lengths of the two parameters are the same. The moisture curve ( ϴ-) is assumed to be spatially uniform, as per the constant parameters shown in Table 2 .

Three types of statistical spatial structures will be generated for the pair of random fields {LnKs, }, and one type of spatial structure for the pair {LnKs, Ln }. The four types of randomly heterogeneous samples (single realizations) are as follows:

• First type: perfectly layered "1D " random field, with X >> Z , or equivalently, X → ∞ and Z finite, in a domain of size 10 m discretized into 1001 nodes. (Ks is lognormal, is Gaussian).

• Second type: imperfectly stratified "2D" anisotropic field with X > Z , in a domain of size 25 m × 50 m, discretized into 501 × 501 grid points. (Ks is lognormal, is Gaussian).

• Third type: "1D and 2D" anisotropic fields with perfect and imperfect random stratification. (Both Ks and are taken to be lognormal random fields).

• Fourth type: 2D isotropic field ( X = Z ), generated in a domain of size 10 m × 10 m, discretized into 1001 × 1001 grid points. (Ks is lognormal, is Gaussian).

In each case, the 2-point covariance function of the single realization soil sample was calculated numerically, then plotted and compared with the theoretically prescribed exponential function. The results were in agreement for a reasonable range of separation distances.

Once the unsaturated flow is computed and analyzed on each single realization sample, transient tracer transport simulations are performed and analyzed in terms of concentrations and spatial moments. These solute plume transport simulations are implemented for 2D anisotropic and 2D isotropic soil samples.

Flow simulation results

Flow through perfectly stratified "1D " soils (

= ∞) One-dimensional (1D) parameter fields, LnK (z) and Gardner's parameter (z), are generated in space for the statistical parameters shown in Table 1 , over a grid of 1001 nodes. Though the two parameters are assumed to be uncorrelated to each other, each parameter is spatially correlated. The two-point correlation length is Z = 0.10 m, which gives a spatial resolution Z / ΔZ of 1/10 on the numerical grid. A replicate of the generated random field LnK (z) is shown in Fig. 1 The computed (empirical) auto-covariances of the generated fields match quite well the theoretical covariances plotted versus the separation distance between the two-points ( "also known as lag "): see Fig. 2 for that of the saturated log conductivity random field (LnKs); the covariance structure is similar for the random field (not shown). In Fig. 2 , the numerically estimated covariance function, calculated for a single replicate of LnKs, is compared to the theoretically prescribed covariance functions (exponentially decreasing with separation distance).

Two principal flow directions will be considered below: flow perpendicular to bedding and flow parallel to bedding. In the first case (perpendicular) the strata are horizontal and flow is vertical, aligned with gravity. In the second case, flow parallel to bedding is obtained by rotating the stratified soil sample by 90°while keeping gravity vertical. In this way, the flow direction remains aligned with the direction of gravity in both the cases. The resulting set of flow fields is used to upscale the hydraulic conductivity against the mean suction. Anisotropy effects are studied and compared with theoretical results on suction-dependent effective conductivity.

Flow experiments perpendicular to stratification : For the generated 1D soil samples of 1001 nodes (each with mesh size Δz = 0.01 m), a set of numerical simulations are carried out with a constant flux at the top of the domain, time marching till the flow reaches a steady state. The steady pressure fields obtained for a series of experiments with different flux rates are shown in Fig. 3 . During our simulations, we come across some cases where internal ponding has occurred locally in some parts of the sample (this is for higher flux rates, not shown here). Since stochastic theory cannot be applied stricto sensuto to partially saturated flow, the high flux cases (q ≥ 0.5 m/d), are not considered for upscaling in this paper. It has been observed that in such cases, it took a long time to reach steady-state, compared to the cases with lower fluxes where saturation was never reached.

Plots of pressure head profiles (not shown here) indicate the stabilization of the pressure head as the constant infiltration rate continues for some time (an example of transient pressure profile evolution is shown in Fig. 13 further below). It has been observed that for some values of infiltration flux, the steady-state is reached very fast, and for some other values, it takes a longer simulation time. The pressure head variation range is broader for lower infiltration fluxes and is narrower for higher fluxes.

Based on the steady flow simulations, the upscaled conductivity curve is plotted point by point vs. mean suction , in terms of K( )/K g and Ln{K( )/Kg} as well. The point-by-point curves are then compared graphically with the spectral perturbation results of Yeh et al., 1985a, Yeh et al., 1985b & Yeh et al., 1985 c) and with the Power Average Model of Ababou (1991a[START_REF] Ababou | Probabilistic Analysis of Unsaturated Conductivity: Anisotropy, Crossing Point, and Upper Envelope: Section 6.3, Chapter 6 (Stochastic Analysis of Flow and Transport)[END_REF].

The comparisons are shown for flow perpendicular to perfect stratification (1D) in Fig. 4 a and b, exhibiting a good match between theoretical and numerical upscaling. The best-fitted curve for the Power Average Model is obtained for with the power p i = -0.0484 for the Ln(K/Kg) fit, and p i = -0.0634 for the (K/Kg) fit. It is seen from this figure that the numerical K( ) curve matches with the harmonic mean at low suctions (wet range), as could be anticipated, but remarkably, the curve moves away from the harmonic mean and closer to geometric mean at higher suctions (dry range).

Note: choosing a linear fit of Ln(K( )/Kg) favors low suctions (wet range of moistures), while choosing instead a power fit of K( )/Kg favors the higher suctions (dryer range of moistures); this will appear more clearly in the next case below (parallel flow).

Flow experiments parallel to stratification : Recall that the case of flow parallel to bedding ( "parallel flow ") is obtained by rotating the strata from horizontal to vertical direction. In this way, "parallel flow " is vertical and aligned with gravity. As anticipated, the parallel flow reaches the steady state very fast, and the resulting pressure and moisture fields are not variable along the flow direction.

In comparison, when the flow is perpendicular to perfect horizontal stratification, the flow field is variable across the direction of flow, in that case, it takes more time to reach the steady-state. The pressure head contours for "parallel flow " will be a set of vertical straight lines, as shown in Fig. 5 . Comparing Fig. 5 (parallel flow) to Fig. 3 (for perpendicular flow), we can clearly see the different structures of the pressure fields: there is a drastic variation of pressure vertically in the case of flow perpendicular to bedding ( Fig. 3 ). The flow parallel to bedding ( Fig. 5 ) reaches steady state faster, and also the pressure head does not vary across the depth since the field is homogeneous in the direction of flow.

Again, the series of simulations performed for "parallel flow " is used to obtain the upscaled K( ) relation point-by-point. The numerical points K( ) are compared with the spectral perturbation theory of Yeh et al., 1985a, Yeh et al., 1985b, Yeh et al., 1985 ), and with the best fitted Power Averaging Model (PAM). These comparisons are shown in Fig. 6 (a,b), along with three special curves, the arithmetic, harmonic and geometric mean curves (these three curves are also computed with the PAM formula, by inserting exponents + 1, 0, -1 respectively).

The best-fitted power average parameter is p i = -0.1749 for the Ln(K/Kg) fit, and p i = + 0.9640 for the (K/Kg) fit. Indeed, it can be seen that the numerical points K( ) coincide with the arithmetic mean curves for low suctions (wet range), but deviate somewhat from it (in favor of the geometric mean) for higher suctions (dry range). This observation confirms that, when estimating the effective K curve, choosing a linear fit of Ln(K( )/Kg) implies favoring low suctions (wet range of moistures) while choosing instead a power fit of K( )/Kg implies favoring more the higher suctions (dryer range of moistures).

Flow through imperfectly stratified field (

>> ; = 1.6266) Here, the case of flow in an imperfectly stratified 2D medium is explored. The hydraulic properties of the soil remain the same as in Tables 1 and2 , with the horizontal correlation length chosen much greater than the vertical one ( x >> z ). This yields horizontally elongated strips of soil, or imperfect layers (hence generalizing the 1D case of perfect layers). The generated random field of saturated hydraulic conductivity is shown in Fig. 7 , the Gardner's scaling parameter is also generated and looks the same; thus not shown here. Numerical flow experiments are performed in two directions, perpendicular and parallel to stratification, by rotating the field by 90 0 , keeping the main flow direction vertical, as explained earlier.

Flow Perpendicular to imperfect stratification : Fig. 8 a andb show the steady-state suction fields for flow perpendicular and parallel to stratification, respectively (flow is vertical and aligned with gravity in both the cases). The effect of stratification is clearly seen in the distribution of suction, as it follows the pattern of stratification. In both the cases, the mean flow direction is vertical, and the streamlines (not shown here) are vertical on average. Fig. 9 shows the comparison of numerical results obtained for flow perpendicular to imperfect stratification with the Power Average Model and with the spectral perturbation theory: both match quite well the point-by-point numerical conductivity curve (numerically upscaled). The figure plots, in particular, the upscaled conductivity curve obtained from PAM with best-fitted power averaging parameter "p i "yield the result as p i = -0.0015 for the Ln(K( )/Kg) linear fit, and p i = -0.864 for the K( )/Kg power fit.

As anticipated, the effective unsaturated conductivity, in this case, does not really correspond to a harmonic mean conductivity curve. Some more tests would be needed to explore the mean behavior of the imperfectly stratified soil near saturation (near zero suctions). Note, however, that the mean suction, mean moisture, and effective conductivity concepts, all become less relevant when the flow regime becomes partially saturated. Our simulations demonstrate that this partial saturation occurs, indeed, when flux rates are increased to higher ranges, q ≈ 1 m/d and above.

Flow parallel to imperfect stratification : Fig. 10 shows the comparison of numerical results with the PAM and with Yeh's spectral perturbation results, which both match quite well the numerically upscaled point-by-point curves. The figures plot the best fit of the upscaled power average conductivity curve: the resulting best fitted power is p i = + +0.0636 for the linear fit of Ln{K( )/Kg}, and p i = -0.0777 for the power fit of K( )/Kg. These results, for this case of "parallel " flow, indicate that the effective conductivity is not the arithmetic mean conductivity curve: it is in fact closer to a geometric mean conductivity curve. This is probably due to the imperfect stratification (as opposed to perfect layers) but also, to the nonlinear nature of unsaturated flow.

More suction points could be tested in the wet range of suctions; however, see the previous remarks on the difficulty of interpreting the mean behavior in the wet range (especially if the soil becomes partially saturated). The dependence of suction variance on the mean suction is shown in Fig. 11 , based on the flow simulation results. This is analyzed just below in relation with theoretical results.

Variance of the suction field (x,z) :

The dependence of suction variance on mean suction is shown in Fig. 11 a and b, based on the flow simulation results, for flow perpendicular and parallel to stratification, respectively. It can be seen from these plots that suction variance 2 is a non-monotonic function of mean suction ⟨ ⟩ (also denoted " " for convenience). Suction variance first decreases towards a minimum around ≈ 0.60 m, and then increases again for larger suctions. This matches with the results of Ababou et al. (1988 ) and [START_REF] Yeh | One dimensional steady state infiltration in heterogeneous soils[END_REF], who observed that the suction variance can approach zero and then increase. Also, according to Yeh et al., 1985a, Yeh et al., 1985 b), [START_REF] Harter | Stochastic analysis of solute transport in heterogeneous, variably saturated soils[END_REF], and many other investigators, suction variance increases with suction for sufficiently high suctions (sufficiently away from the state of total saturation) as occurs here. This behavior can be compared qualitatively with that observed the-oretically for steady flow by Yeh et al., 1985a, Yeh et al., 1985 b), also summarized in Gelhar (1993) (Chap. 4.5: Unsaturated flow in heteroge-neous soils: Steady flow results, pp. 166-175). Thus, for 3D random soil structure, with zero cross-correlation between LnKs and (as we have here), the suction variance is theoretically of the form (given below us-ing our notations): Gelhar's 1993 Eq. [4.5.21]) where g( •, •) is a function plotted numerically in Gelhar 1993 (Fig. 4.38, p. 168), decreasing with (

2 = ( 2 + 2 2 ) 2 ( , ∕ ) (cf.
) and increasing with ( / ). In comparison, for the imperfectly stratified case at hand, and given the statistics of our parameters (shown earlier in Table 1 ), we have here:

LnKs = 0 . 771 ; = 1 . 493 -1 ; ⟨ ⟩ = = 8 . 133 -1 ; , = ( 2 . 0 , 0 . 20 )
Taking now the case of vertical flow orthogonal to bedding, we set = 1.6266 and / = 10 in Gelhar's Eq.[4.5.21] above. For these values, Gelhar's function "g " is roughly equal to 0.5, and we obtain theoretically for our case:

2 ≈ 0 . 02 × 0 . 5944 + 2 . 2290 2 ........ 2 (29)
This theoretical result would yield a suction variance 2 ≈ 0.0279 and standard deviation ≈ 0.167 m for ≈ 0.60 m, 2 ≈ 0.0337 and ≈ 0.184 m for ≈ 0.70 m, and 2 ≈ 0.0761 and ≈ 0.276 m for ≈ 1.20 m. At first sight, these values do not match our experimental values of suction variance in Fig. 11 ; however, the fit of 2 ( ) may be off by only a multiplicative constant (the function "g "), perhaps due to the fact that Gelhar and Yeh's results are for an infinite domain, while we are using here finite soil samples (single realizations).

Randomly heterogeneous samples with both Ks and as lognormal random fields

In all the previous cases, the flow simulations involved a LogNormal random field for Ks and a Gaussian random field for the parameter (Gardner's parameter). In contrast, in this section, the parameter is considered as LogNormal rather than Gaussian, which is, in fact, more realistic. The relations between Gaussian and Log-Normal random fields at single points in space will be addressed, briefly, further below. Fig. 12 shows the generated random field parameters of heterogeneous samples on which unsaturated flow experiments are conducted.

This new set of experiments with LogNormal is justified by the fact that the parameters like Ks and must be positive, therefore they should follow a probability law that constrains them to remain positive. This is the case of the LogNormal law. Thus we take here both Ks and to be Log-Normal; or equivalently, we take both LnKs and Ln to be Gaussian random fields.

The relations between Gaussian and Log-Normal random fields at single points in space are now briefly summarized. For a positive random parameter like , it can be shown that, for moderate variability, the Gaussian and LogNormal distributions are approximately equivalent. Also, in that case, the degree of variability of the parameter can be measured equivalently from the Coeffificent of Variation C or from the standard deviation Ln . This is because of the approximate equivalence formula C = /m ≈ Ln . Other relations can be developed to relate explicitly the 2-point spatial covariance of the random field (x) to that of its logarithm Ln (x).

Fig. 13 shows a set of multiple transient flow experiments (converging to steady state) for the randomly stratified sample, with vertical flow orthogonal to strata induced by unsaturated gravitational boundary conditions. The figure shows intermediate transient pressure profiles h(z,t) during the sequential multi-simulation procedure. Unsaturated zones have negative pressure (h(z,t) < 0). Internal ponding phenomena with h(z,t) > 0 can be seen to occur for this particular experiment. The final steady-state used in this study for K( Ψ) upscaling still has some ponded zones, but much less so than during the transient phase. Interestingly, looking at the steady state pressure profile h(z) in . Figure 13 , we can see that the soil sample has less saturated zones ( h > 0) than unsaturated zones ( h < 0). However, the pressure h(z) is higher in positive zones compared to negative zones (in absolute value). As a result, the mean pressure ⟨h ⟩ is positive (straight line), even though there are more unsaturated zones than saturated zones.

Simulations with such saturated pockets, called internal ponding, were presented in [START_REF] Soraganvi | Analysis and upscaling of unsaturated flow through randomly heterogeneous soil[END_REF] . The phenomenon was analyzed in the form of a ponding ratio, which increases gradually as the unsaturated field evolves towards saturation. We have also developed upscaled conductivity curves for such cases (which are not presented here in this work).

As will be seen below, it is found that the perfectly stratified unsaturated soils can behave like saturated media, with arithmetically averaged K( ) for parallel flow, and harmonically averaged K( ) for perpendicular flow. However in some cases, depending on the flow regime and on layer thickness, this classical behavior does not hold good.

The numerically upscaled K ii ( ) curves : As mentioned earlier, the numerically upscaled K ii ( Ψ) curves are constructed point-by-point (K n , Ψ n ). In some cases, under relatively "wet" flow regimes, internal saturation and ponding could occur and remain stable at steady state, which is an interesting phenomenon and not covered by the various upscaling theories known to us.We choose to show here an example of a numerically constructed curve where the "saturation/ponding" phenomenon did occur at low values of mean suction shown in Fig. 14 .

It is interesting to see how such curves compare to the PAM theory, and to find out the best-fitted value of the power average exponent (p i ) in different cases. For that purpose, we choose to focus on the optimal fit of the Power Average exponent in terms of ln(K/Kg) rather than (K/Kg).

Firstly, we choose to fit the single parameter "p " based on Eqn 30 , which we prefer for theoretical reasons This choice is made in spite of the fact that the best fit p i was usually better with (A, B, C) than with the single (p), as it should.

Secondly, we have found differences in the p-fitted in terms of K( Ψ) vs. ln K( Ψ). The ln(K/Kg) fit was finally chosen because it provides a better-fitted curve over a broader range of suctions. This choice is also based on the fact that, considering all the test cases, the behavior of the fitted exponent (p) obtained with ln K( Ψ) seemed more consistent than that obtained with K( Ψ).

Optimization results: the fitted power average exponent p : Table . 3 summarizes the best-fitted values of the PAM exponent (p) obtained for the vertical flow experiments in the randomly stratified samples with mean flow orthogonal to strata. Other results are not shown and not discussed here. The different columns of Table 3 correspond to different values of the mean scaling factor and, equivalently, of the capillary ratio . The saturated conductivity Ks(z) is a spatially correlated lognormal random field. The scaling factor (z) is perfectly cross-correlated to Ks(z) in the first column, or constant in the other columns.

Behavior of the power exponent vs. capillary / geometric scales : As discussed earlier the numerical experiments on heterogeneous samples were used to obtain upscaled unsaturated K ii ( Ψ) curves described point by point. Parameter LnKs(x,z) is randomly stratified with vertical correlation scale z = 0.10 m and standard deviation (LnKs) = 0.771. Parameter is constant here, and is equal to 8.133 m -1 (this yields a capillary length scale 1/ = 0.13 m). The Power Average exponent "p i " was then identified by a best fit procedure. The results are now used to confirm a previously formulated conjecture concerning the behavior of the Power Average exponent "p i " vs. capillary effects.

The numerically upscaled curve is presented in Fig. 15 in dimensionless variables: dimensionless conductivity K/Kg where Kg ≡ exp( ⟨LnKs ⟩> ), and dimensionless suction Ψ ≡ . ⟨ ⟩ , or Ψ ≡ ⟨ ⟩ ⟨ ⟩ more generally if is variable. Each point in Fig. 15 represents the steady-state solution of a flow simulation experiment. The rectangle in the plot indicates the results obtained in the "wet" flow regime with internal ponding in some regions, resulting in a negative mean suction (positive mean pressure) even though the sample is not saturated everywhere.

At this point, we seek essentially a formulation of the (fitted) exponent "p i " in relation to the capillary/geometry scale ratio which characterizes the heterogeneous unsaturated medium. Accordingly, the fitted exponent "p i " was studied as a function of a dimensionless capillary ratio defined earlier (we write here " ii " because of the dependence on flow direction "i"):

= HE ≈ HE ∕ CAP (30)
This capillary ratio, ii , can be interpreted as the ratio of a geometric scale of heterogeneity (characteristic layer thickness) versus a capillary length scale CAP = 1/ ⟨ ⟩ . The capillary ratio is directional because the geometric heterogeneity length scale is directional ( HET . characterizes heterogeneity along the direction of mean flow). Following these ideas, the results of Table 3 are plotted in Fig. 15 , as p OPT = f ( ii ). Furthermore, the numerical relation p = f ( ) is rather well fitted by the following analytical function (subscript "i " indicates principal direction, e.g. horizontal for i = 1, vertical for i = 3):

= ≈ 2 1∕3 1 + 1∕3 -1 (31)
Fig. 13. Transient numerical simulations: evolution of vertical pressure profiles h(z,t), along the mean flow direction with flow orthogonal to layers. Each graph shows the profiles h(z) between an "initial" time t 0 and a "max" time t MAX for one "iteration" of time marching. Each point in Fig. 15 represents a numerical flow experiment on the randomly stratified sample with flow direction orthogonal to the strata. The curve shows results for as constant in space. The isolated point "o" corresponds to as random and perfectly cross-correlated to Ks. The ordinate value p = 0 corresponds to the geometric mean K( Ψ) curve, the lower "bound" p = -1 to the harmonic mean K( Ψ) curve, and the upper "bound" p = + +1 to the arithmetic mean K( Ψ) curve.

It is important to note that the behavior of the upscaled unsaturated permeability curve depends not only on flow direction and layer thickness but also on CAP = 1/ . For ≈ 1 ( HET ≈ CAP = 1/ ) we obtain an exponent p ≈ 0, and therefore, the upscaled permeability curve is the geometric mean curve, even though the medium is stratified. The harmonic mean curve appears only as an asymptotic case for ii → 0, and our smallest negative exponent was, in fact, p ≈ -0.64 (between geometric and harmonic mean). And finally, the arithmetic mean curve appears as an asymptotic case for ii → ∞ (our largest positive exponent was, in fact, p ≈ + +0.34).

These results on the effective K( Ψ) curve can be explained as follows: (i) if ii ≪ 1, then HET ≪ CAP , and any capillary region of diameter CAP will contain many layers, resulting into harmonic mean behavior; (ii) if

ii ≫ 1, then HET ≫ CAP , and any capillary region of diameter CAP Fig. 15. Optimal power averaging exponent "p", fitted in terms of lnK, and plotted as a function of the capillary ratio . Each point represents a numerical flow experiment on the randomly stratified sample (here, flow is orthogonal to strata). The green curve shows results for constant in space. The isolated point "o" corresponds to random and perfectly cross-correlated to Ks. The ordinate value p = 0 corresponds to the geometric mean K( Ψ) curve; the lower "bound" p = -1 corresponds to the harmonic mean K( Ψ) curve, and the upper "bound" p = + 1 corresponds to the arithmetic mean K( Ψ) curve. Fig. 16. Generated statistically isotropic 2D random field of Ln(Ks(x,z)); Gardner's parameter (x,z) [m -1 ] also follows the same trend and hence not shown here. Note: these two random fields are statistically independent here (zero cross-correlation).

will contain essentially a homogeneous sub-layer, hence the arithmetic mean behavior.

Flow through statistically isotropic 2D media (

= ) In this case, the numerical simulations are performed on a 2D sample of statistically isotropic random soil, with isotropic auto-covariance function. The soil sample is of size 10 m × 10 m, with a numerical grid of size 1000 × 1000 nodes. The statistical parameters given in Table 1 are used to generate 2D isotropic random fields for two parameters: saturated hydraulic conductivity Ks(x,z) [LogNormal] and Gardner's parameter (x,z) [Gaussian].

Fig. 16 shows the random field of LnKs and (x,z) also looks the same and therefore not shown here. Auto covariance functions of the generated random fields are found to be isotropic and exponential, coinciding with the specified auto-covariances. A correlation length of 1.0 m is used in both the directions, with 100 nodes per correlation length along each direction, offering a good statistical resolution in discrete space (see statistical criteria for numerical simulations in [START_REF] Ababou | Numerical simulation of three dimensional saturated flow in randomly heterogeneous porous media[END_REF]). Gravity drainage numerical experiments are carried out for various fluxes, till the flow reaches a steady-state for each prescribed flux. Due to statistical isotropy (i.e., the geometrical isotropy of the random soil structure), the mean suction averaged over depth is nearly equal to the mean suction obtained by averaging over the whole 2D domain (the discrepancy is only 10 -2 m).

Accordingly, a series of numerical infiltration simulations are performed, for several infiltration rates (5E-1, 1E-1, 5E-2, 1E-3, 5E-4, 11E-4, 5E-5, and 1E-5 m/d), and the corresponding steady-state pressure field is obtained in each case. Fig. 18 shows the best-fitted curve of upscaled conductivity from power averaging theory (PAM), with a fitted power "p i " equal to -0.0484 for the best Ln(K/Kg) fit, and -0.1517 for the best K/Kg fit. Thus, as anticipated, the numerical results, in this case, are quite close to the geometric mean curve (p i = 0) although the fitted exponent is slightly negative (p i ≈ -0.05). The dimensionless capillary length scale is z = 8.133 >> 1, which should indeed yield a result close to the geometric mean curve (p i = 0).

In summary, the effective conductivity curve is close to the geometric mean curve, although somewhat smaller, leaning slightly towards the harmonic mean curve. Though the numerical effective conductivity values in Fig. 18 (a,b) match quite well to those obtained analytically, it is seen that the theoretical results diverge from the numerical conductivity values for higher suctions. This occurred also in other cases, 1D as well as 2D. This limitation has been rightly pointed out by many researchers ( Ababou (1991a , Mantoglou andGelhar, 1987 ) and [START_REF] Khaleel | Upscaled flow and transport properties for heterogeneous unsaturated media[END_REF] ), namely, that analytical upscaling can only be applied to relatively low suction values.

On the other hand, the theoretical results are not expected either to hold good near saturation (near zero suctions), the reason being that the theories use the exponential Gardner conductivity curve, which is easy to handle due to its log-linearity, but which cannot be valid near saturation (not valid below some low suction value like the bubbling or air entry suction). Another shortcoming is the limit on soil variability: the analytical and spectral perturbation results are limited to low variability, and may not hold for standard deviations larger than unity in terms of the log-parameters of the soil. In addition, though in numerical upscaling we can increase the variance a bit, a very high variance has been found to cause convergence problems and large computational time.

Overall, concerning the set of results obtained in this section for the upscaled conductivity curves in different cases, it is interesting to note that in all three cases considered, the effective conductivity curves ob-tained by numerical simulations match relatively well with the analyti-cal results given by Eqs.5 and 6 given earlier from spectral perturbation theory ( Yeh et al., 1985a, Yeh et al., 1985b Gelhar 1993) , and matched even better the best fitted Power Average Model of Eqs. (8a), (b), (c) ( Ababou (1991 a), [START_REF] Ababou | Probabilistic Analysis of Unsaturated Conductivity: Anisotropy, Crossing Point, and Upper Envelope: Section 6.3, Chapter 6 (Stochastic Analysis of Flow and Transport)[END_REF], [START_REF] Bagtzoglou | Effective hydraulic property calculations for unsaturated, fractured rock with semi-analytical and direct numerical techniques: review and applications[END_REF]) ; and also better than Classical means.

Another remark is that the variance of tension (suction) for vertical flow parallel to bedding is smaller than the variance for vertical flow perpendicular to bedding. This suggests more tortuous flow paths (comparatively) in the case of vertical flow perpendicular to horizontal stratification.

Transport simulation results

Numerical experiments of tracer transport

In this section, the 2D steady flows are used as input "parameters " for the ADE transport equation. The "flow parameters " are the Darcy velocity vector field q(x,z) and the moisture content field (x,z). The tracer velocity vector v(x,z), or pore velocity, is obtained from: v(x,z) = q (x,z)/ (x,z). The solute is released as a slug input from the top of the domain during a very short time. For all the cases considered, initial conditions and input conditions of the solute are kept the same (as described earlier).

At each sampling time, the 2D concentration plume is averaged over the horizontal width of the flow domain to obtain vertical concentration profiles as a function of depth. Recall that the flow was always vertical, with either horizontal bedding or vertical bedding; the transverse averaging of concentration is performed in both cases, resulting in vertical profiles C(z,t) in both cases.

Analyses of space-time concentration fields C(x,z,t)

The general procedure for transport simulations is to use the stochastic unsaturated flow simulations of the previous section in order to generate transient simulations of unsaturated transport on the previously obtained 2D flow fields in the (x,z) plane. The steady suction fields (x) analyzed in the previous section are re-used here. The steady moisture content field (x) is deduced from (x) using the van Genuchten moisture-suction relationship: (x) = ( (x)). The steady pore velocity field v(x) is computed from unsaturated Darcy's law, q = -K(h,x) ∇ ( h + z ), or equivalently, q = + K ( ,x) ∇ ( -z), and from the velocity/flux relation, v(x) = q (x)/ (x), where v is the pore water velocity ( "tracer velocity "), while q is the so-called Darcy velocity or areal flux density; suction is the negative of pressure head h ( = -h), and is volumetric water content or "moisture ".

Firstly, starting with the case of statistically anisotropic 2D random soil, the solute plume movement is simulated for vertical flow perpendicular to horizontal stratification. The concentration field is shown in Fig. 19 a three days after the release of the slug from the top surface, and Fig. 19 b shows the length averaged concentration distribution plotted at different times, across the depth. Due to horizontal stratification, the solute plume moves down like a strip, whose thickness is related to the duration of the slug and to the dispersion process. At later times, the plume develops a smooth wave-like form with little lateral spreading, as can be seen in Fig. 19 a andb. Because the horizontal strata are perpendicular to the direction of flow, there is a possibility that the stratification is obstructing vertical flow, and obstructing also the vertical spreading of the plume during its descent.

Solute plume transport is analyzed for 2D anisotropic heterogeneous soil samples, for the case of steady vertical flow parallel to bedding also (the strata are now vertical). The effect of heterogeneity is more significant in the present case (vertical flow parallel to bedding), compared to the case of vertical flow perpendicular to bedding. In the present case, the vertical stratification is parallel to the direction of solute movement. The solute plume tends to disperse in the form of "fingers " at low permeability zones, resulting in more dispersion due to steep velocity gradients. The solute plume seems to have a somewhat "erratic " behavior, which increases with time as it continues to spread laterally.

Secondly, we use a statistically isotropic sample of 2D soil, which was studied in the previous section in terms of unsaturated flow. A steady-state suction field of mean suction 0.307 m was previously obtained from the unsaturated flow simulation ( Fig. 17 a). The corresponding moisture and velocity fields are used as input for the transient simulation of unsaturated solute transport. Fig. 21 a (Left) shows the corresponding solute plume at t = 1 day, i.e., one day after the release of the solute at the top of the domain. The concentration values are averaged in the horizontal direction (normal to the vertical flow), and plotted as a function of depth at various times. Fig. 21 (b) shows these transversely averaged concentration profiles as a function of depth, at different times after the release of the conservative tracer at the top boundary. In comparison with the stratified random medium, the solute shows less dispersion but erratic movements here in the statistically isotropic random medium. For example, some discontinuous patches of concentration occur locally, which can be seen in Fig. 21 a.

Comparing now the different cases, it is observed that the lateral spreading of the solute was more significant in the anisotropic case with flow parallel to bedding [Fig. 20] compared to the isotropic case [Fig. 21] and compared to the case of flow perpendicular to bedding [Fig. 19]. For the isotropic case and the anisotropic perpendicular flow case, the width averaged solute plume appears to be Gaussian, even at large times. This indicates the maintenance of the Fickian dispersion regime even at late times. On the other hand, in the case of "parallel flow ", the concentration profiles seem to lose their initial bell curve shape as time grows, which indicates that the solute plume does not disperse according to Fick's law at later times.

Spatial moments of the solute plume vs. time, and macro-dispersion

The concentration field shows considerable variability in space and time, which was illustrated in the previous section. This variability can be attributed to the heterogeneity of the soil, represented by spatially correlated (auto-correlated) random field parameters. Concentration plume variability is accompanied by large scale dispersion or macrodispersion. This can be assessed through the analysis of the spatial moments of the solute plume vs. time.

For the two-dimensional (2D) flow domains (x,z) considered here, we have calculated the spatial moments of the unsaturated solute plume, as a function of time. It has been observed earlier that the variability of the concentration field is much more important in the direction of the mean flow than in the direction perpendicular to flow. This occurs when the medium is stratified horizontally, orthogonally to the mean flow. Thus, we decide here to analyze the plume moments only in the direction of the mean flow (z).

The zero-order, first-order, and second-order spatial moments of a concentration plume were defined earlier in Eq. ( 20) , ( 21) & ( 22) . The calculated first-order moments of the solute plume, Mz(t), are shown graphically vs. time for all the three cases considered as follows: (i) anisotropic field with imperfect stratification: the case of horizontal strata perpendicular to flow ( Fig. 22a ); (ii) anisotropic field with imperfect stratification: the case of vertical strata parallel to flow ( Fig. 23 a); and (iii) isotropic soil structure ( Fig. 24 a).

It can be seen from the figures that the mean vertical displacements Z(t) = Mz(t)/Mo(t) behaves similarly in all the three cases, that is, linearly with time. Since the Lagrangian velocity of the center of mass is defined as U(t) = dZ/dt, and Z(t) is linear in time, the Lagrangian U(t) is therefore constant during the movement of the plume, that is, U(t) ≈ Uo constant. Recall also that the numerically simulated Eulerian flow and A(t) are Lagrangian versions of the vertical macro-dispersion (D) and macro-dispersivity length (A). As previously noted, these quantities would be constant in time if the macro-dispersion process behaved like Fickian diffusion (here we focus on longitudinal, vertical dispersion).

The calculated (fitted) constant transport coefficients, macrodispersion (Do) and macro-dispersivity (Ao), are given in Table 4 .

As can be seen in Table 4 , for flow perpendicular to stratification (horizontal strata), the longitudinal macro-dispersion coefficient is D o = 0.0755 m 2 /d, and the longitudinal macro-dispersivity length scale is A = 0.021 m. For the flow parallel to stratification (vertical strata), the values are D o = 0.8827 m 2 /d and A = 0.2412 m. Thus, vertical macro-dispersivity is much higher for vertical flow parallel to strata, compared with vertical flow perpendicular to strata. This explains (or this is consistent with) the observed stronger vertical spreading of the plume in the flow direction, parallel to stratification. For the statistically isotropic random soil (which is not stratified), the longitudinal vertical macro-dispersion and macro-dispersivity coefficients are even smaller: Do = 0.006 m 2 /d and A o == 0.001 m. This result indicates that there is much less vertical spreading of the pollutant in a non-stratified soil, compared to stratified soils.

For completeness, we note that stochastic solutions for the macrodispersivity through stratified media are given quasi-analytically by [START_REF] Gelhar | Three-dimensional stochastic analysis of macrodispersion in aquifers[END_REF] for saturated flow. In the present study, macrodispersion is more complex: the flow is unsaturated, with spatially variable velocity vector and spatially variable moisture content, through the randomly heterogeneous soil.

Summary and conclusions

Using a single realization approach, upscaling of unsaturated flow and transport properties of the randomly heterogeneous medium is studied by numerical simulations. The effects of perfect layering and imperfect layering (statistical anisotropy) on flow and transport are analyzed. To study the anisotropy dependence, the mean flow direction is taken as either perpendicular or parallel to the layering (perfect or imperfect). For conductivity upscaling, the local conductivity-suction curve of the heterogeneous soil is assumed to follow Gardner's exponential conductivity-suction model, with two random field parameters: saturated conductivity Ks, and the slope parameter . The unsaturated steady flow experiments proposed in this paper are performed under conditions such that the flow is gravitational on average, with < d /dz > = 0 (zero mean suction gradient). However, this is only an average condition. Since the soil is heterogeneous, there occur small scale fluctuations of the suction gradient (implying that d /dz ≠ 0 almost everywhere locally). These local capillary fluctuations influence the mean suction of the profile, and they generate on average some capillary effects on the effective conductivity curve, which are difficult to predict analytically but are seemingly well captured by the semi-empirical Power Average Model (and by the theoretical spectral perturbation solutions with some limitations).

The heterogeneous parameter of Gardner's exponential conductivity curve plays an important role. It can be interpreted as a pore size distribution parameter, whose inverse represents a capillary length scale (1/ = CAP ). In the analyses of results, this capillary scale is compared to layer thickness or correlation length ( Z ) to produce a dimensionless capillary number Z . The cases Z ~1, Z < 1, and Z ⟩> 1 are explored.

It was found that perfectly stratified unsaturated soils can (sometimes) behave like saturated media, with arithmetically averaged K( ) in parallel flow and harmonically averaged K( ) in the perpendicular flow; however in other cases, depending on flow regime and layer thickness, this classical behavior does not hold good.

The numerically upscaled principal conductivity curves match quite well the PAM, better than the classical means (Arithmetic, Geometric, Harmonic), and also somewhat better than the curves obtained from stochastic spectral perturbation theory. It is observed that the upscaled principal components K ii ( ), obtained numerically and with the PAM along directions "i " perpendicular/parallel to perfect stratification coincide with the harmonic/arithmetic mean curves at low suctions (i.e., near saturation), but deviate from it and come closer to the geometric mean at higher suctions (lower saturation). This motivates the use of the PAM, and appears suitable for generation of approximate upscaled conductivity curves to represent the average behavior of heterogeneous unsaturated soils, especially for instance, the mesh-scale or block-scale conductivity curves in large scale simulation codes.

Overall, it was shown that the behavior of effective unsaturated conductivity (nonlinear and anisotropic) can be captured parametrically via a probabilistic nonlinear 'Power Average Model', where the product Z plays a direct role. The latter model was compared to linearized spectral perturbation theory: the two models are in a way "complementary "; they are not exclusive of each other, and they can be made to coincide totally in some cases.

Numerical simulations of unsaturated solute plume migration, advective and dispersive, were also analyzed, based on the single realization approach. The steady-state unsaturated flow fields computed in the first part were used for that purpose (with the mean flow direction taken vertical in all the cases). The simulated solute plumes exhibit different types of spreading depending on soil structure (isotropic, horizontally stratified, and vertically stratified). In the case of statistically isotropic soil (no strata), lateral and longitudinal spreading is less significant (compared to the other cases), and there will be discontinuous concentration patches in the flow domain.

In the case of statistically anisotropic soils, the solute plume moves down as a smooth strip in case of vertical flow perpendicular to stratification, and the concentration distribution is Gaussian even at later times, indicating Fickian macro-dispersion. In the case of vertical flow parallel to stratification, the longitudinal spreading is more significant. Also, the solute plume is not Gaussian at later times, showing that the macro-dispersion processes diverge from Fick's law as time grows. The upscaled macro-dispersivity length scale (in meters) is found to be higher for vertical flow parallel to stratification, compared to vertical flow perpendicular to the stratification. On the other hand, for the structurally isotropic soil (without stratification), the macro-dispersivity length scale is much smaller compared to two previous results.

The results obtained with flow and transport justify the use of single realization simulations, provided that the size of the realization is much larger than the scale of fluctuations due to heterogeneity. Admittedly, the relevant fluctuation scales to be considered for satisfying this "single realization criterion " may depend on the objectives of the anal-ysis, and they may involve a combination of the fluctuation scales of material properties (which are given) and of the resulting random fields such as suction gradient (which are less well known). At any rate, in the absence of other information, this single realization approach can be implemented on a domain much larger than the known fluctuation scales of material properties, as we have done here. With this provision, the single realization approach is more convenient and less tedious nu-merically than Monte-Carlo simulations. In the latter case, the multiple simulations would deal with uncertainty through an ensemble of repli-cates, while the present single realization approach deals with a large (representative) single heterogeneous sample.

In summary, the effects of heterogeneity on flow and transport, and on effective upscaled coefficients, were studied by performing single re-alization simulations of flow and transport on randomly heterogeneous, auto-correlated soil samples, with various spatial structures. In all cases, the hydraulic properties of the coarse sand from the Upper Hanford site (Washington state) were used for the mean hydraulic properties of the random soil samples. Statistically isotropic and anisotropic soils, with steady flow perpendicular and parallel to stratification, were consid-ered for the numerical flow experiments, from which effective conduc-tivity curves were analyzed. Based on the calculated flow fields, two-dimensional transient simulations of tracer transport were then per-formed, from which spatial moments of the unsaturated solute plume were calculated, leading to macro-dispersion analyses.

  , positive in the unsaturated zone [m] m V = (1 -1/n V ) Van-Genuchten / Mualem parameters (dimensionless exponents) V i tracer velocity in direction "i " [m/s] U Z mean effective vertical tracer velocity of the solute plume [m/s] Z vertical height from the bottom of the soil sample [m] Δx , Δz mesh size in the x and z direction, respectively [m]

Fig. 1 .

 1 Fig. 1. Replicate of the 1D random field parameter Ln(Ks).

Fig. 2 .

 2 Fig. 2. Two-point auto-covariance function of the random field parameter Ln(Ks), plotted vs. separation distance.

Fig. 3 .

 3 Fig. 3. Steady pressure heads plotted across depth. Each vertical pressure profile h(z) corresponds to a different flux rate imposed at the top of the soil column.

Fig. 4 .

 4 Fig. 4. Effective conductivity vs. normalized suction: numerical results compared with spectral perturbation results and with best-fitted power averaging expression (PAM), for flow perpendicular to perfect stratification. (a) Ln(K/Kg) best fit p i = -0.0484. (b) (K/Kg) best fit p i = -0.0634.

Fig. 5 .

 5 Fig. 5. Pressure head contours at steady state, for 1D flow, parallel to perfect strata, at mean pressure ⟨h ⟩ = -0.35 m.

Fig. 6 .

 6 Fig. 6. Effective conductivity vs. normalized suction: numerical results compared with spectral perturbation results and with best-fitted power averaging expression (PAM), for flow parallel to perfect stratification. (a) Ln(K/Kg) best fit p i = -0.1749. (b) (K/Kg) best fit p i = 0.9640.

Fig. 7 .

 7 Fig. 7. Replicate of 2D anisotropic random field parameter Log-saturated hydraulic conductivity Ln K S (x,z).

Fig. 8 .

 8 Fig. 8. Steady state suction fields for two flow configurations: a) flow perpendicular, and b) flow parallel to imperfect strata.

Fig. 9 .

 9 Fig. 9. . Effective conductivity vs. normalized suction: numerical results compared with spectral perturbation results and with best-fitted power averaging expression (PAM), for flow perpendicular to imperfect stratification. (a) Ln(K/Kg) best fit p i = -0.0015 (b) (K/Kg) best fit p i = -0.864.

Fig. 10 .

 10 Fig. 10. Effective conductivity vs. normalized suction: numerical results compared with spectral perturbation results and with best-fitted power averaging expression (PAM), for flow parallel to imperfect stratification. (a) Left: Ln(K/Kg) best fit p i = + +0.0636. (b) Right: (K/Kg) best fit p i = -0.0777.

Fig. 11 .

 11 Fig. 11. Suction variance 2 ( ) [m 2 ], obtained from numerical simulations for the case of perfectly stratified anisotropic soil. a) the case of flow perpendicular to bedding. b) the case of flow parallel to bedding.

Fig. 12 .

 12 Fig. 12. Random permeability field K S ( ⃗ x ) [m/day]: (a) perfectly stratified medium, and (b) imperfectly stratified medium. Comparisons of vertical profiles of K S (z) and (z) through the 2D medium at two horizontal positions: (c) at x = 2 m, and (d) at x = 25 m. Table 3 Optimal values of p obtained by fitting the Power Average model to the numerically upscaled K( Ψ) curves [ln(K) fit]. Note: = HET = HET / CAP where CAP = 1/ , and HET is the correlation scale in the direction of flow (thus HET = X if flow // Ox, or HET = Z if flow // Oz as here). P 1D ⊥(flow//Oz) 1D ⊥(flow//Oz) 1D ⊥(flow//Oz) 1D ⊥(flow//Oz) 1D ⊥(flow//Oz) 1D ⊥(flow//Oz)

Fig. 14 .

 14 Fig. 14. Numerically upscaled conductivity curve Kzz( Ψ) for 1D vertical flow orthogonal to strata.

Fig. 17 .

 17 Fig. 17. Suction field obtained on the statistically isotropic 2D soil sample. (a) The steady-state suction field (x,z). (b) Suction variance plotted versus mean suction.

Fig.. 18 .

 18 Fig.. 18. Numerical results compared with Spectral Perturbation theory for flow through statistically isotropic soil samples. (a) Ln(K/Kg) best fit p i = -0.0484. (b) (K/Kg) best fit p i = -0.1517.

  Fig. 17 a shows such a field for a mean suction of 0.307 m, obtained with a flux rate of 0.50 m/d. Fig. 17 b shows the dependence of suction variance on mean suction. This plot is obtained from several numerical experiments, with the resulting suction fields having different mean values.

  Fig. 20 a shows snapshots of the solute plume at five days, after the release of the solute slug at the top of the flow domain. Fig. 20 b shows the length averaged concentrations at different times for this case.

Fig. 19 .

 19 Fig. 19. (a) Left: Two-dimensional solute plume at t = 3 days, for flow perpendicular to bedding, with mean suction 0.30 m. (b) Width averaged concentration profiles at different times, for flow perpendicular to the bedding (anisotropic, imperfect stratification).

Fig. 20 .

 20 Fig. 20. Time evolution of two dimensional solute plume for the case of vertical flow parallel to bedding (mean suction 0.30 m). (a) Top Left: concentration plume C(x,z,t) at t = 5 days. (b) Right: horizontally averaged concentration profiles as a function of depth, C(z,t), shown at three different times ( t = 1 day, t = 5 days, t = 9 days).

Fig. 21 .

 21 Fig. 21. Time evolution of solute plume for flow through a statistically isotropic soil. (a) Relative Concentration distribution after 1 day. (b) Width averaged value of the same at vertical mid-section of the domain at different times

Fig. 22 .

 22 Fig. 22. (a) First moment Mz(t) and (b) Second-moment Mzz(t) about the center of mass for flow perpendicular to stratification (vertical flow, horizontal strata).

Fig. 23 .

 23 Fig. 23. (a) First moment Mz(t) and (b) Second-moment Mzz(t) about the center of mass for flow to parallel to stratification (vertical flow, vertical strata).

Fig. 24 .

 24 Fig. 24. (a) First Moment Mz(t) and (b): Second Moment Mzz(t) about the center of mass Mz(t) for flow in statistically isotropic soil (vertical flow, no stratification).

Table 1

 1 Statistical Parameters of Gardner's exponential conductivity-suction curve (note: the dimensionless number is defined in terms of the vertical scale Z ).

				Horizontal	Vertical
			Standard	Correl. length	Correl. length
	Parameters:	Mean	deviation	X (m)	Z (m)	Δx (m)	Δz (m)	Lx (m)	Lz (m)
	1. Perfectly stratified "1D " random field; X >> Z ; Z = 0.8133
	LnK S ( * )	0.253	0.771	∞	0.10	-	0.01	-	10
	(1/m)	8.133	1.493	∞	0.10	-	0.01	-	10
	2. Imperfectly stratified "2D " random field; X > Z ; z = 1.6266
	LnK S ( * )	0.253	0.771	2.0	0.20	0.05	0.10	25	50
	(1/m)	8.133	1.493	2.0	0.20	0.05	0.10	25	50
	3. Isotropic "2D " random field; X >> Z ; z = 8.133
	LnK S ( * )	0.253	0.771	1.0	1.0	0.01	0.01	10	10
	(1/m)	8.133	1.493	1.0	1.0	0.01	0.01	10	10
	* where Ks is in (m/d) units.

Table 2

 2 Parameters of the moisture retention curve, based on the van Genuchten's moisturesuction model (assumed constant in space).

	s ( m 3 / m 3 )	r ( m 3 / m 3 )	( -1 )	(dimensionless)
	Saturated moisture	Residual moisture	Scaling factor (1/m)	Exponent
	0.397	0.027	4.306	1.82212

Table 4

 4 Macro-dispersion coefficient (D) and macro-dispersivity length scale (A) for solute transport through different kinds of randomly heterogeneous soil samples.

	Upscaled parameters:	Isotropic soil	Anisotropic soil
	macro-dispersion &
	macro-dispersivity
			Perpendicular flow	Parallel flow
	D [m 2 /d]	0.006 [m 2 /d]	0.0755 [m 2 /d]	0.8827 [m 2 /d]
	A [m]	0.001 [m]	0.021 [m]	0.2412 [m]

Acknowledgments

The authors gratefully acknowledge the useful comments and hints by the anonymous reviewers.

Funding

This research did not receive any specific grant from funding agen-cies in the public, commercial, or not-for-profit sectors.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.advwatres.2020.103655.