Effects of Wall Roughness on Coupled Flow and Heat Transport in Fractured Media

TL;DR

This paper tackles non-Fickian heat transport in fractured media by coupling a time-domain random walk (TDRW) for fracture advection–diffusion with a semi-analytical matrix-diffusion exchange model. Matrix trapping times are drawn from the Lévy–Smirnov distribution, derived from first-passage diffusion in a semi-infinite medium, and fracture–matrix heat transfer is computed via a nonlocal convolution (Duhamel) kernel that yields the characteristic $t^{-1/2}$ decay in long-time heat flux. Monte Carlo simulations across varying aperture heterogeneity, correlation length, and Péclet number reveal a transition from early-time superdiffusive (or ballistic) transport to late-time diffusion- or subdiffusion-dominated regimes, with memory effects persisting through the interface. The framework provides physically grounded, computationally efficient predictions of thermal transport in complex fractured systems, with implications for geothermal energy, thermal storage, and engineered heat exchange in low-permeability environments.

Abstract

Heat transfer in fractured media is governed by the interplay between advective transport along rough-walled fractures and conductive transport, both within the fractures and in the surrounding low-permeability matrix. Flow localization induced by aperture heterogeneity, combined with matrix conduction, gives rise to anomalous thermal behavior. To capture these effects, we develop a stochastic modeling framework that couples a time-domain random walk (TDRW) representation of advective and conductive transport in the fractures with a semi-analytical model of conductive heat exchange with the matrix. Matrix trapping times follow a Lévy-Smirnov distribution derived from first-passage theory, capturing the heavy-tailed dynamics typical of fractured systems. Heat flux at the fracture-matrix interface is computed via a nonlocal convolution integral based on Duhamel's principle, accounting for thermal memory effects. The model is validated against analytical benchmarks and finite-element simulations. Monte Carlo simulations over stochastic aperture fields quantify the influence of fracture closure, correlation length, and Péclet number. Results reveal a transition from superdiffusive to subdiffusive regimes, driven by the competition between advective transport along preferential paths, dispersion induced by aperture variability, and matrix-driven heat conduction. In the long-time regime, heat exchange exhibits a characteristic $t^{-1/2}$ decay. At early times, limited thermal penetration into the matrix leads to weaker interfacial fluxes, underscoring the role of matrix thermal inertia. The proposed framework enables physically consistent and computationally efficient simulations of thermal transport in complex fractured systems, with implications for geothermal energy, subsurface thermal storage, and engineered heat exchange in low-permeability environments.

Figures (7)

• Figure 1: (a) Synthetic fracture aperture fields for increasing correlation ratios $L / L_\mathrm{c} = \{2,\,16,\,64\}$ (left to right), all with the same fracture closure $\sigma_a / \langle a\rangle = 0.8$. (b) Dimensionless velocity magnitude field, $\log_{10}(u/\langle u\rangle)$, for the case $L / L_\mathrm{c} = 16$, showing the flow structure under the applied boundary conditions. (c) Transverse fracture profile illustrating wall roughness, aperture geometry, matrix domain, and the initial and boundary conditions for flow and transport. Panels (b) and (c) also indicate the transport mechanisms considered in both fracture and matrix. All synthetic aperture fields were generated using $L_\mathrm{c} = 0.1~\textrm{m}$, $H = 0.8$, $a_\mathrm{m} = 1~\textrm{mm}$, and $\sigma_a / \langle a \rangle = 0.8$.
• Figure 2: Snapshots of the thermal anomaly $\Delta T_{\textrm{f}}$ at successive times for two levels of fracture closure, $\sigma_a/\langle a \rangle = 0.5$ (left) and $\sigma_a/\langle a \rangle = 1.0$ (right). Each row corresponds to a different time instant: $t = \{10^3\,\textrm{s},\,10^4\,\textrm{s},\,10^5\,\textrm{s}\}$, from top to bottom. Simulations were conducted with a mean aperture $\langle a \rangle = 10^{-3}~\mathrm{m}$ and aperture correlation length $L_\mathrm{c} = 10^{-1}~\mathrm{m}$.
• Figure 3: Panel (a) compares the thermal anomaly fields obtained from COMSOL Multiphysics® simulations (first row) and from the TDRW particle tracking scheme (second row). Snapshots at different times, $t = \{50~\textrm{s},\,500~\textrm{s},\,1000~\textrm{s}\}$, are shown from left to right. The first three columns correspond to the parallel plate geometry, while the fourth to sixth columns refer to the heterogeneous case with aperture variability $\sigma_a/\langle a \rangle = 0.21$ and correlation ratio $L/L_\textrm{c} = 1$. Results are shown for a fixed Péclet number $\hbox{Pe} = 51$. Panel (b) displays spatial temperature profiles at a given time (left) and breakthrough curves at the fracture outlet (right) for different Péclet numbers, $\hbox{Pe} = \{10,\,50,\,100\}$. The analytical solution \ref{['eq:CDF_arrival_times_pp']} (black solid line with markers) is superimposed on the TDRW predictions (solid colored lines). Validation was conducted on a fracture realizations with length $L = 5~\mathrm{m}$, other model parameters are listed in Table \ref{['Tab1']}.
• Figure 4: Temporal evolution of the ensemble averages (solid lines) and percentile bands (dashed lines) of the logarithm of the mean displacements, $\log_{10}{\mathcal{M}(t)}$, obtained from the Monte Carlo simulations described in Table \ref{['Tab1']}. Each column, from top to bottom, corresponds to increasing fracture closures, $\sigma_a/\langle a\rangle = \{0.2,\,0.6,\,1.0\}$, reflecting growing aperture variability. Rows correspond to increasing ratios $L/L_\textrm{c} = \{2,\,16,\,64\}$, indicating progressively larger fracture sizes relative to the correlation length, and hence increasing statistical ergodicity of the aperture field. Solid colored lines denote different Péclet numbers: $\hbox{Pe} = 10$ (blue), $\hbox{Pe} = 50$ (red), and $\hbox{Pe} = 100$ (green). Black solid lines are included as visual references for diffusive ($\propto t^{0.5}$) and ballistic ($\propto t$) scaling. Displacement and time are in meters and seconds, respectively. Simulation parameters are summarized in Table \ref{['Tab1']}.
• Figure 5: Temporal evolution of the ensemble averages (solid lines) and percentile bands (dashed lines) of the logarithm of the mean displacements, $\log_{10}{\mathcal{M}(t)}$, obtained from the Monte Carlo simulations summarized in Table \ref{['Tab1']}. Each column (top to bottom) corresponds to increasing fracture closures, $\sigma_a/\langle a\rangle = \{0.2,\,0.6,\,1.0\}$, reflecting growing aperture variability. Rows indicate increasing ratios $L/L_\textrm{c} = \{2,\,16,\,64\}$, corresponding to progressively larger fracture sizes relative to the correlation length, and hence greater statistical ergodicity of the aperture field. Solid colored lines represent different Péclet numbers: $\hbox{Pe} = 10$ (blue), $\hbox{Pe} = 50$ (red), and $\hbox{Pe} = 100$ (green). Black solid lines serve as visual references for diffusive ($\propto t^{0.5}$) and ballistic ($\propto t$) scaling. Displacement and time are reported in meters and seconds, respectively.