# Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields

Tanguy Le Borgne, Marco Dentz, Jesus Carrera

► To cite this version:

Tanguy Le Borgne, Marco Dentz, Jesus Carrera. Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields. Physical Review Letters, 2008, 101, pp.090601. 10.1103/Phys-RevLett.101.090601. insu-00372030

HAL Id: insu-00372030

https://insu.hal.science/insu-00372030

Submitted on 13 Jan 2022

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
Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields

Tanguy Le Borgne,¹ Marco Dentz,² and Jesus Carrera³

¹Geosciences Rennes, UMR 6118, CNRS, Université de Rennes 1, Rennes, France

²Department of Geotechnical Engineering and Geosciences, Technical University of Catalonia (UPC), Barcelona, Spain

³Institute of Environmental Analysis and Water Studies (IDAEA), CSIC, Barcelona, Spain

(Received 27 May 2008; published 26 August 2008)

We define an effective Lagrangian statistical model in phase space (x, t, v) for describing transport in highly heterogeneous velocity fields with complex spatial organizations. The spatial Markovian nature (and temporal non-Markovian nature) of Lagrangian velocities leads to an effective transport description that turns out to be a correlated continuous time random walk. This model correctly captures the Lagrangian velocity correlation properties and is demonstrated to represent a forward model for predicting transport in highly heterogeneous porous media for different types of velocity organizations.

DOI: 10.1103/PhysRevLett.101.090601

PACS numbers: 05.40.Fb, 05.10.Gg, 05.60.-k, 92.40.Kf

Natural flow fields generally display a complex spatio-temporal organization due to heterogeneous geological structures or coherent turbulent structures at different scales, for example, which leads to non-Fickian transport properties (e.g., [1,2]). New experimental developments have recently allowed for significant progress in the modeling of transport in turbulent flows by detailed analysis of the Lagrangian correlation properties (e.g., [3,4]). The resulting effective Lagrangian statistical models often rely on the assumption that the Lagrangian velocities taken at equidistant times describe a Markov process, which leads to effective random walk models that are Markovian in time (e.g., [1,5]). For transport in highly heterogeneous porous media, the direct measurement of Lagrangian velocities is still beyond our capabilities, although some progress has been made recently [6]. Such analysis has been performed numerically for heterogeneous conductivity fields with conductivity values that fluctuate spatially over 10 orders of magnitude [7,8]. This type of medium displays a complex velocity field organization that is characterized by the localization of flow in connected high velocity channels (Fig. 1). The flow organization is manifested by special Lagrangian velocity correlation properties, such as the dependence of the correlation on the local velocity and the possibility for particles to experience abrupt changes in velocities, i.e., large Lagrangian accelerations. As demonstrated in [8], this implies that velocities at equidistant positions along a particle trajectory form a Markov process, which means that here Lagrangian velocities are Markovian in space and not in time. This fundamental property leads naturally to a representation of effective particle motion as a random walk in space-time,

$$x^{(n+1)} = x^{(n)} + \Delta x \quad (1a)$$

$$t^{(n+1)} = t^{(n)} + \frac{\Delta x}{v^{(n)}}, \quad (1b)$$

where the series of successive particle longitudinal veloc-

ities {$v^{(n)}$}$_{n=0}^{\infty}$ at equal longitudinal distances $\Delta x$ along the particle trajectory form a Markov chain in space, whose transitions are characterized by the conditional probability density $r(v, x|v', x')$. The probability for a particle to make a transition from $v'$ at travel distance $x'$ to $v$ at travel distance $x = x' + \Delta x$ is given by (e.g., [9])

$$r(v, x|v', x')dv = \langle \delta(v - v^{(n+1)}) \rangle_{v^{(n)}=v'} dv, \quad (2)$$

where the angular brackets denote the average overall realization of {$v^{(n)}$}$_{n=0}^{\infty}$. For stationary velocity fields, transition distributions are stationary in space, i.e., $r(v, x|v', x') = r(v, x - x'|v')$.

The random walk (1) is a Markov process in the phase space (x, t, v). It can be viewed as a correlated continuous time random walk (CTRW). Unlike classical CTRW mod-

<figure><img src="image_2.png" /><figcaption>FIG. 1 (color online). Example of velocity field realization for the multilognormal conductivity field: spatial distribution of the absolute value of the longitudinal velocity m/s.</figcaption></figure>
els (e.g., [2,10–12]), here the successive particle velocities are not independent, but depend on the velocity at the previous step. Correlated CTRW models have been recently studied in the literature (e.g., [13,14]). While these studies postulate some correlation, here the correlated CTRW approach derives from the Markovian nature of the spatial (Lagrangian) velocity transitions [8]. The proposed Lagrangian statistical model is shown to provide accurate predictions for transport in the preasymptotic non-Fickian transport regime for different velocity organizations. It describes a wide range of transport behaviors, which can be related directly to the Lagrangian velocity correlation properties and thus to the velocity field organization. The effect of including velocity correlation is studied by comparison with the predictions of the corresponding classical CTRW model.

The particle distribution in phase space after n steps of the correlated CTRW is defined by (e.g., [9])

$$f_n(x, t, v) = \langle \delta(x - x^{(n)}) \delta(t - t^{(n)}) \delta(v - v^{(n)}) \rangle, \quad (3)$$

where angular brackets denote the ensemble average over all realizations of the Markov chain {$v^{(n)}$}$_{n=0}^{\infty}$. The Chapman-Kolmogorov equation for $f_n(x, t, v)$ reads as

$$f_{n+1}(x, t, v) = \int dx' \int_0^t dt' \times \int_0^\infty dv' f(x, t, v|x', t', v') f_n(x', t', v'), \quad (4)$$

where the transition probability is defined by (e.g., [9])

$$f(x, t, v|x', t', v') = \langle \delta(x - x^{(n+1)}) |_{x^{(n)}=x'} \delta(t - t^{(n+1)}) |_{t^{(n)}=t'} \delta(v - v^{(n+1)}) |_{v^{(n)}=v'} \rangle. \quad (5)$$

Using (1) and (2) gives

$$f(x, t, v|x', t', v') = \delta(x - x' - \Delta x) \times \delta(t - t' - \Delta x/v') r(v, \Delta x|v'). \quad (6)$$

Inserting (6) into (4), we obtain

$$f_{n+1}(x, t, v) = \int_0^\infty dv' r(v, \Delta x|v') \times f_n(x - \Delta x, t - \Delta x/v', v'). \quad (7)$$

We define the probability per time to just arrive at $(x, v)$, $R(x, v, t)$, by summing $f_n(x, t, v)$ over all n [15],

$$R(x, t, v) = \sum_{n=0}^{\infty} f_n(x, t, v). \quad (8)$$

Summing (7) over n gives for $R(x, t, v)$,

$$R(x, t, v) = f_0(x, t, v) + \int_0^\infty dv' r(v, \Delta x|v') R(x - \Delta x, t - \Delta x/v', v'), \quad (9)$$

where $f_0(x, t, v)$ is the initial distribution of the particles in phase space. We assume here that once a particle arrives at

the phase space position $(x, v)$, it remains at this position for a waiting time $\tau = \Delta x/v$. Thus, the probability $P(x, t, v)$ for a particle to be at the phase space position $(x, v)$ at time $t$ is given by the probability per time to arrive at $(x, v)$ at some earlier time $t' < t$ and to remain there until time $t$,

$$P(x, t, v) = \int_0^t dt' \Theta(\Delta x/v - t') R(x, t - t', v), \quad (10)$$

where $\Theta(t)$ is the Heaviside step function. The particle distribution in space is obtained from $P(x, t, v)$ by integration over all velocities,

$$c(x, t) = \int_0^\infty dv P(x, t, v). \quad (11)$$

For independent successive velocities, i.e., $r(v, x|v', x') = p(v)$, where $p(v)$ is the Lagrangian velocity distribution, the probability distribution in phase space can be decoupled, i.e., $R(x, t, v) = R(x, t)p(v)$ and $P(x, t, v) = P(x, t)p(v)$. Replacing these expressions in Eqs. (9) and (10), the classical CTRW formulation is recovered [15].

In the following, we demonstrate that the proposed Lagrangian statistical model (1) and (2) describes correctly the effective movement of particles in heterogeneous porous media, which are characterized by a spatially varying conductivity $K(x)$. To this end, we compare the solutions of the local scale flow and transport problem to the effective Lagrangian statistical model presented above.

Incompressible fluid flow through a porous medium is described by the Darcy equation. The particular media we focus on are ubiquitously used for modeling solute transport in geological media. The first type of field (Fig. 1) is a multilognormally distributed random conductivity field with a broad range of permeabilities (variance of the log-conductivity $\sigma_{\ln K}^2 = 9$) and a Gaussian correlation function. The second type of conductivity field contains preferentially connected high conductivity zones [8,16], which induces a greater localization of high velocity zones. The two conductivity fields have the same point distributions and similar two point correlation functions, characterized by the correlation length $\lambda$ (here, we set $\lambda = 8$ pixels).

Particle motion in a single realization is described by advection in the flow through the medium and local diffusion. The numerical method for solving flow and transport is described in detail in [8]. The Peclet number, defined by the geometric mean conductivity $K_g$, pixel size $d$, and diffusion coefficient $D$ as $Pe = \frac{dK_g}{D}$, here is set to $Pe = 100$. The mean hydraulic gradient is set to one. The effective (average) transport behavior is obtained by averaging over 100 conductivity field realizations.

The large disorder variance implies the existence of a preasymptotic non-Fickian regime that is characterized by long tails in the first passage time distributions and a preasymptotic non-Fickian scaling of the second centered moment of the particle positions in the direction of the mean flow $\sigma_x^2(t)$ (Figs. 2 and 3). For such high degree of
<figure><img src="image_3.png" /><figcaption>FIG. 2. Comparison of the numerical simulations for the multi-lognormal field with the predictions given by the Lagrangian statistical model and the standard CTRW model. (a) First passage time distribution at a distance of 100 elements (12.5λ) from the inlet (b) variance of the spatial concentration distribution σ<sub>x</sub><sup>2</sup>(t) as a function of time.</figcaption></figure>

heterogeneity, very large grids are needed to assess the full preasymptotic regime [17]. The simulations presented in this Letter address the first part of the preasymptotic regime, for which we find a transition between ballistic motion [σ<sub>x</sub><sup>2</sup>(t) ∝ t<sup>2</sup>] and anomalous transport [σ<sub>x</sub><sup>2</sup>(t) ∝ t<sup>1.3</sup>]. The presence of connected high conductivity patterns, Fig. 3, implies significantly earlier particle breakthrough (about 1 order of magnitude in Fig. 3) and stronger tailing than for the multilognormal field.

The definition of the Lagrangian statistical model (1) as an effective transport model for such media requires the characterization of the spatial velocity transition densities. Ideally, this could be performed analytically. Here, we define it numerically. Notice that solving the full flow and transport on the whole field can be very computationally demanding (we have done it solely for comparison with our effective model). On the other hand, the definition of the spatial velocity transition densities only requires computing transition probabilities over one half of the

<figure><img src="image_4.png" /><figcaption>FIG. 3. Comparison of the numerical simulations for the connected field with the predictions given by the Lagrangian statistical model and the standard CTRW model. (a) First passage time distribution at a distance of 100 elements (12.5λ) from the inlet (b) variance of the spatial concentration distribution σ<sub>x</sub><sup>2</sup>(t) as a function of time.</figcaption></figure>

permeability correlation length [8], whose computational effort is comparatively negligible. To evaluate the spatial Lagrangian velocity transition probabilities (2), we discretize the longitudinal velocity distribution into n classes {C<sub>i</sub>}<sub>1≤i≤n=50</sub> following a logarithmic discretization of the velocity scale. The probability for a particle to make a transition from a longitudinal velocity v' ∈ [v<sub>j</sub>, v<sub>j+1</sub>] at location x to a longitudinal velocity v ∈ [v<sub>j</sub>, v<sub>j+1</sub>] at location x + Δx is given by

$$T_{ij}(\Delta x) = \int_{v_i}^{v_{i+1}} dv \int_{v_j}^{v_{j+1}} dv' r(v, \Delta x|v') p(v')/P_j, \quad (12)$$

where P<sub>j</sub> = ∫<sub>v<sub>j</sub></sub><sup>v<sub>j+1</sub></sup> dv'p(v') is the probability of class C<sub>j</sub>.

T(Δx) represents a velocity transition probability matrix for an incremental length Δx. Note that the longitudinal Lagrangian velocity v that we consider here includes advective and diffusive motions. The specific methodology for determining T(Δx) is given in [8]. For the spatial
increment, we use $\Delta x = \lambda/2$, for which the spatial Markov property is verified [8]. The velocity transition matrix quantifies the dependence of the velocity correlation on the local velocity, such as the increased persistence of high velocities due to connected high velocity channels. It also characterizes the non-negligible probabilities for particles to perform large jumps in the velocity space, such as the transition from a very low velocity to a high velocity when a particle escapes from a stagnation zone to join a high velocity channel (Fig. 1). These properties are not quantified in classical measures of correlation such as the Lagrangian velocity autocorrelation function. By construction and due to the stationarity of the velocity fields considered here, the Lagrangian velocity distribution is an eigenvector of the velocity transition matrix, i.e., $P_i = \sum_{j=0}^{n} T_{ij}(\Delta x)P_j$.

The Lagrangian velocity distribution $p(v)$ and the transition matrix $T(\Delta x)$ together with the spatial Markovian property completely define the Lagrangian statistical model for $\{v^{(n)}\}_{n=0}^{\infty}$ in (1). Using this effective description, we make predictions of the transport behavior over a large range of temporal and spatial scales. The equations of motion (1) of a particle are solved numerically using random walk particle tracking in phase space $(x, t, v)$, which allows for efficient transport simulations. We compare the predictions of this effective random walk model with the numerical random walk simulations of transport through the fully resolved two-dimensional velocity fields. To probe the role of correlation, we also compare the transport behavior resulting from the correlated CTRW with the predictions of a classical CTRW model defined by (1) for the transition probability $r(v, x|v', x') = p(v)$, i.e., independent identically distributed velocities. The model predictions are obtained without fitting of the model parameters.

The proposed Lagrangian statistical model is found to provide very good predictions of first passage time distributions at all distances, including the non-Fickian late time tailing and the early time particle arrival as illustrated in Figs. 2(a) and 3(a). It also predicts the scaling of the second centered moment $\sigma_x^2(t)$ at all times [Figs. 2(b) and 3(b)]. In particular, the scaling $\sigma_x^2(t) \propto t^{1.3}$ in the preasymptotic regime is correctly predicted.

Conversely, the standard CTRW model is found to overestimate the early arrival times by about 1 order of magnitude and to underpredict slightly the largest arrival times [Figs. 2(a) and 3(a)]. It underestimates the second centered moment $\sigma_x^2(t)$ by 1 order of magnitude at large times [Figs. 2(b) and 3(b)]. It predicts a slightly different scaling $\sigma_x^2(t) \propto t^{1.4}$, which is determined by the scaling of the velocity distribution in the low velocity part, $p(v) \propto v^{0.6}$. The scaling predicted by the classical CTRW model is entirely determined by the velocity distribution while the observed scaling $\sigma_x^2(t) \propto t^{1.3}$ is due to both the scaling of

the velocity distribution and the Lagrangian correlation properties.

As a conclusion, the good agreement of the proposed Lagrangian statistical model and the fully resolved numerical simulations demonstrates that the spatial (Lagrangian) velocity transition densities contain all the transport relevant information about the complex organization of the velocity fields studied here. The fact that the model turns out to be a correlated CTRW model is remarkable given the broad use of CTRWs in the description of anomalous transport (e.g., [2, 10, 11, 12]). Note that CTRW is often based on a phenomenological understanding of the small scale behavior. The approach that we propose provides an explicit link between small scale characterization (as represented by the Lagrangian velocity statistics) and large scale dynamics. Here, the (correlated) CTRW representation as a (coupled) random walk in position, time, and velocity space is a direct result of the spatial Markovianity of the Lagrangian velocity transitions. This finding may be a significant step towards the quantification of effective transport in random media and geological flows.

Tanguy Le Borgne acknowledges the financial support of European FP6 Marie Curie IEF Action No. 023858. Marco Dentz acknowledges the financial support of the Spanish MEC program Ramon y Cajal.

[1] P. S. Berloff and J. C. McWilliams, J. Phys. Oceanogr. **32**, 797 (2002).

[2] B. Berkowitz, A. Cortis, M. Dentz, and H. Scher, Rev. Geophys. **44**, RG2003 (2006).

[3] N. Mordant, A. M. Crawford, and E. Bodenschatz, Phys. Rev. Lett. **93**, 214501 (2004).

[4] H. Xu, N. Ouellette, and E. Bodenschatz, Phys. Rev. Lett. **96**, 114503 (2006).

[5] C. Beck, Phys. Rev. Lett. **98**, 064502 (2007).

[6] N. Kleinfelter, M. Moroni, and J. H. Cushman, Phys. Rev. E **72**, 056306 (2005).

[7] T. Le Borgne, J. R. de Dreuzy, P. Davy, and O. Bour, Water Resour. Res. **43**, W02419 (2007).

[8] T. Le Borgne, M. Dentz, and J. Carrera, Phys. Rev. E **78**, 026308 (2008).

[9] H. Risken, *The Fokker-Planck Equation* (Springer, Heidelberg, New York, 1996).

[10] M. F. Shlesinger, B. J. West, and J. Klafter, Phys. Rev. Lett. **58**, 1100 (1987).

[11] B. Berkowitz and H. Scher, Phys. Rev. Lett. **79**, 4038 (1997).

[12] R. Metzler and J. Klafter, Phys. Rep. **339**, 1 (2000).

[13] M. Montero and J. Masoliver, Phys. Rev. E **76**, 061115 (2007).

[14] H. Isliker, arXiv:0710.5152.

[15] H. Scher and M. Lax, Phys. Rev. B **7**, 4491 (1973).

[16] B. Zinn and C. F. Harvey, Water Resour. Res. **39**, 1051 (2003).

[17] J.-R. De Dreuzy, A. Beaudoin, and J. Erhel, Water Resour. Res. **43**, W10439 (2007).