- Email: [email protected]
Effective ADE models for first-order mobile–immobile solute transport: Limits on validity and modeling implications
Advances in Water Resources 86 (2015) 184–192
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Effective ADE models for first-order mobile–immobile solute transport: Limits on validity and modeling implications Scott K. Hansen∗ Department of Earth and Planetary Sciences, Weizmann Institute of Science, 234 Herzl St., Rehovot 76100, Israel
a r t i c l e
i n f o
Article history: Received 6 February 2015 Revised 16 September 2015 Accepted 18 September 2015 Available online 28 September 2015 Keywords: Green’s function Sorption Mobile–immobile Dispersion Anomalous transport
a b s t r a c t Quasi-1D mobile–immobile transport processes which have exponentially distributed random waiting times in both mobile and immobile states are common in hydrologic models (for example, of transport subject to kinetic sorption). The central limit theorem implies that eventually such transport will be expressible with an effective ADE (i.e. a generalization of the common retardation factor approach with an added Fickian dispersion coefficient accounting for the effect of trapping). Previous works have determined formulae for the value of this coefficient based on the transport properties. However, the time until convergence to Gaussian behavior has not previously been quantified. To this end, exact Green’s functions characterizing the transport at all times are derived for the case of pure advection. The Green’s functions are expressed in terms of three dimensionless parameters, representing location, time, and capacity coefficient. In the pre-Gaussian regime, a parametric study characterizing concentration profile asymmetry as a function of the capacity coefficient is performed. Next, heuristics are presented in terms of the dimensionless parameters for the time until the effective ADE adequately reflects reality. For strongly retarded solute, the time until effective ADE validity is found inversely proportional to release (e.g., desorption) rate. The nature of the effective dispersion coefficient is examined, and the possibility of large trapping-driven dispersion even in cases where batch experiments would detect negligible trapping is demonstrated. Collectively, these results call into question reliance on retardation factors derived from batch experiments for many practical transport modeling efforts; knowledge of both the trapping and release kinetics appears essential. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction 1.1. Mathematical treatment of mobile–immobile processes In many subsurface solute transport scenarios, there may be two spatially coextensive domains that are out of equilibrium with each other, each having its own local concentration. Commonly this is seen in models where solute advects only when it is in one, “mobile”, state (i.e. domain), with concentration c(x, t), but can also sometimes be trapped in an “immobile” state, which has its own concentration, cim (x, t), and from which it is eventually released. These models may be expressed with the following set of equations:
∂c ∂c (x, t ) + im (x, t ) = F {c}(x, t ) ∂t ∂t ∂ cim (x, t ) = G{c, cim }(x, t ), ∂t
(1)
∗ Present address: Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. Tel.: (505) 664 0712. E-mail address: [email protected], [email protected]
http://dx.doi.org/10.1016/j.advwatres.2015.09.011 0309-1708/© 2015 Elsevier Ltd. All rights reserved.
where F is a linear differential operator representing some combination of advection, dispersion, and decay, and G is an arbitrary operator. There are a number of different causes of this sort of mobile–immobile non-equilibrium transport: including so-called chemical and physical non-equilibrium, as well as diffusion into porous media [1]. Both the standard chemical and physical nonequilibrium equations, though arising through different conceptual pictures, may be put into an equivalent mathematical form [2], where in our notation we define
G{c, cim } ≡ λc − μcim .
(2)
Here λ represents the probability per unit time that a mobile particle will become immobile, and μ represents the probability per unit time that an immobile particle will become mobile. A variety of dualporosity transport problems can also be modeled in this way [3]. In addition, local equilibrium sorption with a distribution of retardation factors has been shown, when upscaled, to be expressible in the same form [4]. In our analysis, we shall employ an alternative expression for G (which we call G∗ ), previously used by
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
Margolin et al. [5]:
G∗ {c} ≡ λc − λ
t 0
ψim (τ )c(x, t − τ )dτ ,
(3)
which eliminates cim as a dependent variable but introduces time non-locality. As in (2), λ is a spatially homogeneous probability per unit time that a mobile particle will become immobile. Here, ψ im (t) is an arbitrary probability distribution for the length of a single sojourn in the immobile phase. Substituting this into (1) leads to the integrodifferential equation
∂c (x, t ) = F {c}(x, t ) − λc(x, t ) + λ ∂t
t 0
ψim (τ )c(x, t − τ )dτ .
(4)
Similar forms have been used by Schumer et al. [6] to define what the authors call the fractal mobile–immobile (fMIM) paradigm, by Haggerty et al. [7] in the analysis of late-time BTC tails in the presence of trapping processes, and by Benson and Meerschaert [8] in their development of subordination technique for incorporating anomalous transport into ADE-based numerical models. In the Laplace domain, it is not difficult to show that Eqs. (1) and (2) are equivalent to (4) so long as
ψim (t ) = μe−μt .
(5)
While more general forms of ψ im than (5) correspond to some physical systems, we restrict our attention to first-order mobile–immobile behavior, and so will use (4) and (5) as the basis for all analysis. Note that this model implicitly defines a comparable probability distribution,
ψm (t ) = λe−λt ,
(6)
for the time taken by a single sojourn in the mobile phase (i.e. time interval between release and subsequent capture). This particular model, with exponentially-distributed times in both the sorbed and free states, has found wide application in the literature: applications are found in the sediment transport literature [9,10], as well as in analysis of solute trapping [8,11]. 1.2. Trapping-driven dispersion, breakthrough curve asymmetry, and effective ADE models Commonly, transport for each solute particle, including the random trapping, is independent of all other particles. Where this trapping is a result of diffusion into an immobile zone, this is assured. For kinetic sorption, where there exists a great excess of sorption sites, it is justified as well, although in cases of competitive sorption (i.e. limited sites), this may not be the case. Any random transport process that acts independently on different particles must have a dispersive effect, and indeed: the dispersive effect of kinetic sorption has long been recognized in the chromatography literature, dating back at least to [12]. In the subsurface transport literature, some theoretical studies quantifying the dispersive effects of mobile– immobile systems have been published. Late-time first, second, and third temporal moments of breakthrough curves for advectivedispersive subsurface solute transport subject to first-order chemical and physical non-equilibrium were considered by Valocchi [2]. A similar moment analysis, except deriving spatial moments was later performed by Michalak and Kitanidis [13]. The authors further computed effective late-time velocities and dispersion coefficients. These studies were in 1D, and focused on late time limits in the case when (2) holds, and are in the school of classic mobile–immobile zone models with exponential ψ im . Recently, Uffink et al. [11] employed random walk theory to derive an equivalent PDE to (4), and manipulated the equation to derive the same late-time dispersion coefficient shown in [13]. To our knowledge, they were the first to remark that λ and μ may affect the rate of late-time Gaussian convergence, but did not pursue this systematically, save for providing an approximate time until the effective dispersion coefficient stabilizes.
185
Related to the discussion of trapping-driven dispersion is a body of literature on the so-called local equilibrium assumption (LEA). The LEA essentially refers to conditions under which the trapped and free concentrations are related by the same ratio at all times [14]. In our terminology, we express this as c(x, t ) + cim (x, t ) ≈ Rc(x, t ), for some constant, R, referred to as a retardation factor. By substitution of R into (1), one arrives at the retarded transport equation:
R
∂c (x, t ) = F {c}(x, t ). ∂t
(7)
It is apparent that this form simply rescales time (R can be eliminated by a substitution τ ≡ t/R), and so the dispersive effect of trapping is not considered. Nonetheless, some authors conflate applicability of the local equilibrium assumption with usage of a retarded transport equation (e.g. [15]). There are some circumstances under which this is reasonable, namely those in which F is an advection–dispersion operator representing a stronger source of dispersion, compared to which trapping-driven dispersion is negligible. The LEA has been widely studied, with aforementioned paper [2] aiming to validate it by comparing moments, at late time, of exact solutions for chemical and physical non-equilibrium with a retarded advection–dispersion equation. A different approach, seeking directly to find conditions under which local equilibrium is nearly satisfied at a point and relating that to parameters in a system like (1) has also been presented [16]. Later, Wallach examined the domain of validity of the LEA through perturbation theory, treating (7) as an end member in a perturbation expansion for the exact solution which contains trapping-driven dispersion [17]. All of these authors, and many others (see the literature review in [16]) identified conditions under which the dispersive effect of sorption was comparatively negligible. However, since Michalak and Kitanidis [13] computed the effective, late time dispersion coefficient for classic mobile–immobile systems, there is now a more flexible framework available than the dispersion-ignoring LEA approximation. Mobile–immobile solute transport may be treated at late time by a retardation coefficient and this effective dispersion coefficient, even in cases where trappingdriven dispersion cannot be ignored. We term such an approach an effective ADE model. An effective ADE model still cannot apply in circumstances in which there is a significantly asymmetric transport Green’s function, however, since the Green’s function for the ADE is Gaussian. Further, it is well known that in many trapping-and-release processes, plume asymmetry occurs. Uffink et al. [11] showed how this develops in a system described by (1) and (2). Since consideration of mobile–immobile problems is so common in contaminant hydrogeology, there are naturally many alternatives to the effective ADE approach, including exact solutions for specific geometries. A comprehensive survey is beyond the scope of this work, but classes of techniques include analytic solutions such as those included in the aforementioned van Genuchten paper [3], analytic solutions that treat discrete fractures [18,19], and by multi-rate mass transfer schemes [20]. Numerical approaches that treat multiple interacting continua [21,22] are also useful for specific problems. Lastly, it bears noting that in addition to being indicative of earlytime mobile–immobile behavior for exponential ψ im , asymmetry may be caused by power-law distributed immobile times [23], as well as non-sorbing transport in heterogeneous media [24]. The difficulty of distinguishing between mobile–immobile models and transport in heterogeneous media has been remarked upon by Carrera et al. [25], among others. This implies that CTRW approaches [24] represent another viable approach for modeling mobile–immobile systems. Experimentally, at the column scale, asymmetric breakthrough curves have been examined using CTRW approaches by Deng [26] and Li and Ren [27]. This sort of behavior may arise due to anomalous transport
186
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
or mobile–immobile trapping behavior. The tailing effect of mobile– immobile trapping in columns where transport in the mobile state is anomalous has also been considered experimentally by Rubin et al. [28]. 1.3. Goals for this work
∀x > 0 x=0 ∀t
(8)
to generate Green’s functions for concentration profiles and breakthrough curves. Because of this explicit approach, we are able to make a second contribution, which is to characterize the asymmetry of pre-Gaussian behavior as a function of the governing parameters. The Green’s functions for concentration profiles and breakthrough curves are expressed in terms of three dimensionless parameters: a dimensionless location, dimensionless time, and the capacity coefficient, β = λ/μ. A parametric study is performed characterizing the early-time behavior and asymmetry of concentration profiles and breakthrough curves. The time until an effective ADE model is appropriate is quantified in terms of the dimensionless parameters and a rule, independent of β , for the dimensionless time until effective ADE applicability is formulated (implying that, for strongly retarded substances, the dimensional time is inversely proportional to μ; see below). Finally, the unpredictability of the degree of sorption-driven dispersion given β alone is observed. As a consequence of these results, we conclude that one cannot safely use a retarded ADE model (with a retardation factor computed from an equilibrium batch experiment) for sorbing solute transport without characterizing both λ and μ independently. This warning runs contrary to some current practice, and is a practical contribution complimentary to the new theoretical results. 2. Green’s functions for transport In the Laplace domain, (4) can be solved alongside (8) to yield:
cˆ(x, s)[s + λ − λψˆ im (s)] = F {cˆ}(x, s).
(9)
Since we are interested in first-order mass transfer, we define ψ according to (5). We also choose F ≡ −v∂ /∂ x (the advection operator), so that all spreading can be attributed to the mobile–immobile trapping process. Explicit Green’s function solutions for breakthrough curves, b(x, t), and concentration profiles, p(x, t), are developed in dimensional form in Appendix A. From these, one can compute c(x, t ) = t t 0 b(x, t − τ )c(0, τ )d τ , or c(x, t ) = 0 p(x − y, t )c(y, 0)dy. In the following, we will present them in terms of the following dimensionless groups:
t ≡
μt, μx/v, β ≡ λ/μ.
x ≡
2.1. Green’s functions for breakthrough curves Substituting (10) into (21) yields
In light of the above discussion, we continue the explorations begun by Michalak and Kitanidis [13] and Uffink et al. [11], and characterize quantitatively the time until an effective ADE behavior is seen as a function of transport parameters. We perform the analysis explicitly, directly solving (4) subject to the initial and boundary conditions:
c(x, 0) = 0 c(0, t ) = δ(t ) c(∞, t ) = 0
remain at equilibrium, the instantaneous particle trapping rate and particle release rate must be equal.)
(10)
The first two quantities represent, respectively, a dimensionless location and time. The third, β , is sometimes called the capacity factor. It can be thought of as representing the ratio cim /c in a system that is in local equilibrium, or equivalently as the ratio of expected length of one continuous immobile period to the expected length of one continuous mobile period. (The two are equivalent because in order to
b(x , t , β) =
μ
exp{−t + (1 − β)x }
+
β x
t − x
I {2 1
δ(t − x )
β x (t − x )} .
(11)
The terms in the outer square brackets are dimensionless. It may seem counter-intuitive that b, which may be thought of in terms of “concentration” has units of t −1 , but it is important to remember that the initial condition was defined as a temporal Dirac delta function, which also has the units of t −1 , and that b is really a Green’s function which must be integrated against a source history to produce a physical breakthrough curve. 2.2. Green’s functions for concentration profiles We nondimensionalize again, employing (10) and (26). The Green’s function for the concentration profile is
p(x , t , β) =
μ exp{−t + (1 − β)x } δ(t − x ) v β x (t − x )} + I { 2 β x , 1 t −x
(12)
which has the same mathematical form as the breakthrough curve (11), up to a constant factor. It is clear that this must be so, since either Green’s function calculation must essentially represent the plume due to an instantaneous point source, and so the distinct calculations for (11) and (12) corroborate each other. 3. Pre-asymptotic behavior: the effect of time and β Given an ultimate interest in how long it takes for transport to approach a regime in which the ADE is applicable, it is reasonable to consider the concentration profile Green’s functions and consider how large t must become before they are Gaussian. At the same time, how t and β affect their pre-asymptotic evolution and behavior also bears consideration. 3.1. Effect of changing β The capacity factor, β , relates to the effective retardation factor as R = 1 + β , which can be seen by setting G{c, cim } = 0 in (2), to represent equilibrium conditions. Consequently, it is predictable that in Fig. 1, larger values of β correspond to smaller first moments of the Green’s function curves at any fixed time: i.e. more time trapped means the plume travels a shorter distance. Although the distributions become Gaussian eventually for all values of β , their preasymptotic character depends on whether β is greater or less than unity. For small β , there is a heavy tail to the distribution for small x , and for large β the converse is true. Note that the integrals of Green’s functions tend to decrease with increasing β and t . This comes from the fact that the Green’s functions are for solute that is initially all mobile at time 0, and mass transfer between the trapped and free states (each respectively integrated over all space) asymptotically approaches a steady state in which 1/(1 + β) of the initial mass remains free.
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
187
t ' = 10
t '= 1
.1
1
.5
3 100
10 .3
.3 .5 1
100
3
10
.1
.01
.01
t ' = 50
.01
.1
100 10
3
1
.5
.3
Fig. 1. Plots of exact solution Green’s functions ((11) with prefactor μ = 1) for various fixed values of t , as β increases. Note that the various x -axes do not have the same scales. The vertical arrows in each plot are indicative of Dirac distribution representing probability mass for particles that have never been trapped. The heights of the arrows are not significant.
3.2. Effect of changing t The effect of t , the dimensionless time, is illustrated in Fig. 2. Essentially, each of the subplots in that figure can be viewed as a set of superimposed stop-motion snapshots taken at various t . Nonsorbing solute, which travels at the same rate as the underlying fluid flow, will always have traveled to the location t = x (the location of the Dirac spikes). As solute sorbs and then subsequently desorbs it is spread out behind the non-sorbing front. This leads to the progressive depletion of probability mass from the Dirac spikes into the rest of the probability distribution with increasing time. For larger β , more
depletion from the front appears for a fixed t , although actually only λ determines the degree of depletion by a given dimensional time t. 4. Asymptotic behavior: how long until the effective ADE applies? We saw previously that dispersion is unavoidable when dealing with solute transport undergoing mobile–immobile transport. It can also be shown that the transport Green’s functions (having finite statistical moments) will eventually approach Gaussian distributions. We now turn to the question of how to capture transport with an effective ADE.
188
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
(height 0.68)
=1
= .5
t' t'
t'
1
1
3
t'
5
t' t'
3
7
t' t'
5
t'
9
7
t'
=5
=2
t'
9
1
t'
t'
1
3
t' t'
5
t'
3
t'
7
t'
5
t'
9
7
t'
9
Fig. 2. Plots dimensionless Green’s functions ((11) with prefactor μ = 1) with for various fixed values of β as t increases. Vertical axis represents both probability density (for the dotted curves) and total probability mass for the Dirac spikes that are part of the exact solution. The height of each of the solid vertical arrows represents the corresponding amount probability mass each Dirac spike contains—this represents the mass that has never been trapped since being released at t = 0.
4.1. The effective ADE We wish to rewrite (4), where F ≡ −v∂ /∂ x, in the following form
∂c ∂c ∂ 2c (x , t ) = −v∗ (λ, μ) (x , t ) + D∗ (λ, μ) 2 (x, t ). ∂t ∂x ∂x
(13)
The effective v∗ will equal the advection velocity, v, multiplied by the fraction of time that a particle is free (over a long time period, this fraction converges for all particles, even though the absolute time immobile for all particles does not converge, leading to the trappingdriven dispersion that we focus on in this paper and wish to capture with D∗ ). A long time interval can be treated as consisting of an essen-
tially equal number of mobile and immobile periods, and so the retardation factor, R, can be expressed as the ratio between the expected length of a mobile–immobile cycle and the expected length of a single sojourn in the mobile phase: R = E[ψm + ψim ]/E[ψm ]. Since the distributions ψm and ψim are both exponential with parameters λ and μ, respectively, it follows that R = 1 + λ/μ. (We see here that R has a dual interpretation. It is commonly understood as the ratio between peak arrival times for a conservative vs. a trapped-and-released solute. However, on the analysis seen here, it can also be understood as the local ratio between total concentration, c + cim , and c.) Thus,
−1 λ v (λ, μ) = = v 1 + . R μ ∗
v
(14)
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
189
Fig. 3. Comparison of exact and approximate solutions to the example problem at t = 5 (left) and t = 90 (right; note different scale). Solid curves and the Dirac arrow represent the exact solution (12), and dashed curves represent the ADE approximation using an effective velocity (14) and effective dispersion coefficient, D∗ (15). The vertical dotted line, as elsewhere, indicates a discontinuity in exact solution.
The formula for effective dispersion coefficients first-order mobile– immobile transport was computed by Michalak and Kitanidis [13] (later re-derived by Uffink et al. [11], using different means). From their work, it follows that
D∗ (λ, μ) =
β 1 λμv2 = v2 . 3 μ (λ + μ) (1 + β)3
(15)
In this paper, we only analyze a system with no hydrodynamic dispersion. However, note that [11] also reported a total effective dispersion coefficient in the case of trapping and hydrodynamic dispersion when free. 4.2. Time to convergence The effective ADE model can only match true behavior after a Gaussian character has developed in the concentration profile Green’s function. For illustrative purposes, we will compare solutions of (4) and (13) for the case when λ = 0.6, μ = 1, and v = 1 (which ensure that the dimensionless and dimensional coordinates are numerically identical) and β = 0.6. Solutions at early (t = 5) and late (t = 90) times are shown in Fig. 3, illustrating the convergence of the spatial Green’s functions to Gaussian normality as the solute undergoes multiple sorption–desorption cycles. At early time, the ADE approximation under-predicts peak concentration and significantly lags the exact solution. Furthermore, the approximation predicts significant probability density for values of x less than zero and greater than five, both of which are impossible. At late time, however, the coherence between the two solutions is good, suggesting that trapping-driven retardation and dispersion coefficients may be used in the ADE. To confidently use an effective dispersion coefficient in lieu of the exact analysis, we may quantify the degree of coherence between the concentration profiles generated by the exact and the effective ADE approaches as a function of β and t , so as to avoid using the effective approach in a case where they are not sufficiently coherent. To quantify this, we define thenormalized divergence as the L1 norm of the difference between the exact and approximate solutions divided
by the L1 norm of the exact solution, and evaluate this for a number of different values of β and t . This quantity is bounded between zero (curves totally coincident) and two (curves totally disjoint), provided both curves have the same L1 norm. For concreteness, the normalized divergence in Fig. 3 is 0.51 when t = 5, and 0.09 when t = 90; this latter quantity may seem surprisingly high. We adopt a normalized divergence of 0.1 as a heuristic for applicability of an ADE model. In order to understand the regime of applicability, we performed a numerical parametric study. Since we are interested in comparing concentration profiles which are functions of x , there are only two free dimensionless parameters available, t and β . Thus, we may produce contour plots of normalized divergence as a function of these two parameters. These are shown in Fig. 4. From here, we see that for values of β greater than unity and normalized divergence less than 0.5, the contours are more or less vertical, and independent of further increases in β . From this, and the satisfactory performance seen in Fig. 3. We may develop a heuristic that it is reasonable to employ the effective diffusion coefficient approximation when t > 70 and β > 1, since this implies a normalized divergence less than 0.1. That the 0.1 contour continues vertically off the top of the graph for additional orders of magnitude was verified by generating additional plots, not shown. The fact that as β grows, the dimensionless time until the normalized divergence threshold is reached is unchanged, along with (10), implies that the actual time until this occurs varies inversely with μ. When β becomes smaller, Fig. 4 implies a longer dimensionless time, t , is required before the effective ADE model can be employed. This can seem surprising, as small values of β correspond to retardation factors near unity. However, consider the exponential decay parameters λ and that characterize ψ im (5) and ψm (6). It follows from (10) that as μ grows, it will both decrease β in inverse proportion, and increase t proportionally, and so there is no particular reason to expect an increased μ to correspond to a pathologically increased dimensional time, t, required to achieve nearly Gaussian behavior. However, this is exactly what will happen as λ decreases. This is reasonable, since particles must be sorbed at least once in order to move away from the advection front, and must sorb and desorb numerous
190
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
Fig. 4. Semi-log contour plot of normalized divergence for various values of β and t .
Fig. 5. Semi-log contour plot of normalized divergence for various values of β and β t .
times before approaching the quasi-Gaussian behavior that, according to the central limit theorem, they must, in the limit. As λ declines, the mean time between trapping events increases, and so one must wait for a longer time before quasi-Gaussian behavior manifests itself. For large β the time for a cycle to complete is controlled by μ, and for small β it is controlled by λ. If we define tc as the (essentially β -independent) dimensionless time required for convergence in the case of large β , and tc as its dimensional counterpart, then
implies β of a similar magnitude. In this case, it is possible to determine from Eq. (15) that D∗ ≈ v2 μ/λ2 . While β = λ/μ is large, implying that μ/λ 1, it is still possible that μ/λ2 1. For concreteness, consider an example where v = 1, λ = 10−4 , and μ = 10−6 , implying β = 100. These values of λ and μ represent a system that is very slow to reach equilibrium, which is sometimes observed experimentally [26]. In this case D∗ = 97.06. However, multiplying both λ and μ by 104 keeps β = 100, but yields D∗ = 9.706 × 10−3 , so D∗ varies by a factor of 10,000 while corresponding to the same β . Thus, in the case of large β (i.e. large R), a knowledge of the kinetics is required in order to realistically model the system: knowledge of only the equilibrium value capacity factor is not enough.
tc =
tc . E[ψim ]
(16)
This is to say, tc is the approximate number of trap-and-release cycles for convergence to Gaussian behavior (where the brackets “ ” represent rounding up to the nearest greater integer). For small β , It is easier to visualize the behavior by plotting against an alternative dimensionless time, β t (= λt ) instead of μt. This is shown in Fig. 5. It is apparent that the alternative dimensionless time until validity of the dispersion coefficient approach is essentially independent of the value of β for small β , although the threshold of ADE applicability is quite different. From examination of the 0.1 contour, a heuristic β t > 30 and β < 0.3 is developed for applicability of ADE behavior (or β t > 70 and β ≤ 1, to cover all possibilities). Note for comparison that heuristics are presented in [11] for disappearance of the Dirac spikes (rather than applicability of an effective ADE model): which were that t > 3/λ and t > 3/μ. Placing our heuristics in dimensional form, they amount to t > 30/λ and t > 70/μ. These are thus more conservative than those presented previously, which is reasonable since they correspond to more stringent conditions on the plume. 5. The retardation factor is uninformative about trapping-driven dispersion We now discuss the applicability of the effective ADE approach in the context of three special cases in which it is possible to simplify (15). We will see in all three cases that β (and by extension, R) is uninformative about D∗ , which may be large or small depending on the absolute magnitudes of λ and μ. Special case I: β is large. For some heavy organic molecules sorbing in certain environments, R can be greater than 1000 [29]. This scenario
Special case II: β is small. Interestingly, an exactly analogous ambiguity to that shown for large β exists when β → 0; this is to say for R ≈ 1 (using the relation R = 1 + λ/μ = 1 + β ). It follows from (15) that for small β , D∗ ≈ v2 λ/μ2 . Then by permuting λ and μ in the example shown in special case I, exactly the same conclusions follow here. This means the possibility exists for essentially arbitrarily large trappingdriven dispersion for slowly trapped and released solute even when R ≈ 1. Thus, kinetics must be considered in a case when the solute’s trapping affinity is so small that a batch experiment would not predict it to be a factor in transport. 6. Summary and conclusions We considered trapping-driven dispersion in classic, first-order mobile–immobile solute transport systems under pure advection (hydrodynamic dispersion being implicitly taken to be negligible compared to trapping-driven dispersion). A certain time-non-local transport Eq. (4) was used in the analysis; it proved possible to explicitly solve this transport equation using integral transforms, and to write explicitly the Green’s functions for both concentration profiles and breakthrough curves in terms of just three dimensionless parameters. Exact solution for the Green’s functions allowed explicit characterization of behavior over all times and quantification of the time until an effective ADE model with a retardation factor and trappingdriven dispersion coefficient is applicable. In particular, we established three concrete results regarding, respectively, the pre-Gaussian regime, the time to convergence, and the dispersion in the Gaussian regime:
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192 •
•
The value of the capacity factor, β , was seen to qualitatively affect the pre-Gaussian behavior of both concentration profiles and breakthrough curves in a simple way. Values of β < 1 were associated with tailing opposite to the advection direction and late breakthrough. Values of β > 1, contrarily, were associated with tailing in the advection direction and early breakthrough. Heuristics were found for reasonable usage of an effective ADE λ and the D∗ = λμv2 of [11]: λt > 30 for λ model with R = 1 + μ 3 (λ+μ)
•
μ and μt > 70 for μ ≤ λ, with a conservative criterion being t > 70/min(λ, μ). A corollary of this analysis is that for strongly retarded solute, the time until applicability of the effective ADE is inversely proportional to the desorption rate, μ. Similarly, the value of R in the effective ADE model was found uninformative about the strength of sorption-driven dispersion, and the same value may correspond to large or small D∗ , for different combinations of λ and μ. Interestingly, itwas seen possible to have substantialtrapping-driven dispersion with negligible trapping detectable via a batch experiment (i.e. for R ≈ 1).
One might argue that the last two results are of theoretical import only, since if one were to be in possession of λ and μ, it would be easy enough to use the exact solutions (11) or (12), and one would not need an effective ADE at all. However, we believe they provide a warning about the modeling practice of characterizing dispersion separately from sorption, and then simply placing a retardation factor derived from equilibrium arguments into the ADE. Because an arbitrarily large degree of sorption-driven dispersion may be present without impacting R or the conservative-tracer D, and the time until any Gaussian model is applicable cannot be determined without knowledge of λ and μ, there is no substitute for determining them. A possible method for determination of λ and μ is suggested by our exact Green’s function solution, (11). Fitting this to a breakthrough curve from a short column experiment, and simultaneously determining β and pre-multiplication factor μ from the skewness and total recovery, respectively, is a potential means of accomplishing this. (The presence of pore-scale hydrodynamic dispersion, which we have treated as negligible in the above analysis, would naturally be a complicating factor.) Regardless of the viability of this approach, the characterization of both the capture and the release processes separately is suggested as a potentially important topic for future investigation. A concrete example of the significance of these results lies in the proper modeling of transport of strongly sorbing organic pollutants: Hexachlorinated PCB congeners are known to have retardation factors (and thus β values) on the order of 106 [29], indicating times to ADE applicability that are μ-dependent. Very slow desorption rates have also been observed for these congeners, with mean desorption times on the order of weeks [30], implying correspondingly small μ values. From (15) and (16), respectively, this implies very large D∗ and tc values. It is thus apparent that in similar cases, the mere substitution of R directly into the ADE applicable to untrapped solute is apt to lead to seriously misleading results. Instead, it appears only reasonable to employ (4) with appropriately calibrated λ and μ.
191
A1. Green’s function for transient release To determine Green’s functions for breakthrough curves (i.e. kernels b(x, t) that are convolved with c(0, t) to determine c(x, t)), we proceed directly from (9), writing
−v
λμ ˆ ∂ bˆ b = s+λ− ∂x s+μ = (s + μ) + (λ − μ) −
λμ ˆ b. (s + μ)
(17)
This can be solved according to the initial boundary conditions (8) to yield
x bˆ = exp −
v
(s + μ) + (λ − μ) −
Applying the shift theorem,
b(x, t ) = e
−μt −(λ−μ)x/v
e
L
−1
λμ s ( + μ)
exp −
x
s−
v
.
λμ
(18)
.
s
(19)
Defining the quantities a ≡ x/v and d ≡ λμx/v, the inverse Laplace transform in (19) can be rewritten as
L−1 {e−as (ed/s − 1) +e−as } = L−1 e−as ∗ L−1 {ed/s −1} + L−1 {e−as }
=
δ(t − a) +
d I1 {2 d(t − a)}, t −a (20)
which was resolved using table look-up (see [31], p. 255, #5.66) where I1 is the modified Bessel function of the first kind. Successively substituting the definitions of a and d into (20) and into (19) yields a complete time-domain solution:
b(x, t ) = e−μt e−(λ−μ)x/v
× 2
δ(t − x/v) +
λμx (t − x/v) , v
λμ tv x
−1
I1
(21)
This equation is of the same form as that derived by Giddings and Eyring [12] for their probabilistic chromatography equations. An equivalent equation is also derived for the same scenario, using a different conceptual approach, in [8]. In that paper, the authors are concerned with converting between a given moment in time, t, and its corresponding operational time, u (in our context, the amount of time the particle is free). By identifying the operational time with x/v in (26)—natural, since this is a pure advection system—we recover directly their solution. A2. Green’s function for specified initial concentration To develop the Green’s function, p, for the concentration profile, we again solve the same PDE, but require different initial-boundary conditions. Instead of (8), we employ
Acknowledgment
p(x, 0) = δ(x) c(−∞, t ) = 0 c(∞, t ) = 0
I thank the Azrieli Foundation for a postdoctoral fellowship which in part made this work possible.
We then proceed directly from (9), as before, taking the Laplace transform:
Appendix A. Derivation of Green’s functions The derivation of the Green’s functions in dimensional coordinates for both concentration profiles and breakthrough curves in the case when ψ(t ) = μe−μt is presented.
δ(x) − v
t=0 ∀t ∀t
λμ ∂ pˆ ˆ = s+λ− p. ∂x s+μ
(22)
(23)
We then take the Fourier transform from x ω:
pˆ =
1 . λμ s + λ − s+ μ + viω
(24)
192
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
Dividing top and bottom by v, this can be solved by table look-up as
pˆ =
1
v
exp −
x
v
(s + μ) + (λ − μ) −
λμ (s + μ)
H (x),
(25)
where H is the Heaviside step function, which is unity on the interval [0, ∞). Apart from the initial factor of 1/v, then, on the interval of interest this is exactly (18), which has already been solved. Thus, we may conclude the exact solution for the concentration profile Green’s function is
p(x, t ) =
μ −μt −(λ−μ)x/v λμ I1 e e δ(t − x/v) + t v v −1 x
λμx × 2 (t − x/v) . v
(26)
References [1] Weber WJ Jr, Pennell KD, Dekker TJ, Abriola LM. Sorption and retardation of organic contaminants in subsurface systems: effects on transport and fate. In: Aral MM, editor. Advances in groundwater pollution control and remediation. No. 9. NATO ASI series. Netherlands: Springer; 1996. p. 1–31. ISBN 978-94-0106576-4 978-94-009-0205-3. [2] Valocchi AJ. Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils. Water Resour Res 1985;21(6):808–20. http://dx.doi.org/10.1029/WR021i006p00808. http://onlinelibrary.wiley.com/doi/10.1029/WR021i006p00808/abstract. [3] van Genuchten M. A general approach for modeling solute transport in structured soils. In: Proceedings of the 17th international Congress on the hydrogeology of rocks of low permeability, vol. 17, no. (2). In: Mem Int Assoc Hydrogeol. Tucson, AZ; 1985. p. 513–26. [4] Dentz M, Berkowitz B. Exact effective transport dynamics in a one-dimensional random environment. Phys Rev E 2005;72(3):031110. http://dx.doi.org/10.1103/ PhysRevE.72.031110. http://link.aps.org/doi/10.1103/PhysRevE.72.031110. [5] Margolin G, Dentz M, Berkowitz B. Continuous time random walk and multirate mass transfer modeling of sorption. Chem Phys 2003;295(1):71–80. http://dx.doi.org/10.1016/j.chemphys.2003.08.007. http://www.sciencedirect. com/science/article/pii/S0301010403004245. [6] Schumer R, Benson DA, Meerschaert MM, Baeumer B. Fractal mobile/immobile solute transport. Water Resour Res 2003;39(10):1296. http://dx.doi.org/10.1029/2003WR002141. http://onlinelibrary.wiley.com/doi/ 10.1029/2003WR002141/abstract. [7] Haggerty R, McKenna SA, Meigs LC. On the late-time behavior of tracer test breakthrough curves. Water Resour Res 2000;36(12):3467–79. http:// dx.doi.org/10.1029/2000WR900214. http://onlinelibrary.wiley.com/doi/10.1029/ 2000WR900214/abstract. [8] Benson DA, Meerschaert MM. A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations. Adv Water Res 2009;32(4):532–9. http://dx.doi.org/10.1016/j.advwatres.2009.01.002. http: //www.sciencedirect.com/science/article/pii/S0309170809000128. [9] Barry DA, Sander GC, Jomaa S, Heng BCP, Parlange JY, Lisle IG, et al. Exact solutions of the Hairsine–Rose precipitation-driven erosion model for a uniform grain-sized soil. J Hydrol 2010;389(3–4):399–405. http://dx.doi.org/ 10.1016/j.jhydrol.2010.06.016. http://www.sciencedirect.com/science/article/pii/ S0022169410003513. [10] Lisle IG, Rose CW, Hogarth WL, Hairsine PB, Sander GC, Parlange JY. Stochastic sediment transport in soil erosion. J Hydrol 1998;204(1–4):217–30. http://dx.doi.org/10.1016/S0022-1694(97)00123-6. http://www.sciencedirect. com/science/article/pii/S0022169497001236. [11] Uffink G, Elfeki A, Dekking M, Bruining J, Kraaikamp C. Understanding the nonGaussian nature of linear reactive solute transport in 1D and 2D. Transp Porous Media 2012;91(2):547–71. http://dx.doi.org/10.1007/s11242-011-9859-x. http:// link.springer.com/article/10.1007/s11242-011-9859-x. [12] Giddings JC, Eyring H. A molecular dynamic theory of chromatography. J Phys Chem 1955;59(5):416–21. http://dx.doi.org/10.1021/j150527a009. http://dx.doi. org/10.1021/j150527a009. [13] Michalak AM, Kitanidis PK. Macroscopic behavior and random-walk particle tracking of kinetically sorbing solutes. Water Resour Res 2000;36(8):2133–46. http://dx.doi.org/10.1029/2000WR900109. http://onlinelibrary.wiley.com/doi/ 10.1029/2000WR900109/abstract.
[14] Ptacek C, Gillham RW. Laboratory and field measurements of non-equilibrium transport in the Borden aquifer, Ontario, Canada. J Contam Hydrol 1992;10:119– 58. [15] Lee L, Rao P, Brusseau ML, Ogwada RA. Nonequilibrium sorption of organic contaminants during flow through columns of aquifer materials. Environ Toxicol Chem 1988;7:779–93. http://dx.doi.org/10.1002/etc.5620071001. http:// onlinelibrary.wiley.com/doi/10.1002/etc.5620071001/abstract. [16] Bahr JM, Rubin J. Direct comparison of kinetic and local equilibrium formulations for solute transport affected by surface reactions. Water Resour Res 1987;23(3):438–52. http://dx.doi.org/10.1029/WR023i003p00438. http:// onlinelibrary.wiley.com/doi/10.1029/WR023i003p00438/abstract. [17] Wallach R. A small perturbations solution for nonequilibrium chemical transport through soils with relatively high desorption rate. Water Resour Res 1998;34(1):149–54. http://dx.doi.org/10.1029/97WR02490. http://onlinelibrary. wiley.com/doi/10.1029/97WR02490/abstract. [18] Feenstra S, Cherry JA, Sudicky EA, Haq Z. Matrix diffusion effects on contaminant migration from an injection well in fractured sandstone. Ground Water 1984;22(3):307–16. http://dx.doi.org/10.1111/j.1745-6584.1984.tb01403.x. http: //onlinelibrary.wiley.com/doi/10.1111/j.1745-6584.1984.tb01403.x/abstract. [19] Tang DH, Frind EO, Sudicky EA. Contaminant transport in fractured porous media: analytical solution for a single fracture. Water Resour Res 1981;17(3):555– 64. http://dx.doi.org/10.1029/WR017i003p00555. http://onlinelibrary.wiley.com/ doi/10.1029/WR017i003p00555/abstract. [20] Haggerty R, Gorelick SM. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour Res 1995;31(10):2383–400. http://dx.doi.org/10.1029/95WR10583. http:// onlinelibrary.wiley.com/doi/10.1029/95WR10583/abstract. [21] Silva O., Carrera J. mod_process_mrmt.f90: a Fortran module to include multi-rate mass transfer like processes: user’s guide. 2008. http://www.h2ogeo.upc.edu/English/software/MRMT%20module%20in% 20GroupWeb/MRMT%20module%20in%20GroupWeb/User’sManual.rar. [22] Pruess K, Oldenberg C, Moridis G. TOUGH2 users guide, version 2. Tech. rep. LPNL-43134. Earth Sciences Division, Lawrence Berkeley National Laboratory; 2012. http://esd1.lbl.gov/files/research/projects/tough/documentation/TOUGH2_ V2_Users_Guide.pdf. [23] Haggerty R, Wondzell SM, Johnson MA. Power-law residence time distribution in the hyporheic zone of a 2nd-order mountain stream. Geophys Res Lett 2002;29(13):18-1–4. http://dx.doi.org/10.1029/2002GL014743. http:// onlinelibrary.wiley.com/doi/10.1029/2002GL014743/abstract. [24] Berkowitz B, Cortis A, Dentz M, Scher H. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev Geophys 2006;44(2):RG2003. http://dx.doi.org/10.1029/2005RG000178. http:// onlinelibrary.wiley.com/doi/10.1029/2005RG000178/abstract. [25] Carrera J, Sánchez-Vila X, Benet I, Medina A, Galarza G, Guimerà J. On matrix diffusion: formulations, solution methods and qualitative effects. Hydrogeol J 1998;6(1):178–90. http://dx.doi.org/10.1007/s100400050143. http://link. springer.com/article/10.1007/s100400050143. [26] Deng J, Jiang X, Zhang X, Hu W, Crawford JW. Continuous time random walk model better describes the tailing of atrazine transport in soil. Chemosphere 2008;71(11):2150–7. http://dx.doi.org/10.1016/j.chemosphere.2008.01.001. http: //www.sciencedirect.com/science/article/pii/S0045653508000039. [27] Li N, Ren L. Application of continuous time random walk theory to nonequilibrium transport in soil. J Contam Hydrol 2009;108(3–4):134–51. http://dx.doi.org/10.1016/j.jconhyd.2009.07.002. http://www.sciencedirect.com/ science/article/pii/S0169772209000904. [28] Rubin S, Dror I, Berkowitz B. Experimental and modeling analysis of coupled nonFickian transport and sorption in natural soils. J Contam Hydrol 2012;132:28–36. http://dx.doi.org/10.1016/j.jconhyd.2012.02.005. http://www.sciencedirect.com/ science/article/pii/S0169772212000265. [29] Hansen SK, Kueper BH. Forensic misidentification of Aroclor sources in fractured bedrock due to “chromatographic” polychlorinated biphenyl (PCB) congener separation. Environ Forensics 2011;12(1):25–34. http://dx.doi.org/ 10.1080/15275922.2010.547547. http://www.tandfonline.com/doi/abs/10.1080/ 15275922.2010.547547. [30] Ghosh U, Weber AS, Jensen JN, Smith JR. Relationship between PCB desorption equilibrium, kinetics, and availability during land biotreatment. Environ Sci Technol 2000;34(12):2542–8. http://dx.doi.org/10.1021/es9905389. http://pubs.acs. org/doi/abs/10.1021/es9905389. [31] Oberhettinger F, Badii L. Tables of Laplace transforms. Berlin: Springer-Verlag; 1973. ISBN 978-3-642-65645-3.
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Effective ADE models for first-order mobile–immobile solute transport: Limits on validity and modeling implications Scott K. Hansen∗ Department of Earth and Planetary Sciences, Weizmann Institute of Science, 234 Herzl St., Rehovot 76100, Israel
a r t i c l e
i n f o
Article history: Received 6 February 2015 Revised 16 September 2015 Accepted 18 September 2015 Available online 28 September 2015 Keywords: Green’s function Sorption Mobile–immobile Dispersion Anomalous transport
a b s t r a c t Quasi-1D mobile–immobile transport processes which have exponentially distributed random waiting times in both mobile and immobile states are common in hydrologic models (for example, of transport subject to kinetic sorption). The central limit theorem implies that eventually such transport will be expressible with an effective ADE (i.e. a generalization of the common retardation factor approach with an added Fickian dispersion coefficient accounting for the effect of trapping). Previous works have determined formulae for the value of this coefficient based on the transport properties. However, the time until convergence to Gaussian behavior has not previously been quantified. To this end, exact Green’s functions characterizing the transport at all times are derived for the case of pure advection. The Green’s functions are expressed in terms of three dimensionless parameters, representing location, time, and capacity coefficient. In the pre-Gaussian regime, a parametric study characterizing concentration profile asymmetry as a function of the capacity coefficient is performed. Next, heuristics are presented in terms of the dimensionless parameters for the time until the effective ADE adequately reflects reality. For strongly retarded solute, the time until effective ADE validity is found inversely proportional to release (e.g., desorption) rate. The nature of the effective dispersion coefficient is examined, and the possibility of large trapping-driven dispersion even in cases where batch experiments would detect negligible trapping is demonstrated. Collectively, these results call into question reliance on retardation factors derived from batch experiments for many practical transport modeling efforts; knowledge of both the trapping and release kinetics appears essential. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction 1.1. Mathematical treatment of mobile–immobile processes In many subsurface solute transport scenarios, there may be two spatially coextensive domains that are out of equilibrium with each other, each having its own local concentration. Commonly this is seen in models where solute advects only when it is in one, “mobile”, state (i.e. domain), with concentration c(x, t), but can also sometimes be trapped in an “immobile” state, which has its own concentration, cim (x, t), and from which it is eventually released. These models may be expressed with the following set of equations:
∂c ∂c (x, t ) + im (x, t ) = F {c}(x, t ) ∂t ∂t ∂ cim (x, t ) = G{c, cim }(x, t ), ∂t
(1)
∗ Present address: Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. Tel.: (505) 664 0712. E-mail address: [email protected], [email protected]
http://dx.doi.org/10.1016/j.advwatres.2015.09.011 0309-1708/© 2015 Elsevier Ltd. All rights reserved.
where F is a linear differential operator representing some combination of advection, dispersion, and decay, and G is an arbitrary operator. There are a number of different causes of this sort of mobile–immobile non-equilibrium transport: including so-called chemical and physical non-equilibrium, as well as diffusion into porous media [1]. Both the standard chemical and physical nonequilibrium equations, though arising through different conceptual pictures, may be put into an equivalent mathematical form [2], where in our notation we define
G{c, cim } ≡ λc − μcim .
(2)
Here λ represents the probability per unit time that a mobile particle will become immobile, and μ represents the probability per unit time that an immobile particle will become mobile. A variety of dualporosity transport problems can also be modeled in this way [3]. In addition, local equilibrium sorption with a distribution of retardation factors has been shown, when upscaled, to be expressible in the same form [4]. In our analysis, we shall employ an alternative expression for G (which we call G∗ ), previously used by
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
Margolin et al. [5]:
G∗ {c} ≡ λc − λ
t 0
ψim (τ )c(x, t − τ )dτ ,
(3)
which eliminates cim as a dependent variable but introduces time non-locality. As in (2), λ is a spatially homogeneous probability per unit time that a mobile particle will become immobile. Here, ψ im (t) is an arbitrary probability distribution for the length of a single sojourn in the immobile phase. Substituting this into (1) leads to the integrodifferential equation
∂c (x, t ) = F {c}(x, t ) − λc(x, t ) + λ ∂t
t 0
ψim (τ )c(x, t − τ )dτ .
(4)
Similar forms have been used by Schumer et al. [6] to define what the authors call the fractal mobile–immobile (fMIM) paradigm, by Haggerty et al. [7] in the analysis of late-time BTC tails in the presence of trapping processes, and by Benson and Meerschaert [8] in their development of subordination technique for incorporating anomalous transport into ADE-based numerical models. In the Laplace domain, it is not difficult to show that Eqs. (1) and (2) are equivalent to (4) so long as
ψim (t ) = μe−μt .
(5)
While more general forms of ψ im than (5) correspond to some physical systems, we restrict our attention to first-order mobile–immobile behavior, and so will use (4) and (5) as the basis for all analysis. Note that this model implicitly defines a comparable probability distribution,
ψm (t ) = λe−λt ,
(6)
for the time taken by a single sojourn in the mobile phase (i.e. time interval between release and subsequent capture). This particular model, with exponentially-distributed times in both the sorbed and free states, has found wide application in the literature: applications are found in the sediment transport literature [9,10], as well as in analysis of solute trapping [8,11]. 1.2. Trapping-driven dispersion, breakthrough curve asymmetry, and effective ADE models Commonly, transport for each solute particle, including the random trapping, is independent of all other particles. Where this trapping is a result of diffusion into an immobile zone, this is assured. For kinetic sorption, where there exists a great excess of sorption sites, it is justified as well, although in cases of competitive sorption (i.e. limited sites), this may not be the case. Any random transport process that acts independently on different particles must have a dispersive effect, and indeed: the dispersive effect of kinetic sorption has long been recognized in the chromatography literature, dating back at least to [12]. In the subsurface transport literature, some theoretical studies quantifying the dispersive effects of mobile– immobile systems have been published. Late-time first, second, and third temporal moments of breakthrough curves for advectivedispersive subsurface solute transport subject to first-order chemical and physical non-equilibrium were considered by Valocchi [2]. A similar moment analysis, except deriving spatial moments was later performed by Michalak and Kitanidis [13]. The authors further computed effective late-time velocities and dispersion coefficients. These studies were in 1D, and focused on late time limits in the case when (2) holds, and are in the school of classic mobile–immobile zone models with exponential ψ im . Recently, Uffink et al. [11] employed random walk theory to derive an equivalent PDE to (4), and manipulated the equation to derive the same late-time dispersion coefficient shown in [13]. To our knowledge, they were the first to remark that λ and μ may affect the rate of late-time Gaussian convergence, but did not pursue this systematically, save for providing an approximate time until the effective dispersion coefficient stabilizes.
185
Related to the discussion of trapping-driven dispersion is a body of literature on the so-called local equilibrium assumption (LEA). The LEA essentially refers to conditions under which the trapped and free concentrations are related by the same ratio at all times [14]. In our terminology, we express this as c(x, t ) + cim (x, t ) ≈ Rc(x, t ), for some constant, R, referred to as a retardation factor. By substitution of R into (1), one arrives at the retarded transport equation:
R
∂c (x, t ) = F {c}(x, t ). ∂t
(7)
It is apparent that this form simply rescales time (R can be eliminated by a substitution τ ≡ t/R), and so the dispersive effect of trapping is not considered. Nonetheless, some authors conflate applicability of the local equilibrium assumption with usage of a retarded transport equation (e.g. [15]). There are some circumstances under which this is reasonable, namely those in which F is an advection–dispersion operator representing a stronger source of dispersion, compared to which trapping-driven dispersion is negligible. The LEA has been widely studied, with aforementioned paper [2] aiming to validate it by comparing moments, at late time, of exact solutions for chemical and physical non-equilibrium with a retarded advection–dispersion equation. A different approach, seeking directly to find conditions under which local equilibrium is nearly satisfied at a point and relating that to parameters in a system like (1) has also been presented [16]. Later, Wallach examined the domain of validity of the LEA through perturbation theory, treating (7) as an end member in a perturbation expansion for the exact solution which contains trapping-driven dispersion [17]. All of these authors, and many others (see the literature review in [16]) identified conditions under which the dispersive effect of sorption was comparatively negligible. However, since Michalak and Kitanidis [13] computed the effective, late time dispersion coefficient for classic mobile–immobile systems, there is now a more flexible framework available than the dispersion-ignoring LEA approximation. Mobile–immobile solute transport may be treated at late time by a retardation coefficient and this effective dispersion coefficient, even in cases where trappingdriven dispersion cannot be ignored. We term such an approach an effective ADE model. An effective ADE model still cannot apply in circumstances in which there is a significantly asymmetric transport Green’s function, however, since the Green’s function for the ADE is Gaussian. Further, it is well known that in many trapping-and-release processes, plume asymmetry occurs. Uffink et al. [11] showed how this develops in a system described by (1) and (2). Since consideration of mobile–immobile problems is so common in contaminant hydrogeology, there are naturally many alternatives to the effective ADE approach, including exact solutions for specific geometries. A comprehensive survey is beyond the scope of this work, but classes of techniques include analytic solutions such as those included in the aforementioned van Genuchten paper [3], analytic solutions that treat discrete fractures [18,19], and by multi-rate mass transfer schemes [20]. Numerical approaches that treat multiple interacting continua [21,22] are also useful for specific problems. Lastly, it bears noting that in addition to being indicative of earlytime mobile–immobile behavior for exponential ψ im , asymmetry may be caused by power-law distributed immobile times [23], as well as non-sorbing transport in heterogeneous media [24]. The difficulty of distinguishing between mobile–immobile models and transport in heterogeneous media has been remarked upon by Carrera et al. [25], among others. This implies that CTRW approaches [24] represent another viable approach for modeling mobile–immobile systems. Experimentally, at the column scale, asymmetric breakthrough curves have been examined using CTRW approaches by Deng [26] and Li and Ren [27]. This sort of behavior may arise due to anomalous transport
186
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
or mobile–immobile trapping behavior. The tailing effect of mobile– immobile trapping in columns where transport in the mobile state is anomalous has also been considered experimentally by Rubin et al. [28]. 1.3. Goals for this work
∀x > 0 x=0 ∀t
(8)
to generate Green’s functions for concentration profiles and breakthrough curves. Because of this explicit approach, we are able to make a second contribution, which is to characterize the asymmetry of pre-Gaussian behavior as a function of the governing parameters. The Green’s functions for concentration profiles and breakthrough curves are expressed in terms of three dimensionless parameters: a dimensionless location, dimensionless time, and the capacity coefficient, β = λ/μ. A parametric study is performed characterizing the early-time behavior and asymmetry of concentration profiles and breakthrough curves. The time until an effective ADE model is appropriate is quantified in terms of the dimensionless parameters and a rule, independent of β , for the dimensionless time until effective ADE applicability is formulated (implying that, for strongly retarded substances, the dimensional time is inversely proportional to μ; see below). Finally, the unpredictability of the degree of sorption-driven dispersion given β alone is observed. As a consequence of these results, we conclude that one cannot safely use a retarded ADE model (with a retardation factor computed from an equilibrium batch experiment) for sorbing solute transport without characterizing both λ and μ independently. This warning runs contrary to some current practice, and is a practical contribution complimentary to the new theoretical results. 2. Green’s functions for transport In the Laplace domain, (4) can be solved alongside (8) to yield:
cˆ(x, s)[s + λ − λψˆ im (s)] = F {cˆ}(x, s).
(9)
Since we are interested in first-order mass transfer, we define ψ according to (5). We also choose F ≡ −v∂ /∂ x (the advection operator), so that all spreading can be attributed to the mobile–immobile trapping process. Explicit Green’s function solutions for breakthrough curves, b(x, t), and concentration profiles, p(x, t), are developed in dimensional form in Appendix A. From these, one can compute c(x, t ) = t t 0 b(x, t − τ )c(0, τ )d τ , or c(x, t ) = 0 p(x − y, t )c(y, 0)dy. In the following, we will present them in terms of the following dimensionless groups:
t ≡
μt, μx/v, β ≡ λ/μ.
x ≡
2.1. Green’s functions for breakthrough curves Substituting (10) into (21) yields
In light of the above discussion, we continue the explorations begun by Michalak and Kitanidis [13] and Uffink et al. [11], and characterize quantitatively the time until an effective ADE behavior is seen as a function of transport parameters. We perform the analysis explicitly, directly solving (4) subject to the initial and boundary conditions:
c(x, 0) = 0 c(0, t ) = δ(t ) c(∞, t ) = 0
remain at equilibrium, the instantaneous particle trapping rate and particle release rate must be equal.)
(10)
The first two quantities represent, respectively, a dimensionless location and time. The third, β , is sometimes called the capacity factor. It can be thought of as representing the ratio cim /c in a system that is in local equilibrium, or equivalently as the ratio of expected length of one continuous immobile period to the expected length of one continuous mobile period. (The two are equivalent because in order to
b(x , t , β) =
μ
exp{−t + (1 − β)x }
+
β x
t − x
I {2 1
δ(t − x )
β x (t − x )} .
(11)
The terms in the outer square brackets are dimensionless. It may seem counter-intuitive that b, which may be thought of in terms of “concentration” has units of t −1 , but it is important to remember that the initial condition was defined as a temporal Dirac delta function, which also has the units of t −1 , and that b is really a Green’s function which must be integrated against a source history to produce a physical breakthrough curve. 2.2. Green’s functions for concentration profiles We nondimensionalize again, employing (10) and (26). The Green’s function for the concentration profile is
p(x , t , β) =
μ exp{−t + (1 − β)x } δ(t − x ) v β x (t − x )} + I { 2 β x , 1 t −x
(12)
which has the same mathematical form as the breakthrough curve (11), up to a constant factor. It is clear that this must be so, since either Green’s function calculation must essentially represent the plume due to an instantaneous point source, and so the distinct calculations for (11) and (12) corroborate each other. 3. Pre-asymptotic behavior: the effect of time and β Given an ultimate interest in how long it takes for transport to approach a regime in which the ADE is applicable, it is reasonable to consider the concentration profile Green’s functions and consider how large t must become before they are Gaussian. At the same time, how t and β affect their pre-asymptotic evolution and behavior also bears consideration. 3.1. Effect of changing β The capacity factor, β , relates to the effective retardation factor as R = 1 + β , which can be seen by setting G{c, cim } = 0 in (2), to represent equilibrium conditions. Consequently, it is predictable that in Fig. 1, larger values of β correspond to smaller first moments of the Green’s function curves at any fixed time: i.e. more time trapped means the plume travels a shorter distance. Although the distributions become Gaussian eventually for all values of β , their preasymptotic character depends on whether β is greater or less than unity. For small β , there is a heavy tail to the distribution for small x , and for large β the converse is true. Note that the integrals of Green’s functions tend to decrease with increasing β and t . This comes from the fact that the Green’s functions are for solute that is initially all mobile at time 0, and mass transfer between the trapped and free states (each respectively integrated over all space) asymptotically approaches a steady state in which 1/(1 + β) of the initial mass remains free.
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
187
t ' = 10
t '= 1
.1
1
.5
3 100
10 .3
.3 .5 1
100
3
10
.1
.01
.01
t ' = 50
.01
.1
100 10
3
1
.5
.3
Fig. 1. Plots of exact solution Green’s functions ((11) with prefactor μ = 1) for various fixed values of t , as β increases. Note that the various x -axes do not have the same scales. The vertical arrows in each plot are indicative of Dirac distribution representing probability mass for particles that have never been trapped. The heights of the arrows are not significant.
3.2. Effect of changing t The effect of t , the dimensionless time, is illustrated in Fig. 2. Essentially, each of the subplots in that figure can be viewed as a set of superimposed stop-motion snapshots taken at various t . Nonsorbing solute, which travels at the same rate as the underlying fluid flow, will always have traveled to the location t = x (the location of the Dirac spikes). As solute sorbs and then subsequently desorbs it is spread out behind the non-sorbing front. This leads to the progressive depletion of probability mass from the Dirac spikes into the rest of the probability distribution with increasing time. For larger β , more
depletion from the front appears for a fixed t , although actually only λ determines the degree of depletion by a given dimensional time t. 4. Asymptotic behavior: how long until the effective ADE applies? We saw previously that dispersion is unavoidable when dealing with solute transport undergoing mobile–immobile transport. It can also be shown that the transport Green’s functions (having finite statistical moments) will eventually approach Gaussian distributions. We now turn to the question of how to capture transport with an effective ADE.
188
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
(height 0.68)
=1
= .5
t' t'
t'
1
1
3
t'
5
t' t'
3
7
t' t'
5
t'
9
7
t'
=5
=2
t'
9
1
t'
t'
1
3
t' t'
5
t'
3
t'
7
t'
5
t'
9
7
t'
9
Fig. 2. Plots dimensionless Green’s functions ((11) with prefactor μ = 1) with for various fixed values of β as t increases. Vertical axis represents both probability density (for the dotted curves) and total probability mass for the Dirac spikes that are part of the exact solution. The height of each of the solid vertical arrows represents the corresponding amount probability mass each Dirac spike contains—this represents the mass that has never been trapped since being released at t = 0.
4.1. The effective ADE We wish to rewrite (4), where F ≡ −v∂ /∂ x, in the following form
∂c ∂c ∂ 2c (x , t ) = −v∗ (λ, μ) (x , t ) + D∗ (λ, μ) 2 (x, t ). ∂t ∂x ∂x
(13)
The effective v∗ will equal the advection velocity, v, multiplied by the fraction of time that a particle is free (over a long time period, this fraction converges for all particles, even though the absolute time immobile for all particles does not converge, leading to the trappingdriven dispersion that we focus on in this paper and wish to capture with D∗ ). A long time interval can be treated as consisting of an essen-
tially equal number of mobile and immobile periods, and so the retardation factor, R, can be expressed as the ratio between the expected length of a mobile–immobile cycle and the expected length of a single sojourn in the mobile phase: R = E[ψm + ψim ]/E[ψm ]. Since the distributions ψm and ψim are both exponential with parameters λ and μ, respectively, it follows that R = 1 + λ/μ. (We see here that R has a dual interpretation. It is commonly understood as the ratio between peak arrival times for a conservative vs. a trapped-and-released solute. However, on the analysis seen here, it can also be understood as the local ratio between total concentration, c + cim , and c.) Thus,
−1 λ v (λ, μ) = = v 1 + . R μ ∗
v
(14)
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
189
Fig. 3. Comparison of exact and approximate solutions to the example problem at t = 5 (left) and t = 90 (right; note different scale). Solid curves and the Dirac arrow represent the exact solution (12), and dashed curves represent the ADE approximation using an effective velocity (14) and effective dispersion coefficient, D∗ (15). The vertical dotted line, as elsewhere, indicates a discontinuity in exact solution.
The formula for effective dispersion coefficients first-order mobile– immobile transport was computed by Michalak and Kitanidis [13] (later re-derived by Uffink et al. [11], using different means). From their work, it follows that
D∗ (λ, μ) =
β 1 λμv2 = v2 . 3 μ (λ + μ) (1 + β)3
(15)
In this paper, we only analyze a system with no hydrodynamic dispersion. However, note that [11] also reported a total effective dispersion coefficient in the case of trapping and hydrodynamic dispersion when free. 4.2. Time to convergence The effective ADE model can only match true behavior after a Gaussian character has developed in the concentration profile Green’s function. For illustrative purposes, we will compare solutions of (4) and (13) for the case when λ = 0.6, μ = 1, and v = 1 (which ensure that the dimensionless and dimensional coordinates are numerically identical) and β = 0.6. Solutions at early (t = 5) and late (t = 90) times are shown in Fig. 3, illustrating the convergence of the spatial Green’s functions to Gaussian normality as the solute undergoes multiple sorption–desorption cycles. At early time, the ADE approximation under-predicts peak concentration and significantly lags the exact solution. Furthermore, the approximation predicts significant probability density for values of x less than zero and greater than five, both of which are impossible. At late time, however, the coherence between the two solutions is good, suggesting that trapping-driven retardation and dispersion coefficients may be used in the ADE. To confidently use an effective dispersion coefficient in lieu of the exact analysis, we may quantify the degree of coherence between the concentration profiles generated by the exact and the effective ADE approaches as a function of β and t , so as to avoid using the effective approach in a case where they are not sufficiently coherent. To quantify this, we define thenormalized divergence as the L1 norm of the difference between the exact and approximate solutions divided
by the L1 norm of the exact solution, and evaluate this for a number of different values of β and t . This quantity is bounded between zero (curves totally coincident) and two (curves totally disjoint), provided both curves have the same L1 norm. For concreteness, the normalized divergence in Fig. 3 is 0.51 when t = 5, and 0.09 when t = 90; this latter quantity may seem surprisingly high. We adopt a normalized divergence of 0.1 as a heuristic for applicability of an ADE model. In order to understand the regime of applicability, we performed a numerical parametric study. Since we are interested in comparing concentration profiles which are functions of x , there are only two free dimensionless parameters available, t and β . Thus, we may produce contour plots of normalized divergence as a function of these two parameters. These are shown in Fig. 4. From here, we see that for values of β greater than unity and normalized divergence less than 0.5, the contours are more or less vertical, and independent of further increases in β . From this, and the satisfactory performance seen in Fig. 3. We may develop a heuristic that it is reasonable to employ the effective diffusion coefficient approximation when t > 70 and β > 1, since this implies a normalized divergence less than 0.1. That the 0.1 contour continues vertically off the top of the graph for additional orders of magnitude was verified by generating additional plots, not shown. The fact that as β grows, the dimensionless time until the normalized divergence threshold is reached is unchanged, along with (10), implies that the actual time until this occurs varies inversely with μ. When β becomes smaller, Fig. 4 implies a longer dimensionless time, t , is required before the effective ADE model can be employed. This can seem surprising, as small values of β correspond to retardation factors near unity. However, consider the exponential decay parameters λ and that characterize ψ im (5) and ψm (6). It follows from (10) that as μ grows, it will both decrease β in inverse proportion, and increase t proportionally, and so there is no particular reason to expect an increased μ to correspond to a pathologically increased dimensional time, t, required to achieve nearly Gaussian behavior. However, this is exactly what will happen as λ decreases. This is reasonable, since particles must be sorbed at least once in order to move away from the advection front, and must sorb and desorb numerous
190
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
Fig. 4. Semi-log contour plot of normalized divergence for various values of β and t .
Fig. 5. Semi-log contour plot of normalized divergence for various values of β and β t .
times before approaching the quasi-Gaussian behavior that, according to the central limit theorem, they must, in the limit. As λ declines, the mean time between trapping events increases, and so one must wait for a longer time before quasi-Gaussian behavior manifests itself. For large β the time for a cycle to complete is controlled by μ, and for small β it is controlled by λ. If we define tc as the (essentially β -independent) dimensionless time required for convergence in the case of large β , and tc as its dimensional counterpart, then
implies β of a similar magnitude. In this case, it is possible to determine from Eq. (15) that D∗ ≈ v2 μ/λ2 . While β = λ/μ is large, implying that μ/λ 1, it is still possible that μ/λ2 1. For concreteness, consider an example where v = 1, λ = 10−4 , and μ = 10−6 , implying β = 100. These values of λ and μ represent a system that is very slow to reach equilibrium, which is sometimes observed experimentally [26]. In this case D∗ = 97.06. However, multiplying both λ and μ by 104 keeps β = 100, but yields D∗ = 9.706 × 10−3 , so D∗ varies by a factor of 10,000 while corresponding to the same β . Thus, in the case of large β (i.e. large R), a knowledge of the kinetics is required in order to realistically model the system: knowledge of only the equilibrium value capacity factor is not enough.
tc =
tc . E[ψim ]
(16)
This is to say, tc is the approximate number of trap-and-release cycles for convergence to Gaussian behavior (where the brackets “ ” represent rounding up to the nearest greater integer). For small β , It is easier to visualize the behavior by plotting against an alternative dimensionless time, β t (= λt ) instead of μt. This is shown in Fig. 5. It is apparent that the alternative dimensionless time until validity of the dispersion coefficient approach is essentially independent of the value of β for small β , although the threshold of ADE applicability is quite different. From examination of the 0.1 contour, a heuristic β t > 30 and β < 0.3 is developed for applicability of ADE behavior (or β t > 70 and β ≤ 1, to cover all possibilities). Note for comparison that heuristics are presented in [11] for disappearance of the Dirac spikes (rather than applicability of an effective ADE model): which were that t > 3/λ and t > 3/μ. Placing our heuristics in dimensional form, they amount to t > 30/λ and t > 70/μ. These are thus more conservative than those presented previously, which is reasonable since they correspond to more stringent conditions on the plume. 5. The retardation factor is uninformative about trapping-driven dispersion We now discuss the applicability of the effective ADE approach in the context of three special cases in which it is possible to simplify (15). We will see in all three cases that β (and by extension, R) is uninformative about D∗ , which may be large or small depending on the absolute magnitudes of λ and μ. Special case I: β is large. For some heavy organic molecules sorbing in certain environments, R can be greater than 1000 [29]. This scenario
Special case II: β is small. Interestingly, an exactly analogous ambiguity to that shown for large β exists when β → 0; this is to say for R ≈ 1 (using the relation R = 1 + λ/μ = 1 + β ). It follows from (15) that for small β , D∗ ≈ v2 λ/μ2 . Then by permuting λ and μ in the example shown in special case I, exactly the same conclusions follow here. This means the possibility exists for essentially arbitrarily large trappingdriven dispersion for slowly trapped and released solute even when R ≈ 1. Thus, kinetics must be considered in a case when the solute’s trapping affinity is so small that a batch experiment would not predict it to be a factor in transport. 6. Summary and conclusions We considered trapping-driven dispersion in classic, first-order mobile–immobile solute transport systems under pure advection (hydrodynamic dispersion being implicitly taken to be negligible compared to trapping-driven dispersion). A certain time-non-local transport Eq. (4) was used in the analysis; it proved possible to explicitly solve this transport equation using integral transforms, and to write explicitly the Green’s functions for both concentration profiles and breakthrough curves in terms of just three dimensionless parameters. Exact solution for the Green’s functions allowed explicit characterization of behavior over all times and quantification of the time until an effective ADE model with a retardation factor and trappingdriven dispersion coefficient is applicable. In particular, we established three concrete results regarding, respectively, the pre-Gaussian regime, the time to convergence, and the dispersion in the Gaussian regime:
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192 •
•
The value of the capacity factor, β , was seen to qualitatively affect the pre-Gaussian behavior of both concentration profiles and breakthrough curves in a simple way. Values of β < 1 were associated with tailing opposite to the advection direction and late breakthrough. Values of β > 1, contrarily, were associated with tailing in the advection direction and early breakthrough. Heuristics were found for reasonable usage of an effective ADE λ and the D∗ = λμv2 of [11]: λt > 30 for λ model with R = 1 + μ 3 (λ+μ)
•
μ and μt > 70 for μ ≤ λ, with a conservative criterion being t > 70/min(λ, μ). A corollary of this analysis is that for strongly retarded solute, the time until applicability of the effective ADE is inversely proportional to the desorption rate, μ. Similarly, the value of R in the effective ADE model was found uninformative about the strength of sorption-driven dispersion, and the same value may correspond to large or small D∗ , for different combinations of λ and μ. Interestingly, itwas seen possible to have substantialtrapping-driven dispersion with negligible trapping detectable via a batch experiment (i.e. for R ≈ 1).
One might argue that the last two results are of theoretical import only, since if one were to be in possession of λ and μ, it would be easy enough to use the exact solutions (11) or (12), and one would not need an effective ADE at all. However, we believe they provide a warning about the modeling practice of characterizing dispersion separately from sorption, and then simply placing a retardation factor derived from equilibrium arguments into the ADE. Because an arbitrarily large degree of sorption-driven dispersion may be present without impacting R or the conservative-tracer D, and the time until any Gaussian model is applicable cannot be determined without knowledge of λ and μ, there is no substitute for determining them. A possible method for determination of λ and μ is suggested by our exact Green’s function solution, (11). Fitting this to a breakthrough curve from a short column experiment, and simultaneously determining β and pre-multiplication factor μ from the skewness and total recovery, respectively, is a potential means of accomplishing this. (The presence of pore-scale hydrodynamic dispersion, which we have treated as negligible in the above analysis, would naturally be a complicating factor.) Regardless of the viability of this approach, the characterization of both the capture and the release processes separately is suggested as a potentially important topic for future investigation. A concrete example of the significance of these results lies in the proper modeling of transport of strongly sorbing organic pollutants: Hexachlorinated PCB congeners are known to have retardation factors (and thus β values) on the order of 106 [29], indicating times to ADE applicability that are μ-dependent. Very slow desorption rates have also been observed for these congeners, with mean desorption times on the order of weeks [30], implying correspondingly small μ values. From (15) and (16), respectively, this implies very large D∗ and tc values. It is thus apparent that in similar cases, the mere substitution of R directly into the ADE applicable to untrapped solute is apt to lead to seriously misleading results. Instead, it appears only reasonable to employ (4) with appropriately calibrated λ and μ.
191
A1. Green’s function for transient release To determine Green’s functions for breakthrough curves (i.e. kernels b(x, t) that are convolved with c(0, t) to determine c(x, t)), we proceed directly from (9), writing
−v
λμ ˆ ∂ bˆ b = s+λ− ∂x s+μ = (s + μ) + (λ − μ) −
λμ ˆ b. (s + μ)
(17)
This can be solved according to the initial boundary conditions (8) to yield
x bˆ = exp −
v
(s + μ) + (λ − μ) −
Applying the shift theorem,
b(x, t ) = e
−μt −(λ−μ)x/v
e
L
−1
λμ s ( + μ)
exp −
x
s−
v
.
λμ
(18)
.
s
(19)
Defining the quantities a ≡ x/v and d ≡ λμx/v, the inverse Laplace transform in (19) can be rewritten as
L−1 {e−as (ed/s − 1) +e−as } = L−1 e−as ∗ L−1 {ed/s −1} + L−1 {e−as }
=
δ(t − a) +
d I1 {2 d(t − a)}, t −a (20)
which was resolved using table look-up (see [31], p. 255, #5.66) where I1 is the modified Bessel function of the first kind. Successively substituting the definitions of a and d into (20) and into (19) yields a complete time-domain solution:
b(x, t ) = e−μt e−(λ−μ)x/v
× 2
δ(t − x/v) +
λμx (t − x/v) , v
λμ tv x
−1
I1
(21)
This equation is of the same form as that derived by Giddings and Eyring [12] for their probabilistic chromatography equations. An equivalent equation is also derived for the same scenario, using a different conceptual approach, in [8]. In that paper, the authors are concerned with converting between a given moment in time, t, and its corresponding operational time, u (in our context, the amount of time the particle is free). By identifying the operational time with x/v in (26)—natural, since this is a pure advection system—we recover directly their solution. A2. Green’s function for specified initial concentration To develop the Green’s function, p, for the concentration profile, we again solve the same PDE, but require different initial-boundary conditions. Instead of (8), we employ
Acknowledgment
p(x, 0) = δ(x) c(−∞, t ) = 0 c(∞, t ) = 0
I thank the Azrieli Foundation for a postdoctoral fellowship which in part made this work possible.
We then proceed directly from (9), as before, taking the Laplace transform:
Appendix A. Derivation of Green’s functions The derivation of the Green’s functions in dimensional coordinates for both concentration profiles and breakthrough curves in the case when ψ(t ) = μe−μt is presented.
δ(x) − v
t=0 ∀t ∀t
λμ ∂ pˆ ˆ = s+λ− p. ∂x s+μ
(22)
(23)
We then take the Fourier transform from x ω:
pˆ =
1 . λμ s + λ − s+ μ + viω
(24)
192
S.K. Hansen / Advances in Water Resources 86 (2015) 184–192
Dividing top and bottom by v, this can be solved by table look-up as
pˆ =
1
v
exp −
x
v
(s + μ) + (λ − μ) −
λμ (s + μ)
H (x),
(25)
where H is the Heaviside step function, which is unity on the interval [0, ∞). Apart from the initial factor of 1/v, then, on the interval of interest this is exactly (18), which has already been solved. Thus, we may conclude the exact solution for the concentration profile Green’s function is
p(x, t ) =
μ −μt −(λ−μ)x/v λμ I1 e e δ(t − x/v) + t v v −1 x
λμx × 2 (t − x/v) . v
(26)
References [1] Weber WJ Jr, Pennell KD, Dekker TJ, Abriola LM. Sorption and retardation of organic contaminants in subsurface systems: effects on transport and fate. In: Aral MM, editor. Advances in groundwater pollution control and remediation. No. 9. NATO ASI series. Netherlands: Springer; 1996. p. 1–31. ISBN 978-94-0106576-4 978-94-009-0205-3. [2] Valocchi AJ. Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils. Water Resour Res 1985;21(6):808–20. http://dx.doi.org/10.1029/WR021i006p00808. http://onlinelibrary.wiley.com/doi/10.1029/WR021i006p00808/abstract. [3] van Genuchten M. A general approach for modeling solute transport in structured soils. In: Proceedings of the 17th international Congress on the hydrogeology of rocks of low permeability, vol. 17, no. (2). In: Mem Int Assoc Hydrogeol. Tucson, AZ; 1985. p. 513–26. [4] Dentz M, Berkowitz B. Exact effective transport dynamics in a one-dimensional random environment. Phys Rev E 2005;72(3):031110. http://dx.doi.org/10.1103/ PhysRevE.72.031110. http://link.aps.org/doi/10.1103/PhysRevE.72.031110. [5] Margolin G, Dentz M, Berkowitz B. Continuous time random walk and multirate mass transfer modeling of sorption. Chem Phys 2003;295(1):71–80. http://dx.doi.org/10.1016/j.chemphys.2003.08.007. http://www.sciencedirect. com/science/article/pii/S0301010403004245. [6] Schumer R, Benson DA, Meerschaert MM, Baeumer B. Fractal mobile/immobile solute transport. Water Resour Res 2003;39(10):1296. http://dx.doi.org/10.1029/2003WR002141. http://onlinelibrary.wiley.com/doi/ 10.1029/2003WR002141/abstract. [7] Haggerty R, McKenna SA, Meigs LC. On the late-time behavior of tracer test breakthrough curves. Water Resour Res 2000;36(12):3467–79. http:// dx.doi.org/10.1029/2000WR900214. http://onlinelibrary.wiley.com/doi/10.1029/ 2000WR900214/abstract. [8] Benson DA, Meerschaert MM. A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations. Adv Water Res 2009;32(4):532–9. http://dx.doi.org/10.1016/j.advwatres.2009.01.002. http: //www.sciencedirect.com/science/article/pii/S0309170809000128. [9] Barry DA, Sander GC, Jomaa S, Heng BCP, Parlange JY, Lisle IG, et al. Exact solutions of the Hairsine–Rose precipitation-driven erosion model for a uniform grain-sized soil. J Hydrol 2010;389(3–4):399–405. http://dx.doi.org/ 10.1016/j.jhydrol.2010.06.016. http://www.sciencedirect.com/science/article/pii/ S0022169410003513. [10] Lisle IG, Rose CW, Hogarth WL, Hairsine PB, Sander GC, Parlange JY. Stochastic sediment transport in soil erosion. J Hydrol 1998;204(1–4):217–30. http://dx.doi.org/10.1016/S0022-1694(97)00123-6. http://www.sciencedirect. com/science/article/pii/S0022169497001236. [11] Uffink G, Elfeki A, Dekking M, Bruining J, Kraaikamp C. Understanding the nonGaussian nature of linear reactive solute transport in 1D and 2D. Transp Porous Media 2012;91(2):547–71. http://dx.doi.org/10.1007/s11242-011-9859-x. http:// link.springer.com/article/10.1007/s11242-011-9859-x. [12] Giddings JC, Eyring H. A molecular dynamic theory of chromatography. J Phys Chem 1955;59(5):416–21. http://dx.doi.org/10.1021/j150527a009. http://dx.doi. org/10.1021/j150527a009. [13] Michalak AM, Kitanidis PK. Macroscopic behavior and random-walk particle tracking of kinetically sorbing solutes. Water Resour Res 2000;36(8):2133–46. http://dx.doi.org/10.1029/2000WR900109. http://onlinelibrary.wiley.com/doi/ 10.1029/2000WR900109/abstract.
[14] Ptacek C, Gillham RW. Laboratory and field measurements of non-equilibrium transport in the Borden aquifer, Ontario, Canada. J Contam Hydrol 1992;10:119– 58. [15] Lee L, Rao P, Brusseau ML, Ogwada RA. Nonequilibrium sorption of organic contaminants during flow through columns of aquifer materials. Environ Toxicol Chem 1988;7:779–93. http://dx.doi.org/10.1002/etc.5620071001. http:// onlinelibrary.wiley.com/doi/10.1002/etc.5620071001/abstract. [16] Bahr JM, Rubin J. Direct comparison of kinetic and local equilibrium formulations for solute transport affected by surface reactions. Water Resour Res 1987;23(3):438–52. http://dx.doi.org/10.1029/WR023i003p00438. http:// onlinelibrary.wiley.com/doi/10.1029/WR023i003p00438/abstract. [17] Wallach R. A small perturbations solution for nonequilibrium chemical transport through soils with relatively high desorption rate. Water Resour Res 1998;34(1):149–54. http://dx.doi.org/10.1029/97WR02490. http://onlinelibrary. wiley.com/doi/10.1029/97WR02490/abstract. [18] Feenstra S, Cherry JA, Sudicky EA, Haq Z. Matrix diffusion effects on contaminant migration from an injection well in fractured sandstone. Ground Water 1984;22(3):307–16. http://dx.doi.org/10.1111/j.1745-6584.1984.tb01403.x. http: //onlinelibrary.wiley.com/doi/10.1111/j.1745-6584.1984.tb01403.x/abstract. [19] Tang DH, Frind EO, Sudicky EA. Contaminant transport in fractured porous media: analytical solution for a single fracture. Water Resour Res 1981;17(3):555– 64. http://dx.doi.org/10.1029/WR017i003p00555. http://onlinelibrary.wiley.com/ doi/10.1029/WR017i003p00555/abstract. [20] Haggerty R, Gorelick SM. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour Res 1995;31(10):2383–400. http://dx.doi.org/10.1029/95WR10583. http:// onlinelibrary.wiley.com/doi/10.1029/95WR10583/abstract. [21] Silva O., Carrera J. mod_process_mrmt.f90: a Fortran module to include multi-rate mass transfer like processes: user’s guide. 2008. http://www.h2ogeo.upc.edu/English/software/MRMT%20module%20in% 20GroupWeb/MRMT%20module%20in%20GroupWeb/User’sManual.rar. [22] Pruess K, Oldenberg C, Moridis G. TOUGH2 users guide, version 2. Tech. rep. LPNL-43134. Earth Sciences Division, Lawrence Berkeley National Laboratory; 2012. http://esd1.lbl.gov/files/research/projects/tough/documentation/TOUGH2_ V2_Users_Guide.pdf. [23] Haggerty R, Wondzell SM, Johnson MA. Power-law residence time distribution in the hyporheic zone of a 2nd-order mountain stream. Geophys Res Lett 2002;29(13):18-1–4. http://dx.doi.org/10.1029/2002GL014743. http:// onlinelibrary.wiley.com/doi/10.1029/2002GL014743/abstract. [24] Berkowitz B, Cortis A, Dentz M, Scher H. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev Geophys 2006;44(2):RG2003. http://dx.doi.org/10.1029/2005RG000178. http:// onlinelibrary.wiley.com/doi/10.1029/2005RG000178/abstract. [25] Carrera J, Sánchez-Vila X, Benet I, Medina A, Galarza G, Guimerà J. On matrix diffusion: formulations, solution methods and qualitative effects. Hydrogeol J 1998;6(1):178–90. http://dx.doi.org/10.1007/s100400050143. http://link. springer.com/article/10.1007/s100400050143. [26] Deng J, Jiang X, Zhang X, Hu W, Crawford JW. Continuous time random walk model better describes the tailing of atrazine transport in soil. Chemosphere 2008;71(11):2150–7. http://dx.doi.org/10.1016/j.chemosphere.2008.01.001. http: //www.sciencedirect.com/science/article/pii/S0045653508000039. [27] Li N, Ren L. Application of continuous time random walk theory to nonequilibrium transport in soil. J Contam Hydrol 2009;108(3–4):134–51. http://dx.doi.org/10.1016/j.jconhyd.2009.07.002. http://www.sciencedirect.com/ science/article/pii/S0169772209000904. [28] Rubin S, Dror I, Berkowitz B. Experimental and modeling analysis of coupled nonFickian transport and sorption in natural soils. J Contam Hydrol 2012;132:28–36. http://dx.doi.org/10.1016/j.jconhyd.2012.02.005. http://www.sciencedirect.com/ science/article/pii/S0169772212000265. [29] Hansen SK, Kueper BH. Forensic misidentification of Aroclor sources in fractured bedrock due to “chromatographic” polychlorinated biphenyl (PCB) congener separation. Environ Forensics 2011;12(1):25–34. http://dx.doi.org/ 10.1080/15275922.2010.547547. http://www.tandfonline.com/doi/abs/10.1080/ 15275922.2010.547547. [30] Ghosh U, Weber AS, Jensen JN, Smith JR. Relationship between PCB desorption equilibrium, kinetics, and availability during land biotreatment. Environ Sci Technol 2000;34(12):2542–8. http://dx.doi.org/10.1021/es9905389. http://pubs.acs. org/doi/abs/10.1021/es9905389. [31] Oberhettinger F, Badii L. Tables of Laplace transforms. Berlin: Springer-Verlag; 1973. ISBN 978-3-642-65645-3.