A general real-time formulation for multi-rate mass transfer problems

Many flow and transport phenomena, ranging from delayed storage in pumping tests to tailing in river or aquifer tracer breakthrough curves or slow kinetics in reactive transport, display non-equilibrium (NE) behavior. These phenomena are usually modeled by non-local in time formulations, such as multi-porosity, multiple processes non equilibrium, continuous time random walk, memory functions, integro-differential equations, fractional derivatives or multirate mass transfer (MRMT), among others. We present a MRMT formulation that can be used to represent all these models of non equilibrium. The formulation can be extended to non-linear phenomena. Here, we develop an algorithm for linear mass transfer, which is accurate, computationally inexpensive and easy to implement in existing groundwater or river flow and transport codes. We illustrate this approach by application to published data involving NE groundwater flow and solute transport in rivers and aquifers.


Introduction
Solving flow and solute transport phenomena in natural media requires using variables, such as heads and concentrations, that characterize the state of the system at every point.As such, they are termed state variables.State variables are assumed representative of a small portion of water around Correspondence to: O. Silva (orlando.silva@idaea.csic.es)such point.Using them to model flow and transport implicitly requires assuming local equilibrium.
Even though local equilibrium is assumed by default, nonequilibrium (NE) behavior is frequently observed in flow and transport through water bodies.Symptoms of NE depend on the actual phenomenon under study.In flow problems, NE is typically manifested by a delayed drainage (Gerke and van Genuchten, 1993a, b) or a delayed storage mobilization (Fig. 1a).Head in a water body often rises fast at early times in response to an inflow, which suggests a small storage capacity.However, the rate of head rise may decay at later times, suggesting an increase in storage capacity (Boulton, 1955).In tracer transport problems, breakthrough curves typically display fast and sharply rising limbs (Fig. 1b), which might denote a small volume of displaced water, but then display a long tail at late times, which suggests the opposite (Valocchi, 1985;Carrera, 1993;Cortis and Berkowitz, 2004).In reactive transport problems, NE is manifested, for example, by reaction rates, orders of magnitude slower than what would be expected from laboratory experiments (White and Peterson, 1990).
The common thread of these observations is that they can be attributed to spatial variability.The phenomenon (flow, transport, chemical reaction) appears to affect a small portion of the water body at early times, but a large portion at late times.In fact, actual explanations and models of NE behavior are phenomenon specific, but they typically involve visualizing the water body as consisting of mobile and immobile (or less mobile) regions.In water flow through permeable media, NE has been attributed to delayed storage mobilization, either because of resistance at the aquifer free In both cases, the system reacts with small storage capacity at early times, but much larger at late times, suggesting that some time is needed before the full storage capacity is mobilized.
surface (Boulton, 1955;Neuman and Witherspoon, 1971;Neuman and Tartakovsky, 2009), or at low permeability blocks (Barenblatt et al., 1960;Warren and Root, 1963).It has also been attributed to heterogeneity (Cortis and Knudby, 2006).NE has also been observed in multiphase flow and extensively studied in the petroleum field (e.g., Kazemi et al., 1976;Kazemi and Gilman, 1993).Here, departures from equilibrium are attributed to the permeability contrast between porous blocks, which store most of the fluid on the reservoir, and fractures, which have a low storage capacity but a high permeability.Fluid flow occurs mainly through the fractures, whereas the matrix blocks act as fluid sink/sources (Da Prat, 1982).
Water flow and solute transport in unsaturated fractured rock and soils also display non-equilibrium behavior (Gerke and van Genuchten, 1993a, b;Zhou et al., 2006).Usually, the imbibition rate or fluid flow between the fractures and the matrix domains is driven nonlinearly by the local pressure difference in liquid-phase pressure between the fracture and the matrix blocks (Zimmerman et al., 1995).
Non-equilibrium also has been observed in solute transport through rivers that are influenced by the exchange of water between the river and the underlying hyporheic zone (Fernald et al., 2001;Boano et al., 2007;Marion et al., 2008) or by aggregated dead-zones (Beer and Young, 1983;Lees et al., 1998Lees et al., , 2000;;Davis et al., 2000).Dead zones can be associated to a number of effects, such as reverse flows induced by bends and pools, side pockets, zones between dikes, turbulent eddies, and wakes behind bed irregularities and roughness elements (Deng et al., 2004).Non-equilibrium has also been observed in natural systems like saturated hillslopes constituted by well-connected fractures and a main soil matrix (Zhang et al., 2009).
Mathematical formulations to simulate NE behavior are as diverse as the phenomena they represent.However, most of them can be viewed as non-local in time.This means that storage mobilization does not depend solely on heads (or concentration) at the current time, but also on their history.In practice, this can imply that a sink-source term (e.g.Carrera et al., 1998) or that an additional storage term (e.g.Haggerty and Gorelick, 1995)  -Dual porosity formulations (e.g., Warren and Root, 1963) view the medium as consisting of two overlapping continua (mobile and immobile).The mobile continuum allows lateral exchanges of water and solute mass.The immobile continuum only exchanges mass with the mobile region.Actual solution of the problem can be achieved either by spatial discretization, Laplace transform of semi-analytical solutions.Even though dual porosity formulations can be viewed as rather simple, they represent a basic building block for other formulations.
-Integro-differential formulations (Herrera and Rodarte, 1973;Herrera and Yates, 1977;Carrera et al., 1998) consist of representing mass transfer to the immobile region as the convolution of the state variable (head or concentration) past history and a memory function.They can be solved by semi-analytical methods or Laplace transform (Carrera et al., 1998).
-Multi-rate mass transfer (MRMT) (Roth and Jury, 1993;Haggerty and Gorelick, 1995;Carrera et al., 1998;Haggerty et al., 2000;Wang et al., 2005;Gouze et al., 2008;Benson and Meerschaert, 2009) is a generalization of dual-porosity models, in which mass transfer between a mobile and multiple immobile regions is modeled by diffusive or first-order mass transfer terms.Here, the medium heterogeneity is mapped onto a distribution of retention times in the immobile regions.
-Continuous time random walk (CTRW) (Berkowitz and Scher, 1998;Dentz et al., 2004;Berkowitz et al., 2006;Le Borgne et al., 2008).In this modeling approach the movements of solute particles in a heterogeneous medium are modeled as random walks in space and time.That is, both space and time increments are random variables.The latter accounts for the fact that transport in low permeability zones is slower than in Hydrol.Earth Syst.Sci., 13, 1399-1411, 2009 www.hydrol-earth-syst-sci.net/13/1399/2009/ high permeability regions.Here the medium heterogeneity is mapped onto the distribution of time increments.In general the space and time increments are coupled.
These formulations look different.However, when restricted to non-locality in time (some of them can also represent nonlocality in space), they are not.Dual porosity formulations are a special case of MRMT.Dentz and Berkowitz (2003) showed that the MRMT approach and the uncoupled version of the CTRW approach are formally equivalent and gave an explicit map between them.Haggerty et al. (2000) discuss the equivalence between MRMT and integrodifferential formulations.The equivalence of space and time fractional dynamics and these approaches has been discussed by, e.g., Metzler and Klafter (2000), Berkowitz et al. (2006), Schumer et al. (2003), Benson and Meerschaert (2009) and Zhang et al. (2009).In essence, a fractional order time derivative is equivalent to a power-law memory function or a power-law distribution of retention times (MRMT) or time increments (CTRW).
In hindsight, the equivalence of such a diverse set of formulations sounds evident.In essence, they can all be viewed as expressing linear mass transfer between mobile and immobile regions.That this exchange is characterized by a distribution of residence times or time increments, by a memory function or by a fractional order derivative does not look critical.However, a broad set of characterizations hinders comparison and synthesis, thus slowing scientific development.A unified approach is needed.
An additional problem comes from the mathematical difficulty of some of these formulations.Integrodifferential equations or Laplace transform are extremely efficient, but not easy to grasp by many.Worse, their applicability is restricted to linear problems.Therefore, they are not immediate to translate to non-linear phenomena, such as multicomponent reactive transport (Donado et al., 2009;Willmann et al., 2009).
The objective of this work is to propose a unified formulation that is computationally efficient and can be extended to non-linear phenomena.We adopt the MRMT approach for two reasons.First, it is intuitive and does not require mathematical tools beyond basic calculus.Second, it allows localization.As discussed above, the domain is assumed to consist of a mobile continuum and several overlapping immobile continua.These exchange mass linearly with the mobile region.In this way, the state of immobile zones can be characterized by heads (or concentrations).That is, the traditional single state variable is substituted by several state variables (as many as immobile regions).Effectively, these work as local state variables, that can be used to model non-linear processes.

Governing equations
Non-local in time formulations can be used to enrich the behavior of either the flow or transport equations, or both.In either case, they can be viewed in two complementary fashions: (i) as a continuum of delayed storage terms, in which case these equations represent the total mass balance in both mobile and immobile regions, or (ii) as a continuum of sink/source terms, which act as linear mass exchange terms between mobile and immobile zones.In practice, the continuum is substituted by a discrete number of terms.Using the first view, the flow equation becomes where ] represents a sink/source (recharge/extraction), and h im,j [L] and S im,j [L −1 ] are head and specific storage coefficients of the j th immobile zone, respectively.Water storage in each immobile region is fed by a linear exchange with the mobile domain where σ im,j [L 2 L −3 ] is the specific surface of the j th immobile region, L im,j [L] its distance from the mobile zone and K im,j [L T −1 ] its hydraulic conductivity.
Analogously, the solute transport equation expresses the solute mass balance per unit volume of aquifer is the mobile porosity (volume of pores per unit aquifer volume), and R m [−] is the mobile zone retardation factor.Similarly, −3 ] and R im,j [−] are the concentration, porosity and retardation factor of the j th immobile zone.As in flow phenomena, mass balance in the j th immobile region is given by where φ im,j [L 3 L −3 ] is its porosity (volume of pores per unit volume of immobile region), D im,j [L 2 T −1 ] is a molecular diffusion coefficient in the j th immobile region.Equation (4) can be somewhat simplified by writing φ im,j as a O. Silva et al.: A real-time formulation for multi-rate mass transfer problems function of φ im,j (see e.g., Carrera et al., 1998).In practice, the physical meaning of the parameters in Eqs. ( 2) and ( 4) is somewhat approximative.What is important is that (1) N immobile regions are considered, (2) their state is characterized by a state variable (u im ), and (3) they exchange mass with the mobile region proportionally to the difference in state variables.Bearing this in mind, Eqs.(1-4) can be written in general as where u=h, for flow or u = c, for solute transport.The β j [−] coefficients are called capacity coefficients (R im,j φ im,j for transport or S im,j for flow) to account for the distribution of mass in the immobile phases; β [−] is the capacity coefficient of the mobile phase (R m φ m for transport or S m for flow); and α j [T −1 ] is a first-order mass transfer rate coefficient.The right-hand side of Eqs. ( 1) and ( 3) is designated generically by the operator L u .Denoting F j the j th term of the sum in Eq. ( 5), the governing equations for the immobile zones are given as where F is the total exchange between mobile and immobile regions.

Numerical formulation for MRMT
Spatial and time discretization of either the flow or transport equations under local equilibrium assumptions (i.e., without MRMT) leads to a linear system of equations where u m =u k+1 m −u k m , t=t k+1 −t k is the time step, θ is a weighting factor and the superscripts stand for the time in which the variable is evaluated.
Accounting for MRMT can be achieved in two ways: (a) using an appropriate mesh with nodes representing the immobile zones (e.g., Neuman et al., 1982), or (b) by eliminating the unknown in the immobile region as an explicit state variable, i.e. expressing u im,j as a function of u m (e.g., Carrera et al., 1998).Here, we have adopted the later approach because: first, it maintains the number of unknowns unchanged and, second, it is actually simpler to implement into existing generic flow and transport simulation codes.
Figure 2 displays a schematic representation of a hypothetical numerical mesh that includes both the mobile and immobile domains.We assume that each node m of the mobile zone is connected to all adjacent nodes of the mesh and to all the immobile blocks.Node im,j of the immobile region is only connected to node m.Geometrically, node im,j overlaps with node m.We show below that the variable u at node im,j (i.e., u im,j ) can be solved explicitly as a function of u m .Therefore, node im,j needs not be an "uncertain" node, but can be considered as a zero-D node that contributes to the mobile region, and is only updated after mobile state variables have been computed.
We first solve the N first-order ordinary differential Eq. ( 6) in terms of u m , while assuming that u m varies linearly during each time increment.
That is, u m =u k m + u m t t − t k .This leads to N first-order linear differential equations, whose solution is Combining Eqs. ( 6), ( 7a) and ( 9), the flux F j evaluated at time t k+θ will be Notice that this flux is only a function of u at the previous time step and u m .Therefore, the total mass flux, F k+θ , given by Eq. (7b), can be substituted into Eq.( 8).This leads to a system that is identical to Eq. ( 8), except that the storage matrix and sink/source term are modified according to where u k m,i is the value of u at node i of the mobile region and time t k , u k im,j i is the corresponding value in the j th immobile block, and v i is the volume of cell i in volume integrated formulations (e.g., finite element) and is equal to 1.0 in discretized formulations (e.g., finite differences).Finally, it is necessary to update u im,j at the end of each time step using Eq. ( 9).This approach is quite simple to program and should lead to accurate solutions at a very low computational cost.As with the integro-differential approach, the number of nodes/elements is not altered by the addition of the MRMT terms.
These equations were implemented in a Fortran 90 module called mod process MRMT.f90, which is structured following the coding guidelines and rules proposed by Slooten (2009).
The module defines MRMT objects by means of a specific type of variables, that contain the capacity terms β j and mass transfer coefficients α j , and provides services to solve the equation described here.The module supports some characteristics of objectoriented programming.In fact, the module was designed such that its types and operations are available from outside but the details of the implementation are hidden from the user ("black box" principle).These functionalities facilitate linking other Fortran modules or objects to mod process MRMT.f90.For instance, reactive transport may be included in both the mobile and immobile regions, by properly linking the object-oriented tool CHEP-ROO (Bea et al., 2009) to both any conservative transport code and the module mod process MRMT.f90.The module can be downloaded from http://www.h2ogeo.upc.es/English/English/software.htm#Modprocess MRMT.

Equivalence with other similar approaches
As mentioned in the introduction, a large number of nonlocal in time schemes have been presented by different authors.This section is devoted to discussing their equivalence so as to facilitate using them in the proposed formulation.This equivalence is summarized in Table 1 in terms of coefficients α j and β j of the generalized MRMT model proposed in this work (Eqs.5 and 6).Haggerty et al. (2000) provided a comparison table with different MRMT formulations considering governing equations similar to Eqs. ( 5) and ( 6).The present approach is essentially identical to that of Haggerty and Gorelick (1995).The main difference is that they formulated their equations per unit volume of water.Therefore, their capacity coefficients are equal to the coefficients β j of Eqs. ( 5) and ( 6), but divided by the mobile capacity, β.Denoting β j HG and α j HG the capacity and first-order mass transfer coefficients considered by Haggerty and Gorelick (1995), we have the following equivalence relationship

Equivalence with the classical MRMT approach
We have preferred to use capacity coefficients as defined in Eq. ( 5) to keep the physical meaning and consistence of the governing equations as mass balances per unit volume of aquifer.

Equivalence with the integro-differential approach
Many schemes approximate the effect of the immobile region by a continuous memory function.The governing equations are then solved in the Laplace domain.These solutions can also be approximated by expanding the memory function as a sum of exponentials (Carrera et al., 1998).Each summand can then be solved as explained in Sect.3.For the purposes of comparison with our approach, the important issue is to acknowledge that such approaches are typically defined in terms of overall parameters for the whole immobile region (as opposed to independent α j and β j ).For the works of Carrera et al. (1998) and Salamon et al. (2006), the equivalence is given by: where are characteristic parameters of the entire immobile domain.
Coefficients a j and γ j can be found in the literature (e.g., Haggerty and Gorelick, 1995;Carrera et al., 1998;Haggerty et al., 2000;Salamon et al., 2006) for diffusion into different geometries (layered, cylindrical, spherical and veins) and the standard first-order model.These formulations result from the analytical solution of the diffusion equation.A large number of first-order mass transfer rate coefficients and their distributions estimated from field and laboratory test results can be obtained from the works of Cosler (2004) and Haggerty et al. (2004).
The mass flux, F , into the immobile region in memory function based approaches is given by Table 1.Equivalence between the present formulation and other approaches, on the basis of coefficients α j and β j that appear in Eqs. ( 5) and ( 6).L ] is the porosity of the jth immobile zone.
Flow with delayed storage (e.g., Boulton, 1955) L ] is the specific surface of the jth immobile region; Kim,j Integro-differential (e.g., Carrera et al., 1998) 2 T ] and Lim [L] are characteristic parameters of the entire immobile domain with the same meaning as in solute transport; coefficients aj and j γ are shape factors for diffusion into different geometries.a CTRW (e.g., Dentz and Berkowitz, 2003) k are the expansion coefficients of the Laplace transform of the memory function g * .
Fractional derivatives ( ) ( ) Willmann et al., 2008) g is the slope of the memory function in log-log scale; [t1, tN] is the interval of time on which this function displays a power-law behavior.

Solute transport in rivers
(e.g., Lees et al., 2000) is the stream channel cross-sectional area; and As [L 2 ] is the storage zone crosssectional area.
Low permeability blocks (e.g., Warren and Root, 1963) is the porosity of the matrix blocks; Cim [L T 2 M -1 ] is the total compressibility of matrix blocks.
where g is the memory function and * denotes the convolution product.Carrera et al. (1998) approximate this product using the integro-differential approach of Herrera and Rodarte (1973) and Herrera and Yates (1977).An equivalent alternative is to approximate g by where b(α) [T ] is a density function of first-order rate coefficients.Haggerty et al. (2000) provide explicit expressions for the density and memory functions for various models or geometries.To use the approach of Sect.3, we need to express the memory function as α j β j e −α j t (16) Note that Haggerty et al. (2000) included the factor α j β j on the memory function, unlike Carrera et al. (1998) who placed it on flux F .However, both approaches are equivalent.We calculate the convolution product in Eq. ( 14), truncating the memory function at Nth term and following the same algebraic analysis described in the Appendix 1 of Carrera et al. (1998).Thus, we can express F k+θ as β j 1 − e −α j θ t (17a) The equivalence between our approach and integrodifferential approach becomes evident by comparing Eqs.(17a) and ( 10).Also note that, from Eq. ( 9) we obtain the recursive relationship  10) by imposing the condition The coefficients in Eq. ( 13) result from an infinite series expansion that needs to be truncated.Truncation criteria and expression for the final terms of truncated multi-rate series can be found in Haggerty and Gorelick (1995) and Salamon et al. (2006).In the case of diffusion into different geometries, they proposed the same criteria to evaluate β N .However, while Haggerty and Gorelick (1995) suggested writing α N as the rest of the α j coefficients (i.e., Eq. 13a), Salamon et al. (2006) proposed the following expressions where λ for layers, spheres, cylinders are given by Salamon et al. (2006).Dentz and Berkowitz (2003) found a mathematical equivalence between MRMT and the CTRW model (Berkowitz and Scher, 1998;Dentz et al., 2004;Berkowitz et al., 2006;Margolin et al., 2003;Salamon et al., 2006;).They formulated a CTRW approach which is formally equivalent to the integrodifferential formulation of MRMT presented in this paper.They present a map between the memory function defined in the context of integro-differential formulations, g, and the transition time distribution ψ(t)

Equivalence with CTRW
where τ 0 defines which portion of the medium is mobile or immobile and as such is related to the mobile and immobile volume fractions of the medium (see Dentz and Berkowitz, 2003).The Laplace transform of the memory function, g * , can be expanded into a series in s according to where explicit expressions for the a k are given in Dentz and Berkowitz (2003).For g * (s) given by the Laplace transform of Eq. ( 16), we obtain By comparison of Eqs. ( 22) and ( 23), we obtain relations between the β j and the a k for a given series of rates The latter expression can be inverted (numerically) in order to obtain explicit expressions for the weights β j and thus for the memory function g(t) that simulates the transport behavior in a CTRW.

Equivalence with the fractional derivatives approach
Another method to describe non-localities in time is the fractional-order advection-dispersion equations (fADE) (Benson et al., 2000a, b;Zhang et al., 2009).This approach is quite general and can be characterized by a power law memory function when only the time derivative term in the advection-dispersion equation (ADE) is fractional (Willmann et al., 2008).
On the other hand, Dentz and Berkowitz (2003) proposed the use of the truncated power law memory function, which has become widely used because breakthrough curves often display a power law behavior at late times (see, e.g., Zhang et al., 2007;Willmann et al., 2008).The late time behavior of the breakthrough curve can be related to the memory function (Haggerty et al., 2000).This memory function only requires specifying the slope of the memory function in loglog scale, m g , and the interval of time [t 1 , t N ] on which this function displays a power-law behavior.A practical method to calculate the distribution coefficients β j consists of, first, calculate the α j values assuming they are evenly distributed on a logarithmic scale while fixing α 1 = t −1 N and α N =t −1 1 .Secondly, we obtain a recursive relationship for β j values by approximating the memory function with expressions of successive increasing orders, i.e. β i α i = m g log t j − log t j +1 (25) where t j =α −1 j .This leads to To get the values of β j , we first assign an arbitrary value to β N (e.g., β N =1).Then we apply Eq. ( 26) and finally scale these values imposing the condition

Equivalence with other models and phenomena
There are a lot of specific models, problems and other applications that can be cast as non-local in time formulations.Many of those models rely on single rate first-order mass transfer (e.g., Roth and Jury, 1993).Haggerty and Gorelick (1995) have classified these mass transfer models as "chemical" or "physical" with subdivisions into one-site or two-sites models.They have shown how these models are equivalent and interchangeable.
As an examples, we want to mention the model of flow through low permeability blocks (Warren and Root, 1963), the transient storage model of longitudinal solute transport in rivers and streams (Bencala and Walters, 1983), and transport under linear kinetic adsorption (e.g., Nkedi-Kizza et al., 1983).In essence, all of these models are particular cases of the single first-order mass transfer model, in which there is only one immobile zone exchanging mass with the mobile domain.Warren and Root (1963) formulated a model for the singlephase flow of a slightly compressible liquid through formations containing an intergranular porosity and a fissure porosity.In their approach, the continuity equation includes a term accounting for the flow between the two porous regions.Here, the state variables are the pressure in the fractured zone, u m [M L −1 T −2 ], and the pressure in the matrix region, u im,1 [M L −1 T −2 ].Also, we can identify β and β 1 with φ m C m and φ im C im , where C m and C im [L T 2 M −1 ] are the total compressibilities of the fracture and matrix blocks, respectively.The flow between fracture and matrix regions is controlled by the mass transfer rate coefficient α 1 given by where a [L −2 ] is a shape factor reflecting the geometry of the matrix blocks, k im [L 2 ] is the absolute permeability of this region and µ [M L −1 T −1 ] is the fluid viscosity.
For solute transport in rivers the governing Eqs. ( 5) and ( 6) are based on the conservation of mass principle for the stream and storage zone segments (Bencala and Walters, 1983;Lees et al, 2000;Zhang and Aral, 2004).In this context, u m [M L −3 ] is the in-stream solute concentration, u im,1 [M L −3 ] is the storage zone solute concentration, β represents the stream channel cross-sectional area A [L 2 ], and β 1 is the storage zone cross-sectional area A s [L 2 ].The first-order mass transfer rate coefficient is given by where α s [T −1 ] is the storage zone exchange coefficient.Similar identification of parameters can be done with other phenomena and models.Some of them are summarized in Table 1.

Applications
In order to illustrate the proposed formulation, we apply it to four problems.We first test its applicability to flow problems by simulating delayed yield from storage, which we verify by comparison to Boulton's (1955) analytical solution.We also simulate an actual pumping well test carried out at a nearly field scale (Martinez-Landa et al., 2004).Next, we compare our approach with the integro-differential approach, simulating the hypothetical radially convergent tracer test described by Alcolea et al. (2001).We then apply the present approach to solve two problems of radial flow to a pumping well (Haggerty and Gorelick, 1995).Finally, we show how our formulation can be used to simulate a river tracer experiment (Lees et al., 2000).Simulations were carried out with TRACONF (Carrera et al., 1993), a Fortran program for the simulation of water flow and solute transport through porous media.Boulton (1955) developed an analytical solution for unsteady radial flow allowing delayed yield from storage.This problem can be described by Eq. ( 1) with N=1, assuming q=0.We model a hypothetical pumping test, considering a transmissivity of T =0.01 m 2 s −1 , S m =0.001 and S im =0.1, and a pumping rate of 0.04 π m 3 s −1 .Initial head equals zero.We compare our non-local in time approach with the Boulton's solution for three values of the rate coefficient α=2.5×10 −6 , 10 −5 and 5×10 −5 s −1 .Boulton (1955) referred to α as an empirical constant.We considered a mesh of 208 nodes with a spacing size r increasing geometrically with a factor of 1.08.The integration scheme in time was semi-implicit with variable t.

Example of NE in flow phenomena
Figure 3 displays the evolution of heads, at a distance of R=51.6 m from the well, versus dimensionless time (t D =Tt/S m R 2 ).We can see that the numerical solution (solid line) obtained with the present approach matches the analytical solution (circles) obtained by Boulton (1955).The figure illustrates the influence of the rate coefficient α on the system behavior.When α is small (long response time), the system behaves for a long time as if there was no delayed yield.When α increases the system soon behaves as it is constituted only by one domain with the full storage.This is expected because for very large values of α, the water mass transfer between mobile and immobile zones is so fast that they tend to be at equilibrium.
Figure 4 shows the evolution of measured drawdown at a distance of r=0.043 m from the well.We can see that the present approach (solid line) reproduces the experimental data (circles).The figure also shows that an increase in the mass transfer rate coefficient or in the matrix storage coefficient, causes an increase in the fluxes of water from the fracture to the matrix.As expected, this leads to a delay in the drawdown, as shown in Fig. 4 (dotted and dashed lines).

Example of NE in solute transport
We model the convergent tracer test described by Alcolea et al. (2001).A tracer mass of 7.88 g is injected 8 m away from a well pumping 150 m 3 d −1 .The radius of the pumping well and the aquifer thickness are 0.2 m and 5 m, respectively.Porosity of the mobile domain is φ m =0.1.The immobile zone consists of layers of length L im =0.05 m and porosity φ im =0.045.The diffusion coefficient in the immobile zone was set to D im =0.001 m 2 d −1 .Alcolea et al. (2001) used TRANSIN code, which considers the integro-differential approach to solve matrix diffusion problems.Thus, we have calculated coefficients α's and β's according to Eq. ( 13) with N=50.A uniform mesh with a grid size r=0.005m (93 nodes) and a fully implicit integration scheme in time was used in the simulations.
Figure 5 shows the breakthrough curves at the pumping well simulated with the present approach and with TRANSIN, which compare quite well.The maximum relative error was 1%.illustrating the equivalence between the present approach and the integro-differential formulation (Alcolea et al., 2001).Haggerty and Gorelick (1995) presented a case of radial flow to a pumping well, in the context of PCE removal from the Borden sand aquifer under realistic pumping rates.To solve the governing equations, they expressed the MRMT model in dimensionless form and used a semianalytic method.Here we only give the main characteristics concerning with our approach, as the specifications of the problem are well described in their work (Haggerty and Gorelick, 1995).Two hypothetical case studies were considered: the remediation of a homogeneous aquifer (Borden sand) and the remediation of a hypothetical heterogeneous aquifer with a mixture of mass transfer processes.In both cases we used a mesh of 102 nodes with a spacing size r decreasing in such a way that all cells have the same volume and r 1 =1.39 m.The integration scheme in time was semi-implicit with variable t.

Example of NE in solute transport at field scale
The first case study simulates the cleanup history of a homogeneous aquifer where seven immobile zones are associated to a distribution of grain sizes.The grains are assumed to be spherical; the α j and β j coefficients are given in Table 2 of Haggerty and Gorelick (1995).As in their work, we have considered 50 exponential terms to adequately describe each immobile domain.Figure 6a displays the evolution of the PCE mass fraction remaining through an scenario of 500 d.We can see that the numerical solution (solid line) obtained with the present approach match very well the semi-analytical solution (dashed line) obtained by Haggerty and Gorelick (1995).
In the second numerical experiment, we simulate the removal of PCE from an heterogeneous aquifer.Here, heterogeneity arises from four immobile zones of different geometry (porous grains, grain aggregates, clay layers and clay pods) and two immobile zones characterized by a surface reaction (slow and fast reactions).The mass-transfer parameters are those appearing in Table 3 of Haggerty and Gorelick (1995).Note that these parameters are assumed locally heterogeneous, but they have the same distribution at all points in space.Once again, 50 terms were used to describe each geometry and only a single term for each reaction.Figure 6b shows predictions of the mass fraction of PCE remaining in the aquifer during a remediation scenario of 20 000 d. Again, the evolution of the mass fraction remaining obtained through the present approach (solid line) fit well the results predicted by the semi-analytical solution (dashed line).Therefore, the present model can reproduce the behavior of heterogeneous media characterized by different types and rates of mass transfer.

Example of NE due to transient storage in a river tracer experiment
Here we make use of the equivalence between the present approach and the transient storage model described in Sect.4.5 to model a river tracer experiment (Lees et al., 2000).
In one of their experiments, Lees et al. (2000) injected an amount of 10 kg of sodium chloride into a river upstream, and measured the tracer concentration over time at three sampling stations downstreams.They measured a constant discharge of 251 L s −1 .Figure 7 shows the breakthrough curves (BTC) at two reaches (reach B located at 140 m and reach C at 190 m downstream from the injection point) and the BTC's simulated with the present approach.To fit experimental BTC's we fixed the longitudinal dispersion coefficient at D m =0.64 m 2 s −1 .The estimated model parameters (Eq.28) were (A s /A=0.62;α s =0.0056 s −1 ) for reach B, and (A s /A=0.68;α s =0.0058 s −1 ) for reach C.
This example shows that the proposed formulation can reproduce the behavior of mass transfer between active flow and dead-zones in a river system.However, a certain knowledge of the distributions of the storage zone cross-sectional area, storage zone exchange coefficient, and dispersion coefficient could be used to improve the description of experimental data by the MRMT model.In such a case, the proposed approach could better describe the complex mass transfer processes attributed to the wide spectrum of dead zones in natural rivers and streams (Zhang et al., 2009).
flow and transport in heterogeneous media, chemical transport under kinetic adsorption in soils and groundwater and transport in streams and rivers.All of these approaches are equivalent to and can be formulated in terms of the proposed MRMT formulation.We have developed an easy-to-use numerical implementation of multi-rate mass transfer, which embeds existing formulations for non-locality in time.This numerical method is simple and physically consistent with and mathematically equivalent to other non-local in time flow and transport approaches such as integro-differential, CTRW and fractional derivatives formulations.The present approach avoids the spatial discretization of the immobile domain, because it solves state variables of that zone as explicit functions of the state variables in the mobile domain.
The approach also avoids the need to keep track of the past history of state variables at the mobile domain.That is, the proposed formulation yields an algorithm that is local in time to model non-local phenomena.The latter facilitates (i) physical interpretation of parameters and (ii) incorporation of non-linear phenomena, such as chemical reactions, that require local variables.The algorithm was implemented into a Fortran 90 module that is very efficient and accurate, as it involves an analytical solution for the mass balance equations in the immobile zones.Through a series of examples we have illustrated that the proposed formulation can describe a wide spectrum of phenomena displaying non-equilibrium behavior.Thus, the developed MRMT module may find wide application for the realistic modeling of environmental flow and transport problems.
The decomposition of a complex non-local single continuum domain into a simple multicontinuum domain opens a wealth of possibilities for modeling non-linear phenomena.We expect that the actual procedure for extending the proposed formulation to non-linear processes will be problem dependent.Willmann et al. (2009) or Donado et al. (2009) have done it for multicomponent reactive transport because this problem can be decomposed as (1) linear transport of conservative components and (2) non-linear local chemical computations.In this kind of problems, the localization implicit in MRMT is essential.Moreover, in this kind of problems, the immobile regions decomposition obtained for conservative (linear) transport remains unaltered when simulating reactive (non-linear) transport.However, this may not remain true in general.The MRMT formulation can be easily extended for modeling non-linear phenomena such as multiphase flow or rainfall-runoff, but whether the immobile regions decomposition can be left unaltered remains to be proven.
The module has been tested by comparison with published solutions and is publicly available at http://www.h2ogeo.upc.es/English/English/software.htm#Modprocess MRMT.

Fig. 1 .
Fig. 1.Typical response to (a) a constant injection test in a double porosity formation and (b) a tracer test breakthrough curve.In both cases, the system reacts with small storage capacity at early times, but much larger at late times, suggesting that some time is needed before the full storage capacity is mobilized.
are added to the mass balance equations.It is the form of such terms what sets different non-local formulations apart.Their number is too long to list, but the most widely used have been:

Fig. 2 .
Fig. 2. Hypothetical numerical discretization of the mobile and immobile domains.
[L T ] is a molecular diffusion coefficient in the jth immobile region; Rim,j[-]  is the retardation factor of the jth immobile zone; Lim,j [L] is the distance from the mobile zone to the jth immobile zone; is the hydraulic conductivity in the jth immobile region; Sim,j [L -1 ] is the specific storage coefficient of the jth immobile zone; Lim,j [L] is the distance from the mobile zone to the jth immobile zone.and first-order mass transfer coefficients used byHaggerty and Gorelick (1995); β is the mobile capacity.

Fig. 5 .
Fig.5.Breakthrough curve of tracer in a radial convergent test, illustrating the equivalence between the present approach and the integro-differential formulation(Alcolea et al., 2001).