Advective Trapping in the Flow Through Composite Heterogeneous Porous Media

We study the mechanisms of advective trapping in composite porous media that consist of circular inclusions of distributed hydraulic conductivity embedded in a high conductivity matrix. Advective trapping occurs when solute enters low velocity regions in the media. Transport is analyzed in terms of breakthrough curves measured at the outlet of the system. The curve’s peak behavior depends on the volume fraction occupied by the inclusions, while the tail behavior depends on the distribution of hydraulic conductivity values. In order to quantify the observed behaviors, we derive two equivalent upscaled transport models. First, we derive a Lagrangian trapping model using the continuous-time random walk framework that is parameterized in terms of volume fraction and the distribution of conductivities in the inclusions. Second, we establish a non-local partial differential equation for the mobile solute concentration by volume averaging of the microscale transport equation. We show the equivalence between the two models as well as (first-order) multirate mass transfer models. The upscaled approach parameterized by medium and flow properties captures all features of the observed solute breakthrough curves and sheds new light on the modeling of advective trapping in heterogeneous media.


Introduction
Transport of solutes in the subsurface is a relevant topic for many environmental engineering problems, such as managing groundwater quality, storage of substances in the subsurface, energy production or subsurface remediation. Modeling of the transport process is needed to make predictions, to design environmental or industrial operations, or to interpret observations made at a particular site. One of the main difficulties for models is the heterogeneous structure of the subsurface. Hydraulic and transport parameter contrasts can be very high. This is in particular true for fractured media, but also holds for porous media. This poses a challenge for models, as the parameter structure is typically not known, nor is it feasible to resolve it in a model. For this reason, upscaled models are used, where the heterogeneous structure is only partly, or not at all, resolved.
For transport of non-reacting substances in heterogeneous porous media, advection and diffusion or dispersion are the relevant local-scale transport processes and the mass balance for a substance leads to the advection-dispersion equation (ADE). Upscaled models for such transport processes are often formulated in a way that they keep the structure of an ADE, while the (constant) effective transport parameters are defined according to the distribution of medium heterogeneity (Brenner and Edwards 1993). In particular, for highly heterogeneous media, such approaches have been questioned. Heterogeneous structures can lead to a strong variability of transport times, so that spatial concentration distributions or breakthrough curves are poorly represented by those generated with an ADE characterized by constant upscaled transport parameters. As an example, breakthrough curves can show a strong tailing, which is not reproducible with ADE-based solutions. In the last decades, alternative modeling approaches have been developed to account for non-Fickian transport behavior (Frippiat and Holeyman 2008;Dentz et al. 2011;Fiori et al. 2015;Lu et al. 2018;Pedretti and Bianchi 2018), such as multirate mass transfer (MRMT) and memory function models (Haggerty and Gorelick 1995;Carrera et al. 1998;Haggerty et al. 2000), continuous-time random walk (CTRW) models (Berkowitz et al. 2006;Cvetkovic et al. 2014;Noetinger et al. 2016;Comolli et al. 2016) or fractional derivative models (Cushman and Ginn 2000;Benson et al. 2000Benson et al. , 2013Meerschaert and Sikorskii 2011). The parameterization of such models and how to relate them to structural information of the subsurface or material properties have been a field of intense research (e.g., Zinn et al. 2004;Fiori et al. 2007;Willmann et al. 2008;Edery et al. 2014;Zhang et al. 2014;Cvetkovic et al. 2014;Jankovic et al. 2017;Comolli et al. 2019).
In this paper, we focus on the upscaling of transport in composite porous media with sharp material contrast that consist of low-conductivity inclusions embedded in a highly conductive background. This scenario is motivated by fractured porous media, in which the inclusions represent the low conductivity matrix and the connected high conductivity background the fractures. Non-Fickian transport due to diffusive mass exchange between fast transport paths and stagnant regions has been upscaled in the framework of matrix-diffusion or multirate mass transfer approaches (Maloszewski and Zuber 1985;Haggerty and Gorelick 1995;Carrera et al. 1998;Cvetkovic et al. 1999;Gouze et al. 2008;Zhang et al. 2014;Hyman et al. 2019). In these approaches, the impact of mass transfer between a mobile and immobile continua is encoded in memory functions that can in principle be related to the geometry and material properties of the stagnant regions. The MRMT approach was used by Willmann et al. (2008) for the effective modeling of the combined effect of advective-dispersive transport in highly heterogeneous media in terms of empirical memory functions. However, the mechanisms and parameterization of large-scale transport due to advective mass transfer between a highly conductive background and low conductivity inclusions are still open issues. Rubin (1995); Eames and Bush (1999) and  derived macrodispersion coefficients for composite porous media with a bimodal parameter distribution and isotropic and anisotropic inclusions. Zinn et al. (2004) studied experimentally the mechanisms of diffusive and advective trapping in bimodal composite media consisting of circular inclusions with different conductivity contrasts. These authors observed the breakdown of the macrodispersion concept for media with intermediate and high parameter contrasts, and breakthrough curve tailing due to diffusion and slow advection within the inclusions. Hidalgo et al. (2020) derived a memory function model for bimodal composite media that accounts for advective trapping in the inclusions and fast transport in the background.  proposed multi-indicator media consisting of inclusions of different sizes and hydraulic conductivities as models for the conductivity distributions in aquifer and derived expressions for longitudinal dispersivities. Based on this conceptual medium representation, Fiori et al. 2006Fiori et al. , 2007 developed a semi-analytical Lagrangian model for the distribution of solute travel times in highly heterogeneous porous media, which where implemented in a time-domain random walk framework (Cvetkovic et al. 2014) as the multi-indicator self-consistent approximation (MIM-SCA) model. These authors found anomalous transport and strong breakthrough curve tailing due to broad distributions of hydraulic conductivity and explored the relation to CTRW models. The CTRW framework was used by Tyukhova et al. (2016) to model strong tailing of solute breakthrough curves in composite media with random conductivity distributions. These approaches are based on Lagrangian transport formulations in terms of distributions of solute travel times that can be related to the conductivity distribution, while MRMT approaches provide non-local partial differential equations for the mobile and total solute concentrations based on memory functions whose relation to advective mass transfer and medium heterogeneity is elusive.
Here, we study transport in composite media characterized by sharp conductivity contrast between the inclusions and background media. The inclusions are characterized by a distribution of hydraulic conductivities that are much lower than the background conductivity. We generalize the approach presented in Hidalgo et al. (2020) for a bimodal medium. While the medium is similar to the one considered by Fiori et al. (2006), here the main source of heterogeneity comes from the conductivity contrast between the matrix and inclusions. Our upscaling strategy is similar to the one adopted by Fiori et al. (2006) in the sense that large-scale transport is quantified in terms of (advective) transition times over characteristic medium length scales. However, based on the strong contrast between inclusions and background, our upscaling approach combines the mobile-immobile mass transfer and CTRW frameworks. In this approach, mass transfer between the background and the inclusions is modeled as a Poissonian trapping process with trapping times given by the velocity and size of the inclusions. Based on this, we establish an upscaled Lagrangian trapping model from first principles, from which we derive a non-local partial differential equation for the mobile solute concentration. The memory kernel and trapping rates are fully constrained by the medium structure and conductivity distribution, which yields a predictive large-scale transport model. Furthermore, we provide an alternative derivation of the large-scale model using volume averaging and show the equivalence to (first-order) MRMT models.
Section 2 defines the flow and transport model, the model medium and the numerical study cases. Section 3 provides the derivation of the upscaled transport model and Sect. 4 on the comparison of the upscaled model with detailed numerical simulations. Section 5 provides conclusions and outlook for the application of the derived upscaling approach for more general heterogeneity scenarios and multiphase flow conditions.

Background and Methodology
We consider a 2D medium of size L x × L y formed by a homogeneous material of conductivity K m in which disconnected randomly arranged inclusions of circular shape and radius r 0 are embedded. The inclusion arrangement is generated by randomly choosing the centers coordinates from an uniform distribution and discarding inclusions that overlap previously existing ones until the desired volume fraction is covered. The inclusions are homogeneous, but their conductivity is randomly distributed. We will consider log-normal and a gamma distributions. The distribution parameters are chosen so that the mean conductivity of the inclusions is always smaller than the conductivity of the matrix (see Fig. 1).

Flow and Transport Model
We solve the incompressible steady-state flow equation where is the Darcy velocity, K( ) is the hydraulic conductivity and h is the piezometric head. A constant flow rate v 0 is prescribed on the left boundary and zero head on the right one so that the mean flow is aligned with the x axis. The top and bottom boundaries are periodic.
We consider purely advective transport governed by is assumed to be constant and set equal to one. Initially, the solute is uniformly distributed along a line, i.e., c( , t = 0) = c 0 (x) . Zero gradient concentration is considered on the left and right boundaries and periodic boundary conditions on the top and bottom ones.
Transport is solved in a Lagrangian framework. The equation of motion for the position (t; ) of a fluid particle is with (t = 0; ) = . Initially, the particles are uniformly distributed along the left boundary.

Numerical Simulations
Flow and transport problems are solved numerically. Flow (1) is solved using a two-point flux finite volume scheme, which gives the velocity on the cells' sides. The domain is discretized using a regular mesh of square cells of side Δ = r 0 ∕30 . The mesh size was chosen to reproduce the circular shape of the inclusions with enough accuracy. Advective transport (3, 4) is solved using the semi-analytical method of Pollock (1988). The initial solute distribution is represented by N p = 10 6 particles uniformly distributed along the left boundary, which is equivalent to a flux-weighted injection because flow at the inlet is uniform and the streamlines parallel to the x direction. This is achieved by adding a buffer of length ⌈2r 0 ⌉ between the boundary and the first inclusion. The simulation runs until all particles leave the domain and the arrival times at the right boundary are recorded subtracting the travel time through the buffer.
We present the results using dimensionless magnitudes. Length is scaled using the domain height L y , velocity with the prescribed velocity v 0 = q 0 ∕ and conductivity with the background conductivity. The characteristic time is t c = L x ∕v 0 . In dimensionless units, the domain used for the simulation has dimensions 20 × 1 , and the inclusion size is r 0 = 0.0359 . The contrast between the background conductivity and the mean conductivity of the inclusions is ⟨K i ⟩∕K m = 0.01 , where ⟨K i ⟩ is the mean conductivity of the inclusions. We consider two distributions for K i , log-normal and gamma The log-normal distribution is characterized by = 0.01, 2 = 2.3 and the gamma distribution by k = 0.8 and = 0.02638 . The distribution parameters are chosen so that the conductivity fields have equivalent log-statistics regardless of the distribution. The chosen variance corresponds to a mildly heterogeneous distribution of conductivity of the inclusions. However, the main source of heterogeneity in the systems comes from the contrast between inclusions and the matrix. The inclusions volume fraction is 0.1 ≤ ≤ 0.4 . An ensemble of five realizations was considered for all the reported cases.

Transport Upscaling
The upscaling of transport is based on the following rationale. The flow structure is organized into a connected high flow background and embedded low flow inclusions. As particles move through this structure, they follow tortuous flow paths in the background along which they encounter inclusions. Motion in the background is modeled here by translation with the mean background velocity v m , while the impact of flow variability is accounted for by the dispersion coefficient D m . The number of inclusions encountered after a travel distance is modeled as a Poissonian random variable characterized by the spatial rate , The upscaling approach relies on the representation of trapping as a Poisson process, which is valid for uniform random media characterized by a finite mean distance between inclusions. The random structure of the composite media under consideration here is indeed characterized by approximately exponential distributions of distances between inclusions.
Each time a particle gets advectively trapped inside an inclusion, it travels over the small average distance 0 = r 0 ∕2 during a large time = 0 ∕v e , where v e ≪ v m is the average speed inside the inclusion. That is, the particle is essentially retained at the position of the inclusion. The volumetrically sampled Eulerian distribution of average speeds v across the ensemble of inclusions in the following is denoted by p e (v) . The probability to encounter an inclusion of high velocity v h is higher than encountering an inclusion of low velocity v because the stream tube that leads into the high velocity inclusion is wider than the one leading into a low velocity inclusion by a factor v h ∕v as a result of volume conservation. Thus, the distribution of inclusion velocities encountered by a particle is obtained from the Eulerian p e (v) by weighting with the average inclusion speed, where ⟨v e ⟩ denotes the mean of p e (v) . The distribution (t) of residence times = 0 ∕v e can be expressed in terms of p v (v) as Based on these fundamental mechanisms, we derive in the following an upscaled transport model using the CTRW framework.

Model Derivation
In the picture described above, particle motion in the background is given by where (s) is a Gaussian white noise, v m is the mean velocity in the matrix, and D m is the longitudinal dispersion coefficient, which results from velocity fluctuations in the highly conductive matrix. The time the particle spends mobile in the matrix is denoted by s. The actual clock time t describes a compound Poisson process (Feller 1968;Benson and Meerschaert 2009)

and is given by
where n s is the number of trapping events, distributed according to (7) with = v m s the mean distance covered in the matrix after the mobile time s. In order to derive the evolution equation for the particle density c(x, t), we discretize (10) and (11) as where x j ≡ x(s j ) , s j = jΔs , and the j are independent unit Gaussian random variables. Thus, the distribution Λ(x) of transition lengths Δx = x j+1 − x j is given by a Gaussian distribution with mean v m Δs and variance 2D m Δs . The distribution Δ (t) of transition times Δt = t j+1 − t j is given by a compound Poisson distribution, which can be written in Laplace space as The Laplace transform is defined in Abramowitz and Stegun (1972). Laplace transformed quantities in the following are marked by an asterisk, and the Laplace variable is denoted by . Note that (12) describes a continuous-time random walk. Thus, we can write a generalized master equation (Berkowitz et al. 2006) for c(x, t), which reads in Laplace space as A Kramers-Moyal expansion (Risken 1996) of the generalized master equation gives where c 0 (x) is the initial particle distribution. The limit Δs → 0 gives where we used the expansion of expression (13) for small Δs . Furthermore, we defined the memory function where ⟨ ⟩ is the mean residence time in the inclusions. This expression for the memory function reads in time space as It is proportional to the probability that the residence in the inclusion is larger than t. We now define the mobile concentration c m (x, t) through its Laplace transform as With this definition, we can write (16)

Trapping Rate
The trapping rate can be constrained by considering the ratio between the total accumulated mass in the inclusions and in the matrix, which is obtained from (21) for = 0 as where we used that by definition * (0) = 1 . By noting that the ratio between total immobile and mobile masses is equal to the ratios of the volume fractions of inclusions and the volume fraction (1 − ) of the matrix, we obtain where we used that ⟨ ⟩ = 0 ∕⟨v e ⟩ . This implies that = ∕(1 − ).

Velocity Distribution
The distribution of velocities in the inclusions can be estimated based on the analytical expression for the velocity v in in an isolated single inclusion in an infinite domain, which is given by (Wheatcraft and Winterberg 1985;Eames and Bush 1999) where K r = K i ∕K m is the conductivity ration between the inclusion and the matrix, and v m is the background velocity. Since here K r ≪ 1 we can further approximate the velocities in the inclusion as The background velocity v m is estimated from the area fraction covered by the inclusions as The linear relation (27) between the velocity in the inclusion and the conductivity allows to estimate the velocity distribution from the conductivity distribution as This expression provides a direct link between velocity and conductivity distribution.

Results
In the previous section, we have derived a CTRW model described by equations (10) and (11), and we showed that this is equivalent to a memory function model given by equations (22) and (23). In this section, we study the parameterization of these and compare the model data to the results of the direct numerical simulations.

Trapping Model and Velocity Distributions
As outlined above, the input for the transport model is, apart from the parameters describing transport in the background material, the velocity distribution inside of the inclusions p e (v) and the trapping rate . This information is needed for a full parameterization of the model. To be able to estimate the parameters in advance or to have means to measure them by calibration experiments is crucial for the applicability of the model. For this reason, the parameters for different test cases of the setup studied here are shown and shortly discussed in the following. As shown in Fig. 2, the Poisson distribution describes the distribution of trapping events well for both the log-normal and gamma distributions of hydraulic conductivity under different volume fractions. Figure 3 shows that the trapping rate is well estimated by expression (24). Figure 4 compares the analytical estimate (29) for p e (v) with the data from the numerical simulations for the scenarios of lognormal and Gamma distributed conductivities with = 0.4 . The analytical expression provides a solid estimate, also for the other scenarios under consideration (not shown).

Large-Scale Transport
We describe the large-scale transport in terms of particle breakthrough curves measured at the outlet of the domain (Fig. 5). The curves display similar features for both conductivity distribution of the inclusions. The time of the peak arrival is determined by the mean velocity in the matrix, which is proportional to the volume fraction occupied by the inclusions (see (28)). That is, we observe earlier arrival times for higher . The volume fraction is also responsible for the width of the peak. As more volume is covered by the inclusions, streamlines become more tortuous and the number of trapping events increases. The dispersive effect of the heterogeneity is therefore enhanced, and the curves broaden. The tail of the curves depends mainly on the conductivity  Fig. 3 Trapping rate for K log-normally (dots) and gamma (squares) distributed. The trapping rate is proportional to the volume fraction occupied by the inclusions and the ratio between the mean velocity in the inclusions and in the matrix as given by equation (24) distribution. The height of the tail values depends on the volume fraction through the trapping rate . We provide three solution methods of the upscaled transport approach.
(1) A numerical particle tracking approach that implements the upscaled trapping CTRW (see Appendix A.1 for implementation details); (2) a semi-analytical solution based on the explicit Laplace space solution of the non-local partial differential equation for the mobile solute concentration given by (22) and (23)  As expected, all three solution methods provide the same results. They reproduce all salient features of the observed breakthrough curves for different volume fractions of the inclusions (Fig. 6). The peak arrival times and widths are well represented by the mean velocity in the matrix and the dispersion coefficient for impervious inclusions D m = 0.74 r 0 v 0 given by Eames and Bush (1999), where v 0 is the flow rate prescribed at the inlet.
The Poisson distributed trapping events and the velocity distribution in the inclusions are responsible for the late time behavior, that is, the slope of the curves. Only the very early arrival of concentration is not captured by the CTRW model. This corresponds to particles that follow high-velocity streamlines and experience none or very low number of trapping events, whose probability is underestimated by the Poisson distribution (see Fig. 2). This effect is expected to become less important as the size of the domain

Fig. 5
Breakthrough curves measured at the outlet of the domain for (a) K log-normally distributed and (b) K gamma distributed for different volume fraction covered by the inclusions. Results correspond to an ensemble of five realization for each case. The shape of the breakthrough curves is similar for all cases. The slope of the tail of the curves is independent of and controlled by the conductivity distribution. The inclusions have a dispersive effect on the curve, whose peak broadens as increases. However, the time for the arrival of the peak is inversely proportional to .
increases and the probability of not encountering a low conductivity inclusion goes to zero.

Conclusions and Outlook
In this paper, we study the upscaling of advective transport in heterogeneous porous media with a background-inclusion pattern, characterized by circular inclusions with a distribution of conductivity values. We derived a continuous-time random walk model and a memory function model for advective transport in such media. The models were applied to test cases of log-normal and gamma distributed conductivities in the inclusions. We consider a two-dimensional medium, but the upscaling approach can be applied also to three dimensions because the dominant transport mechanisms, retention of particles in low conductivity inclusions and fast transport in the background are the same. A crucial issue of upscaled models is the parameterization. We sought here parametrizations in terms of transport-independent attributes based on material properties and heterogeneity structure only. This means, the models can be parameterized a priori, that is, before the transport experiment or simulation. The upscaled models derived here have as input the trapping time distribution, which can directly be related to the distribution of the velocities inside of the inclusions, as well as the mean trapping rate. For media under consideration here, which are characterized by a strong conductivity contrast between background and inclusions, the velocity distribution is directly related to the conductivity distribution and the mean velocity in the background, which in turn is given in terms of the volume fraction of inclusions.
It needs to be noted that the representation of heterogeneity as circular inclusions of constant size is an idealization. Real media are characterized by a distribution of length scales, which can be captured by introducing a distribution of inclusion sizes. Also, anisotropy as discussed in Eames and Bush (1999) may be a relevant parameter for heterogeneity that needs to be addressed. The presented model provides a good approximation for transport in media characterized by a narrow distribution of length scales. The proposed upscaling approach can be used also for more complex medium geometries as long as there is a large enough contrast between background and inclusions. Also, it needs to be mentioned that upscaled models for advective transport are only a part in the whole picture for transport in highly heterogeneous media. Upscaled models for diffusive mass exchange and for advective trapping can be considered as the two extreme cases of trapping mechanisms. In reality, it is expected that both processes play a role for the mass exchange between fast transport paths and slow regions. How to combine this is an open question.
As an outlook, the upscaling strategy presented here can be also applied to other transport processes. Memory function models have been used for predictions of immiscible displacement in highly heterogeneous media [e.g., Kazemi et al. 1992;Arbogast 1992;Di Donato et al. 2007;Tecklenburg et al. 2013;Spooner et al. 2021]. The modeling principles are similar to those developed for solute transport, but memory functions usually have to be approximated. Although such models are mostly discussed for flow processes for which mass exchange between fast flow paths and stagnant zones can clearly be categorized as capillary driven (Schmid and Geiger 2012;Tecklenburg et al. 2016), it is clear that mass exchange due to buoyancy and viscous flow is also relevant. For mass exchange due to buoyancy, memory function models or equivalent multi-porosity models have been developed (Di Donato et al. 2007). The case of viscous flow in highly heterogeneous media could be considered as a third mass exchange process of fluids that could be upscaled to a memory function model. This case is relevant for the sweep efficiency during oil recovery. Often, the process is well represented as piston-type displacement. The linearized problem would be an advection problem. With this argumentation, an upscaled model for advective transport can be seen as a blueprint for the equivalent upscaled model for viscous dominated two-phase flow in heterogeneous media, restricted for the case of piston-type immiscible displacement. The upscaled models that are derived here for advective transport could be extended to such cases in a straightforward manner.

Appendix B Memory Function Model
In this Appendix, we first provide an alternative derivation of the memory function model (23) that allows to constrain the trapping rate . Then, we give the Laplace space solution for the breakthrough curves and describe the implementation of the memory function model as a first-order multirate mass transfer model.

B.1 Derivation by Volume Averaging
The upscaled transport model (23) with (22) can be derived with the method of volume averaging (Whitaker 1999). For this purpose, the domain is divided into mobile (background) and immobile (inclusions) regions with the volumes V m and V im and nonaveraged variables c m and c im . The transport equation (3) in the Eulerian framework is Using the spatial averaging theorem (Whitaker 1999), we obtain The volume V refers here to the total averaging volume containing both mobile and immobile parts and the angular brackets to the superficial average, where an integral over the full averaging volume is carried out. Notations and definition can be found in the textbook by Whitaker (1999). Reformulating to intrinsic variables, where the average is only carried out over the volume of the mobile (or immobile) domain, this reads where the angular brackets with the index m refer to the volume average over the mobile domain only. The flux term over the inner boundaries can be reformulated using continuity of flux over the boundaries Using and 1 − for the volume fractions of the immobile and the mobile part and c m and c im for the intrinsic volume averages of the concentrations in the mobile and immobile part, one obtains The intrinsic volume average of the mass flux density is here not derived, but one gets an advective part with the intrinsic average of the velocity (here aligned in x-direction) and a macrodispersion term due to the tortuous velocity field (Whitaker 1999). For a line source, this leads to the transport equation (23). Note that the velocity field in the mobile domain alone is not divergence free at the internal interfaces. We make here the assumption that we have a large enough averaging volume such that the divergence in the average is zero.
Assuming that the concentration in the mobile domain is uniformly distributed over the inflow boundary of an inclusion, the intrinsic average of the immobile concentration can be interpreted as the time convolution of the averaged Green's function and the concentration in the mobile domain, which leads to the memory function formulation given by (22).

B.2 Solution in Laplace Space
The breakthrough curve at position x is given by solving (20) in a semi-infinite domain (Dentz and Berkowitz 2003), where (21) is used to substitute c im , with initial conditions and boundary conditions In Laplace space, the solution is where * ( ) is the Laplace transform of the memory function given by (B19) and (B20). The solution in real space is obtained by calculating numerically the inverse Laplace transform of (B13) using the de Hoog algorithm (de Hoog et al. 1982).

B.3 First-Order Multirate Mass Transfer Model
In order to solve the memory function model (23), we formulate it as a first-order multirate mass transfer model, which approximates the memory function (t) as an exponential series (Haggerty and Gorelick 1995;Silva et al. 2009) Thus, we can write the immobile concentration as where c (i) im (x, t) is defined by Thus, the memory function model (23) can be written as the following system of equations because (B16) is the solution of (B17b). The system (B17b) can be solved in a straightforward way by grid-based methods. Here, we employ an explicit in time finite difference method in which the dispersion term is discretized with a second-order scheme and the advection with an upwind scheme.

B.4 Memory Function
The memory function is defined by (17)  where Φ(z) is the cumulative unit Gaussian distribution (zero mean and unit variance), is the mean of ln(v) , and 2 is the variance of ln(v). The exponential expansion of the memory function (B14) used in the multirate mass transfer model has positive frequency coefficients a i which are evenly distributed in a logscale within the chosen frequency range. The amplitude coefficients b i are calculated using a non-linear least squares fitting algorithm. The frequency interval and number of terms are chosen so that the memory function is well reproduced (see Table 1 for the cases in Fig. 6