Stochastic Modeling and Upscaling of Hydrodynamic Transport in Geological Fractures

TL;DR

This work addresses advection-dominated transport in geological fractures, where multiscale aperture heterogeneity leads to anomalous transport such as BTC tailing. The authors generate self-affine, wall-matched fracture apertures with a finite correlation length $L_c$ and prescribed relative closure $\sigma_a/\langle a\rangle$, solve depth-averaged Stokes flow via the lubrication (Reynolds) equation, and track solute transport using two particle-based upscaling approaches: time-domain random walk (TDRW) and a one-dimensional CTRW based on an Ornstein–Uhlenbeck process for Lagrangian velocities. They demonstrate that velocity PDFs are insensitive to $L_c$ but strongly shaped by fracture closure, with low-velocity tails characterized by $\alpha$ (e.g., $\alpha=0.4$ for $\sigma_a/\langle a\rangle=0.75$ and $\alpha=1$ for $\sigma_a/\langle a\rangle=0.25$). The mean plume position scales linearly in time, while the displacement variance shows ballistic early-time behavior and a late-time regime $\mathcal{V}(t)\propto t^{2-\alpha}$ or $\mathcal{V}(t)\propto t\ln t$ depending on $\alpha$, with BTC tails $\mathcal{F}(t,L)\propto t^{-1-\alpha}$ (uniform) and $\mathcal{F}(t,L)\propto t^{-2-\alpha}$ (flux-weighted). The upscaled CTRW model, requiring only $f_e(u)$, $\chi$, and $\ell_c$, accurately reproduces full TDRW results and enables analytical transport scalings, providing an efficient framework for uncertainty quantification and large-scale prediction in fractured media.

Abstract

Characterizing hydrodynamic transport in fractured rocks is essential for carbon storage and geothermal energy production. Multiscale heterogeneities lead to anomalous solute transport, with breakthrough-curve (BTC) tailing and nonlinear growth of plume moments. We study purely advective transport in synthetic fractures with prescribed relative closure $ σ_a/\langle a \rangle $ and correlation length $ L_c $. For each geometry we generate multiple realizations and solve steady, depth-averaged Stokes flow under the lubrication approximation. Flow heterogeneity persists up to $ L_c $. The ensemble-averaged velocity PDFs are insensitive to $ L_c $ but strongly affected by $ σ_a/\langle a \rangle $, particularly their low-velocity power-law scaling. A time-domain random walk (TDRW) yields plume moments and outlet BTCs: the mean longitudinal position grows linearly in time, while the variance shows early ballistic scaling and a late-time regime controlled by the low-velocity power law with exponent $ α$, which depends on $ σ_a/\langle a \rangle $. BTC properties, including peak broadening and tail scaling, are likewise governed by $ α$. We further model advection with a one-dimensional continuous-time random walk (CTRW) that uses only the velocity PDF, flow tortuosity, and $ L_c $. CTRW results closely match TDRW and enable analytical predictions of asymptotic transport scalings.

Figures (8)

• Figure 1: Overview of the geometries and associated spectra from synthetic geological fractures. Figure \ref{['Fig1']}(a) shows one longitudinal profile of the fracture, with the representation of the walls and the definition of the geometrical fields that define the fracture geometry. Figure \ref{['Fig1']}(b) presents the power spectrum for the cases $L/L_\mathrm{c} = 2^5$ in red and $L/L_\mathrm{c} = 2^3$ in green. Figure \ref{['Fig1']}(c) shows four aperture field realizations arranged from left to right, each with its corresponding local aperture probability density function (PDF) displayed below. The first two realizations correspond to $L/L_\mathrm{c} = 2^3$ with closure values of 0.25 and 0.75, respectively, while the last two correspond to $L/L_\mathrm{c} = 2^5$, also with closures of 0.25 and 0.75. In all cases, the local aperture PDFs are approximately Gaussian, with a cutoff at zero when contact regions are present in the fracture plane.
• Figure 2: (a) Representation of the domain partitioning with boundary conditions. (b) Finite volume scheme 5-point stencil: pressure is defined at the centre of each control volume, while the local aperture is estimated along the edge of the cells by arithmetic averaging.
• Figure 3: Maps of fracture apertures (left column) and the corresponding Eulerian velocity magnitude (middle and right columns) at two different times, with $10^6$ superimposed flux-weighted injected particles at the indicated times $t$, for two synthetic fractures with different correlation lengths, $L/L_\mathrm{c} = 2^3$ and $2^5$, for $\sigma_a/\langle a\rangle=0.75$, $L_\mathrm{c}=0.1~\textrm{m}$, and $\langle a\rangle=0.001~\textrm{m}$. Contact zones are depicted in black.
• Figure 4: Flow chart of the numerical modeling workflow. From geometry generation and flow simulation to particle transport and upscaling, all steps are embedded within the Monte Carlo framework.
• Figure 5: Probability density functions (PDFs) of Eulerian (blue) and $s$-Lagrangian (yellow) velocities for the four Monte Carlo realizations: (a) MC1, (b) MC2, (c) MC3, and (d) MC4. The corresponding parameters used for aperture field generation are listed in Table \ref{['tab1']}. Trend lines emphasize the scaling behavior of the low-velocity tails: the Eulerian PDF scales as $u^{\alpha - 1}$ (dashed line), while the $s$-Lagrangian PDF scales as $u^{\alpha}$ (dash-dotted line), in agreement with theoretical predictions for transport in heterogeneous flow fields. The shaded areas represent the confidence interval between the 5th and 95th percentiles.