Upscaling sorption coefficient for numerical simulation ...

Upscaling sorption coefficient for numerical simulation in physically and chemically heterogeneous porous media 1 Jianlin Fu 2 J. Jaime G´ omez-Hern´ andez

3

Department of Hydraulic and Environmental Engineering Universitat Politècnica de València Camino de Vera, s/n, 46022, Valencia, Spain

Part of this paper appeared in: J.J. Gómez-Hernández, Jianlin Fu, and D. Fernandez-Garcia, 2006. Upscaling retardation factors in 2-D porous media, Calibration and Reliability in Groundwater Modeling: From Uncertainty to Decision Making, IAHS Publication 304, 130-136.

1 The

final version for personal research use is available from fu jianlin [email protected]. author. Currently with Department of Energy Resources Engineering, Stanford University, 367 Panama Street, Stanford, CA 94305, USA. E-mail: fu jianlin [email protected] 3 E-mail: [email protected] 2 Corresponding

Contents 1 Introduction

1

2 Mathematical Formulation of the Problem 2.1 Spatially Variable Distribution Coefficient . 2.2 Reactive Transport Problem . . . . . . . . . 2.3 Upscaled Retardation Factor . . . . . . . . . 2.4 Upscaled versus Effective Retardation Factor

. . . .

2 2 3 4 5

. . . .

7 7 8 9 9

. . . . .

11 11 12 15 17 17

. . . .

. . . .

3 Numerical Method for Upscaling 3.1 Joint Generation of K(x) and Kd (x) . . . . . . 3.2 Computation of Equivalent Retardation Factor . 3.3 Computation of p-norm . . . . . . . . . . . . . . 3.4 Configuration of Numerical Simulations . . . . . 4 Results 4.1 Effect 4.2 Effect 4.3 Effect 4.4 Effect 4.5 Effect

of of of of of

Correlation Length Variance . . . . . . Cross-Correlation . Block Size . . . . . Local Dispersion .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

5 Summary

19

References

21

1

Abstract Even with a simple linear isotherm sorption model, the parameters controlling sorption under field conditions can only be approximated by the values derived from batch experiments. First, the measurement scale and conditions are very different from those at the modeling scale. Second, the parameters are heterogeneous in space and, at most, there is only some sparse information about them at a few locations within an aquifer. For these two reasons, there is a need to consider how to treat the heterogeneity of parameters that control sorption, i.e., retardation factor or distribution coefficient, and a need to establish upscaling rules to transfer the information about parameters at the measurement scale to those at the scale of modeling grid blocks. This paper presents some rules of thumb for upscaling distribution coefficient by accounting for different heterogeneous structures of both hydraulic conductivity and distribution coefficient. Exhaustive numerical simulations are carried out by combining different heterogeneity patterns of hydraulic conductivity and distribution coefficient, the cross-correlation between them, the overall degree of variability, and the block size. It is shown that upscaled retardation factor is a function of the correlation length and variance of the natural log of distribution coefficient and hydraulic conductivity, as well as the cross-correlation between both spatial variables. Geometric mean is a good approximate for upscaled retardation factor under certain conditions, e.g., large variances and small correlation lengths of hydraulic conductivity and distribution coefficient. Keywords: Upscaling rules, chemical heterogeneity, spatial structure, stochastic simulation, reactive transport, distribution coefficient, retardation factor

1

Introduction

Sorption processes arising from radionuclides transport in the groundwater system can greatly modify travel times of dissolved contaminants and, consequently, constitute an important factor affecting the long-term safety assessment of deep geological disposal of highactivity and long-life radioactive wastes. The most commonly used model of sorption relies on the concept of distribution coefficient (Kd ) that relates the distribution of solutes between solid and liquid phases according to a linear isotherm. Sorption may cause solute retention, which is of special importance, for instance, in the context of natural attenuation of contaminants from hydrocarbon releases, because hydrocarbons tend to adsorb onto the solid grains. Sorption together with natural biodegradation makes most of the releases from underground tanks go unnoticed since contaminant plumes do not travel much away from the release location before being biodegraded. Another context in which sorption is especially important is the design of deep geological underground repositories for the disposal of radionuclide wastes, since proper modeling of sorption can make the difference between accepting and rejecting a candidate site. As a consequence, sorption has been recognized as an important mechanism that should be taken into account for modeling the process and fate of contaminants in subsurface porous media. Nowadays, risk analysis and decision making are increasingly dependent on the development of complex numerical models (Gómez-Hern´ andez, 2006). The numerical models for the spatiotemporal prediction of solute migration generally use a discretization of the space with grid blocks much larger than the size of samples taken to the laboratory for the determination of distribution coefficients. At the same time, the distribution coefficient is not homogeneous in space. For these two reasons: disparity of scales between observation and modeling, and spatial heterogeneity, it is necessary to build up some upscaling rules that allow to define equivalent distribution coefficients or retardation factors over areas compatible with the size of the numerical model blocks in a proper manner such that the behavior of the model is as approximate as possible to that of the real, heterogeneous, and sparsely sampled aquifer. This study attempts to find these rules by examining some potential controlling factors for upscaling retardation factor (i.e., R, a parameter that is closely related to Kd ) through synthetic simulations. An upscaled retardation factor is first computed for a coarse block by solving an imposed reactive transport problem on that block subject to physical and chemical heterogeneities. Then, the p-norm of the upscaled retardation factor is found by assuming that the block value only depends on the components of individual cells within the block. By observing the p distribution under various controlling factors through exhaustive numerical experiments, one can obtain an ad hoc upscaling rule on how the block value can be computed from its constitutive components. The exhaustive numerical experiments combine different heterogeneity patterns of hydraulic conductivity and distribution coefficient, the cross-correlation between them, and the overall degree of variability. The effects of the block size (through the travel distance) and the local dispersion coefficient are also examined. The paper is organized as follows: the next section gives some details on the mathematical 1

formulation of the problem for the upscaling of the retardation factor; then, the upscaling method and the design of the numerical simulations are described, followed by some detailed analysis of the numerical results; we conclude the paper by a summary on some of the best upscaling rules. Preliminary results on this study were presented in a conference paper (Gómez-Hern´ andez et al., 2006).

2

Mathematical Formulation of the Problem

The steady state flow of incompressible fluid in saturated porous media without external sources or sinks is described by, ∇ · (K · ∇h) = 0,

(1)

where K = (Kij ), i, j ∈ (x, y, z), is the hydraulic conductivity tensor [LT −1 ], and h is the piezometric head [L]. Following Darcy’s law, fluid velocity is calculated by, v=

q 1 = − K · ∇h, φ φ

(2)

where v = (vx , vy , vz )T denotes the mean pore fluid velocity vector [LT −1 ]; q = (qx , qy , qz )T is the Darcian flux vector [LT −1 ]; and φ is the porosity [-]. These two equations are often numerically solved using the finite-difference or the finite-element method in an aquifer with given hydraulic conductivity distribution and boundary conditions. The governing equation of solute transport subject to advection, dispersion and sorption in a physically and chemically heterogeneous aquifer can be expressed as, ∂c ρb ∂s 1 1 + + ∇ · (qc) − ∇ · (D · ∇c) = 0, ∂t φ ∂t φ φ

(3)

where ρb is the bulk density of matrix [M L−3 ]; c is the aqueous-phase concentration [M L−3 ]; s is the nonaqueous-phase concentration of sorbed solutes, [- or M M −1 ]; and D = (Dij ), i, j ∈ (x, y, z), denotes the local-scale dispersion tensor which is closely related to the velocity vector q, [L2 T −1 ]. For the solution of the transport equation, it is necessary to establish first a model relating c and s.

2.1

Spatially Variable Distribution Coefficient

The sorption process occurring at the representative elemental volume (REV) scale is represented by a reversible linear equilibrium isotherm in which the amount of solutes sorbed onto the solid is proportional to the concentration of dissolved solutes (Chrysikopoulos et al., 1990, 1992; Robin et al., 1991). The constant of proportionality is known as the distribution coefficient or partition coefficient, [M −1 L3 ], i.e.,

2

s(x) ≥ 0, c(x)

Kd (x) =

(4)

where x = (x, y, z) represents the spatial location. This parameter is important since the majority of contaminants tends to interact with the porous matrix and their process and fate in the subsurface formations are thus compounded. Numerous laboratory and field observations suggest that distribution coefficient is spatially variable in natural aquifers (Roberts et al., 1986; Mackay et al., 1986; Robin et al., 1991; Brusseau and Srivastava, 1999). Moreover, potential solute plumes in groundwater system will have a wide range of specific distribution coefficients. This variation can be complex and may be caused by other spatially variable parameters such as temperature, pressure, pH, organic matter contents, sorbent sorptive capacity, sorbent surface charge, etc. Therefore, distribution coefficient is often modeled by a random space function, e.g., Dagan (1989), Bellin et al. (1993), Bosma et al. (1993a), Burr et al. (1994), Hu et al. (1995), and Cvetkovic et al. (1998). Different solutes within distinct porous media may reveal quite different values of distribution coefficient. There is not yet a generally accepted statement about the cross-correlation between distribution coefficient and hydraulic conductivity. Robin et al. (1991) reported a weakly negative correlation between K and strontium ion Kd at the Borden site. On the other hand, a significantly positive correlation between lnK and lnKd for tetrachloroethene at the Borden site was obtained by Allen-King et al. (1998). The relationship in a real field may most likely appear to range from weakly negative to mildly positive (Robin et al., 1991), but there is no a priori reason to exclude perfectly positive or negative correlation (Bosma et al., 1993b). In this work, we assume that the spatial variability of lnKd can be modeled by a multiGaussian stochastic process and lnKd may potentially correlate to log hydraulic conductivity in any way. Such multiGaussian models can be effectively created at the measurement scale and honor all the measurement data, if available.

2.2

Reactive Transport Problem

Under the linear isotherm condition, the spatially variable local retardation factor is related to the distribution coefficient through, R(x) = 1 +

ρb Kd (x), φ(x)

(5)

where φ(x) is the spatially variable porosity. The definition of R(x) in Eq. (5) is a quite general expression, meaning that R(x) is independent of the flow and transport conditions but only depends on the physical characteristics of the porous media. The influence of ρb and φ on retardation is considered to be small in comparison with that of Kd which can vary over several orders of magnitude for some solutes (Robin et al., 1991). Combining Eqs. (3), (4) and (5) leads to a reactive transport equation,

3

∂c + ∇ · (φvc) = ∇ · (D · ∇c). (6) ∂t If the local-scale dispersion is set to zero, the transport equation (6) reduces to an advectiondominant transport problem. If the retardation factor is set to one, it reduces to a conservative or passive transport problem. In all cases, the transport problem can be efficiently simulated in a numerical method with the aid of a random-walk particle-tracking scheme under the Lagrangian framework (Tompson, 1993; Wen and Gómez-Hern´ andez, 1996a; Michalak and Kitanidis, 2000; Fu, 2008; Fu and Gómez-Hern´ andez, 2008, 2009b, 2012; Fu et al., 2010, 2011a, 2011b). The travel times of reactive and nonreactive particles at any spatial position x can be obtained from the numerical simulation, i.e., tr (x) and tp (x), respectively. The corresponding breakthrough curves (BTCs), tr (x) and tp (x), can be easily built for the reactive and nonreactive solute plumes at a control plane x. As a result, a dimensionless travel-time-based retardation factor at any spatial position x may be defined by the ratio of travel times of retarded solutes tr (x) to those of nonreactive conservative or passive solutes tp (x), i.e., Rφ

R(x) =

tr (x) ≥ 1. tp (x)

(7)

Or alternatively, the dimensionless retardation factor may be defined as the ratio of the average interstitial fluid velocity to the propagation velocity of solute in the Eulerian framework (Chrysikopoulos et al., 1990), i.e., R(x) =

v(x) ≥ 1, u(x)

(8)

where v(x) and u(x) represent the the fluid velocity and the solute velocity, respectively. The definition of R(x) in Eqs. (7) or (8) is strictly determined by the flow and transport conditions. The retardation factor reduces the magnitude of the solute velocity and dispersion parameters if R > 1.

2.3

Upscaled Retardation Factor

However, a numerical model that has been generated at the scale at which Kd (x) or R(x) have been measured poses a serious challenge for the numerical simulators capable of solving Eqs. (1) and (3) or (6). An upscaling procedure that can transform the size of the high-resolution models accounting for the heterogeneity observed at measurement scale into coarser ones is needed. Similarly to the definition of block hydraulic conductivity by Journel et al. (1986), an upscaled retardation factor for the coarse block model can be defined by a power average of the fine-scale model values within the block, ·

1 RV = V

¸1/p

Z V

Rωp (x)dx 4

;

ω ∈ V,

(9)

where V represents the block support; Rω is the retardation factor at the cell ω within the block V ; and x is the position vector that sweeps inside the block V . However, since there is no general rule to deduce the value of the p-norm, we have to resort to prior numerical simulations to try to find those rules. We will perform numerical simulations at the fine and coarse grids and we will analyze the controlling factors of RV as a function of spatial heterogeneity as described by some geostatistical parameters. To the best knowledge of the authors, upscaling of retardation factor for coarse-scale numerical simulation has received little attention, although the issue of upscaling has been well explored and developed in the past, e.g., Wen and Gómez-Hern´ andez (1996b) and Renard and de Marsily (1997). This study mainly focuses on the properties of block retardation factors and tries to seek some general rules governing the upscaling of retardation factor via numerical simulations.

2.4

Upscaled versus Effective Retardation Factor

It is important to distinguish between the focus of this paper, i.e., upscaled retardation factors, and the concept of effective retardation factor, which has been the subject of numerous analysis in the past (Tompson, 1993; Burr et al., 1994; Metzger et al., 1996; Espinoza and Valocchi, 1997; Rajaram, 1997; Reichle et al., 1998; Fernandez-Garcia et al., 2005; Dai et al., 2007, 2009; Liu et al., 2008; Deng et al., 2010). When we refer to an upscaled retardation factor, we mean the single block value that will produce the same time delay between passive and sorbed contaminants in their transit within the block as the heterogeneous distribution of retardation factors within that block. The effective retardation factor is linked to the stochastic formulation of Eq. (6) and is defined as Ref f (x) =

htr (x)i , htp (x)i

(10)

where tr (x) and tp (x) are the travel times of the reactive and passive solutes, respectively, and h·i indicates expected value computed through the ensemble. Notice that this value is not defined for a block, but for a given location x. The first problem with this definition is that since a particle released in x will take different paths in the different realizations the end point for which the retardation factor must be computed is undefined. For this reason, an average effective retardation factor is defined considering transport up to a control plane that all particle must traverse, located at a distance d downgradient from the release point: ¯ ˆ ef f (x) = htr (d)i ; R ht¯p (d)i

(11)

where tr (d) denotes the BTC at control plane d for the reactive transport problem; tp (d) represents the BTC at the control plane d for the passive case; and t¯ indicates the arithmetic mean of the travel times at control plane d. Only under some conditions, the effective retardation factor could be assimilated to a block retardation factor for blocks much larger 5

than the correlation lengths of the parameters controlling both flow and transport. That is why, in the literature of effective retardation factor, they always talk about macro-scale values. Field findings show an apparent scale- or time-dependency of the macro-scale retardation factor. In the large-scale field tracer test performed at the Canadian Forces Base Borden aquifer (Roberts et al., 1986; Sudicky, 1986), it was observed that the effective retardation factors of several organic sorptive tracers increase with time, even though laboratory column experiments revealed an equilibrium process governed by a linear isotherm. Burr et al. (1994) numerically examined the behavior of two different formulations of the effective retardation factors, calculated from velocity and arrival time of non-reactive and reactive solutes. They concluded that the joint variability of hydraulic conductivity and distribution coefficient produces a scale-dependent increase in the effective retardation factor from the harmonic average to the arithmetic average. Miralles-Wilhelm and Gelhar (1996) derived an expression of effective retardation factor showing that the scale-dependence is caused by the cross-correlated retardation factor and log hydraulic conductivity field and by the autocorrelation of the local retardation factor field. Moreover, at large time, the effective value is equal to the arithmetic mean. Rajaram (1997) developed theoretical expressions to illustrate that, if there is no correlation between hydraulic conductivity and distribution coefficient, effective retardation factor is initially equal to the harmonic mean of ensemble and increases with time approaching the arithmetic mean at large times. Several mechanisms have been proposed to explain the time dependence of the macroscale retardation factor. The spatial variability of local retardation factors coupled with that of hydraulic conductivity may be one of the dominant mechanisms leading to the scaledependent retardation. Roberts et al. (1986), Goltz and Roberts (1987) and Brusseau and Rao (1989) have attributed this phenomenon to a physical non-equilibrium sorption mechanism. On the other hand, Kabala and Sposito (1991) and Miralles-Wilhelm and Gelhar (1996) have shown that the scale-dependence of effective retardation factor can evolve from spatial variation of retardation factor even under local linear adsorption model. Other authors (e.g., Dagan, 1989; Kabala and Sposito, 1991) attribute this phenomenon to the solute degradation. When the contaminant undergoes the first-order solute decay with a spatially varying decay coefficient, it is observed that effective retardation factor is also a function of time. However, it should be stressed that the upscaled retardation factor RV discussed in this paper is essentially different from the macro-scale effective retardation factor Ref f (x). The upscaled retardation factor RV is simply based on the spatial distribution of the local Kd (x) through Eq. (5) and, in essence, is a physical property of porous media (see Eq. (9)). However, just like its counterpart KV (the block upscaled hydraulic conductivity), the computation of RV may also be affected by the boundary conditions applied to the coarse block for its calculation. By contrast, Ref f (x) not only depends on Kd (x) but also on the flow and transport conditions. The computation of RV is independent of the spatial position of the released particles but Ref f (x) is strictly dependent on the particle’s position from the 6

released source. In summary, a block retardation factor is a value assigned to a well-defined heterogeneous block to replace, for numerical computation purposes, the within block heterogeneity by a single homogeneous value that produces the same retardation to the solute travelling through the block as the heterogeneous distribution. An effective retardation factor is a property of the stochastic model for a given release location and a given control plane, computed through the ensemble of realizations, which could be assimilated to a block value only for very large blocks. Once the difference is clarified both in concept and computational method, it is not at all surprising to observe that the upscaled block retardation factor RV as described in this paper behaves totally different from the macro-scale effective Ref f (x). We also hope that this clarification may help the users to make a correct decision on choosing RV or Ref f (x) for numerical modeling purposes.

3 3.1

Numerical Method for Upscaling Joint Generation of K(x) and Kd (x)

Hydraulic conductivity field K(x) is modeled as a lognormally distributed random space function, with Y (x) = ln K(x) having an exponential semivariogram specified by, ½ · ¸¾ 3r 2 γY (r) = σY 1 − exp − , λY where r is the two-point separation distance, σY2 is the logconductivity variance, and λY is the correlation length. Let’s redefine the value of the retardation factor as R(x) = 1+Kd (x), which is equivalent to assuming that ρ/φ(x) is one in Eq. (4), or that Kd (x) is a surrogate parameter that incorporates density and the variability of porosity. Then, consider that Kd (x) can be also modeled through a lognormally distributed random space function, with W (x) = ln Kd (x) having an exponential semivariogram specified by, ½ · ¸¾ 3r 2 γW (r) = σW 1 − exp − . λW The cross semivariogram between Y (x) and W (x) is similarly modeled by an exponential model, ½ · ¸¾ 3r 2 γY W (r) = σY W 1 − exp − . λY W The sill of the cross semivariogram directly p reflects the magnitude of the correlation be2 2 ∈ [−1, 1]. Once the semivariogram and tween Y (x) and W (x), i.e., ρY W = ±σY W / σY2 σW cross semivariogram matrices are specified, the sequential simulation algorithm (GómezHern´ andez and Journel, 1993) can be used for the joint generation of realizations of Y (x) 7

and W (x). Note that λ and σ 2 are two important hyperparameters to describe the spatially heterogeneous structure of the porous media. Upscaling of model parameters generally calls for taking into account such heterogeneity explicitly (e.g., Leung and Srinivasan, 2011). The analysis following will try to establish rules for the upscaling of the retardation factor as a function of these spatial structure parameters.

3.2

Computation of Equivalent Retardation Factor

The main idea of calculating the equivalent retardation factor is to seek individual values, RV , that are capable of reproducing the key features of the reactive mass flux BTC of a random retardation factor field, R(x), subject to identical flow and transport boundary conditions within the same physically heterogeneous media. Since it is impossible to find a single value to match the entire curve, three equivalent values are separately calculated at (5%) (50%) (95%) different matching points, i.e., RV , RV , and RV , which stand for the early, median and late arrived particles, respectively. It should be emphasized that the selection of the most appropriate part of BTC to calculate the block retardation factor is a very important step for the correct application of upscaled retardation factors in engineering environmental operations, e.g., to monitor the extent and degree of groundwater contamination from a known source (Fu and Gómez(5%) Hern´ andez, 2009b; Fu et al., 2011b). RV better reproduces the earliest part of the reference BTC and represents the fastest particles arriving to the control plane as needed, for example, for the design of a radioactive underground repository. The earliest arrivals in the BTC follow the fastest pathways between the release source and the control plane, which are dominated by preferential flow and reactive transport paths. Failing to account for this will yield a too conservative conclusion in risk analysis in that the real arrival time may be too much (50%) faster than the true ones (G´ omez-Hern´ andez and Wen, 1998). RV reproduces the median part of the reference BTC and it is generally linked with the peak arrival time. Public officials assessing health risks associated with contaminant exposure in a drinking water (95%) supply system may be most concerned with this parameter (Lemke et al., 2004). RV reproduces the tail part of the reference BTC and denotes the slowest particles arriving to the control plane as needed for mass removal calculations in remediation problems. The late travel times reflect a more integral behavior, including the effects of flow and reactive transport barriers. An aquifer remediation design without considering such feature may fail because the resident contaminants will be removed more slowly than expected (Wagner and Gorelick, 1989). The procedure to calculate the three block values for a given physically heterogeneous K(x) is the following. First, the BTC tr (x) at a control plane x in the heterogeneous field is obtained by solving the reactive solute transport problem with variable retardation factor R(x). Then, the same problem is solved with a homogeneous R0 value equal to 1 and its BTC (β) (β) (β) determined tc (x). The RV would be equal to the ratio tr (x)/tc (x), for β = 5%, 50%, (β) and 95%. Fig. 1 shows the curves computed with homogeneous values for RV , β = 5%, 8

50%, and 95% for two different blocks. Note that the curves match the three control points in tr (x). Also note that the three curves are not simply shifted versions of one another, since as time passes the dispersion of the BTC curves increase with time.

3.3

Computation of p-norm

In some cases the block retardation factor may be more influenced by low values within the averaging block domain, while in other cases high values are more influential. The power-norm of the values within the block, R(x), allows easily to shift the focus among low, intermediate or high values. Following Journel et al. (1986) and Gómez-Hern´ andez and Gorelick (1989) who introduced the use of a power average to calculate equivalent hydraulic conductivities, we can similarly obtain the p-norm of the block retardation factor as, "

n

1X p RV = R n i=1 i

#1/p ,

where n is the number of grid nodes within the block. Note that p > 1 corresponds to an average mostly influenced by the largest values; p = 1 corresponds to the arithmetic average where all components receive an equal weight; p = 0 corresponds to the geometric average; p = −1 corresponds to the harmonic average, and for p < −1 the power average is mostly influenced by the smallest values (Journel et al., 1999). Since R(x) is always larger than one, the power average formula is strictly monotonous. Once the block retardation value RV has been computed following the procedure delineated above, a unique p value can be determined relating the R(x) within the block and RV . In addition, it is worth to anticipate that the relative values of the dispersion of the two curves tr (x) and tc (x) (recomputed for the constant RV obtained by the p-norm) determines (5%) the sign of the p-norm. If σt2r > σt2c (such as is the case in Fig. 1(A)), the p-norm of RV (95%) will be smaller than 0 but the p-norm of RV will be larger than 0 . If σt2r < σt2c , (such as (5%) (95%) is the case in Fig. 1(B)), the p-norm of RV will be larger than 0 but the p-norm of RV will be smaller than 0.

3.4

Configuration of Numerical Simulations

A synthetic two-dimensional isotropic aquifer is chosen to simulate the flow and transport of reactive solutes for computation of block retardation factor RV and the corresponding p-norm. The size of the computational domain is 200 by 200 which is discretized into 100 by 100 cells with each cell having a dimension 2 by 2. The correlated Y (x) and W (x) fields are generated by a sequential simulation algorithm using the public domain program GCOSIM3D (Gómez-Hern´ andez and Journel, 1993). Both the expected values of ln K(x) and ln Kd (x) are set to zero with an isotropic exponential covariance function model, i.e., µY (x) = µW (x) = 0 such that µK(x) = 1 and µR(x) = 2; and λx = λy for both K(x) and R(x) fields. The 9

1

(A) Computation of R V (ρYW=-1.0)

0.9 0.8

Frequency

0.7 0.6 0.5 0.4 (95%) = 2.03 RV R (50%) =1.95 V R (5%) =1.88 V reference

0.3 0.2 0.1 0 1800

2000

2200

2400

2600

2800

3000

Time

1

(B) Computation of R V (ρYW=1.0)

0.9 0.8

Frequency

0.7 0.6 0.5 0.4 (95%) = 2.17 RV R (50%) =2.22 V R (5%) =2.29 V reference

0.3 0.2 0.1 0

2200

2400

2600

2800

3000

3200

Time

2 Figure 1: Computation of three upscaled retardation factors (λY,W = 1, σY,W = 0.5): (A) ρY W = −1.0; (B) ρY W = 1.0. The reference solid curves denotes the base reactive BTC tr with a spatially variable retardation factor field. The other three curves represent three cases with different constant retardation factors tc . Note that their spreadings are different.

10

geochemical parameter R(x) of the heterogeneous aquifer is upscaled to a homogeneous constant RV with the dimension of the single block equal to 200 by 200. The upscaled ratio, therefore, is of 1002 : 1. Totally 9 × 5 × 9 = 405 scenarios of random combinations of Y (x) and W (x) are exhaustively simulated to form the distribution pattern of the upscaled retardation factors: nine different correlation lengths (λY,W = 1, 2, 5, 10, 20, 30, 40, 50, 100 2 [cells]), five different variances (σY,W = 0.1, 0.5, 1, 2, 6), and nine kinds of cross-correlation between ln K(x) and ln Kd (x), i.e., from negative-, non-, to positive-correlation (ρY W = −1.0, −0.75, −0.5, −0.25, 0.0, 0.25, 0.5, 0.75, 1.0). Note that Y (x) remains unchanged when the cross-correlated W (x) fields with different ρY W are generated. The numerical experiment is designed as a steady-state, natural-gradient, uniform flow in a square field with prescribed heads on two opposite sides and impermeable boundaries on the other two faces. Five-point block-centered finite-difference method is employed to solve the steady flow problems (Eqs. (1) and (2)). There is ten unit-free head drop from the left boundary to the right one along the flow direction such that the imposed mean hydraulic gradient is 0.05 along the principal flow direction, while the top and bottom boundaries are considered as impervious. A constant-displacement random-walk particle-tracking algorithm is employed to solve the conservative and reactive solute transport problems (Eq. (6)), which proves more computationally efficient than the constant-time-step scheme (Wen and GómezHern´ andez, 1996a). For each realization, 40,000 particles released from the left boundary are tracked down until they exit the computational domain. The particles are released by distributing them uniformly along the left boundary, which corresponds to the case of a non-point source. The particles travel through most of the domain when moving along the fluid flow, therefore, sampling the entire K(x) and R(x) fields, hence, we can say that their overall spatial variability has been captured (Fu et al., 2011a). This configuration is also very common in computing the equivalent block value for upscaling hydraulic conductivity (Wen and Gómez-Hern´ andez, 1996b; Fu and Gómez-Hern´ andez, 2009a; Fu et al., 2010, 2011a, 2011b). In addition, in order to examine the effects of the block size (through the change of the travel distance) and of local dispersion, several sets of additional experiments that share the simulation configurations described above are conducted. In the case of the existence of local dispersion, the longitudinal (αL ) and transverse (αT ) dispersivity coefficients are equal to 0.1 and 0.01, respectively.

4 4.1

Results Effect of Correlation Length

Correlation range (or correlation length) is a measure of the distance between two points at which they cease to be correlated [Dagan, 1989]. It reflects the continuity of a formation along a certain direction. Fig. 2 shows how the power factor p is a function of the correlation 2 2 length. An identical trend for both σY,W = 0.1 and σY,W = 1 is that, as the correlation 11

length decreases, p-norm tends to converge to 0 which corresponds to the geometric mean of the heterogeneous field. This suggests that it is proper to use the geometry mean as representative of the random R field to simulate flow and transport of reactive solutes as long as the medium has a short correlation length compared to the size of the coarse grid block. A short correlation length is preferred in applying such geometric mean rule probably because a small λ weakens the spatial dependence between points so that the influence from the adjacent coarse blocks decreases correspondingly and the geometric averaging that imposes on a block is less affected by its neighbors. (5%) It also can be observed from Fig. 2 that, for RV , an estimate by the geometric mean (95%) will introduce a positive error, i.e., the real value is overestimated. For RV , the geometric (50%) mean will underestimate the real block value. For RV , the geometric mean is the best approximate. Furthermore, the higher the degree of heterogeneity is, the more precise the estimate through the geometric mean is. The fluctuation of p-norm with a variance equal to 0.1 is more severe than that with variance equal to 1. Although, it is also true that with smaller heterogeneity, the absolute value of the different means displays less variability than for large heterogeneity (at the limit, for zero variability, all p-norms would be the same). In addition, it is worthwhile noting that the upscaled value is regularly distributed on both sides around the geometric mean, whose exact position is controlled both by the correlation coefficient and the matching position. A general pattern is that the perfectly negative and positive cases are located on the top and bottom, respectively, while the mild one resides between the two limits.

4.2

Effect of Variance

Variance is a measure of the spreading of a random variable around its mean. The heterogeneity of both hydraulic conductivity and retardation factor in terms of variance has a significant effect on block retardation factors. Fig. 3 shows that the block retardation factor converges to the geometric mean of the random field as the variance increases. The trend is more apparent as the correlation length increases. For instance, the p-norm distribution with a correlation length equal to 10 produces a faster convergence than that with a correlation length equal to 1. Note that the scales of the vertical axis may be different between columns. The reason that a large variance results in the geometric mean as being the p-norm to be applied for the upscaling of the retardation factor might be explained as follows. From Eq. (8), the magnitude of the reactive solute transport velocity u(x) is reduced compared to that of the conservative case v(x) since R(x) ≥ 1, which means that the fluctuation or the degree of heterogeneity of velocity field σu2 decreases due to the introduction of the retardation factor. In other words, the degree of heterogeneity of retardation factor σR2 serves to offset that of 2 velocity σu2 resulted from the fluctuation of hydraulic conductivity σln K . The fluctuation of concentration or travel time (σc2 or σt2 ) is also reduced correspondingly. Moreover, a higher degree of heterogeneity in retardation factor yields a more compact distribution in concentration or travel time (or equivalently, a smoother velocity field or a lower degree of 12

(5%)

15

p norm of R V

2

p norm of R V

5

0 ρYW=-1.0 ρYW=-0.75 ρYW=-0.5 ρYW=-0.25 ρYW=0 ρYW=0.25 ρYW=0.5 ρYW=0.75 ρYW=1.0

-5

-10

0

25

50

75

2


-2

-4

-6

100

0

25

correlation length

(50%)

15

p norm of R V

2


-5

-10

p norm of R V

100

2

(σY,W=1)

0

25

50

75

2


-2

-4

-6

100

0

25

correlation length

(95%)

15

p norm of R V

50

75

100

correlation length

2 (σY,W=0.1)

(95%)

6

p norm of R V

2 (σY,W=1)

4

power factor (p)

10

power factor (p)

75

4

power factor (p)

power factor (p)

(50%)

6

(σY,W=0.1)

5

5


-5

-10

-15

50

correlation length

10

-15

2

(σY,W=1)

4

power factor (p)

power factor (p)

10

-15

(5%)

6

(σY,W=0.1)

0

25

50

75

2


-2

-4

-6

100

correlation length

0

25

50

75

100

correlation length

Figure 2: Effect of correlation length λY,W on upscaled retardation factor: The left column 2 displays the p norm distribution of a small variance (σY,W = 0.1) and the right column shows 2 that of a larger variance (σY,W = 1). The first, the second, and the third row display the (5%) (50%) (95%) p-norm distribution of RV , RV , and RV , respectively. 13

(5%)

p norm of R V

3

6

2

4

1 0 ρYW=-1.0 ρYW=-0.75 ρYW=-0.5 ρYW=-0.25 ρYW=0 ρYW=0.25 ρYW=0.5 ρYW=0.75 ρYW=1.0

-1 -2 -3 -4

0

1

2

3

4

(5%)

8

(λY,W=1)

power factor (p)

power factor (p)

4

5

p norm of R V


-2 -4 -6 -8

6

0

1

2

3

variance

(50%)

p norm of R V

3

2

2


-1 -2 -3

0

1

2

3

4

(50%)

4

(λY,W=1)

3

-4

5

p norm of R V

(95%)


-2 -3 -4

6

0

1

2

3

2

4


-1 -2 -3

1

2

3

4

(95%)

8

(λY,W=1)

6

0

4

5

6

variance

3

-4

6

(λY,W=10)

-1

power factor (p)

power factor (p)

p norm of R V

5

1

variance

4

4

variance

power factor (p)

power factor (p)

4

(λY,W=10)

5

p norm of R V

(λY,W=10)


-2 -4 -6 -8

6

variance

0

1

2

3

4

5

6

variance

2 Figure 3: Effect of variance σY,W on upscaled retardation factor: The left column displays the p norm distribution of λY,W = 1 and the right column shows that of λY,W = 10.

14

heterogeneity in hydraulic conductivity). The block retardation factors are less affected by a narrower travel time distribution which corresponds to a larger variance in retardation factor. Therefore, a large variance of retardation factor makes the geometric mean be the best estimator. Furthermore, the cross-correlation between lnK and lnKd also has an influence. For (5%) RV , the geometric mean will overestimate the upscaled value in the case of a perfectly negative correlation but underestimate the upscaled value in the case of a perfectly positive (95%) one. On the other hand, for RV , the results is just the opposite, the geometric mean will underestimate the upscaled value in the case of a perfectly negative correlation but overestimate the upscaled value in a perfectly positive case.

4.3

Effect of Cross-Correlation

The cross-correlation between ln Kd and ln K describes the relationship between the adsorption or desorption of solutes and the physical properties of the porous media. For an aquifer, a positive correlation reveals that the adsorption is proportional to the velocity of solutes while a negative correlation displays an inverse-proportional relationship. In other words, if ρY W > 0, the retardation effect is enhanced as velocity increases; if ρY W < 0, the retardation effect is reduced as velocity increases. The upscaled retardation factors are highly affected by the cross-correlation between log hydraulic conductivity and log distribution coefficient. For a strongly negative correlation, e.g., ρY W = −0.75, the block equivalent retardation factors (5%) (95%) for early times RV are smaller than the ones corresponding to late time arrivals RV . For a strongly positive correlation, e.g., ρY W = 1, the block equivalent retardation factors for early times are larger than those associated with late arrival times (Fig. 4). It is also noted that for a weakly negative correlation, e.g., ρY W = −0.25, the block retardation factors share the same trend as for strongly negative correlation coefficients. Therefore, the block RV values based on the early breakthrough times and their counterparts for the late breakthrough times make remarkable difference, which implies that an upscaling scheme should be adjusted to the specific engineering aim. For example, the p-norm estimate with a goal to remove contaminants should be larger than the one required for monitoring the earliest arrival of contaminants if a negative cross-correlation is observed in the field. Otherwise, if the cross-correlation coefficient is positive, the former should be smaller than the latter. The reasoning on the effect of cross-correlation on the upscaled retardation factors might be explained as follows. From Eq. (8), one can easily find a proportional relationship for the retarded solute velocity, u(x) =

v(x) K(x) ∝ , R(x) Kd (x)

(12)

which holds simply because K(x) relates to v(x) by Eqs. (1) and (2) and R(x) to Kd (x) by Eq. (5). From this relation, one can find that, for a given field ln K(x), a positively crosscorrelated ln Kd (x) will reduce the extremes of u(x) such that its variance σu2 decreases; a 15

p norm of R V (σY,W=0.1; ρYW=-0.75) 2

15

2

4

power factor (p)

power factor (p)

10

5

0

-5 R (5%) V (50%) RV R (95%) V

-10

-15

p norm of R V (σY,W=1.0; ρYW=-0.75)

6

0

25

50

75

2

0

-2 R (5%) V (50%) RV R (95%) V

-4

-6

100

0

25

correlation length

2 p norm of R V (σY,W=0.1; ρYW=-0.25)

15

power factor (p)

power factor (p)

0

-5 R R R

-10

0

25

50

(5%) V (50%) V (95%) V

75

2

0

-2 R (5%) V (50%) RV R (95%) V

-4

-6

100

0

25

correlation length

50

75

100

correlation length

2 p norm of R V (σY,W=0.1; ρYW=1.0)

15

2 p norm of R V (σY,W=1.0; ρYW=1.0)

6

4

power factor (p)

10

power factor (p)

100

4

5

5

0

-5 R R R

-10

-15

75

2 p norm of R V (σY,W=1.0; ρYW=-0.25)

6

10

-15

50

correlation length

0

25

50

75

(5%) V (50%) V (95%) V

2

0

-2 R (5%) V (50%) RV R (95%) V

-4

-6

100

correlation length

0

25

50

75

100

correlation length

Figure 4: Effect of cross-correlation ρY W on upscaled retardation factor: The left column 2 2 displays the p distribution of σY,W = 0.1 and the right column shows that of σY,W = 1.

16

negatively cross-correlated ln Kd (x) will enlarge the extreme values of u(x) such that its variance σu2 increases. On the other hand, from Eq. (3), the travel times are inversely proportional to the retarded solute velocity, tr (x) ∝

1 . u(x)

(13)

Therefore, a positive cross-correlation (ρY W > 0) that reduces the extremes of the transport velocity of reactive solutes u(x) will narrow the spreading of the travel times tr (x) and even (5%) make σt2r smaller than σt2c . The p-norm of the computed RV , accordingly, is larger than (95%) 0 while that of RV is less than 0. The block equivalent retardation factors for early times are larger than those associated with late arrival times. By contrast, a negative crosscorrelation (ρY W < 0) that increases the extremes of u(x) will enlarge the spreading of the travel times tr (x) and even make σt2r much larger than σt2c . The p-norm of the computed (5%) (95%) RV , accordingly, is less than 0 while that of RV is larger than 0. The block equivalent retardation factors for early times are smaller than those associated with late arrival times. Note that tc (x) is independent of the spatial variation of ln Kd (x).

4.4

Effect of Block Size

Considering that the time- or scale-dependence of the effective retardation factor is observed by the laboratory and field experiments and is confirmed through both theoretic closures and numerical simulations, an experiment is designed in this section to check the effect of six different travel distances on upscaled retardation factors at the same matching points, i.e., (5%) (50%) (95%) RV , RV , and RV . The experiment proceeds like this: we set various control planes at different positions, i.e., 50%, 60%, 70%, 80%, and 90% of the whole domain downstream in the main flow direction. The same flow and transport simulations are performed as before. The results are compared with those of the control plane at the end of the random field. This experiment essentially examines the effect of the block size on the upscaling of retardation factor, because different travel distances correspond to different block sizes. Fig. 5 shows that both RV and p have no evident uniform increase or decrease at different distances in the case of a short (smaller than 10) correlation length. Whether for the early, the median or the late arrival time, they share similar patterns. It suggests that block retardation factor is not affected by travel distance once the block large compared to the correlation length of the random field. In other words, it is scale-invariant or independent of the block size if the correlation length is short enough.

4.5

Effect of Local Dispersion

The effect of the local dispersion on flow and transport of solute is often considered to be limited and secondary compared to that of hydraulic conductivity. In order to evaluate the effect of the local dispersion coefficient on the upscaled retardation factor, we compute the 17

9

(5%)

p distribution of R V

3.5

(λY,W=10)

(5%)

RV

distribution (λY,W=10)

6

0

ρ=-0.75; σ2=0.1 2 ρ=-0.75; σ =0.5 ρ=-0.75; σ2=1 2 ρ=0; σ =0.1 ρ=0; σ2=0.5 2 ρ=0; σ =1 2 ρ=1; σ =0.1 2 ρ=1; σ =0.5 2 ρ=1; σ =1

-3

-6

-9

equivalent R

power factor (p)

3 3

50

60

70

80

90

100

110

2.5


2

1.5

120

50

60

70

control plane

9

(50%)


80

90

100

110

120

control plane

3.5

(λY,W=10)

(50%)

RV


6

0


-3

-6

-9

equivalent R

power factor (p)

3 3

50

60

70

80

90

100

110

2.5


2

1.5

120

50

60

70

control plane

9

(95%)


80

90

100

110

120

control plane

3.5

(λY,W=10)

(95%)

RV


6

0

ρ=-0.75; σ2=0.1 ρ=-0.75; σ2=0.5 ρ=-0.75; σ2=1 ρ=0; σ2=0.1 ρ=0; σ2=0.5 ρ=0; σ2=1 2 ρ=1; σ =0.1 ρ=1; σ2=0.5 ρ=1; σ2=1

-3

-6

-9

equivalent R

power factor (p)

3 3

50

60

70

80

90

100

110

2.5


2

1.5

120

control plane

50

60

70

80

90

100

110

120

control plane

Figure 5: Effect of the block size (travel distance) on upscaled retardation factor in heterogeneous velocity fields: The left column displays the p norm distribution and the right column shows the RV distribution.

18

block value using the same parameter configurations but adding the longitudinal (αL ) and transverse (αT ) dispersivity coefficients equal to 0.1 and 0.01, respectively. Fig. 6 shows the difference of upscaled retardation factors between the cases without the local dispersion and those with the dispersion. Generally there are some differences but with no big changes. The only significant difference occur for those fields with a variance equal to 0.1 and perfectly positive cross-correlation. This may be due to the fact that a small perturbation will induce a big difference of p-norm for the smaller variance case as analyzed before.

5

Summary

A set of exhaustive 2D flow and reactive transport simulations in heterogeneous, (cross)correlated hydraulic conductivity and retardation factor fields is conducted to study upscaled retardation factors. Several main findings from this work might be summarized as follows. The upscaled retardation factor is a function of the correlation length of both hydraulic conductivity and distribution coefficient fields, their variance, and the cross-correlation between them. The cross-correlation plays an important role in the upscaling process. For a negative correlation, upscaled retardation factor for early times is smaller than that for late times; for a positive correlation, the block value for early times is larger than that for late times. Block retardation factors tend to converge to the geometric mean as the variance increases, which is more apparent for larger correlation lengths. The block size (or the travel distance) has no obvious influence on the upscaled retardation factors if it is large with respect to the correlation length. The local dispersion is not dominantly influential factors for upscaling. An upscaling scheme should be tailored to a specific engineering aim since the front and tail of a solute plume may display quite different block retardation factors. But in general, the geometric mean of a random field is a good approximate for the equivalent retardation factor under certain conditions, e.g., a short correlation length and a relatively large variance. Note that the upscaling rules built in this work are based on 2D multi-Gaussian fields. The 3D problem might display big differences from the 2D case, but for the 1D case the rules may be quite different. The 1D problem is especially interesting in practice: for example, it is useful for the cases involving streamline tracking and helpful to understand the flow and transport properties in the case of the existence of fractures which usually forms a peculiar path contrasted to the surrounding rock matrix. An extension to the 1D and 3D problems is expected in a future work. Acknowledgements We thank the Spanish Ministry of Science and Innovation for the financial support through projects CGL2011-23295 and ESCOPA REN2002-02428. We also thank support from the Spanish Nuclear Waste Management Company (ENRESA). A doctoral fellowship and an 19

2 p distribution (σY,W=1; ρY,W=-0.75)

6

2 R V distribution (σY,W=1; ρY,W=-0.75)

3.5

4

equivalente R

power factor (p)

3

without α

2

with α 0

-2

-6

2

R (5%) V (50%) RV R (95%) V

-4

0

25

50

75

without α with α

2.5

1.5

100

R (5%) V (50%) RV R (95%) V

0

25

correlation length

2 p distribution (σY,W=1; ρY,W=-0.25)

6

50

75

100

correlation length

2 R V distribution (σY,W=1; ρY,W=-0.25)

3.5

4

without α

equivalente R

power factor (p)

3 2

with α 0

-2

-6

0

25

50

75

with α

2

R (5%) V (50%) RV R (95%) V

-4

without α 2.5

1.5

100

R (5%) V (50%) RV R (95%) V

0

25

correlation length

p distribution (σY,W=1; ρY,W=1.0) 2

6

50

75

100

correlation length

R V distribution (σY,W=1; ρY,W=1.0) 2

3.5

4

equivalente R

power factor (p)

3 2

without α

0

with α

-2

-6

0

25

50

75

without α with α

2

R (5%) V (50%) RV R (95%) V

-4

2.5

R (5%) V (50%) RV R (95%) V 1.5

100

correlation length

0

25

50

75

100

correlation length

Figure 6: Effect of local dispersion coefficient on upscaled retardation factor: The left column displays the p norm distribution and the right column shows the RV distribution. Solid lines display the cases without local dispersion and dashed lines correspond to the cases with local dispersion.

20

extra travel grant awarded to the first author by the Universitat Politècnica de València, Spain, is gratefully acknowledged. We are grateful to Dr. Daniel Fernàndez-Garcia for his comments and help.

References [1] Allen-King, R.M., R.M. Halket, D.R. Gaylord, and M.J.L. Robin, 1998. Characterizing the heterogeneity and correlation of perchloroethene sorption and hydraulic conductivity using facies-based approach, Water Resources Research, 34(3), 385-396. [2] Bellin, A., A. Rinaldo, W.J.P. Bosma, S.E.A.T.M. van der Zee, and Y. Rubin, 1993. Linear equilibrium adsorbing solute transport in physically and chemically heterogeneous porous formations, 1. Analytical solutions, Water Resources Research, 29(12), 4019-4030. [3] Bosma, W.J.P., and S.E.A.T.M. van der Zee, 1993a. Transport of reacting solute in a one-dimensional chemically heterogeneous porous medium, Water Resources Research, 29(1), 117-131. [4] Bosma, W.J.P., A. Bellin, S.E.A.T.M. van der Zee, and A. Rinaldo, 1993b. Linear equilibrium adsorbing solute transport in physically and chemically heterogeneous porous formations, 2. Numerical results, Water Resources Research, 29(12), 4031-4043. [5] Brusseau, M.L., and P.S.C. Rao, 1989. Sorption nonideality during organic contaminant transport in porous medium, Crit. Rev. Environ. Control, 19(1), 33-99. [6] Brusseau, M.L., and R. Srivastava, 1999. Nonideal transport of reactive solutes in heterogeneous porous media: 4. Analysis of the Cape Cod natural-gradient field experiment, Water Resources Research, 35(4), 1113-1125. [7] Burr, D.T., E.A. Sudicky, and R.L. Naff, 1994. Non-reactive and reactive solute transport in three-dimensional heterogeneous porous media: Mean displacement, plume spreading, and uncertainty, Water Resources Research, 30(3), 791-815. [8] Cvetkovic, V., G. Dagan, and H. Cheng, 1998. Contaminant transport in aquifers with spatially variable hydraulic and sorption properties, Proc. R. Soc. Lond. A, 454, 21732207. [9] Chrysikopoulos, C.V., P.K. Kitanidis, and P.V. Roberts, 1990. Analysis of onedimensional solute transport through porous media with spatially variable retardation factor, Water Resources Research, 26(3), 437-446.

21

[10] Chrysikopoulos, C.V., P.K. Kitanidis, and P.V. Roberts, 1992. Macrodispersion of sorbing solutes in heterogeneous porous formations with spatially periodic retardation factor and velocity field, Water Resources Research, 28(6), 1517-1530. [11] Dagan, G., 1989. Flow and transport in porous formations, Springer-Verlag, New York. [12] Dai, Z., A. Wolfsberg, Z. Lu, and R. Ritzi, 2007. Representing Aquifer Architecture in Macrodispersivity Models with an Analytical Solution of the Transition Probability Matrix, Geophysical Research Letters, 34, L20406, doi:10.1029/2007GL031608. [13] Dai, Z., A. Wolfsberg, Z. Lu, and H. Deng, 2009. Scale dependence of sorption coefficients for contaminant transport in satured fractured rocks, Geophysical Research Letters, 36, L01403. [14] Deng, H., Z. Dai, A. Wolfsberg, Z. Lu, M. Ye, and P. Reimus, 2010. Upscaling of reactive mass transport in fractured rocks with multimodal reactive mineral facies, Water Resources Research, 46(6), W06501. [15] Espinoza, C., and A.J. Valocchi, 1997. Stochastic analysis of one-dimensional transport of kinetically adsorbing solutes in chemically heterogeneous aquifers, Water Resources Research, 33(11), 2429-2445. [16] Fernandez-Garcia, D., T.H. Illangasekare, and H. Rajaram, 2005. Differences in the scale dependence of dispersivity and retardation factors estimated from forced-gradient and uniform flow tracer tests in three-dimensional physically and chemically heterogeneous porous media, Water Resour. Res., 41, W03012, doi:10.1029/2004WR003125. [17] Fu, J., 2008, A Markov chain Monte Carlo method for inverse stochastic modeling and uncertainty assessment, Unpublished doctoral dissertation, Technical University of Valencia, Spain, pp.158. [18] Fu, J., and J.J. Gómez-Hernández, 2008, Preserving spatial structure for inverse stochastic simulation using blocking Markov chain Monte Carlo method, Inverse Problem in Science and Engineering, 16 (7), 865-884. DOI: 10.1080/17415970802015781 [19] Fu, J., and J.J. Gómez-Hernández, 2009a, A blocking Markov chain Monte Carlo method for inverse stochastic hydrogeological modeling, Mathematical Geosciences, 41(2), 105-128, DOI: 10.1007/s11004-008-9206-0. [20] Fu, J., and J.J. Gómez-Hernández, 2009b, Uncertainty assessment and data worth in groundwater flow and mass transport modeling using a blocking Markov chain Monte Carlo method, Journal of Hydrology, 364, 328-341, DOI: 10.1016/j.jhydrol.2008.11.014. [21] Fu, J., and J.J. Gómez-Hernández, 2012, A gradient-based blocking Markov chain Monte Carlo method for large-scale stochastic inverse modeling and uncertainty quantification, Mathematical Geosciences, submitted. 22

[22] Fu, J., H.A. Tchelepi, and J. Caers, 2010, A multiscale adjoint method to compute sensitivity coefficients for flow in heterogeneous porous media, Advances in Water Resources, 33(6), 698-709, DOI: 10.1016/j.advwatres.2010.04.005. [23] Fu, J., J. Caers, and H.A. Tchelepi, 2011a, A multiscale method for subsurface inverse modeling: single-phase transient flow, Advances in Water Resources, 34(8), 967-979. DOI: 10.1016/j.advwatres.2011.05.001. [24] Fu, J., C.L. Axness, and J.J. Gómez-Hernández, 2011b, Upscaling transmissivity in the near-well region for numerical simulation: A comparison on uncertainty propagation, Engineering Applications of Computational Fluid Mechanics, 5(1), 49-66. [25] Goltz, M.N., and P.V. Roberts, 1987. Using the method of moments to analyze threedimensional diffusion-limited solute transport from temporal and spatial perspectives, Water Resources Research, 23(8), 1575-1585. [26] Gómez-Hernández, J.J., and S.M. Gorelick, 1989. Effective groundwater model parameter values: influence of spatial variability of hydraulic conductivity, leakance and recharge, Water Resources Research, 25(3), 405-419. [27] Gómez-Hernández, J.J., and A.G. Journel, 1993. Joint sequential simulation of multiGaussian fields. In: Soares, A. (Ed.), Geostatistics Tria’s, vol.1, 85-94, Kluwer. [28] Gómez-Hernández, J.J., and X.-H. Wen, 1998. To be or not to be multi-Gaussian? A reflection on stochastic hydrology, Advances in Water Resources, 21(1), 47-61. [29] Gómez-Hernández, J.J., 2006. Complexity, Gound Water, 44(6), 782-785. [30] Gómez-Hernández, J.J., J. Fu, and D. Fernandez-Garcia, 2006. Upscaling retardation factors in 2-D porous media, Calibration and Reliability in Groundwater Modeling: From Uncertainty to Decision Making, Proceeding of ModelCARE’2005, The Hauge, The Netherlands, June 2005, IAHS Publ. 304, 130-136. [31] Hu, X., F.-W. Deng, and J.H. Cushman, 1995. Nonlocal reactive transport with physical and chemical heterogeneity: II. Linear nonequilibrium sorption with random Kd , Water Resources Research, 31(9), 2239-2252. [32] Journel, A.G., 1999. Conditioning geostatistical operations to nonlinear volume averages, Mathematical Geology, 31(8), 931-953. [33] Journel, A.G., C. Deutsch, and A.J. Desbarats, 1986. Power averaging for block effective permeability, SPE 15128. [34] Kabala, Z.J., and G. Sposito, 1991. A stochastic model of reactive solute transport with time-varying velocity in a heterogeneous aquifer, Water Resources Research, 27(3), 341350. 23

[35] Lemke, L.D., W.A. Barrack II, L.M. Abriola, and P. Goovaerts, 2004. Matching solute breakthrough with deterministic and stochastic aquifer models, Ground Water, 42(6), 920-934. [36] Leung, J.Y., and S. Srinivasan, 2011. Scale-up of mass transfer and recovery performance in heterogeneous reservoirs, Journal of Petroleum Science and Engineering, 86-87, 7186. [37] Liu, C., J.M. Zachara, N.P. Qafoku, and Z. Wang, 2008. Scale-dependent desorption of uranium from contaminated subsurface sediments. Water Resources Research, 44, W08413, doi:10.1029/2007WR006478. [38] Mackay, D.M., D.L. Freyberg, P.V. Roberts, and J.A. Cherry, 1986. A Natural Gradient Experiment on Solute Transport in a Sand Aquifer, 1. Approach and Overview of Plume Movement, Water Resources Research, 22(13), 2017-2029. [39] Metzger, D., H. Kinzelbach, and W. Kinzelbach, 1996. Effective dispersion of a solute cloud in a chemically heterogeneous porous medium: Comparison of two ensembleaveraging procedures, Water Resources Research, 32(11), 3311-3319. [40] Michalak, A.M., and P.K. Kitanidis, 2000. Macroscopic behavior and random-walk particle tracking of kinetically sorbing solutes, Water Resources Research, 36(8), 2133-2146. [41] Miralles-Wilhelm, F., and L.W. Gelhar, 1996. Stochastic analysis of sorption macrokinetics in heterogeneous aquifers, Water Resources Research, 32(6), 1541-1549. [42] Rajaram, H., 1997. Time and scale dependent effective retardation factors in heterogeneous aquifers, Advances in Water Resources, 20(4), 217-230. [43] Reichle, D., H. Kinzelbach, and W. Kinzelbach, 1998. Effective parameters in heterogeneous and homogeneous transport models with kinetic sorption, Water Resources Research, 32(11), 3311-3319. [44] Renard, Ph., and G. de Marsily, 1997. Calculating equivalent permeability: a review, Advances in Water Resources, 20(5-6), 253-278. [45] Roberts, P.C., M.N. Goltz, and D.M. Mackay, 1986. A natural gradient experiment on solute transport in a sand aquifer, 3. Retardation estimates and mass balances for organic solutes, Water Resources Research, 22(13), 2047-2058. [46] Robin, M.J.L., E.A. Sudicky, R.W. Gillham, and R.G. Kachanoski, 1991. Spatial variability of strontium distribution coefficients and their correlation with hydraulic conductivity in the Canadian Forces Base Borden aquifer, Water Resources Research, 27(10), 2619-2632.

24

[47] Sudicky, E.A., 1986. A natural gradient tracer experiment on solute transport in a sand and gravel aquifer: spatial variability of hydraulic conductivity and its role in the dispersion process, Water Resources Research, 22(13), 2069-2082. [48] Tompson, A.F.B., 1993. Numerical simulation of chemical migration in physically and chemically heterogeneous porous media, Water Resources Research, 29(11), 3709-3726. [49] Wagner, B.J., and S.M. Gorelick, 1989. Reliable aquifer remediation in the presence of spatial variable hydraulic conductivity: from data to design, Water Resources Research, 25(10), 2211-2225. [50] Wen, X.-H., and J.J. Gómez-Hernández, 1996a. The constant displacement scheme for tracking particles in heterogeneous aquifers, Ground Water, 34(1), 135-142. [51] Wen, X.-H., and J.J. Gómez-Hernández, 1996b. Upscaling hydraulic conductivities in heterogeneous media: An overview, Journal of Hydrology, 183, ix-xxxii.

25