# Transient Dispersion Regimes in Heterogeneous Porous Media: On the Impact of Spatial Heterogeneity in Permeability and Exchange Kinetics in Mobile-Immobile Transport

Laurent Talon, Emma Ollivier-Triquet, Marco Dentz, Daniela Bauer

► To cite this version:

Laurent Talon, Emma Ollivier-Triquet, Marco Dentz, Daniela Bauer. Transient Dispersion Regimes in Heterogeneous Porous Media: On the Impact of Spatial Heterogeneity in Permeability and Exchange Kinetics in Mobile-Immobile Transport. Advances in Water Resources, 2023, 2023, 174, pp.104425. 10.1016/j.advwatres.2023.104425. hal-04071643

HAL Id: hal-04071643

https://ifp.hal.science/hal-04071643

Submitted on 17 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Transient dispersion regimes in heterogeneous porous media: On the impact of spatial heterogeneity in permeability and exchange kinetics in mobile-immobile transport

Laurent Talon<sup>a,*</sup>, Emma Ollivier-Triquet<sup>a,b</sup>, Marco Dentz<sup>c,*</sup>, Daniela Bauer<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Dispersion

Mobile-Immobile model

Transient regimes

Lattice Boltzmann

CTRW

Heterogeneity

## ABSTRACT

Transport processes in the subsurface are coupled with the heterogeneity of the porous structure. Obtaining an accurate description of mass transport in such media at all time scales remains a crucial task, as the subsurface is strongly impacted by human activity. Spreading of pollutants in the transient but also in the asymptotic regime, as well as the time to reach a critical location (e.g. aquifer, well) significantly depend on the underlying heterogeneity of the permeability field. Moreover, in the case of solute retention, transport is also characterized by local exchange kinetics that depend on the local aquifer properties. Consequently, exchange (retention) times are expected to be spatially heterogeneous. In this work we focus on the influence of spatially heterogeneous permeability and exchange times on the transient and asymptotic transport regimes and provide a parametric study. To this goal, we simulate in a first part the transport in a two-dimensional heterogeneous medium under spatially varying permeability and mobile-immobile mass transfer parameters. Equations are solved using a Lattice-Boltzmann two-relaxation-time (TRT) algorithm. We assume the following relation between the local permeability *K* and the local exchange time: $\tau \propto K^\gamma$. Taking into account this relation, we investigate the impact of the Damköhler number (*Da*, ratio of the advection and exchange time scales), the disorder of the permeability field and the value of the exponent of the coupling function (*γ*) on the spatial evolution of the concentration field and breakthrough curve. We show that, depending on the parameters (*Da*, *γ*, etc.), we can observe transient, non-Fickian dispersion, which is characterized by non-exponential tails of the solute breakthrough curves and a non-linear evolution of the spatial variance of the solute distribution.

In a second part, we present a new continuous time random walk (CTRW) model to upscale these transport behaviors. The model is based on a spatial Markov model for particle velocities that couples advective-dispersive transport and heterogeneous mass transfer through a compound Poisson process. The upscaled model can be fully parameterized by the statistics of permeability and the hydraulic gradient (with no fitting parameter), that is, in terms of medium and flow characteristics. The results of the CTRW model fully capture the non-Fickian transient transport regimes both for the breakthrough curves and spatial concentration variance. In the longtime limit, as expected from the central limit theorem, the CTRW model predicts normal, Fickian behavior. Finally we show, that the time to reach the asymptotic regimes in heterogeneous media (e.g. heterogeneous exchange times) is parameter dependent and in average is two orders of magnitude larger than for the respective homogeneous case. To summarize, the coupling between the heterogeneous permeability field and the local mass transfer properties can strongly influence transient and asymptotic transport regimes and potentially explain experimentally observed non-Gaussian behaviors.

## 1. Introduction

Pollutant transport in the subsurface is coupled to the water content and the heterogeneity of the porous medium. As the subsurface

is strongly influenced by human activity, an accurate description of pollutant transport at all time scales still remains a crucial task. Indeed, for example, the knowledge of such transport is required for efficient

<footnote><sup>a</sup> Université Paris-Saclay, CNRS, FAST, 91405, Orsay, France</footnote>

<footnote><sup>b</sup> IFP Energies Nouvelles, 1 & 4 Av. Bois Préau, 92852 Rueil Malmaison, France</footnote>

<footnote><sup>c</sup> Institute of Environmental Assessment and Water Research (IDAEA-CSIC), C. Jordi Girona, 18-26 08034 Barcelona, Spain</footnote>

<footnote>* Corresponding authors. E-mail addresses: talon@fast.u-psud.fr (L. Talon), marco.dentz@csic.es (M. Dentz).</footnote>
water remediation and for providing access to safe drinking water. For many pollutants (nuclear waste, pharmaceuticals, etc.), contamination remains critical at very low concentrations and cannot be ignored. For this type of pollutants, long-lasting low concentrations represent therefore a major issue in aquifer and soil remediation. Further compaction on these transient, long-lasting regimes and their dependence on the porous structure is therefore crucial in the actual environmental context.

Transport in porous media has been treated in the literature for a long time but still remains an actual topic. While much research focused mainly on the asymptotic regimes, further attention should be paid to the transient and often long-lasting pre-asymptotic regimes. Indeed, experimental field and laboratory observations report that solute transport may not obey classical Fickian dispersion. Non-Fickian transport can be attributed to the heterogeneity of the permeability field, but also to sorption phenomena, that may be space dependent.

A typical example is the reactive transport experiment performed at Cape Cod, where reactive tracer exhibit non-Fickian transport behavior (Brusseau and Srivastava, 1999). Despite the fact that the heterogeneity of the permeability field was rather small $\sigma_{\text{logk}}^2 \approx 0.24$, the velocity of the tracer front was decreasing over time and the variance of the plume spread more rapidly than linearly. Brusseau and Srivastava investigated various mechanism to model this anomalous behavior. They tried different model scenarios where permeability was considered to be either homogeneous or heterogeneous, sorption could be linear or not with either instantaneous or rate limited kinetics. Results of the study show that the early time behavior of the experiment could be reasonably well fitted. However, in contrary, the later times could be fitted only by significantly increasing the heterogeneity of the sorption kinetic coefficient in the second half of the domain. The latter result points out the importance of the choice of the model parameters to describe the long time transport behaviors but also the need to carry out a parametric study.

Dynamics of spreading and the spatio-temporal extension of the pollutants in the transient but also in the asymptotic regime as well as the time to reach a critical location (e.g. aquifer, well, housing) strongly depend on the underlying heterogeneity of the porous medium. This dependence has been pointed out by Dai et al. (2004), Ren et al. (2022) by performing transport simulations in porous media characterized by different criteria (anisotropy, correlation length, scales). However, transport also depends on the spatial variability of exchange kinetics between mobile and immobile solute phases. First-order exchange kinetics have been used to quantify diffusive mass transfer between high and low velocity zones, as well as linear kinetic sorption-desorption reactions (Haggerty and Gorelick, 1995), as most pollutants are subject to sorption phenomena.

The impact of heterogeneity on large scale transport has been largely treated in the literature and quantified in the framework of stochastic modeling for the spatially fluctuating exchange parameters and permeability (Rubin, 2003). Numerical simulations (Selroos and Cvetkovic, 1992; Tompson, 1993; Selroos and Cvetkovic, 1994; Burr et al., 1994) were used to determine the ensemble averaged center of mass velocity and dispersion coefficients or temporal moments. Analytical studies have used first-order perturbation theory in the heterogeneity variance to investigate large scale transport under heterogeneous conductivity and homogeneous kinetic mass transfer (Quinodoz and Valocchi, 1993; Dagan and Cvetkovic, 1993; Massabó et al., 2008; Soltanian et al., 2015), under heterogeneous conductivity and mass transfer properties (Cvetkovic and Shapiro, 1990; Cvetkovic et al., 1998; Miralles-Wilhelm and Gelhar, 1996; Rajaram, 1997; Dentz et al., 2000), and under homogeneous conductivity and heterogeneous sorption properties (Metzger et al., 1996; Reichle et al., 1998; Attinger et al., 1999). However, it is important to stress that all these studies focus on the asymptotic stationary regime and aim at the quantification of effective upscaled transport parameters, such as effective dispersion

coefficients and effective retardation coefficients. Note that these transport parameters are only applicable when the average concentration follows a Fickian equation, which is theoretically expected at very long times. However, as we will see later, before reaching this asymptotic limit, transport may exhibit anomalous behavior. In consequence, it suggests that the behaviors observed, even at the longest times in the Cape Cod experiment, correspond to transient regimes. This is particularly supported by the fact that the velocity of the tracer front decreases with time without reaching a constant value.

Mobile-Immobile models (Coats and Smith, 1964; van Genuchten and Wierenga, 1976; Haggerty and Gorelick, 1995) have been used to describe the impact of kinetic mass transfer between mobile and immobile medium portions on the overall solute migration in saturated and partially saturated media. Moreover, as pointed out by several researcher (e.g., Nkedi-Kizza et al., 1984; Valocchi, 1990; Haggerty and Gorelick, 1995) all first-order kinetic mass transfer models can also be applied to linear sorption problems.

A large number of experiments have been performed to characterize exchange times for transport in porous media with different degrees of heterogeneity. Griffioen et al. (1998) provided a review of exchange times of a large number of laboratory experiments performed in different column samples. Low exchange times (5 h and more) have been observed by Rao et al. (1980) in samples consisting of porous spheres and glass beads. In contrast, Smettem (1984) measured tracer transport in heterogeneous soil columns and determined an average exchange time up to 125 days. Exchange times were computed by fitting a Mobile-Immobile model to the experimental data. However, a major difference between the laboratory and the field scale lies in the fact that the latter presents often spatial variations (e.g., Tompson, 1993; Burr et al., 1994; Brusseau, 1994). Consequently, the permeability but also the exchange time might present large spatial heterogeneity. Here, we consider both heterogeneity in permeability and mobile-immobile mass transfer, which manifest in spatial variability of the flow velocity $\bar{v}(\bar{r})$ and the local mass transfer times $\tau(\bar{r})$, where $\bar{r}$ is the position vector. Transport of the mobile and immobile concentrations $C_m$ and $C_{im}$ can then be modeled as (e.g., Cvetkovic and Shapiro, 1990; Burr et al., 1994; Haggerty and Gorelick, 1995; Reichle et al., 1998)

$$ \theta_m \frac{\partial C_m}{\partial T} + \bar{u}(\bar{r}) \cdot \nabla C_m - D_0 \nabla^2 C_m = - \frac{1}{\tau(\bar{r})} (C_m - C_{im}) \quad (1) $$

$$ \theta_{im} \frac{\partial C_{im}}{\partial T} = \frac{1}{\tau(\bar{r})} (C_m - C_{im}). \quad (2) $$

Here, $\theta_m$ represents the (constant) mobile water content and $\theta_{im}$ the (constant) immobile water content. For simplicity, we assume that local scale dispersion is constant and isotropic as quantified by the dispersion coefficient $D_0$. Note that we assume here for simplicity that the trapping and release rates are constant and equal to the inverse characteristic mass transfer time $1/\tau(\bar{r})$ as in Reichle et al. (1998).

Experimental field and laboratory observations report that solute transport may not obey classical Fickian dispersion (Brusseau and Srivastava, 1999; Berkowitz et al., 2006), and may be characterized by late time tailing of solute breakthrough curves and a non-linear evolution of the spatial second centered moments of the solute distribution, which can be quantified by stochastic perturbation theory for moderately heterogeneous media (Dagan, 1989). Such observations cannot be captured by Fickian advection-dispersion models that are based on constant (asymptotic) effective transport parameters. The heterogeneity-induced large scale transport dynamics are in general non-Fickian, and converge to Fickian behaviors only asymptotically (Dentz et al., 2004; Beaudoin and de Dreuzy, 2013; Comolli et al., 2019). Such behaviors can often be traced back to a broad distribution of mass transfer time scales as a result of spatial heterogeneity in the physical and chemical medium properties (Berkowitz et al., 2006; Frippiat and Holeyman, 2008; Neuman and Tartakovsky, 2008; Dentz et al., 2011). This phenomenological picture is reflected in the multirate mass transfer (MRMT) (Haggerty and Gorelick, 1995; Carrera et al.,
1998), continuous time random walk (CTRW) (Berkowitz et al., 2006) and time-domain random walk (TDRW) (Painter and Cvetkovic, 2005; Noetinger et al., 2016) approaches. These modeling approaches are typically characterized by a probability distribution function of mass transfer times, which are encoded in the memory function (MRMT), and the transition time distribution (CTRW and TDRW). The MRMT approach models solute transport through advection and diffusion or dispersion in a connected mobile region and linear mass transfer between the mobile and immobile regions. The mass transfer properties and distribution of residence times in the immobile regions are encoded in the memory function, which relates the mobile and immobile solute concentrations. The memory function of MRMT can in principle be related to the geometry and distribution of immobile regions and the underlying mass transfer processes (Haggerty and Gorelick, 1995; Carrera et al., 1998; Zinn et al., 2004; Gouze et al., 2008; Zhang et al., 2014; Hidalgo et al., 2020), or is fitted from observed breakthrough curves based on a parametric form for the memory function (Willmann et al., 2008). Similarly, the transition time distribution in the CTRW approach is often fitted based on a suitable parametric equation (Dentz et al., 2004; Berkowitz et al., 2006). The CTRW and TDRW frameworks allow to distinguish between local scale mass transfer mechanisms and to relate the transition probabilities to the heterogeneity distribution and local scale processes such as linear kinetic mass transfer (Margolin et al., 2003), spatially variable retardation properties (Dentz and Castro, 2009), heterogeneous advection (Fiori et al., 2007; Cvetkovic et al., 2014; Comolli et al., 2016, 2019), matrix diffusion (Hyman and Dentz, 2021), and sorption and first-order decay (Painter et al., 2008).

In this article, we study transport in heterogeneous porous media that are characterized by spatially variable permeability and spatially variable kinetic mass transfer properties. In the first part, we solve numerically the direct flow and transport model described by Eqs. (1)–(2). We focus on the pre-asymptotic regime where anomalous behaviors are observed, and their dependence on the system parameters.

In the second part, we use the CTRW framework to derive a new upscaled stochastic model, that is able to predict non-Fickian pre-asymptotic transport as well as the asymptotic behavior. Our CTRW approach uses a compound Poisson process to model linear kinetic exchange under spatially heterogeneous mass transfer properties. It is parameterized by the point distributions of hydraulic permeability, the correlation length of log-permeability, local scale dispersion and advective tortuosity, all of which are flow and medium properties, which means that they are transport-independent. The model is predictive without any fitting parameter.

The article is structured as follows. In Section 2, we first detail the Mobile–Immobile model (MIM) for heterogeneous permeability and exchange time, the generation of the permeability field, and the flow equation. Then, in Section 3 we summarize the numerical simulation of the direct problem using the Lattice-Boltzmann method and present a parametric study of the transient regimes. Section 4 presents a new upscaled transport model (CTRW) and its parameterization. The CTRW model is then validated against the results of the direct simulations.

## 2. The mobile–immobile model for heterogeneous structures

In this part, we present in a first step the derivation of a generic scaling law to relate the exchange time $\tau$ and the permeability $K$. Then, we introduce the non-dimensional equations of the Mobile–Immobile model and the dimensionless parameters for heterogeneous porous media. Afterwards, we detail the generation of the porous medium and present the flow equation.

### 2.1. Generic scaling law to relate the exchange time $\tau$ and the permeability $K$

The characteristic exchange time $\tau$ depends on the underlying physical or chemical process and the porous topology. Using a phenomenological argument, one might, for example, expect that the characteristic

diffusion time in the dead-end pore scales as Reichle et al. (1998): $\tau \sim a^2/D_m$, with $a$ being the typical pore size and $D_m$ the molecular diffusion coefficient. Assuming that $K \propto a^2$, it follows that $\tau \propto K$. The time spent in the dead end thus increases with the permeability because the pores are larger. However, this argument holds only if the permeability is modified by a homothetic transformation. For more complex structures, the residence time may be a non-trivial function of the pore scale geometry. For example, percolation disorder is known to lead to anomalous scaling exponents (Koplik et al., 1988; De Gennes, 1983). Also, based on field measurements, Burr et al. (1994) proposed a power-law between the permeability and the exchange time. In order to be more general, we thus correlate the exchange time to the permeability using a power law:

$$ \tau = AK^{\gamma}, \quad (3) $$

where $A$ is a scaling parameter and the exponent $\gamma \in [-1, 1]$. We consider also negative $\gamma$, which corresponds to an increase of residence time when the permeability decreases. Indeed, if the lower permeability regions have dead-ends with more complex shape, one could expect that the residence times are indeed larger.

### 2.2. Non-dimensional equations

The number of variables of the Mobile–Immobile model can be reduced by changing the time scale as $T/\theta_m \to t$. Eqs. (1) and (2) can then be written as

$$ \left\{ \begin{aligned} \frac{\partial C_m}{\partial t} + \phi_D \frac{\partial C_{im}}{\partial t} + \vec{u} \cdot \vec{\nabla} C_m &= D_0 \nabla^2 C_m \\ \frac{\partial C_{im}}{\partial t} &= \frac{1}{\tau \phi_D} (C_m - C_{im}) \end{aligned} \right. \quad (4) $$

with $\phi_D = \frac{\theta_m}{\theta_e}$ being the ratio of the immobile volume over the mobile one. Dispersion is given by its coefficient $D_0 = du + D_m$, where $d$ is the dispersivity and $D_m$ the molecular diffusion. Note that in heterogeneous porous media the velocity $u$ is also heterogeneous. In order to investigate solely the influence of heterogeneous exchange time $\tau$ on the transport, we simplified the problem and kept $D_0 = \text{const.}$ As mentioned above $\tau$ is then assumed to follow Eq. (3).

Using the imposed mean flow rate $\bar{u}$, Eqs. (4) can be non-dimensionalized by introducing a characteristic length scale $\lambda$ and time scale $\lambda/\bar{u}$:

$$ \left\{ \begin{aligned} \frac{\partial C_m}{\partial t' } + \phi_D \frac{\partial C_{im}}{\partial t' } + \vec{u}' \cdot \vec{\nabla}' C_m &= \frac{1}{P_e} \nabla'^2 C_m \\ \frac{\partial C_{im}}{\partial t' } &= \frac{\lambda}{u \phi_D \bar{\tau}} (C_m - C_{im}) \end{aligned} \right. \quad (5) $$

Two dimensionless numbers can be deduced from Eqs. (5):

* $Pe = \frac{\lambda \bar{u}}{D_0}$, the Péclet number, which characterizes the ratio of the diffusive time scale over the convective one.

* $Da = \frac{\lambda}{\bar{u} \bar{\tau} \phi_D}$, the Damköhler number, corresponding to the ratio of the convection and exchange time scales. Here, we introduce the average exchange time $\bar{\tau} = AK_0^{\gamma}$, where $K_0 = \exp(f_0)$ is the harmonic mean of the permeability field (see Eq. (6) in the next section).

### 2.3. Porous media generation

As it is often assumed to describe the permeability field (Gelhar and Axness, 1983; Burr et al., 1994; Rubin, 2003), $K$ is generated using a correlated log-normal probability distribution function (PDF) given by:

$$ \text{pdf}(f = \ln K) = \frac{1}{\sqrt{2\pi}\sigma_f} e^{-\frac{(f-f_0)^2}{2\sigma_f^2}}, \quad (6) $$
where $f_0$ and $\sigma_f$ are the mean and standard deviation of $\ln K$. Moreover, the field $f$ is spatially correlated according to a Gaussian correlation function:

$$ \langle f(\vec{r}_0)f(\vec{r}_0 + \vec{r}) \rangle_{\vec{r}_0} = \sigma_f^2 e^{-\frac{(x^2+y^2)}{\lambda^2}}, \quad (7) $$

where $\langle \cdot \rangle_{\vec{r}_0}$ represents the spatial average over the whole domain. Here, we use the fixed parameters $f_0 = 0$ and $\lambda = 4$ LBM lattice units. The value $\sigma_f$, which quantifies the heterogeneity of the permeability field, has been varied in the range $\sigma_f^2 \in [0, 1.5]$, which represents moderately heterogeneous permeability fields, but includes most values reported in Rubin (2003).

Since exchange time and permeability are related by Eq. (3), the distribution of the exchange times also follows a log-normal distribution:

$$ \text{pdf}(g = \ln \tau) = \frac{1}{\sqrt{2\pi}\sigma_g} e^{-\frac{(g-\sigma_g)^2}{2\sigma_g^2}}, \quad (8) $$

with $g_0 = \gamma f_0 + \ln(A)$ and $\sigma_g = |\gamma|\sigma_f$. The covariance function of $g$ is thus given by

$$ \langle g(\vec{r}_0)g(\vec{r}_0 + \vec{r}) \rangle_{\vec{r}_0} = \gamma^2 \sigma_f^2 e^{-\frac{x^2+y^2}{\lambda^2}}. \quad (9) $$

## 2.4. Flow equation

Flow in heterogeneous porous media is obtained by solving the Darcy–Brinkmann equation (Brinkman, 1947):

$$ \frac{\rho v}{K(\vec{r})} \vec{u} - \vec{\nabla} P + \rho v \nabla^2 \vec{u} = \vec{0}, \quad (10) $$

where $v$ is the kinematic viscosity, $\rho$ the density of the fluid and $P$ the pressure field. The Darcy–Brinkmann equation is generally used to solve flow in very heterogeneous porous media because it allows momentum diffusion between high and low velocity regions. It also provides stability of the numerical scheme by preventing strong spatial variation of the velocity (Talon et al., 2003; Ginzburg et al., 2015). To solve this equation, we used a Two-Relaxation-Time lattice Boltzmann method (TRT-LBM) described in Appendix A. We also provided informations on the homogeneous case and the validation of the LBM-code against the analytical solution.

# 3. Characterization of transient transport regimes by means of detailed Lattice Boltzmann simulations

In this section, we study transport in heterogeneous porous media under mobile-immobile mass transfer. The geological model is characterized by the permeability field parameters $\phi_D$, $f_0$, $\sigma_f$, and $\lambda$, and by the transport parameters $A$, $\gamma$, and $D_0$. However, we focus here mainly on the influence of the variance of the heterogeneity $\sigma_f$ and the Damköhler number $Da$ using Lattice Boltzmann simulations.

## 3.1. Spatial moments and breakthrough curves

Transport is analyzed by means of the evolution of the longitudinal second centered spatial moment of the normalized mobile concentration

$$ c_m = \frac{C_m}{\int d\vec{r}C_m}. \quad (11) $$

The second centered moment of $c_m$ is defined by

$$ \sigma^2(t) = \int d\vec{r} x^2 c_m - \left( \int d\vec{r} x c_m \right)^2 \quad (12) $$

where $x$ is the component of the position vector $\vec{r}$ in streamwise direction. $\sigma(t)$ measures the spatial extension of the mobile solute plume

in longitudinal direction. The rate of growth of $\sigma^2(t)$ is measured by the macroscopic longitudinal dispersion coefficient, which is defined by

$$ D_L^m(t) = \frac{1}{2} \frac{d\sigma^2(t)}{dt}. \quad (13) $$

Furthermore, we consider the mean concentration, that is, the breakthrough curve, slightly before the outlet boundary of the flow domain, at $x = L - 4\lambda$. It is defined by

$$ BC(t) = \frac{1}{H} \int_y c_m(\vec{r}) dy, \quad (14) $$

where $y$ is the component of the position vector perpendicular to the streamwise direction and $H$ is the height of the flow domain.

## 3.2. Transport in heterogeneous porous media

In this section we detail the results of the transport simulations in heterogeneous porous media. Fig. 1 shows typical mobile and immobile concentration fields for two different values of the Damköhler number, $Da = 1/8$ and $Da = 1/40$ with $\sigma_f = 1$. The two mobile concentration fields are very similar, whereas important differences in the immobile concentration fields can be observed. In the case of low $Da$, there is a larger number of depleted spots (dark blue) at the tracer front while other spots of higher concentration can still be observed far behind the front. Both effects result from the disorder in the exchange time $\tau$, where some spots are characterized by large values of $\tau$. These spots require more time for tracer trapping but also for its release.

Fig. 2 displays the normalized second moment of the average concentration profile (Left) and the breakthrough curve (Right) for the parameters $Da = 1/40$, $\sigma_f = 1$, $Pe = 6.4$ and $\gamma = 1$. For small times, the normalized variance follows a power-law of $t$ with an exponent $\alpha_{early}$ larger than 1, characteristic of anomalous transport. For larger times, the evolution of the variance changes into a different anomalous behavior characterized by another power-law exponent $\alpha$. The breakthrough curve exhibits non-exponential tailing which can also be characterized by a power-law of the following type:

$$ BC(t) \propto t^{-\beta}. \quad (15) $$

Due to the fact that the second moment and the breakthrough curve do not always display a power-law behavior, a criterion needs to be defined to establish its existence. In practice, the fit is done by doing a linear regression of the logarithmic derivative $\frac{d\log f}{d\log t}$ over a decade of the abscissa. The constant of this regression gives then the exponent. Based on the slope, a criterion of relevance of the power-law can then be defined. We fixed the following criterion: if the slope is greater than 0.2, the curve is not considered as a power-law. In other words, if the exponent varies more than 0.2 over a decade, it is not considered as constant.

In the following, we investigate the effect of the Damköhler number $Da$, the disorder $\sigma_f$ of the permeability field and the exponent $\gamma$, on the anomalous transport characterized by the power-law exponents $\alpha_{early}$, $\alpha$ and $\beta$.

Effect of the Damköhler number $Da$. In this paragraph, we investigate transport for different values of $Da$ in a system characterized by the following parameters: $\sigma_f = 1$, $Pe = 6.4$ and $\gamma = 1$. Fig. 3 represents the evolution of the normalized second moment (Left), the normalized dispersion coefficient (Middle) and the breakthrough curve (Right) for different $Da$ values. As can be seen, both, the power-law of the breakthrough curve ($\beta$) and the second regime of the spatial moment ($\alpha$) depend on $Da$. This dependency is particularly obvious when plotting $D_L^m(t)/D_0$. The evolution of the exponent $\alpha$ as a function of $Da$ is given in Fig. 4 (Left). For low and high values of $Da$, $\alpha$ tends to one, suggesting a mere dispersive behavior if the exchange time is either very low or very high. For large values of $Da$, the system is instantaneously at equilibrium (e.g. $C_m = C_{lm}$). In contrast, for small values of $Da$, the exchange between the mobile and immobile regions becomes extremely
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Examples of mobile (Top) and immobile (Bottom) concentration profiles in a heterogeneous porous medium. System parameters are: $\sigma_f = 1$, $Pe = 6.4$ and $\gamma = 1$. Left: $Da = 1/8$, Right: $Da = 1/40$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Temporal evolution of the normalized second moment (Left) and the breakthrough curve (Right) for $Da = 1/40$ and $\sigma_f = 1$, $Pe = 6.4$ and $\gamma = 1$. The fitted exponents are: $\alpha_{early} = 1.3$, $\alpha = 1.7$ and $\beta = 5.2$.</figcaption></figure>

slow and its influence on the overall transport becomes negligible. This case corresponds to the passive transport case studied in the literature, where the dispersion regime is initially anomalous ($\alpha_{early} > 1$) due to the permeability heterogeneity, and becomes normal ($\alpha = 1$) after a certain time period.

In the presence of exchanges between mobile and immobile regions, the transition time between the two regimes (power-law independent of $Da$ and power-law controlled by $Da$ given in Fig. 3) is roughly independent of $Da$. It corresponds therefore to the end of the non-reactive regime which dominates at early time. Above this transition, the effect of the exchange between mobile and immobile zones on transport becomes predominant.

Considering the breakthrough curves (Fig. 3, Right) it can be seen that, for $Da = 0$, the curve takes a Gaussian shape, while breakthrough curves with finite values of $Da$ exhibit significant tailing. Tailing is more pronounced for intermediate $Da$ compared to large and small $Da$ because in both cases the influence of the exchange on the transport is not significant. For certain values of $Da$, tails follow a power-law given by $BC(t) \propto t^{-\beta}$. Fitted exponents of the latter power-law are reported in Fig. 4 (Right). It can be seen that the exponent is increasing with $Da$ corresponding to less significant tailing of the concentration field at high $Da$.

**Effect of disorder of the permeability field $\sigma_f$.** We now investigate the influence of $\sigma_f$ on the transport in a system characterized by $Da = 1/8$, $Pe = 6.4$ and $\gamma = 1$. The evolution of the normalized variance of

the concentration profile, the normalized dispersion coefficient and the breakthrough curves are given in Fig. 5. In comparison with Fig. 3, we can state that the early time regime depends on $\sigma_f$. For low heterogeneity ($\sigma_f \le 0.25$), the early time behavior is diffusive as can be seen from the normalized dispersion coefficient in Fig. 5 (Middle). The transport becomes anomalous for higher values of $\sigma_f$. When fitting the early time data, we can state that the power-law exponent increases with heterogeneity (see Fig. 6, Left).

As already reported in the latter section, after the early time regime there is a transition towards the second regime at a time which is relatively independent of $\sigma_f$. Above this transition, the normalized dispersion coefficient follows different power-laws, with an exponent that also depends on $\sigma_f$ (see Fig. 6, Left). For this particular small value of $Da = 1/8$, the early time regime exponent $\alpha_{early}$ is higher than the exponent $\alpha$ of the second regime. The opposite is however possible for some lower values of $Da$ as shown in Fig. 4.

In conclusion, $\alpha_{early}$ increases with $\sigma_f$ but not with $Da$, this stands for the fact that the early time regime is dominated by the heterogeneity of the permeability but not by the exchange time.

Fig. 5 (Right) represents the corresponding breakthrough curves as a function of $\sigma_f$. For very small $\sigma_f$, the tail is similar to a Gaussian distribution. By increasing $\sigma_f$, a retardation tail appears which becomes a power-law for high heterogeneity values ($\sigma_f \ge 1$). Fig. 6 represents the exponent $\beta$ of the power-law fit (with Eq. (15)) of the tail of the breakthrough curve as function of $\sigma_f$. The exponent decreases with increasing heterogeneity and tailing becomes very significant.
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Normalized second moment (Left), normalized dispersion coefficient (Right) and breakthrough curve (Bottom) as function of the normalized time for different Da values. The legend is common to all three figures. Other parameters are set to $\sigma_f = 1$, $\gamma = 1$ and $Pe = 6.4$. At early times, the second moment follows an anomalous dispersion law with $D_L^m(t)/D_0 \propto t^{-1+\alpha_{early}}$. After a certain time, the system reaches either a normal dispersion regime at low and high Da values or another anomalous dispersion regime characterized by $D_L^t(t)/D_0 \propto t^{-1+\alpha}$ for intermediate values of Da. Except for the case Da = 0 (no trapping), the breakthrough curves present a retardation tail. For intermediate values of Da, the tail follows a power-law given by $BC(t) \propto t^{-\beta}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Exponents $\alpha_{early}$ and $\alpha$ obtained by the power-law fit as function of Da for the same parameters as in Fig. 3 (Left). Whereas the exponent $\alpha_{early}$ is independent of Da, the later time exponent $\alpha$ presents a non-monotonic behavior for intermediate Da values. Exponent $\beta$ obtained, when applicable (see text), by the power-law fit of the tail of the breakthrough curve (Right).</figcaption></figure>

Effect of $\gamma$. In this section we investigate the impact of $\gamma$, characterizing the relationship between the permeability and the exchange time: $\tau \propto K^\gamma$. Fig. 7 (Left) represents the evolution of the normalized dispersion coefficient for $\gamma = -1, 0$ and $1$ with $\sigma_f = 1$, $Da = 1/8$ and $Pe = 6.4$. For better comparison we set the logarithmic average of the permeability field $\exp(f_0)$ equal to one so that the average exchange time is independent of $\gamma$. From Fig. 7 (Left) it can be seen, that in the late time regime, dispersion is anomalous for $|\gamma| \neq 0$ but is normal for $\gamma = 0$. The latter is characterized by $D_L^m(t)/D_0 \simeq \text{const}$ and corresponds to case where the permeability is heterogeneous but the exchange time is homogeneous. For the two cases $\gamma = \pm 1$, the evolution is characteristic of anomalous dispersion and the breakthrough curve (Fig. 7, Middle)

presents a power-law tailing effect. We note that both, $D_L^m(t)/D_0$ and the breakthrough curve, are very similar for $\gamma = \pm 1$.

Fig. 7 (Right) displays the evolution of the anomalous exponent $\alpha$ as function of $\gamma$ obtained from the second moment of this data set. We observe that data points are relatively symmetric with respect to $\gamma = 0$ and increase with $|\gamma|$. Due to the fact, that the exchange time distribution is log-normal, exchange time distributions corresponding to $\gamma$ and $-\gamma$ have the same standard deviation and are thus statistically equivalent. For $\gamma > 0$ high permeabilities correspond to high exchange times, whereas in the case of $\gamma < 0$, high permeabilities correspond to low exchange times. The fact that the curves are very similar, indicates therefore that the anomalous second regime originates from the
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Normalized variance (Left) and normalized dispersion coefficient (Right) as function of the normalized time for different values of $\sigma_f$. Other parameters are set to $Da = 1/8$, $\gamma = 1$ and $Pe = 6.4$. Breakthrough curves for the same parameters (Bottom).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Exponents $\alpha_{early}$ and $\alpha$ as function of the $\sigma_f$ for the same parameters as in Fig. 5 (Left). Both, $\alpha_{early}$ and $\alpha$ vary with $\sigma_f$. Exponent $\beta$ of the power-law fit, when applicable (see text), of the tail of the breakthrough curve (Right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Normalized dispersion coefficient for different values of $\gamma$ (Left), corresponding breakthrough curves (Middle) and exponents of the power-law fit of the second moment (Right). Simulations were performed with the following system parameters: $\gamma = -1, 0, 1$ and $\sigma_f = 1$, $Da = 1/8$ and $Pe = 6.4$.</figcaption></figure>
presence of a broad range of characteristic exchange times. However, this does not depend on whether higher exchange times correspond to regions of higher or lower permeability. This behavior will be discussed later by means of the CTRW model.

## 4. Stochastic transport model

To better understand the physical mechanism at work in this problem, we derive a stochastic transport model based on the continuous time random walk (CTRW) approach. The novelty of our approach lies in the combined modeling of linear kinetic exchange processes under spatially heterogeneous mass transfer properties and heterogeneous hydraulic properties. The derived model is fully predictive in the sense that the parameters are fully constrained in terms of the system properties, that is, there are no fitting parameters.

We focus on longitudinal solute transport, that is, in the direction of the mean pressure gradient. The transport of solute particles in a heterogeneous porous medium is determined by the spatial variability of the medium and flow structure. Consequently, particle velocities and thus particle transport evolve on a characteristic length scale. This notion is naturally accounted for by the CTRW approach, which models particle motion as a space time random walk, in which the average spatial step is constant and the time increment varies according to the local flow velocity. In the following, we define the trapping rate $\omega \equiv 1/\tau$ as the inverse trapping time. Thus, the evolution of the particle position $x_i$ and time $t_i$ after i random walk steps can be written as,

$$x_{i+1} = x_i + \frac{\Delta\ell}{\chi} + \sqrt{2D_0 \frac{\Delta\ell}{v_i}} \xi_i, \quad (16)$$

$$t_{i+1} = t_i + \tau_{v_i} + \theta_s(\tau_{v_i}, \omega_i). \quad (17)$$

These equations read as follows: Eq. (16) describes the particle position after i random walk steps. The second term on the right side of Eq. (16) denotes an advective transition over the constant distance $\Delta\ell/\chi$, where $\Delta\ell$ is an increment along a streamline, which is projected onto the streamwise direction by the advective tortuosity $\chi$ (Comolli et al., 2019). The third term denotes the dispersive transition, which is given by a unit Gaussian random variable $\xi$ and the amplitude $\sqrt{2D_0\Delta\ell/v_i}$, which is the typical dispersion length during the advection time $\Delta\ell/v_i$ with $v_i$ the particle speed at step i. Eq. (17) describes the particle time after i steps. The second term denotes the purely advective transition, which is given by

$$\tau_{v_i} = \frac{\Delta\ell}{v_i}. \quad (18)$$

The single point distribution of $\tau_v$ is denoted by $\psi_v(t)$. It is defined below in terms of the distribution of flow speeds. The third term $\theta_s$ on the right side of Eq. (17) denotes the total trapping time, where the particle can be immobilized many times, with rate $\omega_i$, during an advective step of duration $\tau_{v_i}$. The discretization length $\Delta\ell$ should be smaller than the characteristic length scales, that is, the correlation scales of the underlying disorder.

We note that $\theta_s$ depends on $\tau_v$, as the probability to trapped many times increases with the advective time. In the following, we first summarize the modeling of the velocity transitions for a physically heterogeneous medium. Then, we derive the stochastic representation of spatially variable linear kinetic sorption.

### 4.1. Velocity process

The velocity process, defined as the series of particle speeds $\{v_i\}$, is modeled as a stationary Markov process (Dentz et al., 2016; Morales et al., 2017; Comolli et al., 2019) whose steady-state distribution $p_s(v)$ is equal to the flux-weighted Eulerian flow speed, $v_e$, distribution (Dentz et al., 2016)

$$p_s(v) = \frac{v p_e(v)}{\langle v_e \rangle}. \quad (19)$$

<figure><img src="image_11.png" /><figcaption>Fig. 8. Probability distribution function for the normalized Eulerian velocity resulting from the simulation with $\sigma_f = 1$.</figcaption></figure>

The distribution $\psi_v(t)$ of advective transition times is given in terms of $p_s(v)$ and $p_e(v)$ as

$$\psi_v(t) = \frac{\Delta\ell}{t^2} p_s(\Delta\ell/t) = \frac{\Delta\ell^2}{t^3} \langle v_e \rangle p_e(\Delta\ell/t). \quad (20)$$

Relation (19) can be understood as follows. The Eulerian speed distribution $p_e(v)$ is equal to the steady state distribution of isochronically sampled flow speed, which is a direct consequence of the incompressibility of the flow. The particle speeds in the CTRW, on the other hand, are sampled equidistantly along trajectories, that is, low velocities are sampled less often than for isochronic sampling, which is the reason for the flux-weighting in relation (19). This is discussed in detail by Comolli et al. (2019). The speed distribution $p_e(v)$ is approximately log-normal,

$$p_e(v) = \frac{\exp\left[ -\frac{(\ln v - \mu_e)^2}{2\sigma_e^2} \right]}{v\sqrt{2\pi\sigma_e^2}}, \quad (21)$$

and can be seen in Fig. 8. From relation (19), we find that $p_s(v)$ is also log-normal with the log-mean given by $\mu_s = \mu_e + \sigma_e^2$ and the same log-variance $\sigma_e^2$ as $p_e(v)$. This implies that the distribution $\psi_v(t)$ is also a log-normal.

Particle transitions between successive particle speeds are characterized by the conditional probability $p(v, \Delta\ell|v')$, which depends on the space increment $\Delta\ell$. That is, the particle speed $v_{n+1}$ is drawn from $p(v, \Delta\ell|v')$ given the value $v'$ of the particle speed $v_n$ after n steps. Following Morales et al. (2017), we model $p(v, \Delta\ell|v')$ by the conditional log-normal distribution

$$p(v, \Delta\ell|v') = \frac{\exp\left[ -\frac{(\ln v - \mu_s - \ln v'[1 - \exp(-\Delta\ell/\ell_e)])^2}{2\sigma_e^2[1 - \exp(-2\Delta\ell/\ell_e)]} \right]}{v\sqrt{2\pi\sigma_e^2[1 - \exp(-2\Delta\ell/\ell_e)]}}, \quad (22)$$

where $\ell_e$ denotes the correlation length of flow speeds.

The distribution $p_0(v)$ of initial particle speeds $v_0$ depends on the boundary or initial condition. Here, particles are injected uniformly over a line perpendicular to the mean flow velocity. This implies that $p_0(v) = p_e(v)$ (Comolli et al., 2019).

### 4.2. Linear kinetic exchange process under spatially heterogeneous mass transfer properties

In a particle-based framework, linear kinetic mass exchange processes can be modeled as a compound Poisson process (Feller, 1968). Several authors (Margolin et al., 2003; Benson and Meerschaert, 2009; Comolli et al., 2016; Hyman and Dentz, 2021) have used this approach to model multirate mass transfer and matrix-diffusion under
spatially homogeneous linear mass transfer properties. Based on the compound Poisson process, we derive here a novel CTRW approach for heterogeneous mass transfer systems.

The trapping of particles during a random walk step occurs at rate $\omega = 1/\tau$ such that the average number of trapping events during the advection time $\tau_v$ is given by $\omega\tau_v$. The number $n$ of trapping events is then distributed according to the Poisson distribution

$$p_n(\tau_v, \omega) = \frac{(\omega\tau_v)^n \exp(-\omega\tau_v)}{n!} \quad (23)$$

The total trapping time $\theta_s(\tau_v, \omega)$ during a random walk step is given by

$$\theta_s(\tau_v, \omega) = \sum_{j=1}^{n} \theta_j, \quad (24)$$

where $\theta_j$ is the trapping time for the $j$th sorption event. The $\theta_j$ are independent, exponentially distributed according to

$$\psi(\theta) = \omega\phi_D \exp(-\omega\phi_D\theta) \quad (25)$$

with $\phi_D$ being the ratio of the immobile volume over the mobile one. The total trapping time $\theta_s$ describes a compound Poisson distribution (Feller, 1968). Its distribution can be expressed in Laplace space as (see Appendix D)

$$\psi_c^*(s|\tau_v, \omega) = \exp\left[-\omega\tau_v \left(1 - \frac{\omega\phi_D}{\omega\phi_D + s}\right)\right], \quad (26)$$

where $s$ is the Laplace parameter. The Laplace transform is defined in Abramowitz and Stegun (1972). The inverse Laplace transform of expression (26) can be evaluated explicitly and is given by Bateman (1954, p. 244)

$$\psi_c(t|\tau_v, \omega) = \exp(-\omega\tau_v)\delta(t) + \exp[-\omega(\tau_v + t\phi_D)]\omega\sqrt{\frac{\tau_c\phi_D}{t}}I_1\left(2\sqrt{\omega\phi_D t\omega\tau_v}\right), \quad (27)$$

where $I_1(t)$ is the modified Bessel function of the first kind. The first term denotes the probability that no exchange event occurs, in which case the trapping time is zero. Fig. 9 shows the behavior of $\psi_c(t, \tau_v, \omega)$ for $\omega = 1$, $\phi_D = 1$ and different $\tau_v$.

From these expressions, one can deduce several statistical properties. From the moment generating function Eq. (26), one can infer the mean and variance of the total trapping time:

$$\bar{\theta}_s(\tau_v, \omega) = \frac{\tau_v}{\phi_D} \quad (28)$$

$$\overline{(\theta_s - \bar{\theta}_s)^2}(\tau_v, \omega) = \frac{2\tau_v}{\omega\phi_D^2} = \frac{2\tau_v\tau}{\phi_D^2} \quad (29)$$

The first relation shows that the average total trapping time $\theta_s$, depends only on the advective time $\tau_v$ and not on the exchange time $\tau$. This counter-intuitive result can be understood from the fact that $\omega = 1/\tau$ controls both the number of trapping events during the advection time and the duration of each. The two effects compensate each other since a higher $\tau$ means longer trapping times but also fewer trapping events. The second relationship shows that the width of the distribution is impacted by both the advection time and the exchange time. This can also be understood from the fact that the total exchange time during $\tau_v$ is the sum of $n$ trapping times. The variance is thus proportional to $n\tau^2$ with $n = \tau_v/\tau$.

Another important quantity is the total trapping time distribution tail. In Eq. (27) the limit of $t \gg (\omega\phi_D)^{-1}$, $\psi_c(t|\tau_v, \omega)$ becomes also independent of $\omega$ and can be approximated by $\psi_c(t|\tau_v, \omega) = \delta(t - \tau_v/\phi_D)$, which means that $\theta_s(\tau_v, \omega) \approx \tau_v/\phi_D$ for large values. The longest exchange times are therefore controlled by the advection time only. Furthermore, we can infer that the tail of the distribution is a log-normal.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Distribution $\psi_c(t|\tau_v, \omega)$ of exchange times for $\omega = 1$, $\phi_D = 1$ and $\tau_v = 1, 4, 10$.</figcaption></figure>

### 4.3. Model parameterization

Here, we discuss the parameterization of the upscaled model in terms of the statistical medium description, which is essential as our CTRW approach is not based on a fitting procedure.

#### 4.3.1. Linear kinetic mass exchange

First, we consider the exchange process. The model assumes that the exchange rates $\omega_j$ are independent identically distributed random variables. Their distribution $p_\omega(\omega)$ is given in terms of the distribution of permeability $p_K(k)$ through the relation $\omega = 1/\tau = 1/(AK^\gamma)$,

$$p_\omega(\omega) = \frac{\gamma}{A^{1/\gamma}\omega^{1/(\gamma+1)}} p_K[1/(A\omega)^{1/\gamma}]. \quad (30)$$

That is, it is fully constrained here by the point distribution of permeability, and the parameter of the direct model. There is no fitting parameter.

#### 4.3.2. Velocity process

Second, we consider the process of particle speeds, which is defined by the conditional log-normal distribution Eq. (22) and the log-normal distribution Eq. (21) of Eulerian flow speeds. Thus, the speed process is parameterized by the mean $\mu_e$ and variance $\sigma_e^2$ of the Eulerian flow speed, as well as the correlation length $\ell_c$. These parameters can be constrained by perturbation theory in the fluctuations of $\ln K$ for low and moderate heterogeneity characterized by $\sigma_f^2$ lower than or of the order of one. For strong heterogeneity, empirical relations for the Eulerian speed distribution and correlation length need to be used (Gotovac et al., 2009; Hakoun et al., 2019; Comolli et al., 2019).

To determine the model parameters by perturbation theory, we notice, that $\mu_e$ and $\sigma_e^2$ for the log-normal distribution are related to the mean $\bar{v}$ and variance $\sigma_v^2$ of the Eulerian flow speed by

$$\mu_e = \ln \left( \frac{\bar{v}^2}{\sqrt{\bar{v}^2 + \sigma_v^2}} \right) \quad (31)$$

$$\sigma_e^2 = \ln \left( 1 + \frac{\sigma_v^2}{\bar{v}^2} \right) \quad (32)$$

The mean flow speed $\bar{v}$ and its variance $\sigma_v^2$ can be expressed in terms of the mean $f_0$ and variance $\sigma_f^2$ of $f = \ln K$ by Cvetkovic et al. (1996)

$$\bar{v} = \exp(f_0)\sqrt{\bar{v}} \quad (33)$$

$$\sigma_v^2 = \frac{3}{8}\sigma_f^2\bar{v}^2 \quad (34)$$
<figure><img src="image_13.png" /><figcaption>Fig. 10. Second moment as a function of time for different values of Da with $\sigma_f = 1$. Solid line: LB simulation, red symbols: CTRW simulations. The red dashed line represents a linear trend $\sigma^2 \propto t$.</figcaption></figure>

where $\nabla \bar{p}$ is the mean hydraulic gradient. Thus, expressions (31)–(34) constrain the speed distribution in terms of the K distribution. Similarly, we can relate the speed correlation $\ell_c$ to the correlation length of $\lambda$ of $f$ by $\ell_c = \sqrt{\pi 4\lambda/3}$, where $\lambda$ is the correlation length of permeability (Cvetkovic et al., 1996; Comolli et al., 2019), see also Appendix A. The advective tortuosity $\chi$ is given by the ratio between the mean Eulerian speed $\bar{v}_e$ and the mean streamwise velocity $\bar{u}$, $\chi = \bar{v}_e/\bar{u}$. Thus, all parameters of the upscaled model are constrained by flow and medium properties. There are no fitting parameters. Note that the upscaled model does not represent the correlation between permeability and exchange rate expressed in the detailed model by Eq. (3). Here, particle velocity on one hand, and exchange rate on the other are independent.

### 4.4. Transport predictions with the upscaled CTRW model

In this section, we compare the predictions of the upscaled CTRW to the detailed numerical simulations for three transport scenarios characterized by $\sigma_f = 1$ and $Da = 1.25 \times 10^{-1}$, $6.25 \times 10^{-2}$ and $2.5 \times 10^{-2}$. The coupling parameter $\gamma$ between permeability and exchange time is set to $\gamma = 1$. The correlation length of $f = \ln K$ is $\lambda = 4$ in arbitrary units, the mean of $\ln K$ is $f_0 = 0$, and the mean pressure gradient is unity. The upscaled CTRW model is parameterized based on these characteristics as outlined in Section 4.3. The numerical implementation of the upscaled model described here is efficient and allows to achieve simulation times otherwise prohibitive for the direct numerical LB simulations.

**Spatial moments and breakthrough curve.** The position of a mobile particle in the CTRW framework is given by

$$x(t) = x_{n_t} + v_{n_t}(t - t_{n_t}), \qquad (35)$$

where $n_t = \max(n|t_n \le t)$. A particle is considered mobile if $t - t_{n_t} \le \ell/v_{n_t}$. The second centered longitudinal moment is then defined by

$$\sigma^2(t) = \langle x(t)^2 \rangle - \langle x(t) \rangle^2, \qquad (36)$$

where the angular brackets denote the average over all mobile particles at time $t$. The breakthrough curve $BC_L(t)$ is defined in terms of the arrival times of a particle at distance $L$ from the inlet. The arrival $\tau_L$ is given by

$$\tau_L = \min(t|x_{n_t} \ge L). \qquad (37)$$

The breakthrough curve is then given by the distribution of arrival times as

$$BC_L(t) = \langle \delta[t - \tau_L(x)] \rangle. \qquad (38)$$

#### 4.4.1. Pre-asymptotic regime

Fig. 10 shows the spatial variance normalized by the square of the correlation length for different values of Da obtained from the LB simulations (solid line) and from the stochastic CTRW model (symbols). Fig. 11 displays the normalized dispersion coefficient and Fig. 12 the breakthrough curve. It can be seen that, albeit a slight offset at

<figure><img src="image_14.png" /><figcaption>Fig. 11. Normalized dispersion coefficient for different values of Da with $\sigma_f = 1$. Solid line: LB simulation, red symbols: CTRW simulations.</figcaption></figure>

the beginning, all CTRW predictions agree very well with the direct numerical simulations. From these agreements between the LBM and CTRW, several conclusions can be drawn.

We recall that in the CTRW model, the advection rate and exchange time are drawn independently which is a strong assumption in contrast to the direct simulations where permeability and exchange rate are related by Eq. (3). The distributions of K and $\tau$ of the CTRW model are however related through the exponent $\gamma$. The good agreement between the two approaches indicates that the results are independent of the correlation between $\tau$ and K. However, the results are strongly influenced by the distribution shape of the two quantities (especially their width related by $\gamma$). The last two points explain the above mentioned LBM result that the transport was not affected by the sign of $\gamma$ but only by its magnitude. Indeed, distributions with opposite signs of $\gamma$ have the same log-normal shape, but the K and $\tau$ fields are anti-correlated.

As discussed in the introduction, anomalous dispersion is generally related to a broad distribution of transition times used in the CTRW approach. In our model, there are two contributions to this transition time. The first one is the advection time whose distribution width is directly related to the permeability heterogeneity, $\sigma_f$. This contribution is responsible for the earlier anomalous behavior $\alpha_{early}$ (see Fig. 2). The second one is the waiting time due to sorption. This waiting time distribution depends on both the velocity heterogeneity and the exchange time heterogeneity as described by Eq. (27). Physically, this means that in low velocity regions, particles have more time to be trapped and untrapped. Both $\sigma_f$ and $\sigma_g$ then contribute to the second anomalous behavior, which is also very well captured by the model.

Another indication from these model results is more technical. An important issue when doing direct stochastic numerical simulations is the size of the system which must be representative. This problem does not exist for the CTRW as it is a 1D problem. The good agreement of the results of the two approaches thus shows that the size used in the LBM simulations is sufficiently large for this type of heterogeneity.
<figure><img src="image_15.png" /><figcaption>Fig. 12. Breakthrough curves for different values of Da with $\sigma_f = 1$. Solid line: LB simulation, red symbols: CTRW simulations.</figcaption></figure>

### 4.4.2. Asymptotic regime

The stochastic CTRW model is very efficient because it depends only on point-distributions of the statistical medium properties. Once validated, it allows predictions on very large time scales. In particular, it shows that at very large time the variance scales asymptotically linearly with time, i.e., asymptotically dispersion is normal. This is also reflected in the evolution of the normalized dispersion coefficient $D_L^m(t)$ which converges towards a constant asymptotic value (Fig. 11). Consequently, the power-law behavior, characterized by $\alpha$, identified in the previous section is rather a transient than a true power-law regime. In fact, the distribution $\psi_v(t)$ of the advective transition times is log-normal. Similarly, large exchange times can be approximated by $\theta_s(\tau_v, \omega) \approx \tau_v/\phi_D$ as outlined in Section 4.2, that is, the long-time behavior of the exchange time distribution is also log-normal. This implies that transport is asymptotically normal (Berkowitz et al., 2006) with possible power-law cross-over behaviors as shown in Section 3.2.

Finally, it can be noted that the dispersion coefficient given by $D_{eff} \propto \lim_{t \to \infty} D_L^m(t)/D_0$ decreases with Da. Indeed, as for the homogeneous case (see Appendix C), lower values of Da lead to stronger spreading inducing higher values of $D_L^m(t)/D_0$ in the asymptotic regime. From Fig. 11 it can be seen, that convergence times to reach the asymptotic Gaussian regime are very long. Additionally, this time depends moderately on the Damköhler number Da.

In conclusion, our CTRW model allows robust predictions for all cases we considered, and provides an efficient approach for transport predictions under conditions that are prohibitive for direct numerical simulations.

## 5. Conclusion

In this article, we have studied transport in a heterogeneous porous medium with a spatially varying permeability field combined with the exchange between mobile and immobile regions. To model the coupling between the exchange and the local structure, we assumed the following relation between the local permeability $K$ and the local exchange time: $\tau \propto K^\gamma$. We performed LBM simulations in a large two-dimensional system to study the heterogeneity-induced large-scale transport behavior. The transport behavior is upscaled using a new CTRW approach. The impact of spatial variations in flow velocity and exchange times are modeled in terms of a spatial Markov process for the particle speeds and a compound Poisson process for the mass exchange kinetics. The upscaled model is fully determined by the statistics of permeability and exchange parameters. It has no fitting parameter.

We investigated the influence of the Damköhler number Da, the variance $\sigma_f^2$ of $f = \ln K$ and the value of the coupling exponent $\gamma$ on the temporal evolution of the solute distribution and the breakthrough curves. In the short and intermediate time range, we observed transient anomalous regimes which are characterized by power-law tails of the solute breakthrough curves and a power-law evolution of the spatial variance of the solute distribution. Whereas the short time power-law regime, $\alpha_{early}$, is strongly dependent on the heterogeneity of the permeability field, the intermediate regime, $\alpha$, is mostly dependent on both the exchange time and permeability distribution. The latter regime

becomes more important for large distributions of $\tau$. In particular, we observe an increase in the power-law exponent $\alpha$ for a specific range of Da. In the asymptotic limit, the upscaled CTRW model predicts that the transport behavior becomes Fickian. The resulting dispersion coefficient $D_{eff}$ decreases with increasing Da. We find that the time to attain the asymptotic limit depends on Da.

We also observe that the time to reach the asymptotic regimes in heterogeneous porous media (for $\sigma_f = 1$) is two orders of magnitude higher than the respective time in the homogeneous case (see Appendix C). Thus, the difference in convergence time strongly depends on $\sigma_f$. Indeed, for a given average value of Da, it is expected that the convergence time depends on the largest values of the exchange time distribution which increase with $\sigma_f$.

The upscaled CTRW model provides robust predictions of the behavior observed in the direct numerical simulations. It relates the observed transport behaviors to the distribution of permeability, which controls the distributions of the advective times and the exchange times between mobile and immobile regions. An important advantage of using the CTRW model lies in the fact that the asymptotic limit of the transport behavior can be reached in a very short computation time. For instance, the computation times for the CTRW models is on the order of minutes using one core on a 2 GHz Quad-Core Intel Core i5 processor. In comparison, using a similar architecture, the LBM code would require about 2–3 days to compute $D_L(t)$ and about 2 weeks to compute the BC curve.

The good agreement of the CTRW with the direct simulations reveals several aspects of the transport process. First, the non-Gaussian pre-asymptotic regime is determined by a broad distribution of the total trapping times, which results from the coupling between the heterogeneities of the velocity field and the exchange time. This coupling is due to the fact that in low velocity zones, the particles have a greater probability of being trapped many times. An interesting result is that the correlation between the permeability field and the exchange time has no influence on the transport because the two quantities are drawn independently in the CTRW. This is also confirmed by the similarity between the results of $\gamma = -1$ and $\gamma = 1$.

Also, the CTRW predicts that transport is asymptotically Fickian. The evolution of the spatial variance is characterized by a broad pre-asymptotic cross-over behavior that may be characterized by a power-law. Finally, long lasting pre-asymptotic regimes could potentially explain the experimentally observed non-Gaussian behavior. Indeed, the experimental time is usually limited by the size of the sample. We were able to show that even for weakly heterogeneous media, the non-Gaussian pre-asymptotic regime could have a significant duration because of the heterogeneity of the exchange time.

In future works a challenging task will be to simulate transport in more realistic structures. Indeed, many geological structures are more complex and display a higher heterogeneity. Also, investigating transport in anisotropic or stratified structure will provide crucial knowledge for a better comprehension of pollutant transport in more complex subsurfaces. Moreover, in the present work, we used a constant and isotropic dispersion tensor, independent of the local velocity. This can be considered as a strong assumption. Indeed, if low velocity
regions are also characterized by a low dispersion coefficient, the residence time would increase in these regions. This coupling would then also affect the global transport. In a future work, it could therefore be interesting to evaluate the influence of a heterogeneous and velocity-dependent dispersion tensor.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Laurent Talon reports financial support was provided by CNRS. Marco Dentz reports financial support was provided by Spanish Ministry of Science and Innovation. Daniela Bauer reports financial support was provided by IFP Energies nouvelles.

## Data availability

Data will be made available on request.

## Acknowledgments

M.D. acknowledges the support of the Spanish Research Agency (10.13039/501100011033) and the Spanish Ministry of Science and Innovation through the grants CEX2018-000794-S and PID2019-106887GB-C31 (HydroPore).

## Appendix A. Correlation length of particle velocity

The correlation length of particle velocity can be determined from perturbation theory following Cvetkovic et al. (1996) and Dagan (1989). The longitudinal macrodispersion coefficient can be written as

$$D_L^* = \int_0^\infty dt C_v(t), \qquad (A.1)$$

where $C_v(t)$ is the Lagrangian velocity covariance. This expression can also be written as

$$D_L^* = \sigma_v^2 \tau_c, \quad \tau_c = \frac{1}{\sigma_v^2} \int_0^\infty dt C_v(t), \qquad (A.2)$$

where $\sigma_v^2$ is the velocity variance. The correlation time $\tau_c = \ell_c/\bar{u}$, where $\bar{u}$ is the mean flow velocity and $\ell_c$ is the correlation length. Thus, we obtain

$$\ell_c = \frac{D_L^* \bar{u}}{\sigma_v^2}. \qquad (A.3)$$

Perturbation theory gives $\sigma_v^2 = (3/8)\sigma_f^2 \bar{u}^2$. For the Gaussian covariance function employed in this study, the macrodispersion coefficient is

$$D_L^* = \frac{\sqrt{\pi}}{2} \lambda \bar{u} \sigma_f^2. \qquad (A.4)$$

Thus, we obtain for the correlation length $\ell_c$

$$\ell_c = \frac{\sqrt{\pi 4 \lambda}}{3}. \qquad (A.5)$$

## Appendix B. Numerical method: Two-relaxation-time Lattice Boltzmann method

The Darcy–Brinkmann equation as well as the equations of the Mobile–Immobile Model are solved using a Two-Relaxation-Time (TRT) Lattice Boltzmann scheme, which is summarized in the following. For more details, we refer to Ginzburg et al. (2008, 2015) considering the resolution of the velocity field and to Ginzburg and d'Humières (2010) for the Advection–Diffusion equation.

### B.1. Darcy–Brinkmann equation

In a first step, we provide the fundamentals of LB and the numerical scheme to solve the Darcy–Brinkmann equation given by

$$\frac{\rho \nu}{K(\bar{r})} \bar{u} - \bar{\nabla} P + \rho \nu \Delta \bar{u} = \bar{0}, \qquad (B.1)$$

where $\nu$ is the kinematic viscosity, $\rho$ the density of the fluid and $P$ the pressure field.

The principle of the Lattice Boltzmann Method (LBM) lies in the particular discretization of the Boltzmann equation. Both, space (2D in the present work) and particle velocities are discretized on a regular grid. The nodes of the grid are related by the velocity vectors $\vec{c}_q$ with $q \in [0, Q]$ ($Q = 9$ here) and particles can only move with velocities $\vec{c}_q$. The population $f_q(\bar{r}, t)$ is defined as the density of particles having the velocity $\vec{c}_q$ at the position $\bar{r}$. LBM schemes are based on the decomposition of the temporal evolution of the particle density in two specific steps. The first step is the collision, where the density populations $f_q$ meeting at the same position are redistributed according to a relaxation equation: $\bar{f}_q(\bar{r}, t) = f_q(\bar{r}, t) - s(f_q(\bar{r}, t) - f_q^{eq}(\rho, \bar{u}))$, where $\bar{f}_q$ denotes the post-collision populations. $s$ is a relaxation parameter and $f_q^{eq}$ designates an equilibrium state depending on the local macroscopic variables (density, velocity, permeability, etc.).

The second step is the “propagation”. In this step, all populations are displaced on the grid according to their velocity: $f_q(\bar{r} + \vec{c}_q, t + 1) = \bar{f}_q(\bar{r}, t)$.

The principle of the TRT scheme lies in the modification of the collision step by introducing two different relaxation parameters $s^+$ and $s^-$ for the odd and even components of the population densities. The even and odd components are defined as: $f_q^+ = (f_q + f_{\bar{q}})/2$ and $f_q^- = (f_q - f_{\bar{q}})/2$ where $\bar{q}$ denotes the opposite velocity direction to $q$. The collision becomes then:

$$\bar{f}_q(\bar{r}, t) = f_q - s^+ n_q^+ - s^- n_q^- \quad \text{with} \quad q = 0 \dots 4, \qquad (B.2)$$

$$\bar{f}_{\bar{q}}(\bar{r}, t) = f_{\bar{q}} - s^{+} n_{\bar{q}}^{+} + s^{-} n_{\bar{q}}^{-} \quad \text{with} \quad q = 0 \dots 4, \qquad (B.3)$$

with $n_q^\pm = (f_q^\pm - e_q^\pm)$. $e_q^\pm$ are the equilibrium functions which characterize the physical equations to be solved.

In the case of the Darcy–Brinkmann equation the equilibrium functions depend on the following local macroscopic quantities:

* Pressure:

$$P = c_s^2 \rho = c_s^2 \sum_{q \in [0, Q]} f_q,$$

where $c_s$ is the numerical sound speed ($c_s^2 = 1/3$, here).

* Local velocity:

$$\bar{u} = \frac{1}{\rho_0} \frac{2\bar{J}}{2 + \frac{\nu}{K}},$$

with

$$\bar{J} = \sum_{q \in [0, Q]} f_q \vec{c}_q.$$

The equilibrium functions are given by Ginzburg et al. (2015):

$$e_0 = \rho - 2 \sum_{q=1}^{4} e_q^+, \qquad (B.4)$$

$$e_q^+ = t_q P(\rho), \qquad (B.5)$$

$$e_q^- = \rho_0 t_q (\bar{u} \cdot \vec{c}_q - \Lambda^- \frac{\nu}{K} \vec{c}_q \cdot \bar{u}), \qquad (B.6)$$

where $t_q$ are coefficients depending on the direction. Here, we use $t_q = 1/3$ for the first neighbors and $t_q = 1/12$ for the second neighbors (diagonal). $\rho_0$ is a fixed averaged density ($\rho_0 = 1$, here) and $\Lambda^-$ a parameter depending on the odd relaxation time:

$$\Lambda^- = \frac{1}{s^-} - \frac{1}{2}. \qquad (B.7)$$

The scheme requires the determination of the two relaxation parameters $s^+$ and $s^-$. For numerical reasons (precision, stability), it is however more convenient to prescribe the viscosity $\nu$ and a numerical parameter
<figure><img src="image_16.png" /><figcaption>Fig. B.13. Examples of the permeability field (left) and velocity field (right) for the parameter $\sigma_f = 1$, $\lambda = 4$, $K_0 = 1$.</figcaption></figure>

$\Lambda = \Lambda_+ \Lambda_-$ The two-relaxation parameters are then deduced using the relations: $s^+ = \frac{2}{\delta v+1}$ and $s^- = \frac{2}{\frac{2\Lambda}{3\delta}+1}$.

$\Lambda$ essentially controls the numerical error of the scheme (Ginzburg et al., 2008, 2015), it is set to $\Lambda = 0.25$ in this work.

Advection-diffusion-reaction equation (ADRE). A TRT scheme is also used to solve the advection-diffusion-reaction equation:

$$ \frac{\partial C_m}{\partial t} + \vec{u} \cdot \vec{\nabla} C_m = D \Delta C_m + R, \qquad (B.8) $$

where $R = -\frac{1}{2}(C_m - C_{im})$ is the source term. The principle is identical to the TRT scheme used for the resolution of the flow equation. In this case $g_q(\vec{r}, t)$ represent the densities of particles having the velocity $\vec{c}_q$ at the position $\vec{r}$. The collision step is also similar:

$$ \vec{g}_q(\vec{r}, t) = g_q - s^+(g_q^+ - e_q^+) - s^-(g_q^- - e_q^-) \quad \text{with} \quad q = 0 \dots 4 \qquad (B.9) $$

$$ \vec{g}_q(\vec{r}, t) = g_q - s^+(g_q^+ - e_q^+) + s^-(g_q^- - e_q^-) \quad \text{with} \quad q = 0 \dots 4. \qquad (B.10) $$

In presence of a source term, the local concentration is given by:

$$ C_m = \sum_{q \in [0, Q]} g_q + \frac{R}{2}. \qquad (B.11) $$

Here also, it is convenient to introduce $\Lambda^\pm = \frac{1}{2} \pm \frac{1}{2}$. The main difference with the TRT scheme used for the resolution of the velocity field lies in the equilibrium functions $e_q^\pm$:

$$ e_q^+ = C^{eq} \left[ t_q c_s^2 + \frac{1}{4}(u_x^2 - u_y^2)p_q^{xx} + \frac{1}{4}u_x u_y p_q^{xy} + \frac{1}{12}(u_x^2 + u_y^2) \right] \quad \text{for } q = 1 \dots Q \qquad (B.12) $$

$$ e_q^- = t_q C^{eq}(\vec{u}, \vec{c}) \quad \text{for } q = 1 \dots Q \qquad (B.13) $$

$$ e_0^+ = C^{eq} - \sum_{q \in [1, Q]} e_q^+ \qquad (B.14) $$

where we defined $C^{eq} = C_m + \Lambda^+ R$, $p_q^{xx} = c_{qx}^2 - c_{qy}^2$, $p_q^{xy} = c_{qx}c_{qy}$. The diffusion coefficient is then given by $D = c_s^2 \Lambda^-$. The coefficient $c_s^2$ is a constant parameter and the directional coefficients $t_q$ are now equal to 1/4 for the first neighbors and 1/8 for the diagonal links. The numerical scheme thus requires the definition of the parameters $D$, $c_s^2$ and $\Lambda = \Lambda^+ \Lambda^-$. In this study, we used $c_s^2 = 0.2$ and $\Lambda = 0.2$.

Real time and space data are related to non-dimensional and Lattice Boltzmann variables as follows:

$$ t_{adim} = \frac{u_{LB}^2}{D_{LB}} t_{LB} \qquad (B.15) $$

$$ t_{real} = \frac{D_{reel}}{u_{real}^2} t_{adim} \qquad (B.16) $$

Simulation details. We used a porous medium of size 1024 × 512. Flow in x-direction is driven by a body force and periodic boundary conditions are applied in the x- and y-direction. We supposed $K_0 = 1$ in the Brinkmann equation. The boundary condition for the tracer are $c = 0$ at the inlet and $\frac{\partial C_m}{\partial x} = 0$ at the outlet. After reaching the steady state for the flow, we impose a pulse of concentration at a position slightly after the inlet: $x = 10$ px.

Fig. B.13 shows a zoom of the permeability field (Left) and the velocity field (Right). Figures are given for the following parameters: $\sigma_f = 1$, $\lambda = 4$, $K_0 = 1$.

## Appendix C. Homogeneous porous media

In this section, we present results of transport governed by the Mobile Immobile model in a homogeneous porous medium ($\sigma_f = 0$). In this case, analytical solutions of Eqs. (5) have been determined by Goltz and Roberts (1987), showing that transport becomes Gaussian (Fickian) for large time. By comparison, these results allow further comprehension of transport in heterogeneous porous media but also the validation of the Lattice Boltzmann algorithm.

For homogeneous porous media, the governing equations can be non-dimensionalized with $x' = x - \frac{u}{D_0}$ and $t' = t - \frac{u^2}{D_0}$. It follows:

$$ \left\{ \begin{array}{l} \frac{\partial C_m}{\partial t'} + \phi_D \frac{\partial C_{im}}{\partial t'} + \vec{\nabla}' C_m = \nabla^2 C_m \\ \frac{\partial C_{im}}{\partial t'} = \frac{D a}{P e} (C_m - C_{im}). \end{array} \right. \qquad (C.1) $$

The equations are now governed by only two dimensionless parameters, $\phi_D$ and $\Psi = Pe/Da = \frac{u^2 \tau \phi_D}{D_0}$.

Second spatial moment of the concentration profile and breakthrough curve. Fig. C.14 shows the averaged concentration profiles of the mobile and

<figure><img src="image_17.png" /><figcaption>Fig. C.14. Mobile (plain) and total (dashed) concentration profiles for the homogeneous porous medium at different times for $\Psi = Pe/Da = 25$ (Left).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. C.15. Evolution of the normalized dispersion coefficient $D_L^m(t)/D_0$ as function of time (homogeneous porous medium) for different values of $\Psi = Pe/Da$ (Left). For all parameters, transport becomes normal after a specific time, $\tau_{conv}$, which depends on $\Psi$ (Right).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. C.16. Normalized effective dispersion coefficient for a homogeneous porous medium as function of $\Psi$. Crosses correspond to the simulations and the red line stands for the analytical solution given in Eq. (C.3) (Left). Relative error of the dispersion coefficient as function of $\Psi$ (Right).</figcaption></figure>

the total concentration at different times. The corresponding breakthrough curve is given in Fig. C.14. For short times, the concentration profiles have a slightly asymmetrical bell shape, which is due to the initial condition (very sharp Gaussian distribution), but also to the fact that transport is dominated by molecular diffusion as the influence of trapping events on transport is negligible. Then, for slightly larger time intervals, the concentration profiles become asymmetric as the influence of trapping and advection on tracer transport becomes significant. Finally, in the large time limit, profiles become Gaussian again as steady state is reached. A slight shift can be seen between the mobile and the total concentration, that vanishes for large times.

Fig. C.15 (Left) shows the typical curves of the normalized dispersion coefficient $D_L^m(t)/D_0$ as a function of the non-dimensional time $\frac{u^2}{D_0}$ for different values of the dimensionless system parameter $\Psi = Pe/Da = \frac{u^2AK^\tau\phi_D}{D_0} = \frac{u^2\tau\phi_D}{D_0}$. For all $\Psi$, after a characteristic time ($\tau_{conv}$), the normalized dispersion coefficient reaches a constant value corresponding to a Gaussian asymptotic behavior. Here, the asymptotic dispersion coefficient $D_{eff}$ is defined by

$$D_{eff} = \lim_{t \to \infty} D_L^m(t). \qquad (C.2)$$

The asymptotic longitudinal dispersion coefficient $D_{eff}$ has been determined analytically by Goltz and Roberts (1987) as:

$$\frac{D_{eff}}{D_0} = \frac{1}{1 + \phi_D} + \frac{\Psi \phi_D}{(1 + \phi_D)^3}. \qquad (C.3)$$

Fig. C.16 (Left) shows the analytical and numerical values of $D_{eff}/D_0$. Very good accordance between analytical and numerical data was obtained. Fig. C.16 (Right) gives the dependency of the numerical

error on $\Psi$. Even for low values of $\Psi$ the error remains smaller than 2%.

Considering the shape of $D_L^m(t)/D_0$ (Fig. C.15, Left), two different behaviors can be observed. For small $\Psi$ (e.g. low exchange time), $D_L^m(t)/D_0$ decreases to reach a value below one. This can be explained by the fact that for small values of $\Psi$ the first term in Eq. (C.3) dominates. In this case, the equilibrium is very rapidly reached due to the small exchange time. This regime is the well-known retardation effect based on the assumption that $C_m = C_{im}$ in Eqs. (4).

For larger values of $\Psi$, where either the exchange time or the velocity is significant ($\tau u^2 > D_0/\phi_D$), $D_L^m(t)/D_0$ takes values close to one for small times due to the predominance of molecular diffusion. As trapping events start to be effective, $D_L^m(t)/D_0$ increases significantly, due to the competition between advection and retention, to reach a plateau after a certain time $\tau_{conv}$. Once the system has reached a plateau, transport is Gaussian.

The limiting case $\Psi = \infty$ deserves to be discussed because it corresponds to $\tau = \infty$ and thus to the absence of any exchange. In this case, the observed effective dispersion is necessarily constant and equal to $D_0$ but the theoretical dispersion coefficient given by Eq. (C.3) is infinite. This inconsistency is due to the fact that the convergence time $\tau_{conv}$ also diverges to infinity so that the system remains in the pre-asymptotic regime.

Fig. C.15 (Right) presents the normalized convergence time $\frac{u^2\tau_{conv}}{D_0}$ as a function of $\Psi$. For low values of $\Psi$ the convergence time is independent of $\Psi$ as the exchange time between mobile and immobile zones is very short and the equilibrium is reached quasi instantaneously. For larger $\Psi$ the convergence time increases approximately according to a power-law of $\Psi$. The power-law exponent is larger than one, but it is
expected that for very large $Ψ$, the convergence time is proportional to the exchange time ($τ_{conv} ∝ τ$), and therefore $τ_{conv} ∝ Ψ$. The latter is due to the fact that kinetics solely depend on the first order rate term $\frac{\partial C_{im}}{\partial t} = \frac{1}{\psi}(C_m - C_{im})$.

## Appendix D. Compound Poisson process

We derive here the Laplace space Eq. (26) for the distribution of the total trapping time in the compound Poisson process. To this end, we write the definition of the distribution of $θ_s$ as

$$p_c(t|\tau_v, \omega) = \left\langle \delta \left( t - \sum_{j=1}^{n} \theta_j \right) \right\rangle \quad (D.1)$$

$$= \sum_{n=0}^{\infty} \int \dots \int \delta \left( t - \sum_{j=1}^{n} \theta_j \right) p_n(\tau_v, \omega) \psi(\theta_1) \dots \psi(\theta_n) d\theta_1 \dots d\theta_n \quad (D.2)$$

where $n$ is Poisson distributed and the $θ_j$ are exponential distributed and independent. The Laplace transform of Eq. (D.1) is

$$p_c^*(s|\tau_v, \omega) = \int_0^\infty p_c(t|\tau_v, \omega) \exp(-st) dt \quad (D.3)$$

$$= \int \dots \int \sum_{n=0}^{\infty} p_n(\tau_v, \omega) \exp(-s \sum_{j=1}^{n} \theta_j) \psi(\theta_1) \dots \psi(\theta_n) d\theta_1 \dots d\theta_n \quad (D.4)$$

$$= \sum_{n=0}^{\infty} p_n(\tau_v, \omega) \int \dots \int \prod_{j=0}^{n} [\exp(-s\theta_j)\psi(\theta_j) d\theta_j] \quad (D.5)$$

$$= \sum_{n=0}^{\infty} p_n(\tau_v, \omega) \left[ \int \psi(\theta) \exp(-s\theta) d\theta \right]^n \quad (D.6)$$

Thus,

$$p_c^*(s|\tau_v, \omega) = \sum_{n=0}^{\infty} (\omega\tau_v)^n \frac{\exp(-\omega\tau_v)}{n!} \psi^*(s) = \exp\{-\omega\tau_v[1 - \psi^*(s)]\}, \quad (D.7)$$

where we used the exponential sum. We note that the Laplace transform or $ψ(t)$ given by Eq. (25) is

$$\psi^*(s) = \frac{\omega\phi_D}{\omega\phi_D + s} \quad (D.8)$$

Inserting this expression into (D.7) gives Eq. (26).

## References

Abramowitz, M., Stegun, I.A., 1972. Handbook of Mathematical Functions. Dover Publications, New York.

Attinger, S., Dentz, M., Kinzelbach, H., Kinzelbach, W., 1999. Temporal behavior of a solute cloud in a chemically heterogeneous porous medium. J. Fluid Mech. 386, 77–104.

Bateman, H., 1954. Tables of Integral Transforms. Vol. 1, McGraw-Hill Book Company, New York.

Beaudoin, A., de Dreuzy, J.R., 2013. Numerical assessment of 3-d macrodispersion in heterogeneous porous media. Water Resour. Res. 49, 2489–2496.

Benson, D.A., Meerschaert, M.M., 2009. A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations. Adv. Water Resour. 32, 532–539.

Berkowitz, B., Cortis, A., Dentz, M., Scher, H., 2006. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev. Geophys. 44.

Brinkman, H., 1947. A calculation of the viscous forces exerted by a flowing fluid on a dense swarm of particles. Appl. Sci. Res. Sect. A 1, 27–39.

Brusseau, M.L., 1994. Transport of reactive contaminants in heterogeneous porous media. Rev. Geophys. 32, 285–313.

Brusseau, M.L., Srivastava, R., 1999. Nonideal transport of reactive solutes in heterogeneous porous media: 4. analysis of the cape cod natural-gradient field experiment. Water Resour. Res. 35, 1113–1125.

Burr, D., Sudicky, E.A., Naff, R., 1994. Nonreactive and reactive solute transport in three-dimensional heterogeneous porous media: Mean displacement, plume spreading, and uncertainty. Water Resour. Res. 30, 791–815.

Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guimerà, J., 1998. On matrix diffusion: formulations, solution methods, and qualitative effects. Hydrogeol. J. 6, 178–190.

Coats, K., Smith, B., 1964. Dead-end pore volume and dispersion in porous media. SPE J. 1.

Comolli, A., Hakoun, V., Dentz, M., 2019. Mechanisms, upscaling, and prediction of anomalous dispersion in heterogeneous porous media. Water Resour. Res. 55, 8197–8222.

Comolli, A., Hidalgo, J.J., Moussey, C., Dentz, M., 2016. Non-Fickian transport under heterogeneous advection and mobile-immobile mass transfer. Transp. Porous Media 115, 265–289.

Cvetkovic, V., Cheng, H., Wen, X.H., 1996. Analysis of nonlinear effects on tracer migration in heterogeneous aquifers using Lagrangian travel time statistics. Water Resour. Res. 32, 1671–1680.

Cvetkovic, V., Dagan, G., Cheng, H., 1998. Contaminant transport in aquifers with spatially variable hydraulic and sorption properties. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 454, 2173–2207.

Cvetkovic, V., Fiori, A., Dagan, G., 2014. Solute transport in aquifers of arbitrary variability: A time-domain random walk formulation. Water Resour. Res. 50, 5759–5773.

Cvetkovic, V.D., Shapiro, A.M., 1990. Mass arrival of sorptive solute in heterogeneous porous media. Water Resour. Res. 26, 2057–2067.

Dagan, G., 1989. Flow and Transport in Porous Formations. Springer Berlin Heidelberg.

Dagan, G., Cvetkovic, V., 1993. Spatial moments of a kinetically sorbing solute plume in a heterogeneous aquifer. Water Resour. Res. 29, 4053–4061.

Dai, Z., Ritzi Jr., R.W., Huang, C., Rubin, Y.N., Dominic, D.F., 2004. Transport in heterogeneous sediments with multimodal conductivity and hierarchical organization across scales. J. Hydrol. 294, 68–86.

De Gennes, P.G., 1983. Hydrodynamic dispersion in unsaturated porous media. J. Fluid Mech. 136, 189–200.

Dentz, M., Castro, A., 2009. Effective transport dynamics in porous media with heterogeneous retardation properties. Geophys. Res. Lett 36 (L03403).

Dentz, M., Cortis, A., Scher, H., Berkowitz, B., 2004. Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport. Adv. Water Resour. 27, 155–173.

Dentz, M., Kang, P.K., Comolli, A., Le Borgne, T., Lester, D.R., 2016. Continuous time random walks for the evolution of lagrangian velocities. Phys. Rev. Fluids 1.

Dentz, M., Kinzelbach, H., Attinger, S., Kinzelbach, W., 2000. Temporal behavior of a solute cloud in a heterogeneous porous medium, 1, point-like injection. Water Resour. Res. 36, 3591–3604.

Dentz, M., Le Borgne, T., Englert, A., Bijeljic, B., 2011. Mixing, spreading and reaction in heterogeneous media: a brief review. J. Cont. Hydrol. 120–121, 1–17.

Feller, W., 1968. An Introduction to Probability Theory and Its Applications. In: Wiley Series in Probability and Statistics, Vol. 1, Wiley.

Fiori, A., Jankovic, I., Dagan, G., Cvetkovic, V., 2007. Ergodic transport trough aquifers of non-Gaussian log conductivity distribution and occurrence of anomalous behavior. Water Resour. Res. 43, W09407.

Frippiat, C.C., Holeyman, A.E., 2008. A comparative review of upscaling methods for solute transport in heterogeneous porous media. J. Hydrol. 362, 150–176.

Gelhar, L., Axness, C., 1983. Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resour. Res. 19, 161–180.

Ginzburg, I., d'Humières, A., 2010. Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium. J. Stat. Phys. 139, 1090–1143.

Ginzburg, I., Silva, G., Talon, L., 2015. Analysis and improvement of brinkman lattice Boltzmann schemes: Bulk, boundary, interface, similarity and distinctness with finite elements in heterogeneous porous media. Phys. Rev. E 91, 023307.

Ginzburg, I., Verhaeghe, F., d'Humières, D., 2008. Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions. Commun. Comput. Phys. 3, 427–478.

Goltz, M.N., Roberts, P.V., 1987. Using the method of moments to analyze three-dimensional diffusion-limited solute transport from temporal and spatial perspectives. Water Resour. Res. 23, 1575–1585.

Gotovac, H., Cvetkovic, V., Andricevic, R., 2009. Flow and travel time statistics in highly heterogeneous porous media. Water Resour. Res. 45.

Gouze, P., Le Borgne, T., Leprovost, R., Lods, G., Poidras, T., Pezard, P., 2008. Non-fickian dispersion in porous media: 1. multiscale measurements using single-well injection withdrawal tracer tests. Water Resour. Res. 44.

Griffioen, J., Barry, D., Parlange, J.Y., 1998. Interpretation of two-region model parameters. Water Resour. Res. 34, 373–384.

Haggerty, R., Gorelick, S., 1995. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Ressour. Res. 310, 2383–2400.

Hakoun, V., Comolli, A., Dentz, M., 2019. Upscaling and prediction of lagrangian velocity dynamics in heterogeneous porous media. Water Resour. Res. 55, 3976–3996.

Hidalgo, J.J., Neuweiler, I., Dentz, M., 2020. Transport under advective trapping. J. Fluid Mech. 907.

Hyman, J.D., Dentz, M., 2021. Transport upscaling under flow heterogeneity and matrix-diffusion in three-dimensional discrete fracture networks. Adv. Water Resour. 155, 103994.

Koplik, J., Redner, S., Wilkinson, D., 1988. Transport and dispersion in random networks with percolation disorder. Phys. Rev. A 37 (2619).

Margolin, G., Dentz, M., Berkowitz, B., 2003. Continuous time random walk and multirate mass transfer modeling of sorption. Chem. Phys. 295, 71–80.

Massabó, M., Bellin, A., Valocchi, A.J., 2008. Spatial moments analysis of kinetically sorbing solutes in aquifer with bimodal permeability distribution. Water Resour. Res. 44.
Metzger, D., Kinzelbach, H., Kinzelbach, W., 1996. Effective dispersion of a solute cloud in a chemically heterogeneous porous medium: Comparison of two ensemble-averaging procedures. Water Resour. Res. 32, 3311–3319.

Miralles-Wilhelm, F., Gelhar, L.W., 1996. Stochastic analysis of sorption macrokinetics in heterogeneous aquifers. Water Resour. Res. 32, 1541–1549.

Morales, V.L., Dentz, M., Willmann, M., Holzner, M., 2017. Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory. Geophys. Res. Lett. 44, 9361–9371.

Neuman, S.P., Tartakovsky, D.M., 2008. Perspective on theories of anomalous transport in heterogeneous media. Adv. Water Resour. 32, 670–680.

Nkedi-Kizza, P., Biggar, J.W., Selim, H.M., Genuchten, M.T.V., Wierenga, P.J., Davidson, J.M., Nielsen, D.R., 1984. On the equivalence of two conceptual models for describing ion exchange during transport through an aggregated oxisol. Water Resour. Res. 20, 1123–1130.

Noetinger, B., Roubinet, D., Russian, A., Le Borgne, T., Delay, F., Dentz, M., De Dreuzy, J.R., Gouze, P., 2016. Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale. Transp. Porous Media 1–41.

Painter, S., Cvetkovic, V., 2005. Upscaling discrete fracture network simulations: An alternative to continuum transport models. Water Resour. Res. 41.

Painter, S., Cvetkovic, V., Mancillas, J., Pensado, O., 2008. Time domain particle tracking methods for simulating transport with retention and first-order transformation. Water Resour. Res. 44.

Quinodoz, H.A., Valocchi, A.J., 1993. Stochastic analysis of the transport of kinetically sorbing solutes in aquifers with randomly heterogeneous hydraulic conductivity. Water Resour. Res. 29, 3227–3240.

Rajaram, H., 1997. Time and scale dependent effective retardation factors in heterogeneous aquifers. Adv. Water Resour. 20, 217–230.

Rao, P., Rolston, D., Jessup, R., Davidson, J., 1980. Solute transport in aggregated porous media: Theoretical and experimental evaluation. Soil Sci. Am. J. 44, 1139–1146.

Reichle, R., Kinzelbach, W., Kinzelbach, H., 1998. Effective parameters in heterogeneous and homogeneous transport models with kinetic sorption. Water Resour. Res. 34, 583–594.

Ren, W., Ershadnia, R., Wallace, C.D., LaBolle, E.M., Dai, Z., de Barros, F.P., Soltanian, M.R., 2022. Evaluating the effects of multiscale heterogeneous sediments on solute mixing and effective dispersion. Water Resour. Res. 58, e2021WR031886.

Rubin, Y., 2003. Applied Stochastic Hydrogeology. Oxford University Press, New York.

Selroos, J.O., Cvetkovic, V., 1992. Modeling solute advection coupled with sorption kinetics in heterogeneous formations. Water Resour. Res. 28, 1271–1278.

Selroos, J.O., Cvetkovic, V., 1994. Mass flux statistics of kinetically sorbing solute in heterogeneous aquifers: Analytical solution and comparison with simulations. Water Resour. Res. 30, 63–69.

Smettem, K., 1984. Soil-water residence time and solute uptake: 3. mass transfer under simulated winter rainfall conditions in undisturbed soil cores. J. Hydrol. 67, 235–248.

Soltanian, M.R., Ritzi, R., Dai, Z., Huang, C., Dominic, D., 2015. Transport of kinetically sorbing solutes in heterogeneous sediments with multimodal conductivity and hierarchical organization across scales. Stoch. Environ. Res. Risk Assess. 29, 709–726.

Talon, L., Martin, J., Rakotomalala, N., Salin, D., Yortsos, Y., 2003. Lattice BGK simulations of macrodispersion in heterogeneous porous media. Water Resour. Res. 39, 1135–1142.

Tompson, A.F., 1993. Numerical simulation of chemical migration in physically and chemically heterogeneous porous media. Water Resour. Res. 29, 3709–3726.

Valocchi, A.J., 1990. Use of temporal moment analysis to study reactive solute transport in aggregated porous media. Geoderma 46, 233–247.

van Genuchten, M.T., Wierenga, P., 1976. Mass transfer studies in sorbing porous media i. analytical solution. Soil Sci. Am. J. 4, 473–480.

Willmann, M., Carrera, J., Sanchez-Vila, X., 2008. Transport upscaling in heterogeneous aquifers: What physical parameters control memory functions? Water. Resour. Res. 44, W12437.

Zhang, Y., Green, C.T., Baeumer, B., 2014. Linking aquifer spatial properties and non-Fickian transport in mobile-immobile like alluvial settings. J. Hydrol. 512, 315–331.

Zinn, B., Meigs, L.C., Harvey, C.F., Haggerty, R., Peplinski, W.J., Freiherr von Schwörerin, C., 2004. Experimental visualization of solute transport and mass transfer processes in two-dimensional conductivity fields with connected regions of high conductivity. Environ. Sci. Technol. 38, 3916–3926.