Peter K Kang

Marco Dentz

Tanguy Le Borgne

Seunghak Lee

Ruben Juanes

Anomalous transport in disordered fracture networks: spatial Markov model for dispersion with variable injection modes

Keywords: Discrete Fracture Networks, Injection Modes, Anomalous Transport, Stochastic Modelling, Lagrangian Velocity, Time Domain Random Walks, Continuous Time Random Walks, Spatial Markov Model

port through fractured media: the equivalent porous medium approach [START_REF] Neuman | Stochastic theory of fieldscale Fickian dispersion in anisotropic porous media[END_REF][START_REF] Tsang | Tracer transport in a stochastic continuum model of fractured media[END_REF] 7 and the discrete fracture network approach (DFN) [START_REF] Kiraly | Remarques sur la simulation des failles et du réseau karstique par éléments finis dans les modèles d'écoulement[END_REF]9,10,[START_REF] Moreno | Fluid flow and solute transport in a network 527 of channels[END_REF][START_REF] Juanes | A general and efficient formulation 529 of fractures and boundary conditions in the finite element method[END_REF][START_REF] Park | Transport be-532 havior in three-dimensional fracture intersections[END_REF], 14, 8 [START_REF] Martinez-Landa | An analysis of hydraulic conductiv-538 ity scale effects in granite (Full-scale Engineered Barrier Experiment 539 (FEBEX), Grimsel, Switzerland)[END_REF]16,[START_REF] De Dreuzy | Influence of fracture scale het-545 erogeneity on the flow properties of three-dimensional discrete fracture 546 networks (DFN)[END_REF][START_REF] Schmid | Higher order FE-FV method 548 on unstructured grids for transport and two-phase flow with variable 549 viscosity in heterogeneous porous media[END_REF][START_REF] Makedonska | 552 Particle tracking approach for transport in three-dimensional discrete 553 fracture networks[END_REF][START_REF] Hyman | dfnWorks: A discrete fracture network framework 556 for modeling subsurface flow and transport[END_REF][START_REF] Bernabe | Passive advection-dispersion in 559 networks of pipes: Effect of connectivity and relationship to perme-560 ability[END_REF]. The DFN approach explicitly resolves individ-9 ual fractures whereas the equivalent porous medium approach represents 10 the fractured medium as a single continuum by deriving effective parame-11 ters to include the effect of the fractures on the flow and transport. The 12 latter, however, is hampered by the fact that a representative elementary 13 volume may not exist for fractured media [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF]23]. Dual-porosity mod-14 els are in between these two approaches, and conceptualize the fractured-15 porous medium as two overlapping continua, which interact via an exchange 16 term [START_REF] Barenblatt | Basic concepts in the 567 theory of seepage of homogeneous liquids in fissured rocks [strata[END_REF][START_REF] Warren | The behavior of naturally fractured reservoirs, 570[END_REF][START_REF] Kazemi | Numerical 572 simulation of water-oil flow in naturally fractured reservoirs[END_REF][START_REF] Bibby | Mass transport of solutes in dual-porosity media, Water Re-575 sour[END_REF][START_REF] Feenstra | Matrix diffusion effects on 577 contaminant migration from an injection well in fractured sandstone[END_REF][START_REF] Maloszewski | On the theory of tracer experiments in fis-580 sured rocks with a porous matrix[END_REF][START_REF] Pruess | A practical method for modeling fluid and heat flow in 582 fractured porous media[END_REF][START_REF] Arbogast | Derivation of the double porosity 584 model of single phase flow via homogenization theory[END_REF]32]. rock matrix [33,34,35,36]. In practice, however, their application must 23 account for the uncertainty in the subsurface characterization of fractured 24 media, which is still an considerable challenge [37,38,39]. Thus, there is a Developing an upscaled model for transport in fractured media is espe-31 cially challenging due to the emergence of anomalous (non-Fickian) trans-32 port. While particle spreading is often described using a Fickian framework, 33 anomalous transport-characterized by scale-dependent spreading, early ar-34 rivals, long tails, and nonlinear scaling with time of the centered mean 35 square displacement-has been widely observed in porous and fractured 36 media across multiple scales, from pore [42,43,44,45,46] to single frac-37 ture [47,48,49,50] to column [51,52] to field scale [53,54,55,56,57,41,58].

38

The ability to predict anomalous transport is essential because it leads to fun- The mode of particle injection can have a major impact on transport through porous and fractured media [START_REF] Kreft | On the physical meaning of the dispersion equation 737 and its solutions for different initial and boundary conditions[END_REF][START_REF] Sposito | Predicting solute plume evolution in heteroge-740 neous porous formations[END_REF]83,[START_REF] Frampton | Significance of injection modes and het-745 erogeneity on spatial and temporal dispersion of advecting particles in 746 two-dimensional discrete fracture networks[END_REF][START_REF] Janković | Analysis of the impact of injection mode in 749 transport through strongly heterogeneous aquifers[END_REF][START_REF] Le Borgne | Non-Fickian mixing: Temporal evolution of the scalar dissipa-753 tion rate in heterogeneous porous media[END_REF]58,[START_REF] Dagan | Solute plumes mean velocity in aquifer transport: Impact 756 of injection and detection modes[END_REF]. Two generic injection modes are uniform (resident) injection and flux-weighted injection with distinctive physical meanings as discussed in Frampton and Cvetkovic [START_REF] Frampton | Significance of injection modes and het-745 erogeneity on spatial and temporal dispersion of advecting particles in 746 two-dimensional discrete fracture networks[END_REF]. The work by Sposito and Dagan [START_REF] Sposito | Predicting solute plume evolution in heteroge-740 neous porous formations[END_REF] is one of the earliest studies of the impact of different particle injection modes on the time evolution of a solute plume spatial moments. The significance of injection modes on particle transport through discrete fracture networks has been studied for fractured media [START_REF] Frampton | Significance of injection modes and het-745 erogeneity on spatial and temporal dispersion of advecting particles in 746 two-dimensional discrete fracture networks[END_REF]58]. Dagan [START_REF] Dagan | Solute plumes mean velocity in aquifer transport: Impact 756 of injection and detection modes[END_REF] recently clarified the theoretical relation between injection modes and plume mean velocity. Despite recent advances regarding the significance of particle injection modes, the incorporation of injection methods into effective transport modelling is still an open issue [START_REF] Frampton | Significance of injection modes and het-745 erogeneity on spatial and temporal dispersion of advecting particles in 746 two-dimensional discrete fracture networks[END_REF]. The fundamental challenge is that the Lagrangian velocity distribution experienced by tracer particles evolves with time from its initial distribution

A C C E P T E D M A N U S C R I P T
which is dictated by the injection mode to a stationary velocity distribution [73,[START_REF] Frampton | Significance of injection modes and het-745 erogeneity on spatial and temporal dispersion of advecting particles in 746 two-dimensional discrete fracture networks[END_REF]75,[START_REF] Dentz | Con-758 tinuous time random walks for the evolution of Lagrangian velocities[END_REF]. In this paper, we address these fundamental questions, in the context of anomalous transport through disordered DFNs.

The paper proceeds as follows. In the next section, we present the studied random discrete fracture networks, the flow and transport equations and details of the different particle injection rules. In Section 3, we investigate the The link length between nodes i and j is denoted by l ij . The network is characterized by the distribution of link lengths p l (l) and hydraulic conductivity K. The PDF of link lengths here is exponential

p l (l) = exp(-l/ l) l . (1)
Note that the link length and orientation are independent. The characteristic fracture link length is obtained by taking the average of a link length over all the realizations, which gives l ≈ L x /200. A realization of the random discrete fracture network is generated by assigning independent and identically distributed random hydraulic conductivities K ij > 0 to each link between nodes i and j. Therefore, the K ij values in different links are un-

A C C E P T E D M A N U S C R I P T correlated.
The set of all realizations of the spatially random network generated in this way forms a statistical ensemble that is stationary and ergodic.

We assign a lognormal distribution of K values, and study the impact of conductivity heterogeneity on transport by varying the variance of ln(K).

We study log-normal conductivity distributions with four different variances:

σ lnK = 1, 2, 3, 5.
The use of this particular distribution is motivated by the fact that conductivity values in many natural media can be described by a lognormal law [90, 91].

Flow Field

Steady state flow through the network is modeled by Darcy's law [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF] for the fluid flux u ij between nodes i and j, u ij = -K ij (Φ j -Φ i )/l ij , where Φ i and Φ j are the hydraulic heads at nodes i and j. Imposing flux conservation at each node i, j u ij = 0 (the summation is over nearest-neighbor nodes), leads to a linear system of equations, which is solved for the hydraulic heads at the nodes. The fluid flux through a link from node i to j is termed incoming for node i if u ij < 0, and outgoing if u ij > 0. We denote by e ij the unit vector in the direction of the link connecting nodes i and j.

We

p e (u) = i>j l ij δ(u -u ij ) N l . ( 2
)
where N is the number of links in the network. Link length and flow velocities here are independent. Thus, the Eulerian velocity PDF is given by

p e (u) = 1 N i>j δ(u -u ij ). (3)
Even though the underlying conductivity field is uncorrelated, the mass conservation constraint together with heterogeneity leads to the formation of preferential flow paths with increasing network heterogeneity [92,[START_REF] Kang | Anomalous transport 734 on regular fracture networks: impact of conductivity heterogeneity and 735 mixing at fracture intersections[END_REF]. This is illustrated in Figures 2a andb, which show maps of the relative velocity magnitude for high velocities in networks with log-K variances of 1 and 5. As shown in Figures 2c andd, for low heterogeneity most small flux values occur along links perpendicular to the mean flow direction, whereas low flux values do not show directionality for the high heterogeneity case. This indicates that fracture geometry dominates small flux values for low heterogeneity and fracture conductivity dominates small flux values for high heterogeneity. An increase in conductivity heterogeneity leads to a broader Eulerian velocity PDF, with significantly larger probability of having small flux values as illustrated in Figure 3, which shows p e (u) for networks of different heterogeneity strength.

Transport

Once the fluxes at the links have been determined, we simulate transport

N i 0 = N p i 0 , (4)
where N i 0 is the number of particles injected at node i 0 , N p is the total number of injected particles. Flux-weighted injection introduces particles proportional to the total incoming flux Q i 0 at the injection location i 0

N i 0 = N p Q i 0 i 0 Q i 0 . (5)
Uniform injection simulates an initial distribution of tracer particles extended uniformly over a region much larger than the characteristic heterogeneity scale, and flux-weighted injection simulates a constant concentration pulse where the injected mass is proportional to the local injection flux at an inlet boundary that is much larger than the heterogeneity scale. For the uniform injection, the initial velocity distribution is then equal to the distribution of the Eulerian velocities. For the flux-weighted injection, the initial velocity distribution is equal to the flux-weighted Eulerian distribution. In general the initial velocity distribution may be arbitrary and depends on the conditions at the injection location. More detailed discussions can be found in section 4 and section 5.

Injected particles are advected with the flow velocity u ij between nodes.

To focus on the impact of conductivity variability on particle transport, we assume porosity to be constant. This is a reasonable assumption because

A C C E P T E D M A N U S C R I P T
the variability in porosity is significantly smaller than the the variability in conductivity [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF]93].

At the nodes, we apply a complete mixing rule [94,95,96]. Complete mixing assumes that Péclet numbers at the nodes are small enough that particles are well mixed within the node. Thus, the link through which the particle exits a node is chosen randomly with flux-weighted probability. A different node-mixing rule, streamline routing, assumes that Péclet numbers at nodes are large enough that particles essentially follow the streamlines and do not transition between streamlines. The complete mixing and streamline routing rules are two end members. The local Péclet number and the intersection geometry determine the strength of mixing at nodes, which is in general between these two end-members. The impact of the mixing rule on transverse spreading can be significant for regular DFNs with low heterogeneity [97, 80].

However, its impact is much more limited for random DFNs [96]. Since our interest in this study is the longitudinal spreading in random DFNs, we focus on the case of complete mixing. Thus, the particle transition probabilities p ij from node i to node j are given by

p ij = |u ij | k |u ik | , (6)
where the summation is over outgoing links only, and p ij = 0 for incoming links. Particle transitions are determined only by the outgoing flux distribution.

The particle pathways and times are obtained by the following recursion

A C C E P T E D M A N U S C R I P T relations x n+1 = x n + n e n , (7a)
t n+1 = t n + n u n , (7b)
where x n ≡ x in is the particle position after n random walk steps, n ≡ l ini n+1

the particle displacement and e n ≡ e ini n+1 its orientation; the particle velocity at the nth step is denoted by u n ≡ |u ini n+1 |. The particle displacement, orientation and velocity determined by the transition probability p inj from node i n to the neighboring nodes j given by Eq. ( 6). Equations ( 7) describe coarse-grained particle transport for a single realization of the spatial random network. Particle velocities and thus transition times depend on the particle position. The particle position at time t is x(t) = x in t , where

n t = sup(n|t n ≤ t) (8)
denotes the number of steps needed to reach time t. We solve transport in a single disorder realization by particle tracking based on Eq. ( 7) with the two different injection rules (4) and ( 5) at the inlet at x = 0. The particle density in a single realization is

p(x, t) = δ(x -x nt ) , (9)
where the angular brackets denote the average over all injected particles. As injection. The impact of network heterogeneity and injection method can be clearly seen from projected concentration profiles, f τ (ω). Arbitrary injection modes are also studied and discussed in section 5.2.

Average Solute Spreading Behavior

We first study the average solute spreading behavior for the four different levels of conductivity heterogeneity and the two different injection methods described above. We first illustrate the persistent effect of the particle injection method on particle transport, with the two different injection modes.

To investigate the average spreading behavior, we average over all particles and network realizations. The average particle density is given by

P (x, t) = δ(x -x nt ) , (10)
where the overbar denotes the ensemble average over all realizations. We run Monte Carlo particle tracking simulations for 100 realizations for each combination of conductivity heterogeneity and particle injection rule. In each realization, we release 10 4 particles at the inlet (x = 0) with the two different injection methods.

Breakthrough Curves

The average particle spreading behavior is first studied with the first passage time distribution (FPTD) or breakthrough curve (BTC) of particles at a control plane located at x = x c . The FPTD is obtained by averaging over the individual particle arrival times

τ a (x c ) = inf(t n | |x n -x 0 | > x c ) as f (τ, x c ) = δ[τ -τ a (x c )] . (11)
A

C C E P T E D M A N U S C R I P T (e) (d) flux uniform (f) (b) (a) (c) x x σ lnK = 1 , σ lnK = 1 , σ lnK = 5 , flux flux uniform σ lnK = 5 , σ lnK = 5 , σ lnK = 1 , flux f τ ( ω) f τ ( ω)

Centered Mean Square Displacement

We also study longitudinal spreading in terms of the centered mean square displacement (cMSD) of average particle density, P (x, t). For the longitudinal direction (x), the cMSD is given by σ 2

x (t) = [x(t)x(t) ] 2 where • denotes the average over all particles for a given realization. In Figure 6, we show the time evolution of the longitudinal cMSDs. The time axis is normalized with the mean travel time along the characteristic fracture link length, l. For both injection methods, spreading shows a ballistic regime (∼ t 2 ) at early times, which then transitions to a preasymptotic scaling in

A C C E P T E D M A N U S C R I P T
an intermediate regime and finally to a final asymptotic regime. The time evolutions of cMSDs for the two injection cases are notably different as conductivity heterogeneity increases, while the asymptotic late-time scalings are very similar.

The asymptotic power-law scaling can be understood in the framework of a continuous time random walk (CTRW) description of dispersion. At large times the Lagrangian velocity distributions are in their steady states and subsequent particle velocities are independent. Thus, at large times horizontal particle dispersion can be described by the CTRW

x n+1 = x n + 0 , t n+1 = t n + τ n , (12)
with the transition time τ n = 0 /v n . The velocities v n are distributed according to p s (v) which is space Lagrangian velocity PDF, and 0 is a distance along the streamline that is sufficiently large so that subsequent particle velocities may be considered independent. Thus, the distribution of transit times τ n is given in terms of the space Lagrangian and Eulerian velocity PDFs as [START_REF] Dentz | Con-758 tinuous time random walks for the evolution of Lagrangian velocities[END_REF]

ψ(τ ) = 0 τ 2 p s ( 0 /τ ) = 0 τ 3 v p e ( 0 /τ ), ( 13
)
where v is the average Eulerian velocity, see also Section 4. Specifically, for the scaling p e (v) ∝ v α at small velocities, the transit time PDF scales as

ψ(τ ) ∝ τ -1-β , β = 2 + α. ( 14
)
From Figure 3, we estimate for σ ln K = 5 that α ≈ -0.55, which corresponds to β = 1.45. CTRW theory [66, 65] predicts that the cMSD scales as t 3-β , which here implies t 1.55 . This is consistent with the late-time scaling of the cMSD shown in 6 for σ ln K = 5. In summary, both the increase in conductivity heterogeneity and the uniform injection method enhance longitudinal spreading. For low heterogeneity, the two different injection rules do not affect particle spreading significantly.

A C C E P T E

The difference, however, becomes significant as the conductivity heterogeneity increases. Both the magnitude of the cMSD and the super-diffusive scaling behavior are notably different for the two different injection rules at high heterogeneity. We now analyze the Lagrangian particle statistics to understand the underlying physical mechanisms that lead to the observed anomalous particle spreading.

Lagrangian Velocity Statistics

The classical CTRW approach [64, 65]-see Eq. ( 12)-relies on the independence of particle velocities at subsequent steps and thus spatial po- Particle motion is described here by the recursion relations [START_REF] Tsang | Tracer transport in a stochastic continuum model of fractured media[END_REF]. In this framework, we consider two types of Lagrangian velocities. The t(ime)-Lagrangian velocities are measured at a given time t,

v t (t) = u nt , (15)
where n t is defined by [START_REF] Kiraly | Remarques sur la simulation des failles et du réseau karstique par éléments finis dans les modèles d'écoulement[END_REF]. The s(pace)-Lagrangian velocities are measured at a given distance s along the trajectory. The distance s n traveled by a particle along a trajectory after n steps is given by

s n+1 = s n + n . ( 16
)
The number of steps needed to cover the distance s is described by n s = sup(n|s n ≤ s). Thus, the particle velocity at a distance s along a trajectory is given by

v s (s) = u ns . (17)
The PDF of t-Lagrangian velocities sampled along a particle path is given by

p t (v) = lim n→∞ n i=1 τ i δ(v -u i ) n i=1 τ i , (18)

A C C E P T E D M A N U S C R I P T

where we defined the transit time,

τ i = i u i (19)
The PDF of s-Lagrangian velocities sampled along a particle path are defined analogously as

p s (v) = lim n→∞ n i=1 i δ(v -u i ) n i=1 i . ( 20
)
Note the difference with respect to Eq. ( 2), which samples velocities in the network uniformly, while in Eq. 20 velocities are sampled along trajectories.

Using the definition of the transit time τ i in [START_REF] Makedonska | 552 Particle tracking approach for transport in three-dimensional discrete 553 fracture networks[END_REF], the PDFs of the s-and t-Lagrangian velocities are related through flux weighting as [START_REF] Dentz | Con-758 tinuous time random walks for the evolution of Lagrangian velocities[END_REF] p s (v) = vp t (v)

dv vp t (v) . (21)
Furthermore, for flux-preserving flows and under ergodic conditions, the Eulerian and t-Lagrangian velocity PDFs are equal,

p e (v) = p t (v). (22)
Thus, under these conditions, the s-Lagrangian and Eulerian velocity PDFs are related as [START_REF] Dentz | Con-758 tinuous time random walks for the evolution of Lagrangian velocities[END_REF] p

s (v) = vp e (v)
dv vp e (v)

. ( 23
)
This means that the stationary s-Lagrangian velocity PDF can be determined 272 from the Eulerian velocity PDF. It is important to emphasize that the above definitions of the Lagrangian velocity PDFs refer to stationary conditions. We now define the PDFs of tand s-Lagrangian velocities through sampling between particles and network realizations at a given time (t-Lagrangian) or space (s-Lagrangian) velocities

pt (v, t) = δ[v -v(t)] , ps (v, s) = δ[v -v(s)] . (24)
In general, these quantities evolve in time and with distance along the streamline and are sensitive to the injection conditions because evidently for t = 0 and s = 0 both are equal to the PDF of initial particle velocities pt (v, t = 0) = ps (v, s = 0) = p 0 (v), but their respective stationary PDFs are different, namely

p t (v) = lim t→∞ pt (v, t), p s (v) = lim s→∞ ps (v, s). (25)
Let us consider some further consequences of these properties. First, we notice that under (Eulerian) ergodicity the uniform injection condition (4) corresponds to an initial velocity PDF of

p 0 (v) = p e (v) = p t (v), (26)
that is, the initial velocity PDF is equal to the Eulerian and thus t-Lagrangian 277 velocity PDFs. This means that for the uniform injection method, the t-

278

Lagrangian velocity PDF is steady, pt (v, t) = p t (v), while the s-Lagrangian 279 velocity PDF is not. It evolves from its initial distribution p s (v, s = 0) = 280 p e (v) to the steady state distribution (23).

281

A C C E P T E D M A N U S C R I P T

The flux-weighted injection condition, on the other hand, corresponds to the the initial velocity PDF

p 0 (v) = p s (v), (27)
due to relation [START_REF] Bernabe | Passive advection-dispersion in 559 networks of pipes: Effect of connectivity and relationship to perme-560 ability[END_REF]. The initial velocity PDF is equal to the s-Lagrangian velocity PDF. This means that under flux-weighting, the s-velocity PDF is steady, ps (v, s) = p s (v). Under these conditions, the t-Lagrangian velocity PDF pt (v, t) evolves from the initial distribution pt (v, t = 0) = p s (v) towards the asymptotic p t (v) = p e (v), which is equal to the Eulerian velocity PDF.

These are key insights for the qualitative and quantitative understanding of the average transport behavior.

Space-Lagrangian Velocity Statistics

We analyze particle velocities along their projected trajectories in the longitudinal direction. Spatial particle transitions may be characterized by the characteristic fracture link length, l. The Lagrangian velocity v s (s n ) at a distance x n = n l along the projected trajectory is approximated by the average velocity v n ≡ l/τ n where τ n is the transition time for the distance l at step n. In the following, we investigate the statistical characteristics of the s-Lagrangian velocity series {v n }. For the uniform flow conditions under consideration here, the projected distance x n is a measure for the streamwise distance s n , and v n for the s-Lagrangian velocity v s (s n ). Spatial Lagrangian

A C C E P T E D M A N U S C R I P T

We first study the convergence of the s-Lagrangian velocity PDFs towards a stationary distribution and the invariance of ps (v, s) for a stationary (fluxweighted) initial velocity PDF. We consider the two injection conditions ( 4) and ( 5) and record the distribution of particle velocities at a line located at the control point x c . Under ergodic conditions, we expect ps (v, s) to converge towards its steady state distribution (23) for uniform injection and to remain invariant for the flux-weighted injection. Figure 8 shows ps (v, s = 0) and ps (v, s = x c ) for uniform and flux-weighted injection conditions and two different heterogeneity strengths. We clearly observe that ps (v,

x c ) = p s (v)
is invariant for flux-weighted injection. For uniform injection, ps (v, x c ) has already evolved towards its steady limit after x c = 200 l. This is an indication that the flow and transport system is in fact ergodic. Note that in terms of computational efficiency, this observation gives a statistically consistent way of continuing particle trajectories through reinjection at the inlet. If the outlet is located at a position x c large enough so that ps (v, s = x c ) = p s (v), particles are reinjected at the inlet with flux-weighted probability, this means that the velocity statistics are preserved. Furthermore, this method ensures that the domain is large enough to provide ergodic conditions. In the following, we analyze the statistical properties of streamwise velocity transitions with the aim of casting these dynamics in the frame of a Markov model for subsequent particle velocities.

We first consider the distribution ψ τ (t) of transition times along particle trajectories through sampling the transition times along all particle trajectories and among network realizations. To this end, we consider a flux-weighted injection because it guarantees that the s-Lagranagian velocities are station-

A C C E P T E D M

A N U S C R I P T ary. Figure 9 illustrates the PDF of transition times for different variances of ln(K). As σ ln K increases, the transition time PDFs become broader. The transition time closely follows a truncated power-law distribution.

Next we consider two-point velocity statistics to gain insight into the velocity correlations along a streamline. To this end, we consider the velocity auto covariance for a given lag ∆s = ss . As pointed out above, for fluxweighted injection, the streamwise velocities here are stationary and therefore

C s (s -s ) = [v s (s) -v s (s) ][v s (s ) -v s (s ) ] . (28)
In order to increase the statistics, we furthermore sample along streamlines over a distance of 10 2 l. The velocity variance is

σ 2 v = C s (0). The velocity autocorrelation function χ s (s) = C s (s)/σ 2 v . The correlation length scale c is defined by c = ∞ 0 ds χ s (s). ( 29
)
The inset in Figure 9 shows the increase in the velocity correlation length scale with increasing ln(K) variances for a flux-weighted injection case. This can be attributed to the emergence of preferential flow paths, as shown in

σ lnK =1 σ lnK =2 σ lnK =3 σ lnK =5 ψ τ ( t) • t l t/ t l 1 1 2.45 2.85

A C C E P T E D M A N U S C R I P T

study this feature, we characterize the series of s-Lagrangian velocities {v n } in terms of the transition probabilities to go from velocity v m to velocity v m+n . We determine the transition probabilities under flux-weighted particle injection because, as detailed above, under these conditions, the s-Lagrangian velocity is stationary. Thus, the transition probability is only a function of the number n of steps,

r n (v|v ) = δ(v -v m+n ) | vm=v (30)
Numerically, the transition probability is determined by discretizing the s-Lagrangian velocity PDF into N velocity classes C i = (v s,i , v s,i + ∆v s,i ) and recording the probability for each class given the previous velocity class. This procedure gives the transition matrix T n (i|j) from class j to i after n steps such that r n (v|v ) is approximated numerically as

r n (v|v ) = N i,j=1 I C i (v)T n (i|j)I C j (v ) ∆v i , (31)
where the indicator function

I C i (v) is 1 if v ∈ C i

A C C E P T E D M

A N U S C R I P T correlation scale (see inset of figure 9). This leads to a spatial Markov-chain random walk (which we also termed spatial Markov Model) formulation of particle dispersion that is valid for any initial velocity distribution, and thus for any injection protocol. This modelling approach is in line with the timedomain random walk (TDRW) and continuous time random walk (CTRW) approaches discussed in the Introduction and below.

Markovian Velocity Process

Along the lines of Le Borgne et al. [74] and Kang et al. [START_REF] Kang | Spatial Markov model 725 of anomalous transport through random lattice networks[END_REF], we model the velocity series {v n } as a Markov-chain, which is a suitable model to statistically quantify the evolution of the s-Lagrangian velocities based on the existence of a finite correlation length. In this framework, the n-step transition probability r n (v|v ) satisfies the Chapman-Kolmogorov equation [START_REF] Risken | The Fokker-Planck Equation[END_REF]

r n (v n |v 0 ) = dv k r n-k (v n |v k )r k (v k |v 0 ). ( 32a
)
The velocity process is fully characterized in terms of the one-step transition PDF r 1 (v|v ) and the steady state PDF p s (v) of the s-Lagrangian velocity.

Consequently, the evolution of the s-Lagrangian velocity PDF ps (v, s n ) is given by

p s (v, s n ) = dv r 1 (v|v )p s (v , s n ), (32b)
with the arbitrary initial PDF p s (v, s 0 = 0) = p 0 (v). The number of steps to decorrelate this Markov-chain is given by n c = c / l. Figure 11 shows the evolution of the PDF of s-Lagrangian velocities for the uniform injection (4).

Recall that the uniform injection mode corresponds to the initial velocity

A C C E P T E D M A N U S C R I P T
PDF p e (v). Thus, the numerical Monte Carlo simulations are compared to the predictions of (32b) for the initial condition ps (v, s 0 = 0) = p e (v). The transition PDF r 1 (v|v ) is given by [START_REF] Arbogast | Derivation of the double porosity 584 model of single phase flow via homogenization theory[END_REF] with the velocity transition matrix shown in Figure 10. As shown in Figure 11, the prediction of the Markovian velocity model and the Monte Carlo simulation are in excellent agreement, which confirms the validity of the Markov model (32) for the evolution of s-Lagrangian velocities. Velocity transition dynamics are independent of the particular initial conditions and thus allow predicting the evolution of the Lagrangian velocity statistics for any initial velocity PDF and thus for any injection protocol.

As mentioned above, the Markov-chain {v n } is fully characterized by the stationary PDF of the s-Lagrangian velocities and the transition PDF r 1 (v|v ). The behavior of the latter may be characterized by the number of steps n c needed to decorrelate, i.e., the number of steps n c such that r n (v|v )

for n > n c converges to the stationary PDF r n (v|v ) → p s (s). The number of steps for velocities to decorrelate can be quantified by

n c = c l , (33)
The simplest transition PDF that shares these characteristics is [41, 80, 50, 88]

r 1 (v|v ) = aδ(v -v ) + (1 -a)p s (v), (34)
with a = exp(-l/ c ). This transition PDF is thus fully determined by one single parameter n c . Note that the latter increases with the level of heterogeneity, as illustrated in the inset of figure 9. This parameter is estimated here from the simulated Lagrangian velocities. It may also be measured in We consider particle motion along the mean pressure gradient in xdirection, which is described by the stochastic regression

x n+1 = x n + l, t n+1 = t n + l v n . ( 35
)
The velocity transitions are determined from the Markovian velocity process (32). Note that the {v n } process describes equidistant velocity transitions along particle trajectories, while (35) describes particle motion projected on the x-axis. In this sense, (35) approximates the longitudinal travel distance x n with the distance s n along the streamline, which is valid if the tortuosity of the particle trajectories is low. As indicated in Section 3.2, for travel distances 0 larger than c , or equivalently, step numbers n n c ≡ c / l, subsequent velocities may be considered independent and particle dispersion is fully characterized by the recursion relation [START_REF] Juanes | A general and efficient formulation 529 of fractures and boundary conditions in the finite element method[END_REF] and the transition time PDF [START_REF] Park | Transport be-532 havior in three-dimensional fracture intersections[END_REF]. Thus, as shown in Section 3.2, the CTRW of Eq. ( 12) correctly predicts the asymptotic scaling behavior of the centered mean square displacement. This is not necessarily so for the particle breakthrough and the preasymptotic behavior of the cMSD. As seen in Section 3.1, the late time tailing of the BTC depends on the injection mode and thus on the initial velocity PDF. In fact, the slope observed in Figure 5 for uniform particle injection can be understood through the persistence of the initial velocity PDF.

The first random walk steps until decorrelation at n = n c are characterized

Thus, for an initial velocity PDF p 0 (v) = p e (v), the initial transit time PDF is given in terms of the Eulerian velocity PDF, which is characterized by a stronger probability weight towards low velocities than the PDF of the s-Lagrangian velocities, which is given by (23). The space-time random walk (35) together with the Markov model (32b) is very similar to the TDRW approach [70, 71] and can also be seen as a multi-state, or correlated CTRW approach because subsequent particle velocities and thus transition times are represented by a Markov process [START_REF] Haus | Diffusion in regular and disordered lattices[END_REF][START_REF] Montero | Nonindependent continuoustime random walks[END_REF]74,[START_REF] Chechkin | Continuous-time random walk with correlated waiting times[END_REF]. The joint distribution p(x, v, t) of particle position and velocity at a given time t is given by [74]

p(x, v, t) = t 0 dt H( l/v -t )R(x -vt , v, t -t ), (37)
where H(t) is the Heaviside step function; R(x, v, t) is the frequency by which a particle arrives at the phase space position (x, v, t). It satisfies

R(x, v, t) = R 0 (x, v, t) + dv r 1 (v|v )R(x -l, v , t -l/v ), (38)
where R 0 (x, v, t) = p 0 (x, v)δ(t) with p 0 (x, v) = p(x, v, t = 0). Thus, the right side of (37) denotes the probability that a particle arrives at a position x-vt where it assumes the velocity v by which it advances toward the sampling position x. Equations ( 37) and ( 38) can be combined into the Boltzmann equation In fact, as the initial velocity shows the same behavior at low velocities as the Eulerian velocity PDF, the breakthrough tailing is the same as observed in Figure [START_REF] Juanes | A general and efficient formulation 529 of fractures and boundary conditions in the finite element method[END_REF]. We also observed that the reduced-order Markov velocity model can capture important features of the time evolution of cMSDs. This demonstrates that the proposed model can incorporate arbitrary injection modes into the effective modelling framework.

∂p(x, v, t) ∂t = -v ∂p(x, v, t) ∂x - v l p(x, v, t) + dv v l r 1 (v|v )p(x, v , t), (39

Conclusions

This study shows how the interplay between fracture geometrical properties (conductivity distribution and network geometry) and tracer injection modes controls average particle transport via Lagrangian velocity statistics.

The interplay between fracture heterogeneity and tracer injection methods can lead to distinctive anomalous transport behavior. Furthermore, the injection conditions, for example, uniform or flux-weighted, imply different initial velocity distributions, which can have a persistent impact on particle spreading through DFNs. For uniform injection, the s-Lagrangian velocity distribution evolves from an Eulerian velocity distribution initially to a stationary s-Lagrangian distribution. In contrast, for flux-weighted injection, the s-Lagrangian velocity distribution remains stationary. In this study, we investigated the particle transport and the impact of the injection condition for idealized 2D DFN using a Markov velocity model.

•Introduction 1 Flow 3 A 5 2

 135 Injection modes have major impact on anomalous transport in DFNs • Evolution of the Lagrangian velocity distribution is governed by injection modes • Spatial velocity Markov model for variable injection modes • Equivalence between spatial Markov model and Boltzmann equation 1. and transport in fractured geologic media control many impor-2 tant natural and engineered processes, including nuclear waste disposal, ge-C C E P T E D M A N U S C R I P T ologic carbon sequestration, groundwater contamination, managed aquifer 4 recharge, and geothermal production in fractured geologic media [e.g., 1, , 3, 4, 5]. Two dominant approaches exist for simulating flow and trans-6

17

 17

  DFN modelling has advanced significantly in recent years with the in-18 crease in computational power. Current DFN simulators can take into ac-19 count multiple physical mechanisms occurring in complex 3D fracture sys-20 tems. Recent studies also have developed methods to explicitly model ad-21 vection and diffusion through both the discrete fractures and the permeable 22

25 continued

 25 interest in the development of upscaled transport models that can 26 be parameterized with a small number of model parameters. Ideally, these 27 model parameters should have a clear physical interpretation and should be 28 A C C E P T E D M A N U S C R I P T determined by means of field experiments, with the expectation that the 29 model can then be used for predictive purposes [40, 41].

30

 30

39

  damentally different behavior compared with Fickian transport [59, 60, 61]. 40 The continuous time random walk (CTRW) formalism [62, 63] is a frame-41 work to describe anomalous transport through which models particle motion 42 through a random walk in space and time characterized by random space 43 and time increments, which accounts for variable mass transfer rates due 44 to spatial heterogeneity. It has been used to model transport in heteroge-45 neous porous and fractured media [64, 65, 66, 67, 68, 49, 69] and allows 46 incorporating information on flow heterogeneity and medium geometry for 47 large scale transport modelling. Similarly, the time-domain random walk 48 (TDRW) approach [70, 71, 72] models particle motion due to distributed 49 space and time increments, which are derived from particle velocities and 50 their correlations. The analysis of particle motion in heterogeneous flow 51 fields demonstrate that Lagrangian particle velocities exhibit sustained cor-52 relation along their trajectory [73, 70, 74, 75, 76, 77, 78, 79, 45]. Volume 53 A C C E P T E D M A N U S C R I P T conservation induces correlation in the Eulerian velocity field because fluxes must satisfy the divergence-free constraint. This, in turn, induces correlation in the Lagrangian velocity along a particle trajectory. To take into account velocity correlation, Lagrangian models based on Markovian processes have been proposed [70, 71, 74, 76, 77, 78, 45, 41, 80]. Spatial Markov models are based on the observation that successive velocity transitions measured equidistantly along the mean flow direction exhibit Markovianity: a particle's velocity at the next step is fully determined by its current velocity. The spatial Markov model, which accounts for velocity correlation by incorporating this one-step velocity correlation information, has not yet been extended to disordered (unstructured) DFNs.

Figure 1 .

 1 Figure 1. The position vector of node i in the fracture network is denoted by x i .

  study a uniform flow setting characterized by constant mean flow in the positive x-direction parallel to the principal set of factures. No-flow conditions are imposed at the top and bottom boundaries of the domain, and fixed hydraulic head at the left (Φ = 1) and right (Φ = 0) boundaries. The overbar in the following denotes the ensemble average over all network realizations. The one-point statistics of the flow field are characterized by the Eulerian velocity PDF, which is obtained by spatial and ensemble sampling A C C E P T E D M A N U S C R I P T of the velocity magnitudes in the network

  of a passive tracer by particle tracking. Particles are injected along a line

Figure 1 : 2 Figure 2 :

 122 Figure 1: (a) Example of a two-dimensional DFN studied here, with 2000 fractures (1000 fractures for each fracture set). (b) Subsection of a spatially uncorrelated conductivity field between 0 ≤ x ≤ 1 and 0.5 ≤ y ≤ 1. Conductivity values are assigned from a lognormal distribution with σ lnK = 1. Link width is proportional to the conductivity value; only connected links are shown.

Figure 3 :

 3 Figure 3: Eulerian flux probability density functions for four different levels of conductivity heterogeneity. Increase in conductivity heterogeneity significantly increases the probability of small flux values.

189 shown in Figure 4 ,

 4 both network heterogeneity and injection rule have signifi-190 cant impact on particle spreading. An increase in network heterogeneity leads 191 to an increase in longitudinal particle spreading, and the uniform injection 192 rule significantly enhances longitudinal spreading compared to flux-weighted 193 A C C E P T E D M A N U S C R I P T

Figure 4 :

 4 Figure 4: Particle distribution at t = 20t l for a given realization after the instantaneous release of particles at the inlet, x = 0. t l is the median transition time to travel L x /100. (a) The low heterogeneity case (σ lnK = 1) with flux-weighted injection. (b) The low heterogeneity case (σ lnK = 1) with uniform injection. (c) The projected particle distribution in the longitudinal direction for the low heterogeneity case (σ lnK = 1). (d) The high heterogeneity case (σ lnK = 5) with the flux-weighted injection. (e) The high heterogeneity (σ lnK = 5) with the uniform injection.(f) The projected particle distribution in the longitudinal direction for the high heterogeneity case (σ lnK = 5). For the high heterogeneity case, the injection method has significant impact on particle spreading. The uniform injection method leads to more anomalous spreading.

Figure 5 :Figure 5 shows

 55 Figure 5: (a) FPTDs for σ ln K = 1, 2, 3, 5 with flux-weighted injection at x c = 200 l. Increase in conductivity heterogeneity leads to larger dispersion and stronger late-time tailing. (b) FPTDs for σ ln K = 1, 2, 3, 5 with uniform injection (black solid lines). Uniform injection leads to significantly larger dispersion and late-time tailing compared to the flux-weighted injection (red dashed lines). FPTDs are normalized with the peak arrival time.

  S C R I P T The Monte Carlo simulations show that, in the intermediate regime (t/t l approximately between 1 and 100), the longitudinal cMSD increases linearly with time for flux-weighted injection [Figure 6(a)]. For uniform injection, cMSD increases faster than linearly (i.e., superdiffusively) for intermediate to strong heterogeneity in the intermediate regime [Figure 6(b)]. The stronger heterogeneity led to the increase in the late-time temporal scaling for both flux-weighted and uniform injection cases. The Monte Carlo simulations also show that there is no noticeable difference between the uniform injection and the flux-weighted injection for the low heterogeneity case whereas the difference increases as heterogeneity increases [Figure 6(b), inset].

5 Figure 6 :

 56 Figure 6: Time evolution of longitudinal MSDs for σ ln K = 1, 2, 3, 5 obtained from Monte Carlo simulations (solid lines), and the model predictions from the Markov-Chain CTRW (32) and (35) with the full transition PDF (dashed lines). Increase in conductivity heterogeneity leads to higher dispersion, and the Markov-chain CTRW model is able to accurately capture the time evolution of the MSDs for all levels of heterogeneity and injection rules. (a) Flux-weighted injection. (b) Uniform injection. Inset: Comparison between flux-weighted and uniform injection for σ ln K = 1, 5. Impact of injection rule is significant for high conductivity heterogeneity.

Figure 7

 7 illustrates this relation by compar-273 ing the s-Lagrangian velocity PDFs measured from the numerical simulation 274 to the flux-weighted Eulerian velocity PDFs shown in Figure 3.

Figure 7 : 4 . 2 .

 742 Figure 7: s-Lagrangian velocity PDFs for σ lnK = 1 and σ lnK = 5. The measured s-Lagrangian velocity PDF agrees very well with the PDF obtained by transforming the Eulerian velocity PDF using Eq. (23).

Figure 8 :

 8 Figure 8: Lagrangian flux distributions at the inlet and outlet for uniform and fluxweighted injection rules, for σ lnK = 5. Note that Lagrangian flux distributions at the outlet are identical regardless of the injection method. Inset: the same plot for σ lnK = 1.

Figure 2 .

 2 Figure 2. Painter and Cvetkovic [71] and Frampton and Cvetkovic [16] also reported the existence of clear velocity correlation between successive jumps in DFNs and showed that this correlation structure should be captured for effective transport modelling.The existence of a finite correlation length along the particle trajectories indicates that subsequent velocities, when sampled at a distance much larger than the correlation length c , may be considered independent. In order to

Figure 9 :

 9 Figure 9: (a) Lagrangian transition time distributions for σ ln K = 1, 2, 3, 5 with fluxweighted injection. As the network conductivity becomes more heterogeneous, the transition time distribution becomes broader. Inset: the effective correlation length increases with increasing network heterogeneity: 3.8 l, 5.6 l, 8.2 l, 11.5 l. The correlation step (n c ) is computed by integrating velocity autocorrelation function in space.

and 0 otherwise. 338 Figure 10 5 .A

 338105 Figure10shows the one-step transition matrix T 1 (i|j) for equidistant and

Figure 10 :

 10 Figure 10: (a) One-step velocity transition matrix T 1 (i|j) with linear equiprobable binning for N = 50 velocity classes for σ ln K = 1. (b) Velocity transition matrix with linear equiprobable binning for σ ln K = 5. The color-bar shows the logarithmic scale. (c) Velocity transition matrix with logarithmic binning for σ ln K = 1. (d) Velocity transition matrix with logarithmic binning for σ ln K = 5. Increase in conductivity heterogeneity leads to higher probability close to diagonal entries.

  S C R I P T the field from multiscale tracer tests [41]. In the following, we study particle 373 dispersion in the Markovian velocity model for the full transition PDF shown 374 in Figure 10 and the reduced-order Markov model (34).

375 5 . 2 .

 52 Particle Dispersion and Model Predictions 376

Figure 11 :

 11 Figure11: Evolution of the PDF of s-Lagrangian velocities for uniform injection, i.e., for an initial velocity PDF ps (v, s 0 = 0) = p e (v). The symbols denote the data obtained from the direct numerical simulation, the dashed lines show the predictions of (32b) with the transition matrix shown in Figure10. Inset: Evolution of the PDF of s-Lagrangian velocities for a flux-weighted injection. In this case, the initial velocity PDF is identical to the stationary s-Lagrangian velocity PDF.

  S C R I P Tsee Appendix A. This result provides a bridge between the TDRW approach[70, 71] and the correlated CTRW approach.As illustrated in Figure3, the Eulerian velocity PDF can be characterized by the power-law p e (v) ∝ v α . Thus, the first CTRW steps until the decorrelation at n = n c are characterized by the transit time PDF ψ 0 (t) ∝ t -2-α .The corresponding tail of the BTC is f (t, x c ) ∝ t -2-α . The observed value of α = -0.55 explains the tailing of the BTC in Figure5as t -1.45 , which shows the importance of the initial velocity distribution. We also observe decrease in BTC tailing (larger absolute slope) with travel distance as initial velocity distribution converges to stationary Lagrangian velocity distribution and as tracers sample more velocity values. This implies that one needs to be careful when inferring a β from single BTC measurement because the slope can evolve depending on the injection method, velocity PDF and velocity correlation.First, we compare the results obtained from Monte Carlo simulation in the random DFN to the predictions of the Markov-chain CTRW (32) and (35) with the full transition PDF of Figure10. Figure6shows the evolution of the cMSD for different levels of heterogeneity and different injection modes. As expected from the ability of the Markov model to reproduce the evolution of the s-Lagrangian velocity PDF for both uniform and fluxweighted injection conditions, the predictions of particle spreading are in excellent agreement with the direct numerical simulations. In Figure12we compare breakthrough curves obtained from numerical simulations with the predictions by the Markov model for the uniform and flux-weighted injection modes. Again, the impact of the injection mode and thus initial velocityA C C E P T E D M A N U S C R I P TPDF is fully quantified by the Markov model.We now apply the Markov model (32)-(35), i.e., employing a parsimonious parameterization of the velocity transition PDF, with a single parameter n c (equation (33)), which is estimated here from velocity correlations along streamlines (see inset of figure9). We first compare the reduced-order Markov model to the cases of uniform and flux-weighted injection, and conclude that the proposed parsimonious stochastic model provides an excellent agreement with the direct numerical simulations (Figure13). This implies that the simple correlation model (34) can successfully approximate the velocity correlation structure. Hence it appears that high order correlation properties, quantified from the full transition probabilities (figure10), are not needed for accurate transport predictions in the present case. This suggests promising perspective for deriving approximate analytical solutions for this Markov-chain CTRW model[START_REF] Dentz | Con-758 tinuous time random walks for the evolution of Lagrangian velocities[END_REF]. Furthermore, as discussed in [41], the velocity correlation parameter n c can be estimated in the field by combining cross-borehole and push-pull tracer experiments.Finally, we consider the evolution of the particle BTC and the cMSD for arbitrary injection modes. For real systems both flux-weighted and uniform injections are idealizations. A flux-weighted condition simulates a constant concentration pulse where the injected mass is proportional to the local injection flux at an inlet boundary that is extended over a distance much larger than the correlation scale during a given period of time. A uniform injection represents an initial concentration distribution that is uniformly extended over a region far larger than the correlation length. In general, the initial concentration distribution may not be uniform, and the injection boundaryA C C E P T E D M A N U S C R I P Tmay not be sufficiently large, which leads to an arbitrary initial velocity distribution, biased maybe to low or high flux zones, as for example in the MADE experiments, where the solute injection occurred into a low permeability zone[START_REF] Harvey | Rate-limited mass transfer or macrodispersion: Which dominates plume evolution at the macrodispersion experiment (MADE) site?[END_REF]. For demonstration, we study two scenarios representing injection into low and high flux zones: uniform injections into regions of the 20-percentile highest, and 20-percentile lowest velocities. The initial velocity PDF for the low velocity mode shows the power-law behavior which is the characteristic for the Eulerian PDF, and the initial velocity PDF for the high velocity mode shows narrow initial velocity distribution (Figure15). Eventually, the s-Lagrangian PDFs evolve towards the stationary flux-weighted Eulerian PDF as discussed in the previous section.

Figure 15 shows

 15 Figure 15 shows the predictive ability of the effective stochastic model for these different injection conditions. The reduced-order Markov velocity model compares well with the direct Monte Carlo simulation in the random networks. As expected, the BTCs for injection into low velocity regions have a much stronger tailing than for injection into high velocity regions.

Figure 12 :

 12 Figure 12: Particle BTCs from Monte Carlo simulations and the predictions from the Markov-chain CTRW model with the full velocity transition matrix for (a) flux-weighted injection, and (b) uniform injection at x c = 200 l.

Figure 13 :Figure 14 :

 1314 Figure 13: Particle BTCs from Monte Carlo simulations and predictions from the Markovchain CTRW model with the reduced-order velocity transition matrix for (a) flux-weighted injection, and (b) uniform injection at x c = 200 l.

Figure 15 :

 15 Figure 15: (a) Particle BTCs from Monte Carlo simulations for injection into high-velocity regions (solid line) and predictions from Markov-chain CTRW model with the reducedorder velocity transition matrix (dashed line). (b) Corresponding results for injection into low-velocity regions.

  S C R I P T We have presented a spatial Markov model to quantify anomalous transport through DFNs under arbitrary injection conditions. We derive an analytical relation between the stationary Lagrangian and the Eulerian velocity distribution, and formally incorporate the initial velocity distribution into the spatial Markov model. The proposed model accurately reproduces the evolution of the Lagrangian velocity distribution for arbitrary injection modes. This is accomplished with a reduced-order stochastic relaxation model that captures the velocity transition with a single parameter: the effective velocity correlation c . The agreement between model predictions and direct numerical simulations indicates that the simple velocity correlation model can capture the dominant velocity correlation structure in DFNs.

  These findings can be extended to 3D DFNs, for which similar behaviors regarding the injection mode have been found[58]. Also, Frampton and Cvetkovic [16] reported similar velocity correlation structures for 3D DFNs as in 2D, which suggests that a velocity Markov model such as the one presented in this work can be used for the modelling of particle motion in 3D DFNs.A C C E P T E D MA N U S C R I P T

Acknowledgements: PKK and SL acknowledge a grant (16AWMP-B066761-04) from the AWMP Program funded by the Ministry of Land, Infrastructure and Transport of the Korean government and the support from Future Research Program (2E27030) funded by the Korea Institute of Science and Technology (KIST). PKK and RJ acknowledge a MISTI Global Seed Funds award. MD acknowledges the support of the European Research

Appendix A. Boltzmann Equation

The time derivative of (37) gives

Note that R(x, v, t) denotes the probability per time that a particle has the velocity v at the position x. It varies on a time scale of ¯v.