Tomás Aquino

Tanguy Le Borgne

The diffusing-velocity random walk: A spatial-Markov formulation of heterogeneous advection and diffusion

The spatial distribution of a scalar undergoing advection and diffusion is impacted by the velocity variability sampled by tracer particles. In spatially structured velocity fields, such as porous medium flows, Lagrangian velocities along streamlines are often characterized by a well-defined correlation length and can thus be described by spatial-Markov processes. Diffusion, on the other hand, is generally modeled as a temporal process, making it challenging to capture advective and diffusive dynamics in a single framework. We develop a description of transport based on a spatial-Markov velocity process along Lagrangian particle trajectories, incorporating the effect of diffusion as a local averaging process in velocity space. The impact of flow structure on this diffusive averaging is quantified through an effective shear rate. The latter is fully determined by the point statistics of velocity magnitudes together with characteristic longitudinal and transverse lengthscales associated with the flow field. For infinite longitudinal correlation length, our framework recovers Taylor dispersion, and in the absence of diffusion it reduces to a standard spatial-Markov velocity model. We derive dynamical equations governing the evolution of particle position and velocity and obtain scaling laws for the dependence of longitudinal dispersion on Péclet number. These results provide new insights into the role of shear and diffusion on dispersion processes in heterogeneous media.

Introduction

Transport processes in heterogeneous media are determined by the sampling of the underlying heterogeneous flow field through advection and diffusion, leading to rich dynamical behavior and departure from classical Fickian dynamics [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Klages | Anomalous transport: foundations and applications[END_REF]. The evolution of Lagrangian velocities along particle trajectories can be modeled as a stochastic process, taking into account the statistical properties of the underlying flow field and diffusion [START_REF] Pope | Simple models of turbulent flows[END_REF][START_REF] Sund | Effective models for transport in complex heterogeneous hydrologic systems[END_REF]. Such approaches are greatly simplified if the changes in velocity may be conceptualized as a Markov process, that is, if their evolution depends only on the current state and not on past history [START_REF] Meyer | Particle-based transport model with Markovian velocity processes for tracer dispersion in highly heterogeneous porous media[END_REF][START_REF] Meyer | Testing the Markov hypothesis in fluid flows[END_REF]. In spatially structured velocity fields, such as porous medium flows, it has been found that the Lagrangian velocity structure often follows spatial-Markov dynamics (Le [START_REF] Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF][START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]Puyguiraud et al., 2019b). This means that Lagrangian velocities vary little over spatial scales below a characteristic lengthscale corresponding to a well-defined velocity-field correlation length. In such cases, low velocities persist for longer times than high velocities, because Lagrangian particles take longer times to cross a correlation length at lower velocities. In terms of the temporal Lagrangian velocity statistics, this phenomenon may lead to intermittent velocity timeseries and consequent loss of the Markov property in time [START_REF] Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[END_REF][START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF].

Stochastic Lagrangian methods describe transport in terms of random particle displacements and associated transit times [START_REF] Sund | Effective models for transport in complex heterogeneous hydrologic systems[END_REF]. The stochastic character of these models reflects the statistical properties of the underlying heterogeneity, which can be conceptualized in different manners. For 1 example, time domain random walks consider transport in single-medium realizations with prescribed statistical properties [START_REF] Mccarthy | Continuous-Time random walks on random media[END_REF][START_REF] Banton | A new time domain random walk method for solute transport in 1-D heterogeneous media[END_REF][START_REF] Delay | Time domain random walk method to simulate transport by advectiondiffusion and matrix diffusion in fracture networks[END_REF][START_REF] Painter | Time domain particle tracking methods for simulating transport with retention and first-order transformation[END_REF][START_REF] Russian | Time domain random walks for hydrodynamic transport in heterogeneous media[END_REF][START_REF] Aquino | A coupled time domain random walk approach for transport in media characterized by broadly-distributed heterogeneity length scales[END_REF], whereas in a continuous time random walk successive displacements are independent with statistics determined by medium heterogeneity [START_REF] Scher | Stochastic transport in a disordered solid. I. Theory[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF]. In this context, spatial-Markov descriptions often lead to significant simplifications which facilitate analytical treatment, parameterization based on physical, measurable quantities, and efficient numerical simulation. Spatial-Markov models may be seen as a variation on the continuous time random walk representation of Lagrangian particle movement, with a fixed spatial step and a one-step correlation between successive waiting times or particle velocities (Le [START_REF] Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF][START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]. In some cases, simple Markov processes such as Bernoulli relaxation or Ornstein-Uhlenbeck for the spatial evolution of Lagrangian velocities have been shown to capture the key features of purely advective transport, leading to efficient methods for predicting larger-scale transport properties such as longitudinal dispersion (Puyguiraud et al., 2019a;[START_REF] Comolli | Mechanisms, upscaling, and prediction of anomalous dispersion in heterogeneous porous media[END_REF]. In recent years, spatial-Markov models have also been extensively employed to describe conservative transport in porous media [START_REF] Kang | Spatial Markov model of anomalous transport through random lattice networks[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[END_REF], fractured media [START_REF] Kang | Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model[END_REF][START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF], surface flows [START_REF] Sherman | Parameterizing the spatial Markov model from breakthrough curve data alone[END_REF], and inertial and turbulent flows [START_REF] Bolster | Modeling preasymptotic transport in flows with significant inertial and trapping effects -the importance of velocity correlations and a spatial Markov model[END_REF][START_REF] Sund | Pre-asymptotic transport upscaling in inertial and unsteady flows through porous media[END_REF][START_REF] Kim | Anomalous transport through free-flow-porous media interface: Pore-scale simulation and predictive modeling[END_REF], as well as mixing and reaction (Sund et al., 2017b,a;[START_REF] Wright | Upscaling mixing in highly heterogeneous porous media via a spatial Markov model[END_REF][START_REF] Sherman | A spatial Markov model for upscaling transport of adsorbing-desorbing solutes[END_REF]. Despite the popularity and practical success of spatial-Markov methods over the last decade, a mechanistic model of the role of diffusion in this type of framework remains unavailable. In applications, the transition probabilities characterizing the spatial evolution of Lagrangian velocities in the presence of both advection and diffusion are typically parameterized based on small-scale simulations, and the resulting model is then applied to predict large-scale transport.

As a scalar tracer is transported through a heterogeneous medium, the statistics of velocity sampled by the tracer plume evolve in space and time as the underlying heterogeneity is sampled by moving tracer particles [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]Puyguiraud et al., 2019a;[START_REF] Icardi | Probability density function (pdf) models for particle transport in porous media[END_REF]. Under purelyadvective transport, particles experience this variability as they move along streamlines. In the presence of diffusion, a tracer particle is not confined to a single streamline and also experiences the variability across streamlines. In the classical Brownian motion picture, diffusion is modeled as a temporal process, and coupling advective space-Markovian Lagrangian velocity dynamics with diffusion remains an open problem. Previous approaches [START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF] have introduced a heuristic diffusive cutoff at the level of the crossing times. In this work, we provide an explicit construction of a spatial-Markov velocity process accounting for diffusion and heterogeneous advection.

Diffusion across nearby streamlines leads to a local averaging over the transverse spatial structure of the velocity field. Thus, transverse diffusion in real space translates to an averaging effect in velocity space. We quantify the impact of the spatial structure on this averaging process through an effective shear rate, characterized in terms of flow properties. Transport along the longitudinal direction is then described in terms of equidistant spatial steps along particle trajectories, together with transit times according to the velocity process. Due to the nature of the velocity transitions induced by diffusion, which, as we will show, correspond to a dispersive process in velocity space, we name the approach the diffusing-velocity random walk (DVRW). A conceptual illustration of the DVRW approach for advective-diffusive transitions is presented in Fig. 1.

The general outline of the paper is as follows. First, in Section 2, we discuss the description of transport as a spatial-Markov process in general terms. Section 3 is devoted to the formulation of the DVRW approach, incorporating the impact of advective and diffusive velocity transitions. In Section 4, we present some general considerations about Eulerian velocity statistics. These are then employed to relate the effective shear rate to flow characteristics in Section 5. In Section 6, we develop predictions for asymptotic longitudinal dispersion in the presence of both advective and diffusive transitions. Overall conclusions are presented in Section 7, and some supporting technical derivations and numerical validation results may be found in the appendices. As Lagrangian particles are transported, they undergo advective changes in velocity due to variability along streamlines as well as diffusion-induced averaging over, and transitions to, nearby streamlines. The transverse variability of the flow field is encoded in a velocity-dependent effective shear rate α e (v), which is determined by statistical properties of the flow field. This conceptual illustration is two-dimensional and includes the presence of a solid phase (black circles), but the proposed approach is applicable also to three dimensions and arbitrary flow fields characterized by well-defined characteristic lengths along the longitudinal and transverse directions, as developed in detail in the main text.

Transport as a spatial-Markov process

We present, as a starting point, a generic formulation of transport as a spatial-Markov process (Le [START_REF] Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF]. We start from a discrete formulation and then proceed to consider its continuum limit. Finally, we provide general forms for the dynamical equations of some key transport quantities. This general formulation will be adapted to describe the role of advective and diffusive transitions in the sections that follow.

Discrete formulation

The velocity magnitudes V k after k spatial steps of fixed length ∆s along streamlines are assumed to form a Markov chain. This spatial-Markov velocity process is characterized by its transition probabilities. Discretizing velocity magnitudes into classes, the transition probabilities r ij (s) describe the probability that a tracer particle will be in class i at distance s + ∆s, given that it was in class j at distance s.

We consider transport to be advection-dominated along the local flow direction, with diffusion playing a role locally along the transverse direction(s). At a given velocity v, a spatial step is associated with a duration ∆s/v. The time T k after k spatial steps is thus described by the stochastic recursion relation

T k+1 = T k + ∆s V k , (1)
with the initial time T 0 = 0. The distribution of initial velocities V 0 is determined by the initial spatial distribution of scalar, which is mapped onto the initial distribution of velocities according to the Eulerian velocity field. For example, an homogeneous spatial distribution corresponds to velocities distributed according to the Eulerian velocity probability density function (PDF), which will be discussed in detail in Section 4. In what follows, we will characterize the Markov chain describing the evolution of velocity along a particle's trajectory through its transition probabilities, which depend on the statistical properties of the underlying flow field. The evolution of particle longitudinal positions (along the mean flow direction) is then modeled as

X k+1 = X k + k∆s χ , ( 2
)
where χ is the average tortuosity, which can be computed as the spatial average of the magnitude of Eulerian velocities divided by the spatial average of their projection along the mean flow direction [START_REF] Koponen | Tortuous flow in porous media[END_REF]Puyguiraud et al., 2019b). The use of the average tortuosity to relate longitudinal displacements to displacements along streamlines is common in spatial-Markov models, and it has been applied successfully to describe transport at the pore scale (Puyguiraud et al., 2019a,b). We will assume this approximation in the following developments. For completeness and clarity, we first briefly discuss some of its limitations and possible extensions. An important limitation is the inability to account for recirculation zones, which, if present can be reached by diffusion and lead to long retention times before solute can again exit by diffusion. Such effects would be mostly naturally included in the present framework by introducing transition probabilities into an additional zero-velocity state, with distributed retention times with statistics determined by diffusion and the geometry of recirculation zones, in the spirit of mobileimmobile or multi-rate mass transfer models (see e.g. [START_REF] Coats | Dead-end pore volume and dispersion in porous media[END_REF][START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity[END_REF][START_REF] Comolli | Non-Fickian transport under heterogeneous advection and mobile-immobile mass transfer[END_REF]. Such an approach could in addition be used to account for retention times due to sorption and desorption and related effects. Similarly, the averagetortuosity approximation does not explicitly account for local flow reversal, and it is in general not expected to be adequate if a strong variability in tortousity, rather than in velocity, is responsible for the main effects of transport variability across streamlines. At the expense of greater model complexity and difficulty of parameterization, the present approach could in principle be used in conjunction with variable tortuosity, in terms of statistics or mean values conditioned on longitudinal position and/or velocity. Flow reversal, if important, would require allowing for negative displacements in longitudinal position to be associated with steps along streamlines.

An important simplification brought about by the average-tortuosity approximation relates to the description of quantities at fixed longitudinal distance from injection. An important example, which we will discuss below, are breakthrough curves, representing mass flux per unit time at fixed control planes transverse to the mean flow direction. Under this approximation, distributions on a control plane at distance x from injection coincide with distributions at fixed distance s = χx measured along streamlines. Therefore, fixed-s statistics, which are most naturally obtained in a spatial-Markov model, directly determine control-plane statistics. In general, the relationship between fixed-s and fixed-x distributions is more complex and depends also on the statistics of tortuosity.

Continuum limit

As we will see, it is convenient in our formulation to define velocity classes related to diffusive averaging in terms of the discretization ∆s. Using this approach, the continuum limit of ∆s → 0 will also correspond to infinitesimal velocity class sizes and transition times ∆s/V k . Therefore, it will be associated with a continuous stochastic process for the random time needed to travel distance s along times. In this sense, the recursion relation (1) in the limit ∆s → 0 defines a stochastic process

T (s) = s 0 ds V S (s ) , (3)
where V S is the Markov process corresponding to the continuum limit of the Markov chain defined by the transition probabilities introduced above. The change in a quantity q i (s), depending on velocity class i and given distance s, due to all possible velocity transitions over a step ∆s is given by q i (s + ∆s)q i (s) = j =i r ij q j , where r ij is the transition probability from class j to class i and the sum extends over all velocity classes j = i. Let q(v; s) be the associated continuous density, i.e., q i (s) = bi+1 bi dv q(v; s) ≈ ∆v i q(v i ; s), where the velocity class i is defined as comprising

velocities v ∈ [b i , b i+1 [ and ∆v i = b i+1 -b i . Then, in the limit ∆s → 0, ∂q(v; s) ∂s = Lq(v; s), (4)
where L is the continuum operator describing velocity transitions, defined through

Lq(v i ; s) = lim ∆s→0 r ij (s) -δ ij ∆s∆v i . (5)
We are now in a position to develop general forms for the dynamical equations governing transport quantities. The actual dynamics will then depend on the particular form of transition operator L, which embodies the transition probabilities of the spatial-Markov velocity process.

The space-Lagrangian velocity PDF describes the velocity point statistics of Lagrangian particles at fixed distance traveled along streamlines. Note that, in the presence of diffusion, the same particle trajectory spans multiple streamlines. In addition to advective transport along streamlines, diffusion transverse to the local flow direction induces transitions to nearby streamlines, as will be formalized in what follows. As discussed above, under the average-tortuosity approximation adopted here, this presents no further difficulties, as transport in the average flow direction can be directly related to the total distance traveled along (a collection of) streamlines. Using Eq. ( 4), we obtain a master equation for its evolution in space,

∂p S (v; s) ∂s = Lp S (v; s), (6)
with no-flux boundary conditions for velocity in order to conserve probability. For a given initial condition in s, corresponding to the velocity distribution at s = 0, the solution of Eq. ( 6) represents the distribution of velocities at fixed s over an ensemble of trajectories, irrespective of the arrival time. Equation ( 6), together with the initial condition, defines the continuous stochastic velocity process V S .

Other important transport quantities are breakthrough curves (mass flux per unit time at a given distance) and concentration profiles (PDF of positions at a given time). Due to correlations introduced by the velocity process, it is convenient to first consider the joint probability density of velocity and arrival time at fixed distance, defined such that ψ(v, t; s) dvdt is the probability of arriving at a given distance s at a time in [t, t + dt[ and with velocity in [v, v + dv

The boundary and initial conditions are no-flux in velocity as before, along with ψ(v, t; 0) = p S (v; 0)δ(t) and ψ(v, 0; s) = 0. Note that, by definition, we have ∞ 0 dt ψ(v, t; s) = p S (v; s). Indeed, integrating out t in Eq. ( 7) leads to Eq. ( 6) for the space-Lagrangian velocity PDF. The first passage time PDF at distance s is given by φ(t; s) = ∞ 0 dv ψ(v, t; s). Particle positions as a function of time are given by X T (t) = X[S(t)] = S(t)/χ, where S(t) describes the random distance traveled by advection along streamlines. Since S(t) always increases with time, and longitudinal position is approximated using the average tortuosity, each Lagrangian particle in this model crosses a given longitudinal position at most once. The breakthrough curves, normalized to unit total mass, are thus related to the first passage times by f (t; x) = φ(t; χx).

In order to obtain the PDF of particle positions at a given time, we employ the concept of subordination [START_REF] Feller | An introduction to probability theory and its applications[END_REF][START_REF] Meerschaert | Stochastic models for fractional calculus[END_REF], which may be thought of as a stochastic change of independent variable. For example, consider X U (u) describing the random position of a particle as a function of time u spent moving. If the time spent moving as a function of total time t is itself a random variable U (t), the position of the particle as a function of total time is given by the subordinated process X T (t) = X U [U (t)]. The total time T (u) given time u spent moving is called the subordinator, and U (t) is called its conjugate process. Here, the time T (s) as a function of fixed distance s plays the role of the subordinator, and S(t) is its conjugate process at fixed time t. According to the theory of subordination, the PDF of S(t) is then given by [START_REF] Feller | An introduction to probability theory and its applications[END_REF][START_REF] Benson | A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations[END_REF][START_REF] Meerschaert | Stochastic models for fractional calculus[END_REF]

h(s; t) = t 0 dt ∂f (t ; s) ∂s . (8)
Integrating out velocity in Eq. ( 7), and noting that, as reflected by the boundary conditions, conservation of probability leads to the integral of the right hand side vanishing, this may be written as

h(s; t) = ∞ 0 dv v -1 ψ(v, t; s). (9)
By definition, h(s; t) ds is the probability of having traveled a distance in [s, s + ds[ along streamlines by time t. The spatial concentration, normalized to unit mass, is the PDF of particle positions X T (t) = S(t)/χ, and it is therefore given by c(x; t) = χh(χx; t). Finally, consider the time-Lagrangian velocity PDF, describing particle velocities at fixed time rather than distance, that is, the PDF of velocities V T (t) = V S [S(t)] at fixed t. Using the same approach as before leads to

p T (v; t) = ∞ 0 ds v -1 ψ(v, t; s). ( 10
)
We note that equations for multi-point PDFs, corresponding to the joint PDFs of quantities at more than one time or distance, may be obtained through similar procedures, see [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF].

The diffusing-velocity random walk

As outlined in the introduction, velocity transitions may occur due to both transverse diffusion and velocity variability along each streamline. We now proceed to quantify the impact of these two processes within the generic spatial-Markov framework outlined in the previous section. To this end, we first present an adapted derivation of the formulation of advective velocity transitions developed by [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]. We then develop a new approach to quantify the impact of transverse diffusion on velocity transitions. Finally, we combine the two to arrive at the general DVRW framework.

Advective transitions

We consider first advective transitions along streamlines, in the absence of diffusion. Let r A ij be the velocity transition probabilities of the spatial-Markov process associated with advective transitions, i.e., due to changes in velocity along a streamline. In order to characterize the continuum limit ∆s → 0, we make use of the fact that the velocity process is Markov, with some correlation length / / corresponding to the longitudinal correlation length of the velocity field. We write, to first order in ∆s,

r A ij = ∆s / / β ij (1 -δ ij ) + 1 - ∆s / / (1 -β ii ) δ ij , (11)
where δ •• is the Kronecker delta. The first term corresponds to the probability of changing velocity class and the second of staying in the same class. The β ij encode the dependency of the transition probabilities on the velocity classes. To ensure normalization, i.e., i r A ij = 1, we must have

β jj = 1 -i =j β ij .
Recall that r ij -δ ij , seen as a function of velocity class j for each velocity class i, defines an operator describing the change in velocity class due to a transition, see Eq. (5). Equation ( 11) leads, to first order in ∆s, to

r A ij -δ ij = ∆s / / (β ij -δ ij ). ( 12
)
Thus, in the continuum limit, we recover the dynamical equations of Section (2) with a transition operator L = L A describing advective transitions along streamlines. This is in general an integral (as opposed to differential) operator, because the velocity transitions may be long-range in velocity space. Indeed, according to Eq. ( 5), we have, for an arbitrary density q(v) as a function of velocity v,

L A q(v) = -1 / / ∞ 0 dv [β(v|v ) -δ(v -v )] q(v ), ( 13
)
where β is the transition PDF along streamlines (units [β] = 1/V = T /L), i.e., β(v|v ) dv is the transition probability from velocity v to a velocity in [v, v + dv[. It corresponds to the limit

β(v i |v j ) = lim ∆s→0 β ij ∆v i . ( 14
)
In this case, it has been shown that, under ergodicity and incompressibility, the equilibrium space-Lagrangian velocity PDF is the flux-weighted Eulerian PDF (see Section 4), irrespective of the initial condition, and the equilibrium time-Lagrangian velocity PDF is the Eulerian velocity PDF [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]Puyguiraud et al., 2019a,b). Non-stationary transition probabilities can be encoded in sdependent β ij and would lead to s-dependent β(v|v ). We note that, in natural media such as geological structures, non-stationarity is typically most naturally described in terms of distance x along the mean flow direction, which under the average-tortuosity approximation is directly related to s via s = χx.

The details of the dynamics depend on the choice of process governing the transition probabilities between velocity classes and on the underlying Eulerian velocity distribution. As in [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF], we will focus on the case of Bernoulli relaxation. This process is defined by the transition probabilities

r A ij = e -∆s/ / / δ ij + 1 -e -∆s/ / / p ∞ i . (15)
The physical interpretation of this setup is as follows. Velocities persist on the scale of the correlation length / / , and when a transition occurs, the probability of the new velocity being in class i is given by the prescribed equilibrium probability p ∞ i . Thus, the velocity distribution "relaxes" towards the equilibrium distribution on a spatial scale on the order of / / . The exponential probability of persistence in Eq. ( 15) is a consequence of the assumption that the probability of transition per unit length is constant and equal to 1/ / / . In this sense, the Bernoulli process may be seen as the simplest Markov process converging to a given equilibrium distribution on a given scale. For small ∆s, this corresponds, according to Eq. ( 11), to β ij = p ∞ i , the equilibrium probability of velocity class i. The continuum-limit transition PDF is given by β(v i |v j ) = p ∞ (v i ), the equilibrium PDF evaluated at v i , irrespective of v j . As mentioned above, the equilibrium PDF should in this context be taken equal to the flux-weighted Eulerian PDF.

Under the Bernoulli relaxation process, all velocities relax towards the equilibrium distribution at the same spatial rate. Other processes may be used to account for velocity-dependent relaxation while retaining the Markov property and stationarity of the transitions. For example, a process based on the classical Ornstein-Uhlenbeck process describing temporal velocity fluctuations in Brownian motion [START_REF] Uhlenbeck | On the theory of the brownian motion[END_REF] has been successfully employed to capture slower spatial relaxation of low velocities in connection with transport in a heterogeneous porous medium (Puyguiraud et al., 2019a). Conversely, in the presence of preferential high-velocity channels, one expects stronger correlation, and consequently slower spatial relaxation, at high velocities. In general, in the present formulation, velocity-dependent relaxation corresponds to j-dependent β ij . From Eq. ( 12), where the latter always occur in the combination β ij / / / , we see that this leads, in effect, to a velocity-dependent correlation length. In this sense, we may interpret / / / i =j β ij = / / /(1β jj ) as an effective correlation length associated with velocity class j. The probability of leaving a velocity class in a given step is then inversely proportional to the effective correlation length. Note that the β ij are in principle arbitrary, so long as probability is conserved, i β ij = 1, and they lead to the correct equilibrium distribution, the flux-weighted Eulerian PDF. Even when the overall velocity distribution across all particles has reached equilibrium, single-particle velocities remain dynamic, with the spatial velocity series of each particle undergoing purely-advective transport being directly determined by the choice of β ij .

Diffusive transitions

Next, we examine the role of diffusion, disregarding advective transitions for the moment. This corresponds to the limit of infinite correlation length / / → ∞. It is directly applicable to stratified flow, where velocity is constant along each streamline, and velocity transitions are due only to transverse diffusion across streamlines. This will lead us to discrete and continuous formulations of transverse diffusion as a spatial-Markov process.

Discrete formulation

The characteristic transverse length explored by diffusing particles in a time interval ∆t is given by √ 2D∆t, and a spatial displacement of length ∆s at velocity v corresponds to a duration ∆t = ∆s/v. The local averaging of velocities due to transverse diffusion during this time interval depends on the spatial structure of the velocity field. The cornerstone of our approach is the notion of a velocitydependent effective shear, see Fig. 1. We describe the local variation of the flow around regions of velocity v in terms of an effective transverse shear rate magnitude α e (v), so that the range of velocities averaged by diffusion around velocity v is given by ∆v

(v) = α e (v) 2D∆s/v. ( 16
)
The actual transverse velocity gradient magnitude depends in general on position, and a given value of velocity at a randomly chosen spatial location is thus associated with a PDF of possible shear values.

Our approach may be seen as a mean-field formulation associating an average shear rate α e (v) with each velocity value v.

We consider for now that the effective shear rate α e (v) is known as a function of velocity magnitude, and we will derive the consequences for transport. Section 5 will be devoted to relating the effective shear rate to physical properties, in particular the underlying flow statistics. Note that, analogously to above with the transition PDF β characterizing advective transitions, we have assumed that α e (v) does not depend on distance s. A non-stationary model can be formulated using the present approach, but we refrain from exploring it here for simplicity.

The changes in velocity at each step are modeled as a spatial-Markov process, which is characterized by the diffusive transition probabilities r D ij from each velocity class j to each velocity class i over a step ∆s along the flow direction. We discretize velocity magnitudes

v into classes v ∈ [b i , b i+1 [. The width ∆v i = b i+1 -b i of
each class is determined by transverse diffusive averaging according to Eq. ( 16). That is, we set ∆v i = α i 2D∆s/v i , where v i is the average velocity within class i and α i = α e (v i ). A recursive construction valid in the limit of small class widths is given in Appendix B. The class widths ∆v i vanish in the limit ∆s → 0, so that small-∆v i approximations are reasonable for small ∆s. The same is true of the transition times ∆s/v i .

In each step ∆s, diffusion averages over the current velocity class, and induces transitions to the nearest classes according to the transition probabilities

r D ij = r + j δ i,j+1 + r - j δ i,j-1 , (17)
where r ± j are the transition probabilities from class j to class j ± 1, with r + j + r - j = 1. For the j = 0 velocity class, we have r + 0 = 1 and r - 0 = 0, since there is no class below. Conversely, for the highest velocity class, transitions are always to the class below. For the remaining classes, the diffusive transition probabilities are obtained as follows. By construction, during a transition, diffusion homogenizes a transverse length corresponding to a velocity class. The transition is thus associated with a transition time ∆t = ∆s/v j , where v j is the (arithmetic) average velocity in the class. However, for a given ∆t, the amount of distance traveled at a given velocity v is proportional to v. This means that the probability of a particle having velocity v when it finishes the spatial step ∆s is proportional to v. In other words, the velocity distribution within a class after a spatial transition is flux weighted. Imposing a diffusive transition to the velocity class below if the particle has a velocity lower than the class average (and to the class above otherwise) leads to r

- j = vj bj dv v/ bj+1 bj dv v . Thus, r - j = (v 2 j -b 2 j )/(b 2 j+1 -b 2 j ).
Approximating the class average v j by the class center, we have b

j = v j -∆v j /2 and b 2 j+1 -b 2 j = (b j+1 -b j )(b j+1 + b j ) = 2v j ∆v j , so that r ± j = 1 2 ± ∆v j 8v j . ( 18
)

Continuum limit

According to the previous construction, as ∆s → 0, both the velocity class widths ∆v i → 0 and the corresponding transition times ∆s/v i → 0, indicating that this corresponds to a genuine continuous limit. We now examine this limit in detail, and we obtain the continuum stochastic process underlying diffusive transitions, as well as the associated transition operator.

In the continuum limit, Eq. ( 17) leads to the operator L = L D associated with diffusive transitions, see Eq. ( 5). Because diffusive transitions are local, the corresponding operator is differential. For an arbitrary density q,

L D q(v) = ∂ ∂v γ(v) ∂q(v) ∂v -µ(v)q(v) , (19)
see Appendix C for details on the derivation. The spatial velocity diffusivity γ

([γ] = V 2 /L = L/T 2 ) corresponds to the limit γ(v i ) = lim ∆s→0 ∆v 2 i /(2∆s)
and is given by

γ(v) = Dv -1 α e (v) 2 , (20)
and 21) is a spatial velocity drift ([µ] = V /L = 1/T ). The first term in this drift is due to the fact that particles are more likely to transition to higher velocities by diffusion within a class, due to the flux-weighting effect discussed above. The second arises because velocity classes, corresponding to spatial averaging by diffusion over a constant spatial step, have velocity-dependent sizes. Thus, particles spend longer distances at velocities where classes are smaller, giving rise to an effective drift towards these velocities.

µ(v) = v -1 - ∂ ∂v γ(v) 2 (
The size of the velocity class decreases with velocity, because shorter crossing times are associated with smaller diffusion lengths, and increases with effective shear, because higher effective shear corresponds to stronger variation of velocity over the same transverse distance (see Fig. 1). In particular, the space-Lagrangian velocity PDF obeys the master equation ( 6) with L = L D . The boundary conditions in v must ensure conservation of probability, so that we have the no-flux condition γ∂p S /∂vµp S = 0 at the minimum and maximum velocities.

Combining advective and diffusive transitions

We now combine the diffusive and advective transition mechanisms to obtain the complete transition probabilities of the DVRW framework, see Fig. 1. We impose that, if (and only if) a particle does not undergo an advective transition along a streamline, diffusion causes a transition to one of the nearest velocity classes. That is, the transition probabilities in the presence of both advection and diffusion become

r ij = r A ij (1 -δ ij ) + r A jj r + j δ i,j+1 + r - j δ i,j-1 . (22)
Thus, to first order in ∆s, the changes in velocity class are determined by

r ij -δ ij = ∆s / / (β ij -δ ij ) + r + j δ i,j+1 -δ ij + r - j δ i,j-1 . (23)
Due to the Markovian nature of the process, the somewhat artificial requirement that diffusive transitions occur only when an advective transition does not occur is inconsequential in the limit of small ∆s. This can be seen from the fact that r ijδ ij , which defines the operator characterizing the change in velocity classes over ∆s, is composed of the sum of terms associated with purely-advective and purely-diffusive velocity transitions to leading order in ∆s; thus, the transition operator in the continuum limit becomes the sum of the diffusive and advective contributions. This leads to the same continuum-limit dynamical equations as before, with the transition operator now given by L = L A +L D , see Eqs. ( 13) and ( 19), representing the effect of both advective and diffusive transitions. That is, for an arbitrary density q,

Lq(v) = -1 / / ∞ 0 dv [β(v|v ) -δ(v -v )] q(v ) + ∂ ∂v γ(v) ∂q(v) ∂v -µ(v)q(v) . ( 24
)
Note that, if there is no velocity variation along streamlines, or equivalently / / → ∞, then r A ij = δ ij , and we recover the pure-diffusion formulation, valid for stratified flow. Conversely, D = 0 recovers the purely-advective scenario discussed above.

Eulerian velocity statistics

As seen from Eqs. ( 19)-( 21), the effective shear α e (v) plays a key role in quantifying diffusive transitions.

In order to relate it to flow properties, and in particular to velocity statistics, let us first discuss some properties of the Eulerian PDF of velocity magnitudes, defined as the probability of finding a certain velocity magnitude value at a uniformly random spatial location.

Denoting the spatial velocity field magnitude at position x in a domain Ω by v E (x), the Eulerian velocity PDF is then defined as

p E (v) = |Ω| -1 Ω dx δ[v -v E (x)], (25)
where δ(•) is the Dirac delta, and for a d-dimensional spatial domain A, |A| denotes its measure (number of elements, area, or volume, respectively for d = 1, 2, 3). Assuming a smooth, non-constant velocity field, changing variables in the Dirac delta [START_REF] Hörmander | The analysis of linear partial differential operators I: Distribution theory and Fourier analysis[END_REF] leads to

p E (v) = Λ(v) dσ(x) |Ω||∇v E (x)| . (26)
In d spatial dimensions, Λ(v) is the (d -1)-dimensional spatial surface where the velocity field has magnitude v, Λ(v) = {x ∈ Ω : v E (x) = v}, and dσ(x) is the corresponding (d -1)-area element at the point x on Λ(v). The harmonic average of the local shear rate magnitude |∇v E (x)| over this surface is given by

α h (v) =   |Λ(v)| -1 Λ(v) dσ(x) |∇v E (x)|    -1 , ( 27
) leading to p E (v) = |Λ(v)| |Ω|α h (v) , (28)
which shows that the Eulerian PDF is directly related to the harmonic average of the local shear rate given a velocity magnitude. Note that this result is valid for an arbitrary smooth velocity field that is not constant. For a piecewise-smooth velocity field, the result applies piecewise. If the velocity field has fully degenerate maxima or minima, that is, regions of finite volume (in two dimensions, area)

where velocity is constant, they yield Dirac-delta contributions corresponding to that velocity, with a probability mass (coefficient) given by the ratio of the region volume and |Ω|. This can easily be seen from the definition of the Eulerian PDF, Eq. ( 25).

Throughout, an overline will denote the ensemble average (over tracer particles), and V E will stand for a random variable distributed according to p E , corresponding to the velocity at a uniformly random spatial point. Using the previous result, it follows that

α h (V E ) = |Ω| -1 ∞ 0 dv |Λ(v)|. ( 29
)
We introduce also the flux-weighted Eulerian PDF, defined according to

p F (v) = v V E p E (v), (30)
which plays an important role in our formulation.

We now consider the Eulerian PDF of velocities at a fixed distance s along streamlines. Let Ω ⊥ (s) be the (d -1)-dimensional cross-section of Ω at fixed s, and let Λ(v; s) = {x ∈ Ω ⊥ (s) : v E (x) = v} be the (d -2)-surface of constant velocity on Ω ⊥ (s). The gradient of the velocity magnitude transverse to the flow direction is given by

∇ ⊥ v E , where ∇ ⊥ = ∇ -(v/|v| 2 )v • ∇. Adapting the previous derivation, p E (v; s) = |Λ(v; s)| |Ω ⊥ (s)|α h (v; s) , (31)
where α h (v; s) is the harmonic average of |∇ ⊥ v E | over Λ(v; s). Note that the Eulerian velocity PDF at fixed s is not sensitive to gradients along the flow direction; their contribution to the full PDF arises through their effect on the variation of |Λ(v; s)|/α h (v; s) with s. For stratified flows, where velocity is constant along each streamline, p E (v; s) = p E (v) always holds. This equality also holds for more general flows, as long as the point statistics of velocity over any given transverse plane coincide with those of the full domain. We will assume this to be the case in what follows for simplicity.

As an example, consider Poiseuille flow in d = 2 dimensions. The Eulerian velocity field is in this case given by

v E (y) = v M 1 - 2y L 2 , ( 32
)
where L = |Ω ⊥ | is the transverse domain width, y ∈ [-L/2, L/2] is the position in the transverse direction, and v M is the maximum velocity, occurring at y = 0. The full Eulerian PDF is equal to the PDF over the transverse direction, since there is no variability along the longitudinal direction. The absolute value of the gradient of velocity is uniquely determined by the velocity, and we have

α h (v) = 4v M L 1 - v v M . ( 33
)
Since the same absolute value of the gradient occurs at exactly two points (except at the maximum, where it is zero), we have |Λ(v)| = 2. According to Eq. ( 28), the Eulerian PDF is thus given by

p E (v) = 1 2v M 1 -v/v M , (34)
and the average velocity is V E = 2v M /3. Note the square-root divergence near the maximum, which is in agreement with the discussion in Appendix D for d = 1, the effective dimension of variability of this flow field.

Effective shear rate

This section is devoted to linking the effective shear α e (v) rate to flow characteristics, in particular point statistics as encoded in the Eulerian PDF. In what follows, we will achieve this by requiring that the DVRW formulation satisfy two criteria in the limit of longitudinal correlation length / / → ∞ (i.e. when only diffusive transitions are present): (i) reproducing the asymptotic space-Lagrangian velocity PDF as distance s → ∞; and (ii) reproducing the asymptotic Taylor dispersion coefficient as time t → ∞.

Equilibrium velocity PDF

As shown in Section 3.2, the space-Lagrangian velocity PDF for purely-diffusive velocity transitions obeys the master equation ∂p S (v; s)/∂s = L D p S (v; s), with the diffusive transition operator L D given by Eq. ( 19). Thus, the equilibrium PDF solves L D p ∞ S (v) = 0. Integrating this equation using no-flux boundary conditions, and imposing normalization, we obtain

p ∞ S (v) = v α e (v)   ∞ 0 dv v α e (v )   -1 . ( 35
)
According to Taylor dispersion theory, the equilibrium distribution of velocities for a stratified flow after diffusion samples the full transverse variability is the flux-weighted Eulerian distribution, p ∞ S (v) = p F (v). When only velocity transitions by transverse diffusion are considered, which is equivalent to taking the limit of an infinite longitudinal correlation length / / in Eq. ( 24) for the transition operator L, the DVRW becomes equivalent to transport in stratified flow. Asymptotic dispersion should then agree with the Taylor result.

According to Eq. ( 35),

p ∞ S (v) ∝ v/α e (v)
, where the proportionality factor is v-independent and ensures normalization. Obtaining the flux-weighted Eulerian PDF, Eq. ( 30), as the space-Lagrangian equilibrium PDF thus requires α e (v) ∝ 1/p E (v), so that p ∞ S (v) ∝ vp E (v). Therefore, we write for the effective shear rate

α e (v) = α e (V E ) δv p E (v) , (36)
for v ∈]v m , v M [ (and zero otherwise) and with δv = v Mv m the difference between the maximum and minimum velocities. Note that we have fixed the velocity-independent coefficient in terms of the average effective shear rate α e (V E ) = ∞ 0 dv p E (v)α e (v). We will show shortly that this average is determined by asymptotic longitudinal dispersion. It is worth noting also that a physical instance of a velocity field in a finite domain always exhibits a finite maximum velocity. However, theoretical Eulerian PDFs, applying in principle to an infinite domain or an infinite number of domain realizations, may be defined for any positive velocity magnitude. We will later obtain a form of the effective shear rate which does not explicitly depend on the maximum and minimum velocity values and is suitable for direct computation for an arbitrary Eulerian PDF.

In order to provide intuition for the velocity-dependence of the effective shear rate, consider the stratified flow profile shown in Fig. 2, corresponding to a cut in the direction transverse to a twodimensional flow, which we name the M-flow. This synthetic example, where the local shear magnitude

Transverse position y

Eulerian velocity magnitude v E (y)

v m v M δv M δv m Λ M = 4 Λ m = 2 v c ℓ m /2 ℓ M ℓ m /2

Longitudinal dispersion

Next, we turn to the asymptotic behavior of longitudinal dispersion for purely-diffusive velocity transitions, that is, in the limit of infinite longitudinal correlation length / / as before. First, consider V T (t), the Lagrangian velocity process along particle trajectories as a function of travel time t. We have, for particle positions as a function of time,

X T (t) = χ -1 t 0 dt V T (t ), (37)
from which X T (t) = χ -1 t 0 dt V T (t ). The previous results imply that the equilibrium time-Lagrangian velocity PDF coincides with the Eulerian PDF, as expected from Taylor dispersion theory. To see this, consider the joint PDF of velocity and arrival time, Eq. ( 7) with L = L D . Integrating out s and using Eq. ( 10) for the time-Lagrangian velocity PDF leads to

∂p T (v; t) ∂t = p S (v; 0)δ(t) + L D ∞ 0 ds ψ(v, t; s). (38)
For the steady-state, we must thus have ∞ 0 ds ψ(v, t; s) ∝ p F (v), and therefore, using normalization,

p ∞ T (v) = p E (v). ( 39
)
Since, by definition, the time-Lagrangian PDF is the PDF of V T (t), this implies that V T (t) converges to V E , and therefore, asymptotically,

X T (t) = χ -1 V E t. Similarly, calculating X 2 T (t) leads to a longitudinal dispersion σ 2 (t) = χ -2 t 0 dt t 0 dt V T (t )V T (t ), ( 40
)
where for late times V T (t) = V T (t) -V E are the velocity fluctuations about the mean as a function of particle travel time. In a statistically stationary velocity field, the velocity correlations at late times depend only on the time difference,

V T (t )V T (t ) = C v (|t -t |)
. This yields the Green-Kubo relation [START_REF] Kubo | Statistical Physics II: Nonequilibrium Statistical Mechanics[END_REF] for the longitudinal dispersion coefficient,

D / / (t) = 1 2 dσ 2 (t) dt = χ -2 t 0 dt C v (t ). ( 41
)
The integral of the correlation function of velocity fluctuations, when it converges, is given by the product of the velocity autocorrelation and a correlation time. The first is proportional to V E 2 and, for purely-diffusive velocity transitions, the second is of the order of the diffusion time

τ D = L 2 /(2D),
where L is the characteristic width of the domain cross-section. For large t, the integral can be approximated by extending the upper limit to infinity, leading to the asymptotic Taylor dispersion coefficient

D T = η V E 2 L 2 χ 2 D , ( 42
)
where η is a dimensionless coefficient characterizing the impact of the spatial organization of velocities. This is the form of the asymptotic longitudinal dispersion coefficient whenever the only mechanism for sampling velocity variability is diffusive (for example, in stratified flows), and it arises once the full variability has been sampled. While η does not have a simple general form, it can be expressed in terms of the normalized velocity fluctuations within the transverse domain (see [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF] and Appendix E),

ν(y) = v E (y) V E -1, y ∈ Ω ⊥ . (43)
In two dimensions, it is given by

η = L -3 L/2 -L/2 dy   y 0 dy ν(y )   2 , (44)
whereas for an axisymmetric flow we have, with ν(r) in terms of the radial coordinate r = √ y • y,

η = 8 L 4 L/2 0 dr r   r 0 dr r ν(r )   2 . ( 45
)
The longitudinal dispersion coefficient corresponding to the DVRW for purely-diffusive velocity transitions can be computed from the dynamical equation ( 7) for the joint PDF of velocity and arrival time, with the transition operator L = L D given by Eq. ( 19), together with Eq. ( 9) for the concentration PDF. As shown in Appendix F.1, the asymptotic longitudinal dispersion coefficient is given by

D ∞ = Iδv 2 α e (V E ) 2 L 2 V E 2 L 2 χ 2 D , (46)
where the dimensionless coefficient I is a property of the Eulerian velocity PDF,

I = ∞ 0 dv p E (v)[C E (v) -C F (v)] 2 . ( 47
)
Here, we have introduced C E,F (v ) = v 0 dv p E,F (v ) as the Eulerian and flux-weighted Eulerian cumulative distribution functions (CDFs), respectively.

We require that, when only diffusive transitions are present, the DVRW should reproduce the correct Taylor dispersion coefficient, D ∞ = D T . The first fraction in Eq. ( 46) must then equal η, see Table 1: Values of ω, I, and η for different stratified flows. The last column shows the value I/(ηω 2 ) of the ratio of the average effective shear to the harmonic average shear. Poiseuille flow for d = 3 dimensions refers to a paraboloid profile in a cylindrical pipe. Coefficients for the M -flow are given in terms of λ = |Λ(V E )|; the polynomial coefficients for η and I are given by c i = (2368, -5760, 5840, -2880, 740, -96, 5) and d i = (3904, -9600, 10160, -4800, 1100, -120, 5), i = 0, . . . , 6, and

K(λ) = [(6 -λ) 2 -2(4 -λ) 2 ] 2 . Flow ω I η α e (V E )/ α h (V E ) Couette (d = 2) 1 1/30 1/30 = 1 Triangular (d = 2) 2 1/30 1/120 = 1 M -flow (d = 2) 8 6 -λ 6 i=0 c i λ i 120K(λ) 6 i=0 d i λ i 1920K(λ) ≈ 7.74 • 10 -1 -1.02 Poiseuille (d = 2) 2 2/105 1/210 = 1 Poiseuille (d = 3) 8/3 1/30 1/192 ≈ 9.49 • 10 -1
Eq. ( 42). This represents a second constraint on the average effective shear rate, which must be given by

α e (V E ) = I η δv L . (48)
Together with Eq. ( 36), which enforced the first constraint of reproducing the correct asymptotic space-Lagrangian velocity PDF, we obtain

α e (v) = I/η Lp E (v) . ( 49
)

Role of statistical and spatial flow structure

From Eqs. ( 29) and ( 48), we may write

α e (V E ) = I ηω 2 α h (V E ), ω = L δv ∞ 0 dv |Λ(v)|. (50)
It may be surprising to observe that the ratio between the average effective shear and the average harmonic shear may differ from unity. The reason for this lies in the fact that the average harmonic shear rate, as defined by Eq. ( 27), does not provide information about the spatial organization of the velocity field. As we have seen, the velocity dependence of the effective shear rate is fixed by the Eulerian PDF. Different spatial organizations, however, lead to different Taylor dispersion coefficients, as encoded in η, and this information must thus be included in the effective shear rate through its average value. The impact on dispersion due solely to the point statistics of velocity is quantified by the coefficient I. Using the definition of the Eulerian PDF, Eq. ( 25), we can express I as a spatial integral,

I = |Ω ⊥ | -3 Ω ⊥ dy   Ω ⊥ dy ν(y )H[v E (y) -v E (y )]   2 , ( 51
)
where H(•) is the Heaviside step function, meaning the inner integral extends over positions y at which the velocity magnitude is smaller than at y. Details on the derivation may be found in Appendix F.1. The relationship between I and η, i.e., between the impact of the statistical and spatial structures, becomes apparent when one considers certain special cases. Consider the following particular scenario in d = 2 dimensions. First, assume constant |Λ(v)| = |Λ(V E )|. This case, for which we have also ω = |Λ(V E )|, describes a situation where every velocity value appears ω times. Further assume that the velocity gradient magnitude is the same at each of these values, so that the shear rate magnitude may be written as a function of velocity. This implies the velocity profile has a simple spatial structure: the transverse domain is divided into subregions of equal length L/ω, within each of which the velocity magnitude is monotonic. Furthermore, the velocity profile in adjacent regions is symmetric with respect to reflection across the boundary between them. Examples of this type (in d = 2) are Couette flow, Poiseuille flow, and triangular flow (corresponding to the M -flow, see Fig. 2, with δv M = 0). This spatial structure, together with the fact that the integral of the velocity fluctuations in each subregion is null, leads, according to Eq. ( 51), to

I = L -3 L/2 -L/2 dy   ω y 0 dy ν(y )   2 , ( 52
)
where the factor of ω in the inner integral is due to the equal contribution of each subregion. Comparing to Eq. ( 44), this corresponds to I = ω 2 η, so that, from Eq. ( 50), α e (V E ) = α h (V E ). However, breaking either of the previous hypotheses will in general result in α e (V E ) = α h (V E ). For example, for an axisymmetric flow in d = 3, we have

I = 512 L 6 L/2 0 dr r   r 0 dr r ν(r )   2 , (53)
so that, in general, α e (V E ) = α h (V E ) even in this simple scenario, compare Eq. ( 45). For illustration purposes, we computed η, ω, I, and α e (V E )/ α h (V E ) = I/(ηω 2 ) for a number of simple stratified flows; the results are summarized in Table 1. With these results in mind, we identify

⊥ = η I L (54)
as the relevant lengthscale characterizing the variability of the velocity field in the transverse direction. When I/η = ω, which holds in the simple scenario discussed above, this lengthscale coincides with the spatial period over which the Eulerian statistics of the velocity field repeat. Substituting in Eq. ( 49), the effective shear rate becomes α e (v) = 1

⊥ p E (v) . ( 55
)
This is the central result connecting the effective shear rate, which, together with the usual diffusion coefficient, determines the diffusive transitions of the DVRW, to flow properties. This connection was achieved by enforcing that the DVRW predictions for the asymptotic space-Lagrangian velocity distribution and longitudinal dispersion coefficient reduce to known results for purely-diffusive velocity transitions, that is, for infinite velocity correlation length along streamlines.

The characteristic transverse lengthscale embodies the impact of the interplay between the statistical and spatial structures of the flow field in the transverse direction. In the absence of detailed information, it can be estimated based on the characteristic transverse length over which the full variability of the velocity structure is covered. If the spatial structure of the flow is known on a transverse section of the domain, I may be computed according to the corresponding Eulerian PDF, see Eq. ( 47), and η may be obtained by computing the Taylor dispersion coefficient for a stratified flow characterized by the flow field magnitudes on the cross section, see Eq. ( 42). The assumption of stationarity we have made here, according to which the effective shear rate does not depend on distance s, is satisfied when these quantities remain constant across different domain cross-sections. Although we do not explore the resulting non-stationary model here, the approach generalizes to cross-section-dependent quantities.

We numerically validated the formulation of diffusive transitions in terms of the effective shear rate for transport in some instances of the stratified flows from Table 1 by comparing simulations of the discretized Lagrangian formulation of the DVRW to standard particle tracking simulations with a fixed time step. The results are discussed in Appendix G.

Dispersion under heterogeneous advection and diffusion

In order to illustrate the impact of the interplay between advection and diffusion on transport, let us consider longitudinal dispersion under heterogeneous advection and diffusion. Following [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF], we consider as a rich example a gamma distribution of Eulerian velocities (Fig. 3b). The gamma PDF, expressed in terms of the mean velocity V E , is given by

p E (v) = θv V E θ e -θv/V E vΓ(θ) , ( 56
)
where Γ(•) is the gamma function and θ > 0 is a parameter controlling the power-law dependence ∝ v θ-1 at low velocities. This type of PDF, which combines low-velocity power-law behavior with an exponential cutoff at high velocities, has been employed to model the velocity PDFs in porous medium flows both at the pore and Darcy scales [START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]. The flux-weighted Eulerian PDF is again gamma,

p F (v) = θv V E θ e -θv/V E V E Γ(θ) , (57)
with the same exponential cutoff and a low-velocity dependency ∝ v θ .

Pure advection

We first summarize some known results on dispersion in gamma-distributed velocity fields, in the absence of diffusion [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]. Recall that, at late times, the PDF of velocities associated with crossing a fixed distance is p F , and let V F be a random variable distributed according to the latter. Consider that velocity transitions due to changes along streamlines happen roughly after each correlation length / / . The times / / /V F to cross a correlation length then have PDF p / / (t) = p F ( / / /t) / / /t 2 , and their average is τ / / = / / /V F = / / /V E . The large-time tail of this PDF is associated with small velocities, which, for the gamma PDF, are associated with the behavior p

F (v) ∼ (v/V E ) θ /V E . Therefore, p / / (t) ∼ ( / / /V E ) 1+θ t -2-θ .
The stable exponent corresponding to this tailing behavior is 1 + θ [START_REF] Feller | An introduction to probability theory and its applications[END_REF][START_REF] Meerschaert | Stochastic models for fractional calculus[END_REF]. Thus, the crossing times associated with a gamma distribution of Eulerian velocities always have a finite mean, but their variance is infinite for θ 1 (for a general discussion of crossing times with infinite moments, see e.g. [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]). When θ > 1, the crossing times have both a finite mean and variance, and the relevant correlation time in the Green-Kubo formula (41) for dispersion is of the order of τ / / . This leads to a dispersion coefficient ∼ V E 2 τ / / = V E / / . For θ < 1, on the other hand, the mean of the crossing times is finite, but their variance is infinite. In this case, the variability of waiting times about the mean dominates dispersion. The variance of crossing times observed by time t is σ 2 c (t) ∼ t 0 du u 2 p / / (u) ∼ τ 2 / / (t/τ / / ) 1-θ . The relative importance of the fluctuations about the mean with respect to the mean value is σ 2 c (t)/τ 2 / / , and we estimate the relevant correlation time as τ / / σ 2 c (t)/τ 2 / / ∼ τ / / (t/τ / / ) 1-θ . This leads, according to Eq. ( 41), to a dispersion coefficient

∼ V E 2 τ / / (t/τ / / ) 1-θ = V E 2-θ θ / / t 1-θ .
This semi-heuristic argument can be made precise, leading to an exact form for dispersion for an arbitrary Eulerian PDF under the Bernoulli process, see Appendix F.2. In agreement with the scaling arguments above, for gamma-distributed Eulerian velocities and an initial condition according to the Eulerian velocity distribution, we find the late-time (t τ / / ) longitudinal dispersion coefficient

χ 2 D / / (t) =      (1 -θ) -1 θ θ V E / / t/τ / / 1-θ , 0 < θ < 1 V E / / ln e -1 t/τ / / , θ = 1 (θ -1) -1 V E / / , θ > 1 . ( 58
)
The special case θ = 1, for which the gamma distribution reduces to exponential, leads to a logarithmic correction. Note also that the leading coefficient, in the case of θ < 1, is sensitive to the initial condition, although the scalings with time, mean velocity, and correlation length are not. These results correspond to Fickian diffusion for θ > 1, and superdiffusive behavior, that is, dispersion σ 2 (t) growing superlinearly with t, for θ 1.

Advection-diffusion

We now employ the DVRW formulation to quantify asymptotic longitudinal dispersion when both advective and diffusive transitions are present. Diffusion changes the way in which tracer particles sample the velocity field when compared to pure advection, and it is therefore expected to impact dispersion properties (see Fig. 1). The equality of the equilibrium PDFs for advective and diffusive transitions, which holds under conditions of ergodicity and incompressibility, implies that the flux-weighted Eulerian PDF is also the equilibrium PDF in the presence of both types of transition. However, diffusion impacts the temporal correlation properties, which determine dispersion as discussed in Section 5.2. This leads to a complex dependency of dispersion properties on Péclet number [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Sallès | Taylor dispersion in porous media. Determination of the dispersion tensor[END_REF][START_REF] Deng | Longitudinal dispersion coefficient in straight rivers[END_REF][START_REF] Kandhai | Influence of stagnant zones on transient and asymptotic dispersion in macroscopically homogeneous porous media[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF]Puyguiraud et al., 2019b). Previous attempts to quantify this behavior have introduced a cutoff resulting from longitudinal diffusion in the distribution of local transit times arising from the distribution of velocities [START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF]. The DVRW framework allows for deriving a mechanistic description of the interplay between transverse diffusion and advection. This interplay is mediated by the effective shear rate, and it leads to new dispersion laws. When diffusion is present, averaging takes place across nearby velocities. In other words, a Lagrangian particle that would be retained in a low-velocity area can be removed by transverse diffusion into a faster streamline, effectively cutting off transition times and enforcing a loss of velocity correlation. Even at arbitrarily high Péclet number, this effect is important because it is concerned with velocities that are arbitrarily smaller than the Eulerian mean value. According to the DVRW formulation, transverse diffusion over a correlation length averages over a velocity range ∆v / / (v) = α e (v) 2D / / /v. Diffusion thus enforces a minimum average velocity that can be sampled by a Lagrangian particle over a correlation length according to ∆v / / (v min ) = 2v min , so that v min = [Dα 2 e (v min ) / / /2] 1/3 . This minimum velocity corresponds to a maximum transit time τ max = / / /v min , so that

τ max = 2 2 / /
Dα 2 e (v min )

1/3 . ( 59
)
The maximum correlation time decreases with increasing effective shear rate at low velocities, which corresponds to increased velocity variation over the same spatial distance, as τ max ∼ α e (v min ) -2/3 . Note that this result constitutes an implicit equation for τ max , since v min = / / /τ max .

As we have seen, the persistence of low velocities under purely-advective transport may drive non-Fickian dispersion behavior. As discussed above, when the Eulerian PDF of low velocities scales like a power law, p E (v) ∼ (v/V E ) θ /v for θ < 1, the transit times across a correlation length are broadly-distributed and longitudinal dispersion grows asymptotically like D / / (t) ∼ V E 2-θ θ / / t 1-θ . In the presence of diffusion, after a time t = τ max , larger crossing times are cut off due to diffusive averaging at low velocities and dispersion stabilizes at a value

D ∞ = D / / (τ max ) ∼ V E 2-θ 2+θ 3 / / Dα 2 e (v min ) -1-θ 3 . (60)
Using the relationship between the effective shear rate and the Eulerian PDF, Eq. ( 55), and defining the Péclet number Pe = / / V E /D in terms of the longitudinal correlation length, the previous results yield the scalings

v min ∼ Pe 2θ 1+2θ , τ min ∼ Pe -2θ 1+2θ , α e (v min ) ∼ Pe 3θ 1+2θ . (61)
We have used the power-law dependency of the Eulerian PDF for low velocities compared to the mean at v min , which holds for high Péclet number because v min /V E decreases with Pe. We see that the effective shear rate associated with the minimum velocity increases with Pe. For 0 < θ < 1, this impacts the dispersion coefficient through a factor of α e (v min ) -2(1-θ)/3 ∼ Pe -2θ(1-θ)/(1+2θ) , which decreases with Pe.

As expected, diffusive averaging decreases the dispersion coefficient, because it smooths out velocity variability. However, as we have seen, correlation properties along streamlines contribute an additional factor ∼ V E 2-θ , leading overall to an asymptotic dispersion coefficient which grows with Pe,

D ∞ D ∼ Pe 2+θ 1+2θ . (62)
Our approach thus predicts a nontrivial scaling of the dispersion coefficient at high Pe for broadlydistributed crossing times. The presence of a maximum correlation time implies a return to Fickianity at sufficiently late times whenever diffusion is present, and the corresponding dispersion coefficient results from the interplay between diffusion and advective heterogeneity, which is mediated by the effective shear rate. The scaling exponent characterizing the behavior of the dispersion coefficient with Pe decreases with θ. As θ approaches unity, we recover the usual linear Pe dependency characteristic of advection-dominated transport with finite-mean crossing times. On the other hand, as θ approaches zero, the dependency approaches quadratic, as for classical Taylor dispersion.

It should be noted that our scaling predictions differ from previous results. [START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF] employed a heuristic cutoff at the level of the transit times, which they identified with the diffusion time associated with homogenizing a longitudinal correlation length. This led, for 0 < θ < 1, to the prediction of an intermediate scaling of the longitudinal dispersion coefficient related to the stable exponent of the crossing times as Pe 3-(1+θ) = Pe 2-θ , followed by a regime linear in Pe. In contrast, our results suggest that shear plays a key role in the interplay between heterogeneous advection and diffusion, leading to D ∞ ∼ Pe (2+θ)/(1+2θ) . The case presented in [START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF] corresponds to θ = 0.8, leading to a transition between D ∞ ∼ Pe 1.2 and D ∞ ∼ Pe in their model for high Péclet numbers. Our approach predicts instead a single regime in this range of Péclet with an intermediate scaling exponent ≈ 1.08. Note that this scaling might be difficult to distinguish experimentally from a transition between the two exponents 1.2 and 1. The largest differences between the two modelling frameworks would be observed for the smallest θ > 0, where our approach predicts a scaling exponent approaching 2 rather than unity.

In order to characterize the different regimes that arise for late-time longitudinal dispersion in more detail, we proceed to nondimensionalize the dynamical variables. Diffusive transitions are associated with a timescale τ ⊥ = 2 ⊥ /(2D) and a corresponding longitudinal lengthscale L / / = V E τ ⊥ . On the other hand, advective transitions are associated with the correlation length / / , which corresponds to the timescale τ / / = / / /V E . We nondimensionalize velocity as v * = v/V E , distance as s * = s/L / / , and time as t * = t/τ ⊥ . Consider now the transition operator L = L A + L D , Eq. ( 24). In terms of the nondimensional variables, the dimensionless transition operator

L * = L / / L = ζ Pe L * A + L * D /2, with L * A,D = L / / L A,D and ζ = τ ⊥ τ / / = 2 ⊥ 2 2 / / (63)
the ratio of longitudinal and transverse timescales.

This nondimensional form highlights the dominant transition mechanisms occurring under different flow conditions. Diffusive transitions dominate for Pe 1/ζ, and advective transitions for Pe 1/ζ. For Pe 1, diffusion dominates longitudinal dispersion, a scenario which we have disregarded here. According to Taylor dispersion, Eq. ( 42), and the definition of the transverse lengthscale, Eq. ( 54), the late-time longitudinal dispersion coefficient when only diffusive transitions are considered is given by

χ 2 D ∞ D = 2Iζ Pe 2 . ( 64
)
Thus, the effects of diffusion in the longitudinal direction equal those of transverse diffusion on the asymptotic dispersion coefficient when Pe = Pe ⊥ , with

Pe ⊥ = (2Iζ) -1/2 , (65)
for which χ 2 D ∞ /D = 1. This implies that the transverse-diffusion-dominated regime has a width of order 1/ √ ζ in the Péclet scale and disappears for sufficiently small ζ. The details of the behavior of the large-Pe regime and its onset depend on the particular form of the advective transitions and Eulerian PDF. We focus now on the example of a gamma distribution of Eulerian velocities together with a Bernoulli relaxation process of correlation length / / , as above. We have, by direct computation according to Eq. ( 47),

I = Γ(3θ) 27 θ Γ(θ)Γ(1 + θ) 2 . ( 66
)
Computing the maximum correlation time according to Eq. ( 59) with the effective shear rate according to the gamma PDF leads to

τ * max = τ max τ ⊥ = 2 Γ(θ) θ ζ Pe θ 2 1+2θ , (67)
valid for high Pe so that the minimum average velocity v min is small compared to V E . We consider an initial condition according to the Eulerian velocity PDF, corresponding to a transversely homogeneous spatial injection, at the longitudinal origin x = 0. Dispersion in the absence of diffusion is then given by Eq. ( 58). Nondimensionalizing the latter by dividing by the diffusion coefficient and substituting the maximum cutoff time for θ 1, Eq. ( 67), leads to an explicit form of Eq. ( 62) for asymptotic dispersion,

χ 2 D ∞ D =      [2/Γ(1 + θ)] 2(1-θ) 1+2θ [(1 -θ)ζ] -1 (θζ Pe) 2+θ 1+2θ , 0 < θ < 1 Pe ln 4e -3 ζ Pe /3, θ = 1 (θ -1) -1 Pe, θ > 1 , (68)
for t τ / / for θ > 1 and t τ max for θ 1. These results can be used to estimate the characteristic Pe corresponding to the onset of the large-Pe regime more precisely. Equating the two dispersion-regime expressions ( 64) and ( 68), we approximate the critical Pe = Pe / / associated with the onset of the regime dominated by advective transitions as

Pe / / =      θ[2/Γ(θ)] 2(1-θ) 3θ [2ζ θ (1 -θ)I] -1/θ , 0 < θ < 1 e 4 /(4ζ), θ = 1 [2ζ(θ -1)I] -1 , θ > 1 . ( 69
)
As expected from the discussion above, the advection-dominated regime is characterized by an onset inversely proportional to ζ in all cases. Note that for θ = 1 the asymptotic expressions for the intermediate and high Péclet regimes are never equal, due to the divergent behavior of the logarithm for small arguments. Thus, we estimate in this case Pe ⊥ as the value for which the two expressions 20 have a ratio closest to unity. The high-Pe expression is always smaller by a factor of 9e -4 ≈ 0.165 at this point. Thus, for sufficiently late times, we have three possible scaling regimes for the dependence of the dispersion coefficient on Péclet number: (i) Pe Pe ⊥ : longitudinal diffusion dominates; in this case, which we do not treat explicitly, D ∞ /D = 1; (ii) Pe ⊥ Pe Pe / / : transverse diffusion across streamlines dominates velocity transitions, leading to classical Taylor dispersion, D ∞ /D ∼ Pe 2 ; and (iii) Pe Pe / / : advective transitions along streamlines dominate, but transverse diffusion ensures Fickianity, and we have D ∞ /D ∼ Pe ϑ , with an exponent 1 ϑ < 2 that depends on the Eulerian PDF through θ (Eq. ( 68)). If Pe / / Pe ⊥ , there is no intermediate Pe 2 regime.

Figure 3 shows the behavior of the asymptotic longitudinal dispersion coefficient as a function of Péclet number for this setup, for fixed ζ = 10 -1 and θ = 1/2, 1, 3/2. Time is nondimensionalized by τ ⊥ , distance s by L / / , and position x by L / / /χ. A convenient parameterization of the nondimensional problem is obtained by setting D = 1/2, χ = 1, ⊥ = 1, / / = (2ζ) -1/2 , and V E = (ζ/2) 1/2 Pe. The free nondimensional parameters are Pe, ζ, and θ. These set, respectively, the relation between diffusive and advective timescales, the relation between the timescales of these two processes, and the low-velocity scaling of the Eulerian velocity PDF. The asymptotic dispersion coefficient in these units is given by

D * ∞ = χ 2 τ ⊥ D ∞ /L 2 / / , so that χ 2 D ∞ /D = 2D *
∞ in our parameterization. We set the discretization according to the minimum between Eq. ( 112) with ∆v = 5 • 10 -2 V E /θ and ∆s = 5 • 10 -2 / / , and we used 10 3 particles (see Appendix G for a discussion of the convergence properties of the DVRW). Dispersion was computed using the DVRW for 20 logarithmically-spaced asymptotic times between t * = 5 and t * = 10. The corresponding dispersion coefficient was obtained by backward-difference differentiation and averaged over the resulting values to reduce numerical error.

The results are in agreement with the theoretical predictions. For θ = 1/2, 1, 3/2, the regime onsets are estimated as (Pe ⊥ , Pe / / ) ≈ (6, 9 • 10 1 ), (8, 1 • 10 2 ), (10, 2 • 10 2 ), respectively. The three corresponding scaling regimes of dispersion with Péclet number are shown in Fig. 3. Note that, for θ = 3/2, there is a wide intermediate range of Péclet numbers (∼ 10 2 -10 6 ) over which different mechanisms are at play, causing the behavior of dispersion with Pe to deviate from the pure scaling laws derived above.

The pure advection-dominated regime scaling is then reached only at very large Pe ∼ 10 6 , whereas for θ = 1/2, 1 the asymptotic regime provides accurate predictions at much lower Pe ∼ 10 2 , 10 3 .

Conclusions

We have presented a new theoretical framework that couples heterogeneous advection and diffusion into a spatial-Markov stochastic process. Our model, the DVRW, allows for taking advantage of the fact that Lagrangian velocity dynamics in spatially structured flows are Markovian in space, while incorporating the effect of diffusion. Aside from its potential as a simulation technique, the DVRW formulation allows for the analytical development of dynamical equations for key transport quantities, highlighting the role of statistical properties of the flow in their evolution. In particular, we have shown how this approach provides new scaling laws describing dispersion properties at different Péclet numbers, explicitly incorporating the interplay between transverse diffusion and advective shear.

In the DVRW framework, diffusive transitions across nearby streamlines are determined by the following flow properties: (i) point statistics of velocity magnitude, as embodied by the Eulerian PDF;

(ii) average tortuosity; and (iii) a characteristic transverse lengthscale, encoding information about the transverse spatial and statistical structure of Eulerian velocities. These properties are combined into an effective shear rate, which, together with the usual diffusion coefficient, governs diffusive transitions in velocity space. As in previous works, transitions along streamlines are described by transition probabilities, which can typically be modeled with recourse to simple processes parameterized by (i) the Eulerian PDF; (ii) average tortuosity; and (iii) a longitudinal velocity correlation length.

The present work has been mainly concerned with formulating the approach while clarifying its physical meaning and characterizing the necessary parameters. Subsequent work is necessary towards parameterizing and applying the DVRW formulation to predict transport in realistic heterogeneous media, and in particular to test the practical applicability of the predicted scaling laws for dispersion. Since it quantifies single-particle trajectories in the presence of diffusion and heterogeneous advection, the DVRW formulation also opens up other promising avenues for future research, such as upscaling effective reaction dynamics in the presence of flow heterogeneity. again to leading order in ∆s,

ψ i (t; s + ∆s) = ∆v i j 0 ∆v j δ (T k + ∆s/v j -t) δ(V k -v j )δ(V k+1 -v i ), = j 0 r ij (s)ψ j (t -∆s/v j ; s) , (71)
where we have used the Markov property to break the average and identified the average of ∆v i δ(V k+1v i ) given V k = v j as the transition probability from velocity class j to i at distance s. Taylor-expanding ψ j (t -∆s/v j ; s) around t leads to

ψ i (t; s + ∆s) -ψ i (t; s) ∆s∆v i = j 0 r ij -δ ij ∆s∆v i ψ j (t; s) - 1 ∆v i v j ∂ψ i (t; s + ∆s) ∂t . ( 72
)
This yields the dynamical equation ( 7) in the limit ∆s → 0.

B Diffusion-averaged velocity class widths

Consider velocities discretized into classes [b i , b i+1 [. The width ∆v i = b i+1b i of class i represents the range of velocities averaged by diffusion over a spatial step of length ∆s along streamlines. The width of class i, associated with arithmetic average velocity v i , is given by ∆v i = α i 2D∆s/v i , where α i = α e (v i ). In the limit of small class widths, consider the approximation v i ≈ (b i+1 + b i )/2 and α i ≈ α e (b i ) respectively for the class average velocities and effective shear rates. Replacing

v i = b i + ∆v i /2 in ∆v i = α i 2D∆s/v i yields the cubic equation ∆v 2 i (∆v i + 2b i ) = 4Dα 2 i ∆s.
Using the standard cubic solutions and some algebra we arrive at

∆v i = a i (1 + ξ i ) 1/3 -(1 -ξ i ) 1/3 2 + (1 -ξ i ) 1/3 , ( 73
)
where a i = (Dα 2 i ∆s) 1/3 and ξ i = 1 -8b 3 i /(27a 3 i ). Class edges can then be obtained recursively up to a maximum velocity v M according to b j+1 = b i + ∆v i with b 0 = v m .

For 0 ξ 1, the expression for ∆v i is clearly real and positive. At small velocities, for 8b 3 i /(27Dα 2 i ∆s) 1, we obtain ∆v i ≈ (4Dα 2 i ∆s) 1/3 , which in particular holds exactly for class i = 0 if v m = 0. For ξ < 0, it may be convenient to work with the explicitly real expression

∆v i = a i (1 + |ξ|) 1/3 1 -4 sin 2 arctan |ξ| 3 , (74)
which is obtained by expressing the conjugate pair 1 ± |ξ| in terms of their absolute value and phase and simplifying the resulting expression. Note that the argument of the sine is smaller than π/6, so that this expression is also positive. Expanding for |ξ i | ≈ 8b 2 i /(27Dα 2 i ∆s) 1, we obtain ∆v i ≈ 2Dα 2 0 ∆s/b i , which holds at sufficiently high velocities. This situation corresponds to ∆v i 2b i , or v i ≈ b i .

The approach above breaks down if α e (b i ) is zero or divergent. This happens where the Eulerian PDF diverges or is zero at b i , respectively, see Eq. (55). The above construction can be employed separately within velocity intervals where the Eulerian PDF p E (v) is finite and nonzero, except that the approximation α e (v i ) ≈ α e (b i ) cannot be used at the leftmost class within each interval. In this case, we use α i = α e (b i + ∆v i /2), which leads to an implicit equation for ∆v i . This equation can be solved numerically in general, and possibly analytically if the local form of the velocity PDF is known.

C Diffusive transitions in the continuum limit

The transition probabilities associated with diffusive transitions, characterized by Eqs. ( 17) and ( 18

), give r D ij -δ ij = r + j δ i,j+1 -δ ij + r - j δ i,j-1 . (75)
Let q i be a quantity dependent on velocity class i, and let q(v) be the associated continuous density, i.e., q i = bi+1 bi dv q(v) ≈ ∆v i q(v i ). We have

j 0 r D ij -δ ij ∆s∆v i q j = r - i+1 q i+1 -q i + r + i-1 q i-1 ∆s∆v i . (76)
Expanding in Taylor series around v i ,

∆v i±1 ≈ ∆v i 1 ± ∂∆v(v) ∂v vi , (77)
where ∆v(v) = 2Dα e (v) 2 ∆s/v, and

∆v i±1 q(v i±1 ) ≈ ∆v i q(v i ) ± ∆v i ∂∆v(v)q(v) ∂v vi + ∆v 2 i 2 ∂ 2 ∆v(v)q(v) ∂ 2 v vi . (78)
We thus find, to leading order in ∆s,

r D ij -δ ij ∆s∆v i = 1 2∆s ∂∆v(v) ∂v ∂∆v(v)q(v) ∂v + ∆v(v) ∂ 2 ∆v(v)q(v) ∂v 2 -∆v(v) ∂ ∂v q(v) 2v vi , = 1 2∆s ∂ ∂v ∆v(v) 2 ∂q(v) ∂v -q(v) 1 v - ∂ ∂v ∆v(v) 2 2 vi . (79)
This leads to ∂q(v; s) ∂s = L D q(v; s),

with the diffusive transition operator L D given according to Eq. ( 19).

D Spatial extrema and the Eulerian PDF

We present here a simple scaling argument to identify the qualitative impact of a spatial extremum (maximum or minimum) on the Eulerian PDF.

E (v M + ∆v) ∝ |∆x| d-1 |Ω||∇ 2 v||∆x| = |∆v| d/2-1 |Ω||∇ 2 v| d/2 . ( 81
)
This means that a non-degenerate spatial extremum corresponds to a square-root divergence in d = 1, a ∆v-independent contribution in d = 2, and a zero in d = 3. It is interesting to note that the Eulerian PDF associated with Poiseuille flow in a cylindrical pipe is velocity-independent. This can be understood as follows. The probability density for a given velocity is inversely proportional to the corresponding velocity gradient magnitude, which in turn is proportional to 1v/v M . In a cylindrical pipe, the flow profile is a paraboloid, and the higher gradients near the pipe walls are compensated by the fact that, due to the axisymmetry of the profile, the spatial frequency of occurrence of a given velocity is proportional to |Λ(v)| ∝ 1v/v M . This results in a constant probability density across velocities. Couette flow in two spatial dimensions, the constant-shear flow arising between a stationary plate and a plate moving at constant velocity, also has this property, because the velocity gradient is constant and each velocity occurs the same number of times (exactly once). However, this flow presents no non-degenerate spatial maxima.

The average position is thus given by X T (t) = ρ 1 (t) and longitudinal dispersion by the variance σ 2 (t) = ρ 2 (t)ρ 1 (t) 2 . Multiplying Eq. ( 7) by s and integrating in s, we find, for k 1, the recursion relations

v -1 ∂I k (v, t) ∂t = L D I k (v, t) + kI k-1 (v, t). (87)
Integrating in v yields also

dρ k (t) dt = χ -k k ∞ 0 dv I k-1 (v, t). (88)
This gives directly ρ 0 (t) = 1, confirming that the concentration PDF is correctly normalized. It is easy to see by inspection that Eq. ( 87) for k = 1 admits asymptotic solutions of the form I ∞ 0 (v) ∝ p F (v), so that, in order to satisfy ρ 0 (t) = 1, we have

I ∞ 0 (v) = V E p F (v)
. Then, Eq. ( 88) for k = 1 yields the asymptotic mean position ρ 1 (t) = χ -1 V E t. This means that, as expected, the average position grows asymptotically according to the mean velocity (corrected by the tortuosity).

The asymptotic longitudinal dispersion coefficient is obtained through a similar, albeit more involved, approach. We expect the asymptotic form ρ 2 (t) = 2D ∞ t + χ -2 V E 2 t 2 for the second raw moment of position, where D ∞ is the asymptotic longitudinal dispersion coefficient. Therefore, we look for an asymptotic solution to Eq. ( 87) for k = 1 of the form

I 1 (v, t) = a(v) + b(v)t. Satisfying Eq. (88) with k = 2 requires ∞ 0 dv a(v) = D ∞ and ∞ 0 dv b(v) = χ -2 V E 2 .
On the other hand, Eq. ( 87)

requires L D b(v) = 0; we conclude that b(v) = χ -2 V E 2 p F (v) and L D a(v) = χ -2 V E [p E (v) -p F (v)] . (89)
Using the boundary and initial conditions for ψ(v, t; s), integrating this equation yields

a(v) = Ap F (v) + V E v χ 2 Dα e (v) v 0 dv C E (v ) -C F (v ) α e (v ) , (90)
where A is a so-far arbitrary coefficient. Using Eq. ( 36) for the effective shear, we obtain

a(v) = p F (v) A + V E 2 δv 2 χ 2 Dα e (V E ) 2 g(v) , g(v) = v 0 dv p E (v )[C E (v ) -C F (v )]. (91)
Next, we use Eq. ( 86) for ρ 1 (t), which requires ∞ 0 dv a(v)/v = 0 (as well as

∞ 0 dv b(v)/v = χ -2 V E , which is satisfied). This fixes A and leads to a(v) = V E 2 δv 2 p F (v) χ 2 α e (V E ) 2 ∞ 0 dv g(v ) [δ(v -v ) -p E (v )] . ( 92
) Setting ∞ 0 dv a(v) = D ∞ yields D ∞ = IV E 2 δv 2 χ 2 Dα e (V E ) 2 , (93)
where v). Integrating by parts leads to

I = ∞ 0 dv [p F (v) -p E (v)]g (
I = ∞ 0 dv p E (v)[C E (v) -C F (v)] 2 . ( 94
)
The quantity I may also be expressed as a spatial integral. To this end, we express p E (v) as a spatial average according to its definition, Eq. ( 25). Then, the Eulerian CDF is given by

C E (v) = v 0 dv |Ω ⊥ | -1 Ω ⊥ dy δ[v E (y) -v], = |Ω ⊥ | -1 Ω ⊥ dy H[v -v E (y)], (95)
where H(•) is the Heaviside step function. Similarly, the flux-weighted Eulerian CDF reads

C F (v) = |Ω ⊥ | -1 Ω ⊥ dy v E (y) V E H[v -v E (y)]. (96)
Substituting, we find

I = |Ω ⊥ | -3 Ω ⊥ dy ∞ 0 dv δ[v E (y) -v]   Ω ⊥ dy ν(x )H[v -v E (y )]   2 , (97)
where ν(y) = v E (y)/V E -1 are the normalized velocity fluctuations. Performing the integral over v leads directly to Eq. ( 51).

F.2 Purely-advective velocity transitions

The temporal velocity correlation function is given, by definition, by

C v (t, t ) = V T (t)V T (t ) -V T (t) V T (t ), (98)
where V T (t) is the Lagrangian velocity at time t. Following [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF], we introduce the propagator g(v, t) = e -vt/ / / (99)

and the transition time PDFs across a correlation length associated with the velocity PDF p J ,

φ J (t) = -1 / / ∞ 0 dv g(v, t)vp J (v), (100)
where J ∈ {0, E, F } and p 0 (v) = p T (v; 0) is the initial velocity PDF. We have for the Bernoulli process [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF],

V T (t)V T (t ) = ∞ 0 dv V T (t -t |v)vp T (v, t ), (101)
where V T (t|v) is the average velocity at time t + t conditioned on velocity v at time t . Its Laplace transform is given by [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF] ṼT (λ|v) = vg(v, λ)

+ ∞ 0 dv vv 2 p E (v ) / / V E g(v, λ)g(v , λ) 1 -φF (λ) . ( 102
)
where the tilde denotes the Laplace transform (with respect to time) and λ is the associated Laplace variable. Using the definition of φ F (t), this can be rewritten as ṼT (λ|v) = vg(v, λ) 1 -φF (λ) .

(103)

G Numerical validation of diffusive transitions

We consider, for validation and illustration purposes, transport in some instances of the stratified flows from Table 1. We employ a pulse initial condition of unit mass, at the origin along the flow direction, and uniform along the direction transverse to the flow, corresponding to an initial velocity distribution according to the Eulerian PDF. We fix, in arbitrary units, L = 1, D = 5 • 10 -4 , and V E = 1. This corresponds to a diffusion time τ D = L 2 /(2D) = 10 3 to homogenize the transverse direction, an associated advective distance L D = V E τ D = 10 3 , and a Péclet number associated with the width L of V E L/D = 2 • 10 3 . We compare the results of our DVRW, implemented using the discrete Lagrangian formulation, to standard fully-resolved particle tracking random walk (PTRW) simulations with fixed-time-step spatial transitions associated with advection and transverse diffusion. The simulation results are in very good agreement for the breakthrough curves and concentration distributions as well as the velocity PDFs at fixed times and distances, as shown for two-dimensional Poiseuille flow in Figs. 4 and5, respectively. The results in these figures correspond to ∆s = 10 -6 L D for the DVRW and ∆t = 10 -5 τ D for the PTRW simulations, using 10 5 particles in both cases.

A slight mismatch can be observed for concentrations very close to injection, see the right panel of Fig. 4, which occur at early times. This is due to the impact of discretization, which is most noticeable at positions and times corresponding to few transitions. At large distances (compared to L D ), both concentrations and breakthrough curves become approximately Gaussian, in accordance with the central limit theorem. Regarding the velocity PDFs, Fig. 5, note how there is no change in the time-Lagrangian PDF, since the initial condition coincides with the Eulerian PDF and is therefore the equilibrium solution. As for the space-Lagrangian PDF, note the quick evolution towards the equilibrium flux-weighted Eulerian PDF: for the last two distances, corresponding s = 2.6•10 -1 L D and 10L D , the PDF has essentially already converged. Slight mismatches between the simulation results and the theoretical equilibria are due to the sampling discretization when computing the PDFs.

The temporal evolution of the mean tracer position and longitudinal dispersion for four examples of stratified flow are shown in Fig. 6. The discretizations used in this case were ∆s = 10 -6 L D , ∆t = 10 -6 τ D , and 10 4 particles. The results show excellent agreement between PTRW and DVRW simulations.

As usual for particle methods, precision scales approximately with the inverse of the square root of the number of particles. As for the spatial discretization, note that the relevant discretization for the DVRW is at the level of the Eulerian PDF, which is achieved indirectly through discretizing the distance ∆s along streamlines. In order to estimate the discretization necessary to achieve a given precision in velocity space, consider Eq. ( 16) describing diffusive velocity averaging and characterizing velocity class widths in the discretized DVRW. Inverting this relation for ∆s and using Eq. ( 55) for the effective shear rate, we have, for a given precision ∆v, ∆s = min

i {v i p E (v i ) 2 }τ ⊥ ∆v 2 , ( 112
)
where τ ⊥ = 2 ⊥ /(2D) is the diffusion time associated with homogenizing a region of characteristic length ⊥ . If vp E (v) has zeros, the minimum value on the right hand side occurs within ∆v of a zero. Accurate results typically require ∆v/V E 1. The convergence with ∆s of longitudinal dispersion for Poiseuille flow is illustrated in Fig. 7. The order of convergence is found to be approximately 1/3, corresponding to the dependence of velocity class sizes on ∆s 1/3 at small velocities, see Appendix B. Here, fine discretizations are needed for accurate results. However, when velocity transitions along streamlines are included, the discretization requirements become less strict at high Péclet number. Although we do not consider this point further The baseline solution used for comparison is computed using a PTRW simulation with ∆t = 10 -5 . All solutions use 10 5 particles. Linear fitting to the logarithm of the error versus the logarithm of ∆s indicates an order of convergence of 0.31, with 95% confidence bounds of 0.30 and 0.33.

here, we expect it may be possible to develop faster-converging Lagrangian discretizations of the continuous DVRW description. An alternative is to directly solve the Eulerian continuum equations for the PDFs of quantities of interest, where the discretization may be imposed directly at the level of velocity magnitudes.

Figure 1 :

 1 Figure1: Conceptual illustration of the central concepts behind the diffusing-velocity random walk (DVRW). As Lagrangian particles are transported, they undergo advective changes in velocity due to variability along streamlines as well as diffusion-induced averaging over, and transitions to, nearby streamlines. The transverse variability of the flow field is encoded in a velocity-dependent effective shear rate α e (v), which is determined by statistical properties of the flow field. This conceptual illustration is two-dimensional and includes the presence of a solid phase (black circles), but the proposed approach is applicable also to three dimensions and arbitrary flow fields characterized by well-defined characteristic lengths along the longitudinal and transverse directions, as developed in detail in the main text.

  is constant and equal to α, is chosen to highlight the role of the multiplicity (number of spatial occurrences) |Λ(v)| associated with each velocity. Velocity magnitude v occurs |Λ(v)| = Λ M = 4 times for v larger than a critical velocity v c , and |Λ(v)| = Λ m = 2 times for v < v c . This implies that the center and outer parts of the flow cover velocity variations δv M,m = |v M,mv c | over lengths M,m = Λ M,m δv M,m /α, respectively. Similarly, the effective shear rate α M,m e ∝ δv M,m / M,m = α/Λ M,m for velocities above and below the critical velocity v c . Noting that p E (v) ∝ |Λ(v)|/α, see Eq. (28), we see that, indeed, α e (v) ∝ 1/p E (v).

Figure 2 :

 2 Figure 2: Transverse velocity magnitude profile for the M -flow. The key relevant features of this illustrative example are: (i) The multiplicity (number of spatial occurrences) of velocity magnitudes above (below) the critical value v c is Λ M = 4 (Λ m = 2); (ii) The transverse gradient magnitude is constant. It is given by α = δv m /( m /Λ m ) = δv M /( M /Λ M ).

Figure 3 :

 3 Figure 3: Asymptotic longitudinal dispersion coefficient as a function of Péclet number, with velocity transitions along streamlines according to a Bernoulli relaxation process and transitions across streamlines according to transverse diffusion. The three values of θ shown highlight the three possible asymptotic qualitative behaviors discussed in the text. The inset shows corresponding Eulerian velocity gamma PDFs, with low velocities having a probability density ∼ v θ-1 .

Figure 4 :Figure 5 :

 45 Figure 4: Breakthrough curves (left) and concentration profiles (right) for two-dimensional Poiseuille flow, under a pulse initial condition homogeneous along the direction transverse to the flow. Times are shown in units of the diffusion time τ D and positions in units of the advective distance L D .

Figure 6 :

 6 Figure 6: Mean (left) and variance (right) of tracer positions for different stratified flows, under a pulse initial condition homogeneous along the direction transverse to the flow. Times are shown in units of the diffusion time τ D and position in units of the advective distance L D . Variances are normalized in terms of the Taylor dispersion coefficient as indicated.

Figure 7 :

 7 Figure 7: Convergence of the DVRW with discretization for longitudinal dispersion under Poiseuille flow in two dimensions. The convergence properties are quantified through the L 2 relative error norm, computed for the solutions with different ∆s at 1000 equally-spaced times between t = 0 and t = 10τ D .The baseline solution used for comparison is computed using a PTRW simulation with ∆t = 10 -5 . All solutions use 10 5 particles. Linear fitting to the logarithm of the error versus the logarithm of ∆s indicates an order of convergence of 0.31, with 95% confidence bounds of 0.30 and 0.33.

  Near a non-degenerate (point) extremum v M in d dimensions, at a spatial distance ∆x, we have |Λ(v)| ∝ |∆x| d-1 . The associated change in velocity is |∆v| ∝ |∇v||∆x|, and the gradient behaves as |∇v| ∝ |∇ 2 v||∆x|, so that |∆x| ∝ |∆v| 1/2 /|∇ 2 v| 1/2 . Thus, the contribution near a spatial extremum to the Eulerian PDF is p

Acknowledgements

TA is supported by a Marie Sk lodowska Curie Individual Fellowship, under the project ChemicalWalks 792041. TLB gratefully acknowledges funding by the ERC under the project ReactiveFronts 648377.

Declaration of interests

The authors report no conflict of interest.

A Joint PDF of velocity and arrival time

The joint PDF of velocity and arrival time, given distance s traveled along streamlines, is defined as

where the overline denotes an ensemble average over tracer particles and δ(•) is the Dirac delta. Reverting to the discrete version of the dynamics, see Eq. (1), we consider the joint density ψ i (t; s) of arriving at time t and having velocity in class i given distance s. We have

to first order in ∆s. Using Eq. ( 1) to express T k+1 in terms of T k and V k and introducing a partition of unity for the latter being in velocity class j we obtain,

E Taylor dispersion

The Taylor dispersion coefficient is given by D T = ηχ -2 V E 2 L 2 /D. Here, we discuss the dimensionless coefficient η. It is given by η = |Ω ⊥ | -1 Ω ⊥ dy ν(y)ϕ(y), where the dimensionless function ϕ solves [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF])

where n is the unit outward normal to the transverse domain boundary ∂Ω ⊥ . We note that this determines ϕ up to an additive constant, which does not affect η because the spatial integral of the velocity fluctuations is null.

If the flow depends only on a one-dimensional transverse coordinate y, integrating this equation leads, up to an additive constant, to

Thus,

Integrating by parts in the outer integral leads to Eq. ( 44). For axisymmetric flows, depending only on r = √ y • y, we have L/2 0 dr rν(r) = 0. Using ∇ 2 = r -1 d 2 /dr 2 r for the Laplacian in cylindrical coordinates and |Ω ⊥ | = π(L/2) 2 , and following the same procedure as before, we obtain Eq. ( 45).

F Longitudinal dispersion

This appendix is concerned with determining asymptotic longitudinal dispersion. First, we present the derivation for the DVRW formulation under purely-diffusive velocity transitions. Next, we discuss the dispersion coefficient for advective transitions only, for the case of a Bernoulli relaxation process.

F.1 Purely-diffusive velocity transitions

We start by defining

According to Eq. ( 9), the raw moments of concentration are given, for k 0, by

Similarly, we have [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF])

from which

We focus here on the case of an initial condition according to the Eulerian PDF, p 0 (v) = p E (v). Direct computation from Eq. ( 100) shows that 1 -t 0 dt φ F (t ) = φ E (t) / / /V E , from which we conclude that 1 -φF (λ) = ( / / λ/V E ) φE (λ). This shows immediately that ṼT (λ) = V E /λ, so that, inverting the Laplace transform, the average velocity is constant and equal to V E , as expected. Furthermore, as shown in the main text, the time-Lagrangian velocity PDF is in this case stationary, and we have p T (v, t) = p E (v). Therefore, from Eq. ( 101),

We conclude that the velocity correlation function is stationary; for t > t ,

Taking Laplace transforms and using Eq. ( 103) yields

Thus, using the Green-Kubo relation (41) gives

For gamma-distributed Eulerian velocities, Eqs. ( 56) and (57), we take the Laplace transform of Eq. ( 100) and perform the integral to obtain φF

where E θ (x) = ∞ 1 dt exp(-xt)/t θ is an exponential integral. This admits the small-λ expansions φF (λ) ≈

    

1τ / / λ -Γ(-θ)(τ / / λ) 1+θ , θ < 1 1τ / / λln(e γ E τ / / λ)(τ / / λ) 2 , θ = 1 1τ / / λ + [θ/(θ -1)](τ / / λ) 2 , θ > 1

where τ / / = / / /V E and γ E ≈ 0.577 is the Euler-Mascheroni constant. Substituting in Eq. ( 109) and inverting the Laplace transform to leading order in λ leads to Eq. ( 58).