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

TL;DR

This paper introduces a spatial Markov model that captures the evolution of tracer particle velocities in disordered fracture networks, accounting for different injection modes and predicting anomalous transport behavior.

Contribution

It develops a novel Markov-based stochastic model that incorporates injection mode effects into transport predictions in fractured media.

Findings

Model accurately predicts anomalous transport in heterogeneous fracture networks.

Incorporates initial injection mode effects into velocity evolution.

Successfully matches observed transport behaviors across various heterogeneity levels.

Abstract

We investigate tracer transport on random discrete fracture networks that are characterized by the statistics of the fracture geometry and hydraulic conductivity. While it is well known that tracer transport through fractured media can be anomalous and particle injection modes can have major impact on dispersion, the incorporation of injection modes into effective transport modelling has remained an open issue. The fundamental reason behind this challenge is that---even if the Eulerian fluid velocity is steady---the Lagrangian velocity distribution experienced by tracer particles evolves with time from its initial distribution, which is dictated by the injection mode, to a stationary velocity distribution. We quantify this evolution by a Markov model for particle velocities that are equidistantly sampled along trajectories. This stochastic approach allows for the systematic incorporation of the initial velocity distribution and quantifies the interplay between velocity distribution and spatial and temporal correlation. The proposed spatial Markov model is characterized by the initial velocity distribution, which is determined by the particle injection mode, the stationary Lagrangian velocity distribution, which is derived from the Eulerian velocity distribution, and the spatial velocity correlation length, which is related to the characteristic fracture length. This effective model leads to a time-domain random walk for the evolution of particle positions and velocities, whose joint distribution follows a Boltzmann equation. Finally, we demonstrate that the proposed model can successfully predict anomalous transport through discrete fracture networks with different levels of heterogeneity and arbitrary tracer injection modes.

Equations

$$p_l(l)=\frac{\exp (-l/\bar{l})}{\bar{l}}.$$

$$p_e(u)=\frac{{\sum }_{i>j}\overline{l_{ij}\delta (u-u_{ij})}}{N_{\ell }\bar{l}}.$$

$$p_e(u)=\frac{1}{N_{\ell }}\sum _{i>j}\overline{\delta (u-u_{ij})}.$$

$$N_{i_0}=\frac{N_p}{{\sum }_{i_0}},$$

$$N_{i_0}=N_p\frac{Q_{i_0}}{{\sum }_{i_0}{Q}_{i_0}}.$$

$$p_{ij}=\frac{|u_{ij}|}{{\sum }_k|u_{ik}|},$$

$${\mathbf{x}}_{n+1}$$

$$t_{n+1}$$

$$n_t=\sup (n|t_n\le t)$$

$$p(\mathbf{x},t)=⟨\delta (\mathbf{x}-{\mathbf{x}}_{n_t})⟩,$$

$$\overline{P}(\mathbf{x},t)=\overline{⟨\delta (\mathbf{x}-{\mathbf{x}}_{n_t})⟩},$$

$$f(\tau ,x_c)=\overline{⟨\delta [\tau -{\tau }_a(x_c)]⟩}.$$

$$x_{n+1}=x_n+{\ell }_0,$$

$$\psi (\tau )=\frac{{\ell }_0}{{\tau }^2}{p}_s({\ell }_0/\tau )=\frac{{\ell }_0}{{\tau }^3\overline{v}}{p}_e({\ell }_0/\tau ),$$

$$\psi (\tau )\propto {\tau }^{-1-\beta },$$

$$v_t(t)=u_{n_t},$$

$$s_{n+1}=s_n+{\ell }_n.$$

$$v_s(s)=u_{n_s}.$$

$$p_t(v)=\underset{n\to \infty }{\lim }\frac{{\sum }_{i=1}^n{\tau }_i\delta (v-u_i)}{{\sum }_{i=1}^n{\tau }_i},$$

$${\tau }_i=\frac{{\ell }_i}{u_i}$$

$$p_s(v)=\underset{n\to \infty }{\lim }\frac{{\sum }_{i=1}^n{\ell }_i\delta (v-u_i)}{{\sum }_{i=1}^n{\ell }_i}.$$

$$p_s(v)=\frac{v{p}_t(v)}{\int dvv{p}_t(v)}.$$

$$p_e(v)=p_t(v).$$

$$p_s(v)=\frac{v{p}_e(v)}{\int dvv{p}_e(v)}.$$

$${\hat{p}}_t(v,t)=\overline{⟨\delta [v-v(t)]⟩},$$

$$p_t(v)$$

$$p_0(v)=p_e(v)=p_t(v),$$

$$p_0(v)=p_s(v),$$

$$C_s(s-s^{\prime })=\overline{⟨[v_s(s)-⟨v_s(s)⟩][v_s(s^{\prime })-⟨v_s(s^{\prime })⟩]⟩}.$$

$${\ell }_c=\underset{0}{\overset{\infty }{\int }}ds{\chi }_s(s).$$

$$r_n(v|v^{\prime })=⟨\delta (v-v_{m+n})⟩|{}_{v_m=v^{\prime }}$$

Full text

1 Introduction

Flow and transport in fractured geologic media control many important natural and engineered processes, including nuclear waste disposal, geologic carbon sequestration, groundwater contamination, managed aquifer recharge, and geothermal production in fractured geologic media [e.g., 1, 2, 3, 4, 5]. Two dominant approaches exist for simulating flow and transport through fractured media: the equivalent porous medium approach [6, 7] and the discrete fracture network approach (DFN) [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. The DFN approach explicitly resolves individual fractures whereas the equivalent porous medium approach represents the fractured medium as a single continuum by deriving effective parameters to include the effect of the fractures on the flow and transport. The latter, however, is hampered by the fact that a representative elementary volume may not exist for fractured media [22, 23]. Dual-porosity models are in between these two approaches, and conceptualize the fractured-porous medium as two overlapping continua, which interact via an exchange term [24, 25, 26, 27, 28, 29, 30, 31, 32].

DFN modelling has advanced significantly in recent years with the increase in computational power. Current DFN simulators can take into account multiple physical mechanisms occurring in complex 3D fracture systems. Recent studies also have developed methods to explicitly model advection and diffusion through both the discrete fractures and the permeable rock matrix [33, 34, 35, 36]. In practice, however, their application must account for the uncertainty in the subsurface characterization of fractured media, which is still an considerable challenge [37, 38, 39]. Thus, there is a continued interest in the development of upscaled transport models that can be parameterized with a small number of model parameters. Ideally, these model parameters should have a clear physical interpretation and should be determined by means of field experiments, with the expectation that the model can then be used for predictive purposes [40, 41].

Developing an upscaled model for transport in fractured media is especially challenging due to the emergence of anomalous (non-Fickian) transport. While particle spreading is often described using a Fickian framework, anomalous transport—characterized by scale-dependent spreading, early arrivals, long tails, and nonlinear scaling with time of the centered mean square displacement—has been widely observed in porous and fractured media across multiple scales, from pore [42, 43, 44, 45, 46] to single fracture [47, 48, 49, 50] to column [51, 52] to field scale [53, 54, 55, 56, 57, 41, 58]. The ability to predict anomalous transport is essential because it leads to fundamentally different behavior compared with Fickian transport [59, 60, 61].

The continuous time random walk (CTRW) formalism [62, 63] is a framework to describe anomalous transport through which models particle motion through a random walk in space and time characterized by random space and time increments, which accounts for variable mass transfer rates due to spatial heterogeneity. It has been used to model transport in heterogeneous porous and fractured media [64, 65, 66, 67, 68, 49, 69] and allows incorporating information on flow heterogeneity and medium geometry for large scale transport modelling. Similarly, the time-domain random walk (TDRW) approach [70, 71, 72] models particle motion due to distributed space and time increments, which are derived from particle velocities and their correlations. The analysis of particle motion in heterogeneous flow fields demonstrate that Lagrangian particle velocities exhibit sustained correlation along their trajectory [73, 70, 74, 75, 76, 77, 78, 79, 45]. Volume 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.

The mode of particle injection can have a major impact on transport through porous and fractured media [81, 82, 83, 84, 85, 86, 58, 87]. Two generic injection modes are uniform (resident) injection and flux-weighted injection with distinctive physical meanings as discussed in Frampton and Cvetkovic [84]. The work by Sposito and Dagan [82] 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 [84, 58]. Dagan [87] 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 [84]. The fundamental challenge is that the Lagrangian velocity distribution experienced by tracer particles evolves with time from its initial distribution which is dictated by the injection mode to a stationary velocity distribution [73, 84, 75, 88]. 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 emergence of anomalous transport by direct Monte Carlo simulations of flow and particle transport. In Section 4, we analyze Eulerian and Lagrangian velocity statistics to gain insight into the effective particle dynamics and elucidate the key mechanisms that lead to the observed anomalous behavior. In Section 5, we develop a spatial Markov model that is characterized by the initial velocity distribution, probability density function (PDF) of Lagrangian velocities and their transition PDF, which are derived from the Monte Carlo simulations. The proposed model is in excellent agreement with direct Monte Carlo simulations. We then present a parsimonious spatial Markov model that quantifies velocity correlation with a single parameter. The predictive capabilities of this simplified model are demonstrated by comparison to the direct Monte Carlo simulations with arbitrary injection modes. In Section 6, we summarize the main findings and conclusions.

2 Flow and Transport in Discrete Random Fracture Networks

2.1 Random Fracture Networks

We numerically generate random DFNs in two-dimensional rectangular regions, and solve for flow and tracer transport within these networks. The fracture networks are composed of linear fractures embedded in an impermeable rock matrix. The idealized 2D DFN realizations are generated by superimposing two different sets of fractures, which leads to realistic discrete fracture networks [89, 67]. Fracture locations, orientations, lengths and hydraulic conductivities are generated from predefined distributions, which are assumed to be statistically independent: (1) Fracture midpoints are selected randomly over the domain size of $L_{x}\times L_{y}$ where $L_{x}=2$ and $L_{y}=1$; (2) Fracture orientations for two fracture sets are selected randomly from Gaussian distributions, with means and standard deviation of $0^{\circ}\pm 5^{\circ}$ for the first set, and $90^{\circ}\pm 5^{\circ}$ for the second set; (3) Fracture lengths are chosen randomly from exponential distributions with mean $L_{x}/10$ for the horizontal fracture set and mean $L_{y}/10$ for the vertical fracture set; (4) Fracture conductivities are assigned randomly from a predefined log-normal distribution. An example of a random discrete fracture network with $2000$ fractures is shown in Figure 1.

The position vector of node $i$ in the fracture network is denoted by $\mathbf{x}_{i}$. 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

[TABLE]

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 $\bar{l}\approx 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 uncorrelated. 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: $\sigma_{\text{ln}K}=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].

2.2 Flow Field

Steady state flow through the network is modeled by Darcy’s law [22] for the fluid flux $u_{ij}$ between nodes $i$ and $j$, $u_{ij}=-K_{ij}(\Phi_{j}-\Phi_{i})/l_{ij}$, where $\Phi_{i}$ and $\Phi_{j}$ are the hydraulic heads at nodes $i$ and $j$. Imposing flux conservation at each node $i$, $\sum_{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 $\mathbf{e}_{ij}$ the unit vector in the direction of the link connecting nodes $i$ and $j$.

We 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 ($\Phi=1$) and right ($\Phi=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 of the velocity magnitudes in the network

[TABLE]

where $N_{\ell}$ is the number of links in the network. Link length and flow velocities here are independent. Thus, the Eulerian velocity PDF is given by

[TABLE]

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, 80]. This is illustrated in Figures 2a and b, 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 and d, 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.

2.3 Transport

Once the fluxes at the links have been determined, we simulate transport of a passive tracer by particle tracking. Particles are injected along a line at the inlet, $x=0$, with two different injection methods: (1) uniform injection, and (2) flux-weighted injection. Uniform (resident) injection introduces particles uniformly throughout the left boundary; this means that an equal number of particles is injected into each inlet node $i_{0}$,

[TABLE]

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}$

[TABLE]

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 the variability in porosity is significantly smaller than the the variability in conductivity [22, 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

[TABLE]

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 relations

[TABLE]

where $\mathbf{x}_{n}\equiv\mathbf{x}_{i_{n}}$ is the particle position after $n$ random walk steps, $\ell_{n}\equiv l_{i_{n}i_{n+1}}$ the particle displacement and $\mathbf{e}_{n}\equiv\mathbf{e}_{i_{n}i_{n+1}}$ its orientation; the particle velocity at the $n$th step is denoted by $u_{n}\equiv|u_{i_{n}i_{n+1}}|$. The particle displacement, orientation and velocity determined by the transition probability $p_{i_{n}j}$ 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 $\mathbf{x}(t)=\mathbf{x}_{i_{n_{t}}}$, where

[TABLE]

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

[TABLE]

where the angular brackets denote the average over all injected particles. As shown in Figure 4, both network heterogeneity and injection rule have significant impact on particle spreading. An increase in network heterogeneity leads to an increase in longitudinal particle spreading, and the uniform injection rule significantly enhances longitudinal spreading compared to flux-weighted injection. The impact of network heterogeneity and injection method can be clearly seen from projected concentration profiles, $f_{\tau}(\omega)$. Arbitrary injection modes are also studied and discussed in section 5.2.

3 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

[TABLE]

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.

3.1 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 $\tau_{a}(x_{c})=\text{inf}(t_{n}|\;|x_{n}-x_{0}|>x_{c})$ as

[TABLE]

Figure 5 shows FPTDs at the outlet, $f(\tau,x_{c}=200\bar{l})$, for different conductivity heterogeneities and injection rules. Conductivity heterogeneity has a clear impact on the FPTD by enhancing longitudinal spreading. This is so because stronger conductivity heterogeneity leads to broader particle transition time distribution, which in turn leads to enhanced longitudinal spreading. The injection rule also has a significant impact on FPTDs especially for high conductivity heterogeneity. FPTDs between the two different injection rules are similar for $\sigma_{\ln K}=1$, but uniform injection shows significantly stronger tailing for $\sigma_{\ln K}=2,3,5$ [Figure 5(b)]. As conductivity heterogeneity increases, the flux values at the inlet also becomes broader. For flux-weighted injection, most of particles are injected at the nodes with high flux values. However, for uniform injection, particles are uniformly injected across the injection nodes and relatively large number of particles are released at the nodes with low flux values. This leads to notable difference between the two injection rules and the difference grows as the conductivity heterogeneity increases.

3.2 Centered Mean Square Displacement

We also study longitudinal spreading in terms of the centered mean square displacement (cMSD) of average particle density, $\overline{P}(\mathbf{x},t)$. For the longitudinal direction ($x$), the cMSD is given by $\sigma^{2}_{x}(t)=\overline{\langle[x(t)-\langle x(t)\rangle]^{2}\rangle}$ where $\langle\cdot\rangle$ 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, $\bar{l}$. For both injection methods, spreading shows a ballistic regime ($\sim t^{2}$) at early times, which then transitions to a preasymptotic scaling in 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

[TABLE]

with the transition time $\tau_{n}={\ell_{0}}/{v_{n}}$. The velocities $v_{n}$ are distributed according to $p_{s}(v)$ which is space Lagrangian velocity PDF, and $\ell_{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 $\tau_{n}$ is given in terms of the space Lagrangian and Eulerian velocity PDFs as [88]

[TABLE]

where $\overline{v}$ is the average Eulerian velocity, see also Section 4. Specifically, for the scaling $p_{e}(v)\propto v^{\alpha}$ at small velocities, the transit time PDF scales as

[TABLE]

From Figure 3, we estimate for $\sigma_{\ln K}=5$ that $\alpha\approx-0.55$, which corresponds to $\beta=1.45$. CTRW theory [66, 65] predicts that the cMSD scales as $t^{3-\beta}$, which here implies $t^{1.55}$. This is consistent with the late-time scaling of the cMSD shown in 6 for $\sigma_{\ln K}=5$.

The Monte Carlo simulations show that, in the intermediate regime ($t/\overline{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].

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. 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.

4 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 positions. Recent studies, however, have shown that the underlying mechanisms of anomalous transport can be quantified through an analysis of the statistics of Lagrangian particle velocities such as velocity distribution and correlation [74, 76, 98, 77, 78, 45, 80]. In the following, we briefly introduce two viewpoints for analyzing Lagrangian velocities—equidistantly and isochronally along streamlines—and the relation between them [88]. We then proceed to a detailed analysis of the Lagrangian velocity statistics measured equidistantly along streamlines.

4.1 Lagrangian Velocities

Particle motion is described here by the recursion relations (7). In this framework, we consider two types of Lagrangian velocities. The t(ime)–Lagrangian velocities are measured at a given time $t$,

[TABLE]

where $n_{t}$ is defined by (8). 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

[TABLE]

The number of steps needed to cover the distance $s$ is described by $n_{s}=\sup(n|s_{n}\leq s)$. Thus, the particle velocity at a distance $s$ along a trajectory is given by

[TABLE]

The PDF of t–Lagrangian velocities sampled along a particle path is given by

[TABLE]

where we defined the transit time,

[TABLE]

The PDF of s–Lagrangian velocities sampled along a particle path are defined analogously as

[TABLE]

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 $\tau_{i}$ in (19), the PDFs of the s– and t–Lagrangian velocities are related through flux weighting as [88]

[TABLE]

Furthermore, for flux-preserving flows and under ergodic conditions, the Eulerian and t-Lagrangian velocity PDFs are equal,

[TABLE]

Thus, under these conditions, the s–Lagrangian and Eulerian velocity PDFs are related as [88]

[TABLE]

This means that the stationary s-Lagrangian velocity PDF can be determined from the Eulerian velocity PDF. Figure 7 illustrates this relation by comparing the s-Lagrangian velocity PDFs measured from the numerical simulation to the flux-weighted Eulerian velocity PDFs shown in Figure 3.

4.2 Evolution of Lagrangian Velocity Distributions

It is important to emphasize that the above definitions of the Lagrangian velocity PDFs refer to stationary conditions. We now define the PDFs of t– and s–Lagrangian velocities through sampling between particles and network realizations at a given time (t-Lagrangian) or space (s-Lagrangian) velocities

[TABLE]

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 $\hat{p}_{t}(v,t=0)=\hat{p}_{s}(v,s=0)=p_{0}(v)$, but their respective stationary PDFs are different, namely

[TABLE]

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

[TABLE]

that is, the initial velocity PDF is equal to the Eulerian and thus t–Lagrangian velocity PDFs. This means that for the uniform injection method, the t–Lagrangian velocity PDF is steady, $\hat{p}_{t}(v,t)=p_{t}(v)$, while the s–Lagrangian velocity PDF is not. It evolves from its initial distribution $p_{s}(v,s=0)=p_{e}(v)$ to the steady state distribution (23).

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

[TABLE]

due to relation (21). The initial velocity PDF is equal to the s–Lagrangian velocity PDF. This means that under flux-weighting, the s–velocity PDF is steady, $\hat{p}_{s}(v,s)=p_{s}(v)$. Under these conditions, the t–Lagrangian velocity PDF $\hat{p}_{t}(v,t)$ evolves from the initial distribution $\hat{p}_{t}(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.

4.3 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, $\bar{l}$. The Lagrangian velocity $v_{s}(s_{n})$ at a distance $x_{n}=n\bar{l}$ along the projected trajectory is approximated by the average velocity $v_{n}\equiv\bar{l}/\tau_{n}$ where $\tau_{n}$ is the transition time for the distance $\bar{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 velocities have been studied by Cvetkovic et al. [73] and Gotovac et al. [75] for highly heterogeneous porous media and by Frampton and Cvetkovic [16] for 3D DFNs in view of quantifying particle travel time statistics and thus modelling effective particle motion.

We first study the convergence of the s–Lagrangian velocity PDFs towards a stationary distribution and the invariance of $\hat{p}_{s}(v,s)$ for a stationary (flux-weighted) 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 $\hat{p}_{s}(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 $\hat{p}_{s}(v,s=0)$ and $\hat{p}_{s}(v,s=x_{c})$ for uniform and flux-weighted injection conditions and two different heterogeneity strengths. We clearly observe that $\hat{p}_{s}(v,x_{c})=p_{s}(v)$ is invariant for flux-weighted injection. For uniform injection, $\hat{p}_{s}(v,x_{c})$ has already evolved towards its steady limit after $x_{c}=200\bar{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 $\hat{p}_{s}(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 $\psi_{\tau}(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 stationary. Figure 9 illustrates the PDF of transition times for different variances of $\ln(K)$. As $\sigma_{\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 $\Delta s=s-s^{\prime}$. As pointed out above, for flux-weighted injection, the streamwise velocities here are stationary and therefore

[TABLE]

In order to increase the statistics, we furthermore sample along streamlines over a distance of $10^{2}\bar{l}$. The velocity variance is $\sigma^{2}_{v}=C_{s}(0)$. The velocity autocorrelation function $\chi_{s}(s)=C_{s}(s)/\sigma_{v}^{2}$. The correlation length scale $\ell_{c}$ is defined by

[TABLE]

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 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 $\ell_{c}$, may be considered independent. In order to 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,

[TABLE]

Numerically, the transition probability is determined by discretizing the s-Lagrangian velocity PDF into $N$ velocity classes $\mathcal{C}_{i}=(v_{s,i},v_{s,i}+\Delta 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^{\prime})$ is approximated numerically as

[TABLE]

where the indicator function $\mathbb{I}_{\mathcal{C}_{i}}(v)$ is $1$ if $v\in\mathcal{C}_{i}$ and [math] otherwise.

Figure 10 shows the one-step transition matrix $T_{1}(i|j)$ for equidistant and logarithmically equidistant velocity classes for different network heterogeneity. Higher probabilities along the diagonal than in the off-diagonal positions indicate correlation between subsequent steps, which, however, decreases as the number of steps along the particle trajectory increases, as indicated by the existence of a finite correlation scale $\ell_{c}$.

5 Stochastic Particle Motion and Effective Transport Model

In the following, we describe the evolution of the s-Lagrangian velocities by a Markov-chain, which is motivated by the existence of a finite spatial 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 time-domain random walk (TDRW) and continuous time random walk (CTRW) approaches discussed in the Introduction and below.

5.1 Markovian Velocity Process

Along the lines of Le Borgne et al. [74] and Kang et al. [77], 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^{\prime})$ satisfies the Chapman–Kolmogorov equation [99]

[TABLE]

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}=\ell_{c}/\bar{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 PDF $p_{e}(v)$. Thus, the numerical Monte Carlo simulations are compared to the predictions of (32b) for the initial condition $\hat{p}_{s}(v,s_{0}=0)=p_{e}(v)$. The transition PDF $r_{1}(v|v^{\prime})$ is given by (31) 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^{\prime})$. 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^{\prime})$ for $n>n_{c}$ converges to the stationary PDF $r_{n}(v|v^{\prime})\to p_{s}(s)$. The number of steps for velocities to decorrelate can be quantified by

[TABLE]

The simplest transition PDF that shares these characteristics is [41, 80, 50, 88]

[TABLE]

with $a=\exp(-\bar{l}/\ell_{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 the field from multiscale tracer tests [41]. In the following, we study particle dispersion in the Markovian velocity model for the full transition PDF shown in Figure 10 and the reduced-order Markov model (34).

5.2 Particle Dispersion and Model Predictions

We consider particle motion along the mean pressure gradient in $x$–direction, which is described by the stochastic regression

[TABLE]

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 $\ell_{0}$ larger than $\ell_{c}$, or equivalently, step numbers $n\gg n_{c}\equiv\ell_{c}/\bar{l}$, subsequent velocities may be considered independent and particle dispersion is fully characterized by the recursion relation (12) and the transition time PDF (13). 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 by the transit time PDF

[TABLE]

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 [100, 101, 74, 102]. The joint distribution $p(x,v,t)$ of particle position and velocity at a given time $t$ is given by [74]

[TABLE]

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

[TABLE]

where $R_{0}(x,v,t)=p_{0}(x,v)\delta(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^{\prime}$ 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

[TABLE]

see A. This result provides a bridge between the TDRW approach [70, 71] and the correlated CTRW approach.

As illustrated in Figure 3, the Eulerian velocity PDF can be characterized by the power-law $p_{e}(v)\propto v^{\alpha}$. Thus, the first CTRW steps until the decorrelation at $n=n_{c}$ are characterized by the transit time PDF $\psi_{0}(t)\propto t^{-2-\alpha}$. The corresponding tail of the BTC is $f(t,x_{c})\propto t^{-2-\alpha}$. The observed value of $\alpha=-0.55$ explains the tailing of the BTC in Figure 5 as $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 $\beta$ 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 Figure 10. Figure 6 shows 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 flux-weighted injection conditions, the predictions of particle spreading are in excellent agreement with the direct numerical simulations. In Figure 12 we 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 velocity PDF 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 figure 9). 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 (Figure 13). 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 (figure 10), 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 [88]. 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 boundary may 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 [103]. 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 (Figure 15). Eventually, the s-Lagrangian PDFs evolve towards the stationary flux-weighted Eulerian PDF as discussed in the previous section.

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. 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 (12). 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.

6 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.

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 $\ell_{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.

In this study, we investigated the particle transport and the impact of the injection condition for idealized 2D DFN using a Markov velocity model. 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.

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 Council (ERC) through the project MHetScale (617511). TLB acknowledges the support of European Research Council (ERC) through the project ReactiveFronts (648377). RJ acknowledges the support of the US Department of Energy through a DOE Early Career Award (grant DE-SC0009286). The data to reproduce the work can be obtained from the corresponding author.

Appendix A Boltzmann Equation

The time derivative of (37) gives

[TABLE]

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 $\bar{\/}v$. Thus we can approximate (37) as

[TABLE]

Using this approximation and combining (40) with (38) gives for $t>0$

[TABLE]