# Reduced-Order Models Unravel the Joint Impact of Aperture Heterogeneity and Shear-Thinning Rheology on Fracture-Scale Flow Metrics

A. Lenci, Yves Méheust, V. Di Federico, V. Ciriello

<figure><img src="image_1.png" /></figure>

## To cite this version:

A. Lenci, Yves Méheust, V. Di Federico, V. Ciriello. Reduced-Order Models Unravel the Joint Impact of Aperture Heterogeneity and Shear-Thinning Rheology on Fracture-Scale Flow Metrics. Water Resources Research, 2024, 60 (1), pp.e2023WR036026. 10.1029/2023WR036026. insu-04431399

HAL Id: insu-04431399

https://insu.hal.science/insu-04431399

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Reduced-Order Models Unravel the Joint Impact of Aperture Heterogeneity and Shear-Thinning Rheology on Fracture-Scale Flow Metrics

A. Lenci, Yves Méheust, V. Di Federico, V. Ciriello

► To cite this version:

A. Lenci, Yves Méheust, V. Di Federico, V. Ciriello. Reduced-Order Models Unravel the Joint Impact of Aperture Heterogeneity and Shear-Thinning Rheology on Fracture-Scale Flow Metrics. Water Resources Research, 2024, 60 (1), pp.e2023WR036026. 10.1029/2023WR036026. insu-04431399

HAL Id: insu-04431399

https://insu.hal.science/insu-04431399v1

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR036026

**Key Points:**

* We investigate the dependence of the fracture's hydraulic behavior on its geometry, externally imposed pressure gradient, and fluid rheology

* Reduced-order models are used to approximate the model's responses and to explore a wider parameter space with a much-reduced CPU time

* The obtained dependence of the hydraulic properties on the problem's parameters confirms and widely extends our previous findings

**Correspondence to:**

A. Lenci,
alessandro.lenci@unibo.it

**Citation:**

Lenci, A., Méheust, Y., Di Federico, V., & Ciriello, V. (2024). Reduced-order models unravel the joint impact of aperture heterogeneity and shear-thinning rheology on fracture-scale flow metrics. *Water Resources Research*, *60*, e2023WR036026. https://doi.org/10.1029/2023WR036026

Received 8 AUG 2023
Accepted 6 JAN 2024

**Author Contributions:**

**Conceptualization:** A. Lenci, V. Ciriello

**Data curation:** A. Lenci

**Formal analysis:** A. Lenci

**Funding acquisition:** A. Lenci, V. Di Federico

**Investigation:** A. Lenci, Y. Méheust, V. Ciriello

**Methodology:** A. Lenci, Y. Méheust, V. Di Federico, V. Ciriello

**Project Administration:** V. Di Federico

**Resources:** A. Lenci, V. Di Federico

**Software:** A. Lenci

**Supervision:** Y. Méheust, V. Di Federico, V. Ciriello

**Validation:** A. Lenci, V. Ciriello

**Visualization:** A. Lenci

**Writing – original draft:** A. Lenci, V. Ciriello

© 2024. The Authors.

This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_3.png" /></figure>

## Reduced-Order Models Unravel the Joint Impact of Aperture Heterogeneity and Shear-Thinning Rheology on Fracture-Scale Flow Metrics

A. Lenci¹ ID, Y. Méheust² ID, V. Di Federico¹ ID, and V. Ciriello¹ ID

¹Department of Civil, Chemical, Environmental and Materials Engineering, Università di Bologna, Bologna, Italy,

²University of Rennes, CNRS, Géosciences Rennes—UMR 6118, Rennes, France

**Abstract** Subsurface industrial operations often make use of complex engineered fluids. In fractured media, the hydraulic behavior of a geological fracture is affected by the shear-thinning (ST) rheology of such fluids, depending on the imposed pressure gradient. Besides, owing to the stochastic nature of heterogeneities in the fracture's local apertures, obtaining the generic behavior of a fracture under ST flow requires large statistics of fractures whose aperture fields have the same statistical properties, namely their mean value, standard deviation, correlation length $L_c$, and Hurst exponent (which controls the scale invariances at scales smaller than $L_c$). The first such Monte Carlo study was recently proposed by Lenci, Putti, et al. (2022) using a new lubrication-based model relying on a non-linear Reynolds equation. This model describes the flow of an Ellis fluid, with both a low shear rate quasi-Newtonian viscosity plateau and a large shear rate power law ST behavior, as measured for engineered fluids such as biopolymer solutions (e.g., xanthan gum). Here we aim to obtain the dependence of the fracture transmissivity and its dispersion over the statistics, as well as of the flow correlation length, on fracture closure, geometry correlation length, and applied pressure gradient, over a vast volume of the parameter space, and in a simple mathematical form. Employing reduced-order models based on the polynomial chaos expansion theory to this aim, we discuss the properties of the obtained topographies of interest in the parameter space.

### 1. Introduction

Industrial activities in the subsurface commonly involve the injection of engineered working fluids in fractured aquifers or reservoirs. These operations aim at resource recovery (Kamal et al., 2015), geothermal energy exploitation (Pruess, 2006), carbon sequestration (Leung et al., 2014) or soil remediation (Zhang, 2003); reaching depths under the ground level that range from a few meters (e.g., reclamation of chlorinated solvent-contaminated groundwater, see Yang et al. (2017)) up to tens of kilometers (e.g., enhanced geothermal systems, see Zhu et al. (2015)).

Crystalline geological formations are made of fractured rocks, with blocks of low porosity matrix interspersed by interconnected fractures forming networks. Fluids mainly flow through such networks, as fractures provide a preferential path due to the typical high contrast between the fractures and matrix conductance (Bodin et al., 2003). The flow across networks is strongly affected by the heterogeneity in fracture density, lengths (Bour & Davy, 1998) and mean apertures within the network (de Dreuzy et al., 2002), but also by the fracture-scale heterogeneity resulting from spatial variations in the aperture fields of individual fractures (de Dreuzy et al., 2012); the latter can even impact the network-scale topology of the flow. Such typical millimetric variations of local fracture apertures are due to mechanical wear and chemical weathering (Brown, 1987). Incorporating the effect of fracture-scale heterogeneity in network-scale simulations is thus critical to obtain accurate predictions. In the literature, this has been done by meshing all fractures in a fracture network in two dimensions (2D), considering realistic aperture variations within fracture planes, and solving the flow over a single mesh encompassing all fractures (de Dreuzy et al., 2012) or by employing upscaling techniques that rely on the statistical characterization of macroscale flow metrics (MFM s) (Boso & Tartakovsk y, 2018).

At the single fracture scale, theoretical (Ge, 1997), numerical (Thompson & Brown, 1991) and experimental (Méheust & Schmittbuhl, 2000) studies have highlighted the importance of considering the influence of walls roughness to properly assess the macroscale hydraulic behavior (Brown, 1989; Durham & Bonner, 1995; Tsang & Tsang, 1987). The aforementioned spatial variations in fracture apertures induce flow channeling, and in
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR036026

**Key Points:**

* We investigate the dependence of the fracture's hydraulic behavior on its geometry, externally imposed pressure gradient, and fluid rheology

* Reduced-order models are used to approximate the model's responses and to explore a wider parameter space with a much-reduced CPU time

* The obtained dependence of the hydraulic properties on the problem's parameters confirms and widely extends our previous findings

**Correspondence to:**

A. Lenci,
alessandro.lenci@unibo.it

**Citation:**

Lenci, A., Méheust, Y., Di Federico, V., & Ciriello, V. (2024). Reduced-order models unravel the joint impact of aperture heterogeneity and shear-thinning rheology on fracture-scale flow metrics. *Water Resources Research*, *60*, e2023WR036026. https://doi.org/10.1029/2023WR036026

Received 8 AUG 2023
Accepted 6 JAN 2024

**Author Contributions:**

**Conceptualization:** A. Lenci, V. Ciriello

**Data curation:** A. Lenci

**Formal analysis:** A. Lenci

**Funding acquisition:** A. Lenci, V. Di Federico

**Investigation:** A. Lenci, Y. Méheust, V. Ciriello

**Methodology:** A. Lenci, Y. Méheust, V. Di Federico, V. Ciriello

**Project Administration:** V. Di Federico

**Resources:** A. Lenci, V. Di Federico

**Software:** A. Lenci

**Supervision:** Y. Méheust, V. Di Federico, V. Ciriello

**Validation:** A. Lenci, V. Ciriello

**Visualization:** A. Lenci

**Writing – original draft:** A. Lenci, V. Ciriello

© 2024. The Authors.

This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Reduced-Order Models Unravel the Joint Impact of Aperture Heterogeneity and Shear-Thinning Rheology on Fracture-Scale Flow Metrics

A. Lenci¹ ID, Y. Méheust² ID, V. Di Federico¹ ID, and V. Ciriello¹ ID

¹Department of Civil, Chemical, Environmental and Materials Engineering, Università di Bologna, Bologna, Italy,

²University of Rennes, CNRS, Géosciences Rennes—UMR 6118, Rennes, France

**Abstract** Subsurface industrial operations often make use of complex engineered fluids. In fractured media, the hydraulic behavior of a geological fracture is affected by the shear-thinning (ST) rheology of such fluids, depending on the imposed pressure gradient. Besides, owing to the stochastic nature of heterogeneities in the fracture's local apertures, obtaining the generic behavior of a fracture under ST flow requires large statistics of fractures whose aperture fields have the same statistical properties, namely their mean value, standard deviation, correlation length $L_c$, and Hurst exponent (which controls the scale invariances at scales smaller than $L_c$). The first such Monte Carlo study was recently proposed by Lenci, Putti, et al. (2022) using a new lubrication-based model relying on a non-linear Reynolds equation. This model describes the flow of an Ellis fluid, with both a low shear rate quasi-Newtonian viscosity plateau and a large shear rate power law ST behavior, as measured for engineered fluids such as biopolymer solutions (e.g., xanthan gum). Here we aim to obtain the dependence of the fracture transmissivity and its dispersion over the statistics, as well as of the flow correlation length, on fracture closure, geometry correlation length, and applied pressure gradient, over a vast volume of the parameter space, and in a simple mathematical form. Employing reduced-order models based on the polynomial chaos expansion theory to this aim, we discuss the properties of the obtained topographies of interest in the parameter space.

## 1. Introduction

Industrial activities in the subsurface commonly involve the injection of engineered working fluids in fractured aquifers or reservoirs. These operations aim at resource recovery (Kamal et al., 2015), geothermal energy exploitation (Pruess, 2006), carbon sequestration (Leung et al., 2014) or soil remediation (Zhang, 2003); reaching depths under the ground level that range from a few meters (e.g., reclamation of chlorinated solvent-contaminated groundwater, see Yang et al. (2017)) up to tens of kilometers (e.g., enhanced geothermal systems, see Zhu et al. (2015)).

Crystalline geological formations are made of fractured rocks, with blocks of low porosity matrix interspersed by interconnected fractures forming networks. Fluids mainly flow through such networks, as fractures provide a preferential path due to the typical high contrast between the fractures and matrix conductance (Bodin et al., 2003). The flow across networks is strongly affected by the heterogeneity in fracture density, lengths (Bour & Davy, 1998) and mean apertures within the network (de Dreuzy et al., 2002), but also by the fracture-scale heterogeneity resulting from spatial variations in the aperture fields of individual fractures (de Dreuzy et al., 2012); the latter can even impact the network-scale topology of the flow. Such typical millimetric variations of local fracture apertures are due to mechanical wear and chemical weathering (Brown, 1987). Incorporating the effect of fracture-scale heterogeneity in network-scale simulations is thus critical to obtain accurate predictions. In the literature, this has been done by meshing all fractures in a fracture network in two dimensions (2D), considering realistic aperture variations within fracture planes, and solving the flow over a single mesh encompassing all fractures (de Dreuzy et al., 2012) or by employing upscaling techniques that rely on the statistical characterization of macroscale flow metrics (MFM s) (Boso & Tartakovsk y, 2018).

At the single fracture scale, theoretical (Ge, 1997), numerical (Thompson & Brown, 1991) and experimental (Méheust & Schmittbuhl, 2000) studies have highlighted the importance of considering the influence of walls roughness to properly assess the macroscale hydraulic behavior (Brown, 1989; Durham & Bonner, 1995; Tsang & Tsang, 1987). The aforementioned spatial variations in fracture apertures induce flow channeling, and in
**Writing – review & editing:** Y. Méheust,
V. Di Federico, V. Ciriello

particular flow localization in areas of higher conductance, where the fracture aperture is larger. The resulting hydraulic behavior of the (rough) geological fracture of interest is usually compared to the parallel plate model where both walls are assumed to be perfectly smooth planes, parallel to each other and separated by a distance that is identical to the mechanical aperture (distance between the mean planes of the two rough walls) of the rough fracture. Flow channeling induces either a flow-enhancing or flow-inhibiting hydraulic behavior with respect to the corresponding parallel plate model (Méheust & Schmittbuhl, 2000). However, over a large population of fractures with identical statistical properties, flow-inhibition dominates (Méheust & Schmittbuhl, 2001), but the statistical dispersion (i.e., variability) over the fractures' hydraulic behavior depends on their correlation length (Méheust & Schmittbuhl, 2003).

This study focuses both on Newtonian and shear-thinning (ST) fluids. Subsurface industrial activities widely adopt complex fluids, such as foams (Faroughi et al., 2018; Hosseini et al., 2019), clayey muds, polymer solutions, and aqueous and non-aqueous suspensions, which are ST, that is, whose viscosity tends to decrease under increasing shear rates. CO₂-based ST fluids (Jiang, 2021) have also been recently proposed as a water-saving solution to reduce the use of freshwater in subsurface activities. ST fluids have been demonstrated to enhance the cleanup efficiency (Hosseini et al., 2019) and carrying capability of remedial amendments (Tosco et al., 2012; Zhong et al., 2008).

The effects of ST rheology on the transmissivity of variable aperture fractures were first investigated in the past using simplified deterministic (Di Federico, 1997, 2001) or stochastic (Di Federico, 1998) conceptualizations. More recently, the fracture's hydraulic behavior has been shown to be augmented by the interplay between pore-scale heterogeneity and ST rheology. This potentially leads to fracture transmissivities orders of magnitude larger than the Newtonian counterpart (Lenci, Méheust, et al., 2022). The transmissivity is computed from the ratio of the volumetric flow rate to the applied macroscopic pressure gradient. In numerical simulations, the flow rate is obtained from simulating the fluid velocity field. The Reynolds equation, based on the lubrication approximation (Brown, 1987), is often used to model Newtonian flow in geological fractures. This depth-averaged approach relies on the local validity of the Hagen-Poiseuille law, which holds as long as the conditions for Stokes flow are satisfied and spatial variations of the aperture fields are slow.

Recently, an efficient algorithm has been proposed by Lenci, Putti, et al. (2022) to solve a generalized non-linear Reynolds equation that extends the model to ST fluids whose rheology includes a transition from a Newtonian plateau at low shear rates to a power law ST behavior at large shear rates. This numerical strategy has been adopted to produce a Monte Carlo (MC) study of the ST flow through rough fractures (Lenci, Méheust, et al., 2022), overcoming the computation burden that is typical of stochastic analysis of non-linear problems in spatially varying domains and thus allowing to evidence the impact of the fluid's rheology as discussed above. However, exploring the space of controlling parameters for ST flow in geological fractures over a wide range of values of these parameters, using fracture-scale simulations, would inevitably require a vast number of MCs, and thus a tremendous amount of computational time, due to the dimensionality of the parameter space (consisting of at least three dimensions: fracture closure, fracture correlation length, imposed pressure gradient), the nonlinearity of the model, and the variability of the aperture field's topology over different realizations.

In this study, reduced-order models (ROMs) are thus employed to drastically decrease the computational cost of this statistical framework and to obtain the dependence of given quantities of interest (QoIs) describing a fracture's hydraulic behavior on the controlling parameters in a simple mathematical form, rather than using a pure MC approach as previously done by Lenci, Méheust, et al. (2022). The QoIs, whose response surfaces are reconstructed in this study, are the ensemble average and difference between the third and first quartiles (computed over the number of MCs of the aperture field) of a set of MFMs, which are related to the effective transmissivity and flow correlation length.

The polynomial chaos expansion (PCE) technique is used for model reduction since it has demonstrated its effectiveness in different hydrogeological settings. For example, PCE has been used to accelerate variance-based global sensitivity analysis and uncertainty quantification (otherwise prohibitive from a computational point of view) in transport problems in subsurface hydrology (Ciriello & Barros, 2020; Ciriello et al., 2012, 2017; Kim et al., 2022; Oladyshkin et al., 2012) and for a sensitivity-based transport model calibration (Ciriello et al., 2013, 2015). PCE has also been used to compute distribution-based global sensitivity analysis (Ciriello et al., 2019), and to predict the future behavior of water resources under climate change (Focaccia et al., 2021). Merchán-Rivera et al. (2021) used PCE to quantify the uncertainty in surface water-groundwater interaction estimates. The key advantage of
**Writing – review & editing:** Y. Méheust,
V. Di Federico, V. Ciriello

particular flow localization in areas of higher conductance, where the fracture aperture is larger. The resulting hydraulic behavior of the (rough) geological fracture of interest is usually compared to the parallel plate model where both walls are assumed to be perfectly smooth planes, parallel to each other and separated by a distance that is identical to the mechanical aperture (distance between the mean planes of the two rough walls) of the rough fracture. Flow channeling induces either a flow-enhancing or flow-inhibiting hydraulic behavior with respect to the corresponding parallel plate model (Méheust & Schmittbuhl, 2000). However, over a large population of fractures with identical statistical properties, flow-inhibition dominates (Méheust & Schmittbuhl, 2001), but the statistical dispersion (i.e., variability) over the fractures' hydraulic behavior depends on their correlation length (Méheust & Schmittbuhl, 2003).

This study focuses both on Newtonian and shear-thinning (ST) fluids. Subsurface industrial activities widely adopt complex fluids, such as foams (Faroughi et al., 2018; Hosseini et al., 2019), clayey muds, polymer solutions, and aqueous and non-aqueous suspensions, which are ST, that is, whose viscosity tends to decrease under increasing shear rates. CO₂-based ST fluids (Jiang, 2021) have also been recently proposed as a water-saving solution to reduce the use of freshwater in subsurface activities. ST fluids have been demonstrated to enhance the cleanup efficiency (Hosseini et al., 2019) and carrying capability of remedial amendments (Tosco et al., 2012; Zhong et al., 2008).

The effects of ST rheology on the transmissivity of variable aperture fractures were first investigated in the past using simplified deterministic (Di Federico, 1997, 2001) or stochastic (Di Federico, 1998) conceptualizations. More recently, the fracture's hydraulic behavior has been shown to be augmented by the interplay between pore-scale heterogeneity and ST rheology. This potentially leads to fracture transmissivities orders of magnitude larger than the Newtonian counterpart (Lenci, Méheust, et al., 2022). The transmissivity is computed from the ratio of the volumetric flow rate to the applied macroscopic pressure gradient. In numerical simulations, the flow rate is obtained from simulating the fluid velocity field. The Reynolds equation, based on the lubrication approximation (Brown, 1987), is often used to model Newtonian flow in geological fractures. This depth-averaged approach relies on the local validity of the Hagen-Poiseuille law, which holds as long as the conditions for Stokes flow are satisfied and spatial variations of the aperture fields are slow.

Recently, an efficient algorithm has been proposed by Lenci, Putti, et al. (2022) to solve a generalized non-linear Reynolds equation that extends the model to ST fluids whose rheology includes a transition from a Newtonian plateau at low shear rates to a power law ST behavior at large shear rates. This numerical strategy has been adopted to produce a Monte Carlo (MC) study of the ST flow through rough fractures (Lenci, Méheust, et al., 2022), overcoming the computation burden that is typical of stochastic analysis of non-linear problems in spatially varying domains and thus allowing to evidence the impact of the fluid's rheology as discussed above. However, exploring the space of controlling parameters for ST flow in geological fractures over a wide range of values of these parameters, using fracture-scale simulations, would inevitably require a vast number of MCs, and thus a tremendous amount of computational time, due to the dimensionality of the parameter space (consisting of at least three dimensions: fracture closure, fracture correlation length, imposed pressure gradient), the nonlinearity of the model, and the variability of the aperture field's topology over different realizations.

In this study, reduced-order models (ROMs) are thus employed to drastically decrease the computational cost of this statistical framework and to obtain the dependence of given quantities of interest (QoIs) describing a fracture's hydraulic behavior on the controlling parameters in a simple mathematical form, rather than using a pure MC approach as previously done by Lenci, Méheust, et al. (2022). The QoIs, whose response surfaces are reconstructed in this study, are the ensemble average and difference between the third and first quartiles (computed over the number of MCs of the aperture field) of a set of MFMs, which are related to the effective transmissivity and flow correlation length.

The polynomial chaos expansion (PCE) technique is used for model reduction since it has demonstrated its effectiveness in different hydrogeological settings. For example, PCE has been used to accelerate variance-based global sensitivity analysis and uncertainty quantification (otherwise prohibitive from a computational point of view) in transport problems in subsurface hydrology (Ciriello & Barros, 2020; Ciriello et al., 2012, 2017; Kim et al., 2022; Oladyshkin et al., 2012) and for a sensitivity-based transport model calibration (Ciriello et al., 2013, 2015). PCE has also been used to compute distribution-based global sensitivity analysis (Ciriello et al., 2019), and to predict the future behavior of water resources under climate change (Focaccia et al., 2021). Merchán-Rivera et al. (2021) used PCE to quantify the uncertainty in surface water-groundwater interaction estimates. The key advantage of
using PCE is that this technique provides an approximation of the response surface of the QoIs in the form of a polynomial series in the space of variability of the selected governing parameters (Ghanem & Spanos, 1991; Wiener, 1938). To be determined, the PCE requires a negligible number of model runs, thus being advantageous over onerous MCs for mapping out the governing parameter space (Ciriello & Barros, 2020).

Hence, we greatly extend our previous findings on the hydraulic behavior of Newtonian and ST fracture flow, by obtaining a prediction of the generic hydraulic behavior of a geological fracture, under Newtonian and ST flow, which is virtually continuous over the parameter space of practical interest. In addition, we characterize the relationship between the structural correlation length and that of the velocity field; to our knowledge, this has never been investigated for fracture flow, let alone for ST fracture flow. Furthermore, our approach provides both the average behavior over the statistics and the statistical dispersion as a function of the parameters, both for the hydraulic behavior and for the flow correlation length. These findings are only achievable through model reduction, which, to our knowledge, has never been used to address fracture-scale hydraulics.

The paper is organized as follows. In Section 2, the geometry of geological fractures is described and the flow equations, numerical model, and stochastic framework are described; the principle of the ROM is also explained. In Section 3, the results of the analysis are presented and interpreted. We draw conclusions in Section 4.

## 2. Materials and Methods

### 2.1. Problem Statement

*Synthetic geological fractures:* Geological planar fractures consist of two facing rough walls which are self-affine (Power et al., 1987) at all scales (Schmittbuhl et al., 1995). Assuming statistical isotropy of the topographies, this scale invariance manifests in the Fourier transform of the topographies of the upper and the lower wall, $z_u$ and $z_l$ respectively, which scale as (Méheust & Schmittbuhl, 2001)

$$|z_{u||}(\mathbf{k})| \propto |k|^{-1-H}, \quad (1)$$

where $|k| = \sqrt{k_x^2 + k_y^2}$ is the wave number and $H$ is the Hurst exponent that controls the scale invariance. The measured values of $H$ fall in the range [0.5; 0.8] (Boffa et al., 1999; Bouchaud et al., 1990), so that the topographies are planar at large scales, the corresponding asymptotic planes being also the mean planes of the topographies.

Assuming that these two mean planes are parallel to each other, the local fracture aperture is the distance between the two fracture walls, measured along the direction normal to the mean planes. Further assuming that $z_{u||}$ have been defined with a zero mean, the aperture field is defined as

$$a = z_u - z_l + a_m, \quad (2)$$

where the mechanical aperture $a_m$ (Brown, 1987) is the distance between the mean planes of the two walls and is equal to the mean aperture $\langle a \rangle$ if the walls are not in contact. When $a_m$ is set sufficiently low, contact zones result from eventual interpenetrations between the rough walls; in such cases the negative values of $a$ from Equation 2 are set to zero and the mean aperture $\langle a \rangle$ becomes larger than $a_m$. To simplify notations, we shall, in the following, use the notation $\langle a \rangle$ instead of $a_m$; $\langle a \rangle$ is then to be understood as the distance between the mean planes of the fracture walls, and is only strictly equal to the aperture field's mean value when the walls are not in contact.

Fracture walls are generally not entirely uncorrelated with respect to each other: their Fourier modes are identical above a given characteristic length scale which we denote correlation length $L_c$ and assume to be a constitutive parameter of the rock formation, resulting from its history of tectonic activity and weathering. Consequently, the Fourier transform of the aperture field $a$ only possesses the same self-affinity as the walls, expressed by Equation 1, at scales smaller than $L_c$. Depending on the fracture length $L$, this will include the entire relevant scale range (if $L \le L_c$), or not.

A geological fracture's geometry can thus be described from the following statistical parameters: its length $L$, the correlation length $L_c$, its mechanical aperture $a_m$, its standard deviation $\sigma_a$, the Hurst exponent $H$, and a parameter $\omega$ controlling the stochastic nature of the aperture field (e.g., the seed of the random number generator used to generate it numerically). An example of synthetic fracture is represented in Figure 1. The synthetic aperture wall topographies have been generated following the pseudocode described in Lenci, Putti, et al. (2022), which
using PCE is that this technique provides an approximation of the response surface of the QoIs in the form of a polynomial series in the space of variability of the selected governing parameters (Ghanem & Spanos, 1991; Wiener, 1938). To be determined, the PCE requires a negligible number of model runs, thus being advantageous over onerous MCs for mapping out the governing parameter space (Ciriello & Barros, 2020).

Hence, we greatly extend our previous findings on the hydraulic behavior of Newtonian and ST fracture flow, by obtaining a prediction of the generic hydraulic behavior of a geological fracture, under Newtonian and ST flow, which is virtually continuous over the parameter space of practical interest. In addition, we characterize the relationship between the structural correlation length and that of the velocity field; to our knowledge, this has never been investigated for fracture flow, let alone for ST fracture flow. Furthermore, our approach provides both the average behavior over the statistics and the statistical dispersion as a function of the parameters, both for the hydraulic behavior and for the flow correlation length. These findings are only achievable through model reduction, which, to our knowledge, has never been used to address fracture-scale hydraulics.

The paper is organized as follows. In Section 2, the geometry of geological fractures is described and the flow equations, numerical model, and stochastic framework are described; the principle of the ROM is also explained. In Section 3, the results of the analysis are presented and interpreted. We draw conclusions in Section 4.

## 2. Materials and Methods

### 2.1. Problem Statement

*Synthetic geological fractures:* Geological planar fractures consist of two facing rough walls which are self-affine (Power et al., 1987) at all scales (Schmittbuhl et al., 1995). Assuming statistical isotropy of the topographies, this scale invariance manifests in the Fourier transform of the topographies of the upper and the lower wall, $z_u$ and $z_l$ respectively, which scale as (Méheust & Schmittbuhl, 2001)

$$|z_{u||}(\mathbf{k})| \propto |k|^{-1-H}, \quad (1)$$

where $|k| = \sqrt{k_x^2 + k_y^2}$ is the wave number and $H$ is the Hurst exponent that controls the scale invariance. The measured values of $H$ fall in the range [0.5; 0.8] (Boffa et al., 1999; Bouchaud et al., 1990), so that the topographies are planar at large scales, the corresponding asymptotic planes being also the mean planes of the topographies.

Assuming that these two mean planes are parallel to each other, the local fracture aperture is the distance between the two fracture walls, measured along the direction normal to the mean planes. Further assuming that $z_{u||}$ have been defined with a zero mean, the aperture field is defined as

$$a = z_u - z_l + a_m, \quad (2)$$

where the mechanical aperture $a_m$ (Brown, 1987) is the distance between the mean planes of the two walls and is equal to the mean aperture $\langle a \rangle$ if the walls are not in contact. When $a_m$ is set sufficiently low, contact zones result from eventual interpenetrations between the rough walls; in such cases the negative values of $a$ from Equation 2 are set to zero and the mean aperture $\langle a \rangle$ becomes larger than $a_m$. To simplify notations, we shall, in the following, use the notation $\langle a \rangle$ instead of $a_m$; $\langle a \rangle$ is then to be understood as the distance between the mean planes of the fracture walls, and is only strictly equal to the aperture field's mean value when the walls are not in contact.

Fracture walls are generally not entirely uncorrelated with respect to each other: their Fourier modes are identical above a given characteristic length scale which we denote correlation length $L_c$ and assume to be a constitutive parameter of the rock formation, resulting from its history of tectonic activity and weathering. Consequently, the Fourier transform of the aperture field $a$ only possesses the same self-affinity as the walls, expressed by Equation 1, at scales smaller than $L_c$. Depending on the fracture length $L$, this will include the entire relevant scale range (if $L \le L_c$), or not.

A geological fracture's geometry can thus be described from the following statistical parameters: its length $L$, the correlation length $L_c$, its mechanical aperture $a_m$, its standard deviation $\sigma_a$, the Hurst exponent $H$, and a parameter $\omega$ controlling the stochastic nature of the aperture field (e.g., the seed of the random number generator used to generate it numerically). An example of synthetic fracture is represented in Figure 1. The synthetic aperture wall topographies have been generated following the pseudocode described in Lenci, Putti, et al. (2022), which
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Representation of a geological fracture, with definitions of the various relevant geometric quantities and flow boundary conditions. Wall topographies are both generated over a 1,024 × 1,024 nodes mesh with the following parameters: $L/L_c = 8$, $\sigma_d/\langle a \rangle = 0.2$, $L_c = 0.1$ m, $\langle a \rangle = 1$ mm, and $H = 0.8$.</figcaption></figure>

implement the above theoretical geometrical properties. The same pseudocode has been used to generate all the synthetic aperture fields which are considered in this study.

**Fluid rheology:** The ST fluid rheology is modeled with the three-parameter Ellis constitutive law,

$$ \mu = \mu_0 \left[ 1 + \left( \frac{\tau}{\tau_{1/2}} \right)^{\frac{1}{n}-1} \right]^{-1}, \qquad (3) $$

which relates the fluid's dynamics viscosity to the shear stress $\tau$; $\mu_0$ defines the low-shear rate quasi-Newtonian viscosity plateau, $n$ is the ST exponent that defines the power-law trend at large shear rate, and the characteristic shear stress $\tau_{1/2}$ shapes the transition between those two asymptotic behaviors. These rheological parameters can be derived experimentally from rheometric measurements by least-squares regression, or inferred from other constitutive laws as proposed for the Carreau-Yasuda model by Lenci, Putti, et al. (2022).

**Lubrication flow model:** The steady-state isothermal viscous flow of an incompressible Ellis fluid through a geological fracture is modeled using a lubrication-based approach. The local pressure distribution $P$ across the mean fracture plane is governed by the generalized Reynolds equation (Lenci, Putti, et al., 2022):

$$ \nabla \cdot \left[ \left( \frac{a^3}{12} + \frac{n}{2n+1} \left( \frac{a^{2n+1}}{2^{1+n} \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}-1} \| \nabla P \|^{\frac{1}{n}-1} \right) \nabla P \right] = 0, \qquad (4) $$

which holds for Stokes flow conditions (Re $\ll$ 1) as long as the aperture field's gradient remains small ($\|\nabla a\| \ll 1$). For Newtonian fluids, the classical Reynolds equation can be recovered for $n = 1$ or $\tau_{1/2} \to +\infty$. Solving Equation 4 yields the spatial distribution of the pressure ($P$), and the velocity field is then obtained as

$$ \mathbf{u} = \left[ \frac{a^2}{12\mu_0} + \frac{n}{(2n+1)\mu_0} \left( \frac{a^{n+1}}{2^{1+n}\tau_{1/2}^{1-n}} \right)^{\frac{1}{n}-1} \|\nabla P\|^{\frac{1}{n}-1} \right] \nabla P. \qquad (5) $$

**Transmissivity:** The volumetric flow rate $Q$ can then be computed from the velocity field to estimate the apparent transmissivity,

$$ T = \frac{Q\mu_0}{\nabla P}. \qquad (6) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Representation of a geological fracture, with definitions of the various relevant geometric quantities and flow boundary conditions. Wall topographies are both generated over a 1,024 × 1,024 nodes mesh with the following parameters: $L/L_c = 8$, $\sigma_d/\langle a \rangle = 0.2$, $L_c = 0.1$ m, $\langle a \rangle = 1$ mm, and $H = 0.8$.</figcaption></figure>

implement the above theoretical geometrical properties. The same pseudocode has been used to generate all the synthetic aperture fields which are considered in this study.

**Fluid rheology:** The ST fluid rheology is modeled with the three-parameter Ellis constitutive law,

$$ \mu = \mu_0 \left[ 1 + \left( \frac{\tau}{\tau_{1/2}} \right)^{\frac{1}{n}-1} \right]^{-1}, \qquad (3) $$

which relates the fluid's dynamics viscosity to the shear stress $\tau$; $\mu_0$ defines the low-shear rate quasi-Newtonian viscosity plateau, $n$ is the ST exponent that defines the power-law trend at large shear rate, and the characteristic shear stress $\tau_{1/2}$ shapes the transition between those two asymptotic behaviors. These rheological parameters can be derived experimentally from rheometric measurements by least-squares regression, or inferred from other constitutive laws as proposed for the Carreau-Yasuda model by Lenci, Putti, et al. (2022).

**Lubrication flow model:** The steady-state isothermal viscous flow of an incompressible Ellis fluid through a geological fracture is modeled using a lubrication-based approach. The local pressure distribution $P$ across the mean fracture plane is governed by the generalized Reynolds equation (Lenci, Putti, et al., 2022):

$$ \nabla \cdot \left[ \left( \frac{a^3}{12} + \frac{n}{2n+1} \left( \frac{a^{2n+1}}{2^{1+n} \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}-1} \| \nabla P \|^{\frac{1}{n}-1} \right) \nabla P \right] = 0, \qquad (4) $$

which holds for Stokes flow conditions (Re $\ll 1$) as long as the aperture field's gradient remains small ($\|\nabla a\| \ll 1$). For Newtonian fluids, the classical Reynolds equation can be recovered for $n = 1$ or $\tau_{1/2} \to +\infty$. Solving Equation 4 yields the spatial distribution of the pressure ($P$), and the velocity field is then obtained as

$$ \mathbf{u} = \left[ \frac{a^2}{12\mu_0} + \frac{n}{(2n+1)\mu_0} \left( \frac{a^{n+1}}{2^{1+n}\tau_{1/2}^{1-n}} \right)^{\frac{1}{n}-1} \|\nabla P\|^{\frac{1}{n}-1} \right] \nabla P. \qquad (5) $$

**Transmissivity:** The volumetric flow rate $Q$ can then be computed from the velocity field to estimate the apparent transmissivity,

$$ T = \frac{Q\mu_0}{\nabla P}. \qquad (6) $$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Average longitudinal autocorrelation functions of the velocity magnitude for the flow of a Newtonian and a shear-thinning fluid in two fracture realizations, one with $L/L_c = 1$ (a) and one with $L/L_c = 16$ (b). The parameters adopted otherwise for the simulations are: $\sigma_{\text{d}}/\langle a \rangle = 0.5$, $\langle a \rangle = 10^{-3}$ m, $H = 0.8$, and $L_c = 0.1$ m.</figcaption></figure>

Even though for Newtonian flow the transmissivity is an intrinsic property of the fracture's geometry, for ST fluids it depends on both that geometry, the magnitude of the external pressure gradient, and the rheological parameters $n$ and $\tau_{1/2}$.

Hereinafter, the subscript “0” is adopted to denote quantities related to the Newtonian model, while the subscript $\parallel$ stands for quantities estimated for the parallel plate geometry, that is, assuming a uniform aperture field $a = \langle a \rangle$. In this work two different dimensionless formulations of the transmissivity are thus used: $T/T_0$ to study the deviation, due to the ST rheology, from the hydraulic behavior of the fracture under Newtonian flow, and $T/T_{\parallel}$ to analyze the compound effect of fracture heterogeneity and fluid rheology.

In the literature, the transition from the quasi-Newtonian ($T/T_0 \approx 1$) to the ST hydraulic behavior ($T/T_0 > 1$) has been observed to occur when the external pressure gradient reaches the crossover value

$$ \overline{\nabla P_c} = \frac{2\tau_c}{\langle a \rangle}, \qquad (7) $$

where $\tau_c$ is a crossover shear stress obtained from $\tau_c/\tau_{1/2} = \left[ 1 + (\tau_c/\tau_{1/2})^{\frac{1}{n}-1} \right]^{-1}$ (Lenci, Putti, et al., 2022).

**Correlation length of the velocity field:** The finite volume scheme (Lenci, Putti, et al., 2022) provides edge velocities, which can be used to define the cell-centered flow velocity by bilinear interpolation scheme, as proposed by Pollock (1988). Considering a regular partition of the domain, the longitudinal and transverse components of the flow velocity located at the center of each mesh cell are obtained by arithmetic averaging of the edge velocities at opposite cell faces. The correlation length of the velocity magnitude $\ell$ is then analyzed in relation to the correlation length of the aperture field $\lambda$. Both correlation lengths are defined as the first zero-crossing lag of the autocorrelation function (ACF), estimated by averaging the ACFs of all longitudinal profiles (respectively of the aperture and velocity fields) of the fracture, and discarding the contact zones. Figure 2 provides examples of such average ACFs of the velocity magnitude for the flows of a Newtonian fluid and a ST fluid in two fracture realizations, one with $L/L_c = 1$ and one with $L/L_c = 16$.

Note that $\lambda$ is related to $L_c$, which is the upper limit scale for the self-affinity of the aperture field and is imposed when generating the geometry, but the former is proportional to the latter. Indeed, when a random fractal is truncated beyond a given length scale, say $L_c$, the correlation length of the resulting field is proportional to $L_c$, as demonstrated numerically by Cintoli et al. (2005). The proportionality factor between $\lambda$ and $L_c$ depends on $L_c$ but not on the fracture closure. More details are provided in Appendix A.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Average longitudinal autocorrelation functions of the velocity magnitude for the flow of a Newtonian and a shear-thinning fluid in two fracture realizations, one with $L/L_c = 1$ (a) and one with $L/L_c = 16$ (b). The parameters adopted otherwise for the simulations are: $\sigma_{\text{d}}/\langle a \rangle = 0.5$, $\langle a \rangle = 10^{-3}$ m, $H = 0.8$, and $L_c = 0.1$ m.</figcaption></figure>

Even though for Newtonian flow the transmissivity is an intrinsic property of the fracture's geometry, for ST fluids it depends on both that geometry, the magnitude of the external pressure gradient, and the rheological parameters $n$ and $\tau_{1/2}$.

Hereinafter, the subscript “0” is adopted to denote quantities related to the Newtonian model, while the subscript $\parallel$ stands for quantities estimated for the parallel plate geometry, that is, assuming a uniform aperture field $a = \langle a \rangle$. In this work two different dimensionless formulations of the transmissivity are thus used: $T/T_0$ to study the deviation, due to the ST rheology, from the hydraulic behavior of the fracture under Newtonian flow, and $T/T_{\parallel}$ to analyze the compound effect of fracture heterogeneity and fluid rheology.

In the literature, the transition from the quasi-Newtonian ($T/T_0 \approx 1$) to the ST hydraulic behavior ($T/T_0 > 1$) has been observed to occur when the external pressure gradient reaches the crossover value

$$ \overline{\nabla P_c} = \frac{2\tau_c}{\langle a \rangle}, \qquad (7) $$

where $\tau_c$ is a crossover shear stress obtained from $\tau_c/\tau_{1/2} = \left[ 1 + (\tau_c/\tau_{1/2})^{\frac{1}{n}-1} \right]^{-1}$ (Lenci, Putti, et al., 2022).

**Correlation length of the velocity field:** The finite volume scheme (Lenci, Putti, et al., 2022) provides edge velocities, which can be used to define the cell-centered flow velocity by bilinear interpolation scheme, as proposed by Pollock (1988). Considering a regular partition of the domain, the longitudinal and transverse components of the flow velocity located at the center of each mesh cell are obtained by arithmetic averaging of the edge velocities at opposite cell faces. The correlation length of the velocity magnitude $\ell$ is then analyzed in relation to the correlation length of the aperture field $\lambda$. Both correlation lengths are defined as the first zero-crossing lag of the autocorrelation function (ACF), estimated by averaging the ACFs of all longitudinal profiles (respectively of the aperture and velocity fields) of the fracture, and discarding the contact zones. Figure 2 provides examples of such average ACFs of the velocity magnitude for the flows of a Newtonian fluid and a ST fluid in two fracture realizations, one with $L/L_c = 1$ and one with $L/L_c = 16$.

Note that $\lambda$ is related to $L_c$, which is the upper limit scale for the self-affinity of the aperture field and is imposed when generating the geometry, but the former is proportional to the latter. Indeed, when a random fractal is truncated beyond a given length scale, say $L_c$, the correlation length of the resulting field is proportional to $L_c$, as demonstrated numerically by Cintoli et al. (2005). The proportionality factor between $\lambda$ and $L_c$ depends on $L_c$ but not on the fracture closure. More details are provided in Appendix A.
## 2.2. Numerical Method

Stochastic apertures fields defined by the set of parameters ($L, L_c, a_m, \sigma_a, H, \omega$) are generated using a random generator and shaping the Fourier transform into the product of a power law cone (see Equation 1) and of the Fourier transform of a random noise (see Lenci, Putti, et al., 2022; Méheust & Schmittbuhl, 2003).

We adopt a depth-averaged approach to solve Equation 4 over a 2D square domain of the mean fracture plane. Permanent flow occurs under the effect of a longitudinal pressure drop $\Delta P$ imposed through Dirichlet boundary conditions on the eastern and western boundaries and “no-flow-through-the-boundary” conditions imposed on the northern and southern boundaries. The latter are Neumann boundary conditions that stipulate a zero pressure gradient across the boundary; they are well suited to simulating fractures with a finite lateral extent and are particularly relevant to simulating geological fracture flow. The other type of boundary conditions that could be applied on the lateral boundary are periodic boundary conditions, which are much less relevant to real fractures.

Equation 4 is solved on a $2^{10} \times 2^{10}$ grid with the finite volume scheme proposed by Lenci, Putti, et al. (2022), whose implementation relies on an inexact preconditioned conjugate gradient (PCG) Newton-Krylov method combined with a parameter continuation strategy (Lenci, Putti, et al., 2022). The mesh size is selected based on previous studies. In particular, the nonlinearity of the constitutive law results in a higher computational expense as compared to Newtonian (and hence, Darcian) flow, so we choose a convenient value for the grid size ($2^{10} = 1,024$) to keep an adequate level of detail while allowing to perform the required number of simulations in a reasonable time. Note that larger modes of the fracture aperture field dominate smaller modes in controlling the flow heterogeneity; as such, further grid refinement is unnecessary because it wouldn't impact the predicted hydraulic behavior.

Here, we have replaced the PCG method with the aggregation-based algebraic multigrid method (Napov & Notay, 2012; Notay, 2012), which performs better in achieving global convergence. As already assessed by Lenci, Méheust, et al. (2022), this highly efficient scheme allows a stochastic analysis of the phenomenon via MCs, resulting in a reasonable amount of computational time.

## 2.3. Stochastic Framework

The MFMs which we evaluate with the model are: (a) the fracture transmissivity ($T$) and (b) the flow correlation length ($\ell$). In particular, the former is expressed by means of the ratios $T_0/T_{0\parallel}$ and $T/T_{\parallel}$; these are used to measure the effect of the aperture field's heterogeneity on the transmissivity with respect to Newtonian and ST fluids, respectively. The parallel plate model, characterized by parallel planar walls and thus a uniform aperture $a = \langle a \rangle$, is taken as the homogeneous reference case for the transmissivity. As for the flow's characteristic length scale, it is relevant for upscaling transport in heterogenous media (Comolli et al., 2019), where the medium's correlation length is generally used as a surrogate to relate macroscopic solute transport to flow conditions (Dentz et al., 2016). In particular, the flow's correlation length has been adopted as the only parameter present in Markov models for predicting the evolution of the spatially sampled Lagrangian velocities in Bernoulli and Ornstein-Uhlenbeck processes (Dentz et al., 2016; Kang et al., 2017; Morales et al., 2017). However, the incomplete characterization of the flow's correlation length in both porous and fractured media, caused it to be replaced with the medium's correlation length, which is supposed to be of the same order of magnitude. Here, we focus on the ratios between the flow's and medium's correlation lengths for Newtonian and ST fluids, $\ell_0/\lambda$ and $\ell/\lambda$ respectively, to quantify how they differ under different geometrical, rheological, and external forcing conditions. Indeed, if for heterogeneous porous media the two characteristic lengths diverge significantly (Cvetkovic et al., 1996), little has been investigated for fractured media. Finally, to explore the transition from the quasi-Darcian to the non-linear ST regime over the entire parameter space, the ratio between the transmissivities of the two fluids, $T/T_0$, is also analyzed for different flow conditions and degrees of heterogeneity. This allows us to extend the results obtained by Lenci, Méheust, et al. (2022), which were constrained to a limited set of scenarios.

Once defined the set of MFMs in this study, that is, $\{T_0/T_{0\parallel}, T/T_{\parallel}, \ell_0/\lambda, \ell/\lambda, T/T_0\}$, the QoIs are key statistics of these selected metrics, computed over 100 MCs of the aperture field. We focus, in particular, on the ensemble average (hereinafter this will be indicated directly with the name (·) of each MFM for the sake of brevity) and the difference between the third and first quartiles (hereinafter $\Delta q_{(\cdot)}$). Among these QoIs, $T_0/T_{0\parallel}$ and $\ell_0/\lambda$ (and the corresponding quartiles) are functions only of the set of parameters dealing with the heterogeneity:
## 2.2. Numerical Method

Stochastic apertures fields defined by the set of parameters ($L, L_c, a_m, \sigma_a, H, \omega$) are generated using a random generator and shaping the Fourier transform into the product of a power law cone (see Equation 1) and of the Fourier transform of a random noise (see Lenci, Putti, et al., 2022; Méheust & Schmittbuhl, 2003).

We adopt a depth-averaged approach to solve Equation 4 over a 2D square domain of the mean fracture plane. Permanent flow occurs under the effect of a longitudinal pressure drop $\Delta P$ imposed through Dirichlet boundary conditions on the eastern and western boundaries and “no-flow-through-the-boundary” conditions imposed on the northern and southern boundaries. The latter are Neumann boundary conditions that stipulate a zero pressure gradient across the boundary; they are well suited to simulating fractures with a finite lateral extent and are particularly relevant to simulating geological fracture flow. The other type of boundary conditions that could be applied on the lateral boundary are periodic boundary conditions, which are much less relevant to real fractures.

Equation 4 is solved on a $2^{10} \times 2^{10}$ grid with the finite volume scheme proposed by Lenci, Putti, et al. (2022), whose implementation relies on an inexact preconditioned conjugate gradient (PCG) Newton-Krylov method combined with a parameter continuation strategy (Lenci, Putti, et al., 2022). The mesh size is selected based on previous studies. In particular, the nonlinearity of the constitutive law results in a higher computational expense as compared to Newtonian (and hence, Darcian) flow, so we choose a convenient value for the grid size ($2^{10} = 1,024$) to keep an adequate level of detail while allowing to perform the required number of simulations in a reasonable time. Note that larger modes of the fracture aperture field dominate smaller modes in controlling the flow heterogeneity; as such, further grid refinement is unnecessary because it wouldn't impact the predicted hydraulic behavior.

Here, we have replaced the PCG method with the aggregation-based algebraic multigrid method (Napov & Notay, 2012; Notay, 2012), which performs better in achieving global convergence. As already assessed by Lenci, Méheust, et al. (2022), this highly efficient scheme allows a stochastic analysis of the phenomenon via MCs, resulting in a reasonable amount of computational time.

## 2.3. Stochastic Framework

The MFMs which we evaluate with the model are: (a) the fracture transmissivity ($T$) and (b) the flow correlation length ($\ell$). In particular, the former is expressed by means of the ratios $T_0/T_{0\parallel}$ and $T/T_{\parallel}$; these are used to measure the effect of the aperture field's heterogeneity on the transmissivity with respect to Newtonian and ST fluids, respectively. The parallel plate model, characterized by parallel planar walls and thus a uniform aperture $a = \langle a \rangle$, is taken as the homogeneous reference case for the transmissivity. As for the flow's characteristic length scale, it is relevant for upscaling transport in heterogenous media (Comolli et al., 2019), where the medium's correlation length is generally used as a surrogate to relate macroscopic solute transport to flow conditions (Dentz et al., 2016). In particular, the flow's correlation length has been adopted as the only parameter present in Markov models for predicting the evolution of the spatially sampled Lagrangian velocities in Bernoulli and Ornstein-Uhlenbeck processes (Dentz et al., 2016; Kang et al., 2017; Morales et al., 2017). However, the incomplete characterization of the flow's correlation length in both porous and fractured media, caused it to be replaced with the medium's correlation length, which is supposed to be of the same order of magnitude. Here, we focus on the ratios between the flow's and medium's correlation lengths for Newtonian and ST fluids, $\ell_0/\lambda$ and $\ell/\lambda$ respectively, to quantify how they differ under different geometrical, rheological, and external forcing conditions. Indeed, if for heterogeneous porous media the two characteristic lengths diverge significantly (Cvetkovic et al., 1996), little has been investigated for fractured media. Finally, to explore the transition from the quasi-Darcian to the non-linear ST regime over the entire parameter space, the ratio between the transmissivities of the two fluids, $T/T_0$, is also analyzed for different flow conditions and degrees of heterogeneity. This allows us to extend the results obtained by Lenci, Méheust, et al. (2022), which were constrained to a limited set of scenarios.

Once defined the set of MFMs in this study, that is, $\{T_0/T_{0\parallel}, T/T_{\parallel}, \ell_0/\lambda, \ell/\lambda, T/T_0\}$, the QoIs are key statistics of these selected metrics, computed over 100 MCs of the aperture field. We focus, in particular, on the ensemble average (hereinafter this will be indicated directly with the name (·) of each MFM for the sake of brevity) and the difference between the third and first quartiles (hereinafter $\Delta q_{(\cdot)}$). Among these QoIs, $T_0/T_{0\parallel}$ and $\ell_0/\lambda$ (and the corresponding quartiles) are functions only of the set of parameters dealing with the heterogeneity:
p₁ = {L, L_c, σ_a, ⟨a⟩, H}; while the remaining QoIs depend on p₁, √P and p₂ = {τ₁/₂, n}, the latter collecting the fluid rheological parameters.

Among the parameters the QoIs depend on, we set the following:

* L/L_c ∈ {1; 4; 16}, assuming a constant value of the characteristic length used in the fracture generator L_c = 0.1 m to account for different fracture sizes (L ∈ {0.1 m; 0.4 m; 1.6 m}) (Lenci, Méheust, et al., 2022),

* H = 0.8, the typical value for basalt and granite rock (Bouchaud et al., 1990),

* The rheological parameters for a ST fluid (Xanthan gum) to n = 0.36, μ₀ = 4.42 Pa s, and τ₁/₂ = 0.199 Pa, and for a Newtonian fluid the dynamic viscosity to μ = 4.42 Pa s (Lenci, Putti, et al., 2022);

While we span continuously the QoIs over a proper range of variation of the remaining parameters by modeling them as uniformly distributed random variables. Specifically, we consider:

* The fracture closure σ_a/⟨a⟩ ∈ [0; 1], with ⟨a⟩ = 0.001 m,

* The dimensionless pressure gradient √P/√P_c ~ [0.1; 10] to consider an interval ranging from naturally occurring groundwater head gradients to industrial operating conditions (Jung, 1989).

In order to be able to explore such a vast parameter space, instead of relying solely on the MC simulations of the numerical flow model, we apply model reduction as reported in the next section.

## 2.4. Reduced-Order Models

In order to reconstruct the dependence of the QoIs on the selected parameters, that is, the topography of these QoIs in the selected random parameter space, continuously and at a negligible computational cost, ROMs are used according to the PCE theory.

Let Ω denote the generic QoI of this study. Let f be the high-fidelity model (HFM) described in Sections 2.1 and 2.2, which provides Ω as a function of the aforementioned set of parameters (see Section 2.3), M of them being random and collected in vector p: Ω = f(p), with dim(p) = M. If the variance of Ω is finite, the approximation of the HFM provided by the PCE is admitted (Ghanem & Spanos, 1991; Wiener, 1938; Xiu & Karniadakis, 2002) as

$$ \hat{\Omega} = \sum_{a \in \mathbb{N}^M} s_a \Psi_a(p). \qquad (8) $$

In Equation 8, multi-indices a = {a₁, ..., a_M} ∈ N^M are associated with multivariate polynomials Ψ_a of degree |a| = Σ_M l_i, which constitute an orthonormal basis with respect to the joint probability density function (PDF) of p (Xiu & Karniadakis, 2002); while the PCE coefficients s_a are deterministic coordinates of the spectral decomposition (Ghanem & Spanos, 1991). In practice equation 8 is truncated as

$$ \hat{\Omega} = \sum_{i=0}^{P-1} s_i \Psi_i(p), \quad P = \frac{(M+q)!}{M!q!}, \qquad (9) $$

where q is the maximum degree of the expansion, that is, |a| ≤ q for all a ∈ N^M.

The PCE coefficients are computed through a regression-based approach consisting in the minimization of the variance of the residual, ε = |Ω̂ - Ω|, with respect to the s_i; the set of regression points in the random parameter space is provided by the probabilistic collocation method (PCM) (Webster et al., 1996).

In our study,

$$ \Omega \in \left\{ T_0/T_{0\parallel}, T/T_{\parallel}, \ell_0/\lambda, \ell/\lambda, T/T_0, \Delta q_{T/\parallel}, \Delta q_{\lambda/\parallel}, \Delta q_{T_0/\parallel} \right\} \text{ and } p = \left\{ \sigma_a/\langle a \rangle, \sqrt{P}/\sqrt{P_c} \right\}, \qquad (10) $$

with M = 2 and q = 4, thus resulting in a number of HFM runs equal to P = 15 according to the PCM (in place of a full MC series in the parameter space). Since the parameters in p are distributed uniformly, the Legendre polynomials are used as the basis in the ROMs given by Equation 9 (Xiu & Karniadakis, 2002).
p₁ = {L, L_c, σ_a, ⟨a⟩, H}; while the remaining QoIs depend on p₁, √P and p₂ = {τ₁/₂, n}, the latter collecting the fluid rheological parameters.

Among the parameters the QoIs depend on, we set the following:

* L/L_c ∈ {1; 4; 16}, assuming a constant value of the characteristic length used in the fracture generator L_c = 0.1 m to account for different fracture sizes (L ∈ {0.1 m; 0.4 m; 1.6 m}) (Lenci, Méheust, et al., 2022),

* H = 0.8, the typical value for basalt and granite rock (Bouchaud et al., 1990),

* The rheological parameters for a ST fluid (Xanthan gum) to n = 0.36, μ₀ = 4.42 Pa s, and τ₁/₂ = 0.199 Pa, and for a Newtonian fluid the dynamic viscosity to μ = 4.42 Pa s (Lenci, Putti, et al., 2022);

While we span continuously the QoIs over a proper range of variation of the remaining parameters by modeling them as uniformly distributed random variables. Specifically, we consider:

* The fracture closure σ_a/⟨a⟩ ∈ [0; 1], with ⟨a⟩ = 0.001 m,

* The dimensionless pressure gradient √P/√P_c ~ [0.1; 10] to consider an interval ranging from naturally occurring groundwater head gradients to industrial operating conditions (Jung, 1989).

In order to be able to explore such a vast parameter space, instead of relying solely on the MC simulations of the numerical flow model, we apply model reduction as reported in the next section.

## 2.4. Reduced-Order Models

In order to reconstruct the dependence of the QoIs on the selected parameters, that is, the topography of these QoIs in the selected random parameter space, continuously and at a negligible computational cost, ROMs are used according to the PCE theory.

Let Ω denote the generic QoI of this study. Let f be the high-fidelity model (HFM) described in Sections 2.1 and 2.2, which provides Ω as a function of the aforementioned set of parameters (see Section 2.3), M of them being random and collected in vector p: Ω = f(p), with dim(p) = M. If the variance of Ω is finite, the approximation of the HFM provided by the PCE is admitted (Ghanem & Spanos, 1991; Wiener, 1938; Xiu & Karniadakis, 2002) as

$$ \hat{\Omega} = \sum_{a \in \mathbb{N}^M} s_a \Psi_a(p). \qquad (8) $$

In Equation 8, multi-indices a = {a₁, ..., a_M} ∈ N^M are associated with multivariate polynomials Ψ_a of degree |a| = Σ_M l_i, which constitute an orthonormal basis with respect to the joint probability density function (PDF) of p (Xiu & Karniadakis, 2002); while the PCE coefficients s_a are deterministic coordinates of the spectral decomposition (Ghanem & Spanos, 1991). In practice equation 8 is truncated as

$$ \hat{\Omega} = \sum_{i=0}^{P-1} s_i \Psi_i(p), \quad P = \frac{(M+q)!}{M!q!}, \qquad (9) $$

where q is the maximum degree of the expansion, that is, |a| ≤ q for all a ∈ N^M.

The PCE coefficients are computed through a regression-based approach consisting in the minimization of the variance of the residual, ε = |Ω̂ - Ω|, with respect to the s_i; the set of regression points in the random parameter space is provided by the probabilistic collocation method (PCM) (Webster et al., 1996).

In our study,

$$ \Omega \in \left\{ T_0/T_{0\parallel}, T/T_{\parallel}, \ell_0/\lambda, \ell/\lambda, T/T_0, \Delta q_{T/\parallel}, \Delta q_{\lambda/\parallel}, \Delta q_{T_0/\parallel} \right\} \text{ and } p = \left\{ \sigma_a/\langle a \rangle, \sqrt{P}/\sqrt{P_c} \right\}, \qquad (10) $$

with M = 2 and q = 4, thus resulting in a number of HFM runs equal to P = 15 according to the PCM (in place of a full MC series in the parameter space). Since the parameters in p are distributed uniformly, the Legendre polynomials are used as the basis in the ROMs given by Equation 9 (Xiu & Karniadakis, 2002).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Impact of aperture heterogeneity on the ratio $T_0/T_{0\parallel}$ and normalized flow correlation length along the direction of the external pressure gradient, $\ell_0/\lambda$—panels (a, b) respectively. Parameters that describe the aperture heterogeneity are the fracture closure $\sigma_a/\langle a \rangle$ (x-axes) and the ratio $L/L_c$. Different colors refer to different values of $L/L_c$: yellow is 1, orange is 4, and 16 for purple. Solid lines are the polynomial chaos expansion estimations, circles denote data from the MC simulations (among which black edges denote ensemble averages), and bars indicate the first and third quantiles.</figcaption></figure>

# 3. Results and Discussion

## 3.1. Newtonian Flow

We first focus the analysis on the QoIs for Newtonian fluids, $T_0/T_{0\parallel}$ and $\ell_0/\lambda$, depicted in the left and right panels of Figure 3, respectively. Figure 3 shows, for different correlation lengths ($L/L_c \in \{1; 4; 16\}$) distinguishable by color variations, the MC ensembles provided by the PCM and adopted to generate the predictions of QoIs (continuous lines) as a function of $\sigma_a/\langle a \rangle$ in order to calibrate the PCE. In particular, each MC ensemble is provided with the first and third quartile (bars), 100 single realizations (dots with colored edges), and the relative average (circle symbols with black edges). In Figure 3a, single realizations exhibit two opposite hydraulic behaviors referring to the parallel plate model transmissivity: flow-inhibiting if $T_0/T_{0\parallel} < 1$, or flow-enhancing if $T_0/T_{0\parallel} > 1$. When the fluid rheology is Newtonian, the transmissivity is an implicit feature of the fracture geometry, hence these behaviors are solely affected by the morphology of the aperture field. The flow-inhibiting behavior is dominant in the statistics and promoted by higher fracture closure, as it induces major viscous energy losses due to roughness. Obstructions, such as contact areas, hinder fluid motion and increase flow path tortuosity. On the other hand, the presence of areas with large aperture aligned parallel to the external pressure gradient is associated with a flow-enhancing behavior, although this geometrical feature is not favored statistically (Lenci, Putti, et al., 2022; Méheust & Schmittbuhl, 2001). The effect of the fracture closure on the transmissivity is similar among different values of $L/L_c$. However, small fractures (i.e., $L \approx L_c$), present a higher statistical dispersion around the ensemble average with several outliers displaying prominent flow-enhancing behavior.

Figure 3b shows the same plots as Figure 3a but for the ratio $\ell_0/\lambda$ of the flow correlation length to the geometrical correlation length. To our knowledge, $\ell$ (including $\ell_0$ for Newtonian flow) has hardly been investigated in the literature. Neither the ensemble average nor the width of the confidence interval are observed to be significantly sensitive to fracture closure, hence $L/L_c$ is the most significant parameter that substantially affects $\ell_0/\lambda$. In particular, the higher the ratio $L/L_c$, the higher the ensemble average and the statistical dispersion.

## 3.2. Shear-Thinning Flow

**Comparison to Newtonian transmissivity:** We now extend the analysis to complex fluids that may develop a ST behavior. Figure 4 shows the dependency on fracture closure and dimensionless pressure gradient of the PCE-based predictions of QoIs: the ensemble average of the ratio of the transmissivity to the corresponding Newtonian transmissivity, $T/T_0$, and the corresponding amplitude of the confidence interval, $\Delta q_T/T_0$. The former is adopted to quantify the departure of the ST hydraulic behavior from the Newtonian transmissivity, while
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Impact of aperture heterogeneity on the ratio $T_0/T_{0\parallel}$ and normalized flow correlation length along the direction of the external pressure gradient, $\ell_0/\lambda$—panels (a, b) respectively. Parameters that describe the aperture heterogeneity are the fracture closure $\sigma_a/\langle a \rangle$ (x-axes) and the ratio $L/L_c$. Different colors refer to different values of $L/L_c$: yellow is 1, orange is 4, and 16 for purple. Solid lines are the polynomial chaos expansion estimations, circles denote data from the MC simulations (among which black edges denote ensemble averages), and bars indicate the first and third quantiles.</figcaption></figure>

# 3. Results and Discussion

## 3.1. Newtonian Flow

We first focus the analysis on the QoIs for Newtonian fluids, $T_0/T_{0\parallel}$ and $\ell_0/\lambda$, depicted in the left and right panels of Figure 3, respectively. Figure 3 shows, for different correlation lengths ($L/L_c \in \{1; 4; 16\}$) distinguishable by color variations, the MC ensembles provided by the PCM and adopted to generate the predictions of QoIs (continuous lines) as a function of $\sigma_a/\langle a \rangle$ in order to calibrate the PCE. In particular, each MC ensemble is provided with the first and third quartile (bars), 100 single realizations (dots with colored edges), and the relative average (circle symbols with black edges). In Figure 3a, single realizations exhibit two opposite hydraulic behaviors referring to the parallel plate model transmissivity: flow-inhibiting if $T_0/T_{0\parallel} < 1$, or flow-enhancing if $T_0/T_{0\parallel} > 1$. When the fluid rheology is Newtonian, the transmissivity is an implicit feature of the fracture geometry, hence these behaviors are solely affected by the morphology of the aperture field. The flow-inhibiting behavior is dominant in the statistics and promoted by higher fracture closure, as it induces major viscous energy losses due to roughness. Obstructions, such as contact areas, hinder fluid motion and increase flow path tortuosity. On the other hand, the presence of areas with large aperture aligned parallel to the external pressure gradient is associated with a flow-enhancing behavior, although this geometrical feature is not favored statistically (Lenci, Putti, et al., 2022; Méheust & Schmittbuhl, 2001). The effect of the fracture closure on the transmissivity is similar among different values of $L/L_c$. However, small fractures (i.e., $L \approx L_c$), present a higher statistical dispersion around the ensemble average with several outliers displaying prominent flow-enhancing behavior.

Figure 3b shows the same plots as Figure 3a but for the ratio $\ell_0/\lambda$ of the flow correlation length to the geometrical correlation length. To our knowledge, $\ell$ (including $\ell_0$ for Newtonian flow) has hardly been investigated in the literature. Neither the ensemble average nor the width of the confidence interval are observed to be significantly sensitive to fracture closure, hence $L/L_c$ is the most significant parameter that substantially affects $\ell_0/\lambda$. In particular, the higher the ratio $L/L_c$, the higher the ensemble average and the statistical dispersion.

## 3.2. Shear-Thinning Flow

**Comparison to Newtonian transmissivity:** We now extend the analysis to complex fluids that may develop a ST behavior. Figure 4 shows the dependency on fracture closure and dimensionless pressure gradient of the PCE-based predictions of QoIs: the ensemble average of the ratio of the transmissivity to the corresponding Newtonian transmissivity, $T/T_0$, and the corresponding amplitude of the confidence interval, $\Delta q_T/T_0$. The former is adopted to quantify the departure of the ST hydraulic behavior from the Newtonian transmissivity, while
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Dependence of the fracture transmissivity normalized by the corresponding transmissivity for Newtonian flow, *T/T*₀, on both the dimensionless external pressure gradient $\nabla P/\nabla P_c$ and the fracture closure $\sigma_d/\langle a \rangle$, produced for three different *L/L*<sub>c</sub> ratios: 1, 4, and 16, illustrated respectively in panels (a-c). For each of these panels, the inset within the left side figure shows the topography of *T/T*₀ in the ($\nabla P/\nabla P_c$, $\sigma_d/\langle a \rangle$) parameter space, while the main plots are 1D profiles extracted from this topography at given closures, the black dashed line denoting the parallel plate model transmissivity *T/T*₀ and the solid black line the asymptotic power law trend; the right hand side smaller figure is a plot similar to the aforementioned inset, but showing the amplitude of the confidence interval.</figcaption></figure>

the latter measures the statistical dispersion around the ensemble average. The panels report the model's response across the parameter space for three values of the correlation length ($L/L_c \in \{1; 4; 16\}$), illustrated respectively in panels 4a-b-c. Each panel shows a color map of the dependence of $T/T_0$ on $\log(\nabla P/\nabla P_c)$ and the fracture closure $\sigma_d/\langle a \rangle$, a similar color map for the corresponding confidence interval amplitude, and plots of $\ln(T/T_0)$ as a function of $\log(\overline{\nabla P}/\overline{\nabla P_c})$ at several fracture closures ($\sigma_a/\langle a \rangle \in \{0.1; 0.5; 0.7; 0.9\}$).

Under low external pressure gradients, the Ellis fluid exhibits a quasi-Darcian behavior ($T \propto \overline{\nabla P}$), as the transmissivity is comparable with the Newtonian transmissivity ($T \approx T_0$). This results from low shear rates across the flow domain remaining under the crossover value below which the Ellis fluid does not develop a relevant ST behavior ($\mu \approx \mu_0$, low-shear rate plateau). Consequently, the hydraulic behavior is controlled solely by medium heterogeneity, as for Newtonian fluids (Figure 3). Increasing external pressure gradients triggers a transition to a non-linear regime, which is present in all panels of Figure 4 and exhibits a power law trend of slope $\frac{1}{n} - 1$ when the dimensionless pressure gradient is sufficiently higher than 1. The power law trend of $T/T_0$ is weakly affected by the fracture closure and the correlation length, although heterogeneity contributes to the regulation of the transition between the two regimes. In particular, the establishment of the non-linear regime due to the fluid's ST behavior, as well as the statistical dispersion around the ensemble average, are promoted by higher fracture closures and lower values of $L/L_c$, while $\overline{\nabla P}/\overline{\nabla P_c}$ is the control parameter of the transition. Indeed, for a given $\overline{\nabla P}/\overline{\nabla P_c}$, higher closures mean a higher degree of flow channeling, which results in a wider range of hydraulic behaviors (from flow-inhibiting to flow-enhancing).

**Comparison to parallel plate transmissivity:** Figure 5 depicts the PCE-predicted behavior of the ratio of the fracture transmissivity *T* to that of the corresponding parallel plate geometry, *T*<sub>‖</sub>, and of the amplitude of the relative confidence interval $\Delta q_{\overline{T}/T_0}$. The figure's layout is otherwise identical to that of Figure 4. Again, for low dimensionless pressure gradients the hydraulic behavior tends to the Newtonian transmissivity (see Figure 3a), as *T* → *T*₀ and *T*<sub>‖</sub> → *T*<sub>0‖</sub>. For sufficiently high dimensionless pressure gradient, indicatively $\overline{\nabla P}/\overline{\nabla P_c} > 3$ (a transition roughly occurs between 0.1—quasi-Newtonian behavior—and 4), almost horizontal contour lines are seen in the $T/T_0$ topography, indicating that the fracture closure's impact on the transmissivity scales equivalently in the rough-walled fracture and the corresponding parallel plate geometry. This is all the truer as the correlation length $L_c$ is larger. Moreover, as depicted in the maps on the right-hand side, the fracture closure does not significantly impact the statistical dispersion of $T/T_0$ when the flow is strongly ST. Furthermore, when the ST hydraulic behavior is developed, the observed relation between the transmissivity and

the fracture closure is non-monotonic, for any value of $L/L_c$. The ST nature of the fluids mitigates the viscous energy losses induced by the higher heterogeneity due to fracture closure, increasing the statistical likelihood of flow-enhancing hydraulic behavior. As the aperture field becomes more heterogeneous, both the transmissivity and statistical dispersion increase, with many realizations characterized by a transmissivity that is orders of magnitude higher than that of the corresponding parallel plate case. This can be explained by the formation of channels where the fluid shear-thins due to high shear rates, so that the resulting apparent viscosity decreases and velocities increase, resulting in preferential flow paths of very high conductance due to minor viscous energy losses. The higher statistical dispersion observed when fracture closure increases and $L/L_c = 1$ depends on the more chaotic nature of the medium, which then affects the flow by generating larger channels.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dependence of the fracture transmissivity normalized by the corresponding transmissivity for Newtonian flow, *T/T*₀, on both the dimensionless external pressure gradient ∇*P*/∇*P*<sub>c</sub> and the fracture closure *σ*<sub>d</sub>/*⟨*a**⟩*, produced for three different *L/L*<sub>c</sub> ratios: 1, 4, and 16, illustrated respectively in panels (a-c). For each of these panels, the inset within the left side figure shows the topography of *T/T*₀ in the (∇*P*/∇*P*<sub>c</sub>, *σ*<sub>d</sub>/*⟨*a**⟩*) parameter space, while the main plots are 1D profiles extracted from this topography at given closures, the black dashed line denoting the parallel plate model transmissivity *T/T*₀ and the solid black line the asymptotic power law trend; the right hand side smaller figure is a plot similar to the aforementioned inset, but showing the amplitude of the confidence interval.</figcaption></figure>

the latter measures the statistical dispersion around the ensemble average. The panels report the model's response across the parameter space for three values of the correlation length (*L*/*L*<sub>c</sub> ∈ {1; 4; 16}), illustrated respectively in panels 4a-b-c. Each panel shows a color map of the dependence of *T/T*₀ on log(∇*P*/∇*P*<sub>c</sub>) and the fracture closure *σ*<sub>d</sub>/*⟨*a**⟩*, a similar color map for the corresponding confidence interval amplitude, and plots of ln(*T*/*T*₀) as a function of log(√*P*/√*P*<sub>c</sub>) at several fracture closures (*σ*<sub>a</sub>/*⟨*a**⟩* ∈ {0.1; 0.5; 0.7; 0.9}*).

Under low external pressure gradients, the Ellis fluid exhibits a quasi-Darcian behavior (T ~ √P), as the transmissivity is comparable with the Newtonian transmissivity (T ≈ T₀). This results from low shear rates across the flow domain remaining under the crossover value below which the Ellis fluid does not develop a relevant ST behavior (μ ≈ μ₀, low-shear rate plateau). Consequently, the hydraulic behavior is controlled solely by medium heterogeneity, as for Newtonian fluids (Figure 3). Increasing external pressure gradients triggers a transition to a non-linear regime, which is present in all panels of Figure 4 and exhibits a power law trend of slope 1/n - 1 when the dimensionless pressure gradient is sufficiently higher than 1. The power law trend of *T/T*₀ is weakly affected by the fracture closure and the correlation length, although heterogeneity contributes to the regulation of the transition between the two regimes. In particular, the establishment of the non-linear regime due to the fluid's ST behavior, as well as the statistical dispersion around the ensemble average, are promoted by higher fracture closures and lower values of *L/L*<sub>c</sub>, while √P/√P<sub>c</sub> is the control parameter of the transition. Indeed, for a given √P/√P<sub>c</sub>, higher closures mean a higher degree of flow channeling, which results in a wider range of hydraulic behaviors (from flow-inhibiting to flow-enhancing).

**Comparison to parallel plate transmissivity:** Figure 5 depicts the PCE-predicted behavior of the ratio of the fracture transmissivity *T* to that of the corresponding parallel plate geometry, *T*<sub>‖</sub>, and of the amplitude of the relative confidence interval Δ*q* *T*/*T*<sub>‖</sub>. The figure's layout is otherwise identical to that of Figure 4. Again, for low dimensionless pressure gradients the hydraulic behavior tends to the Newtonian transmissivity (see Figure 3a), as *T* → *T*₀ and *T*<sub>‖</sub> → *T*<sub>0‖</sub>. For sufficiently high dimensionless pressure gradient, indicatively √*P*/√*P*<sub>c</sub> > 3 (a transition roughly occurs between 0.1—quasi-Newtonian behavior—and 4), almost horizontal contour lines are seen in the *T/T*<sub>‖</sub> topography, indicating that the fracture closure's impact on the transmissivity scales equivalently in the rough-walled fracture and the corresponding parallel plate geometry. This is all the truer as the correlation length *L*<sub>c</sub> is larger. Moreover, as depicted in the maps on the right-hand side, the fracture closure does not significantly impact the statistical dispersion of *T/T*<sub>‖</sub> when the flow is strongly ST. Furthermore, when the ST hydraulic behavior is developed, the observed relation between the transmissivity and

the fracture closure is non-monotonic, for any value of *L/L*<sub>c</sub>. The ST nature of the fluids mitigates the viscous energy losses induced by the higher heterogeneity due to fracture closure, increasing the statistical likelihood of flow-enhancing hydraulic behavior. As the aperture field becomes more heterogeneous, both the transmissivity and statistical dispersion increase, with many realizations characterized by a transmissivity that is orders of magnitude higher than that of the corresponding parallel plate case. This can be explained by the formation of channels where the fluid shear-thins due to high shear rates, so that the resulting apparent viscosity decreases and velocities increase, resulting in preferential flow paths of very high conductance due to minor viscous energy losses. The higher statistical dispersion observed when fracture closure increases and *L/L*<sub>c</sub> = 1 depends on the more chaotic nature of the medium, which then affects the flow by generating larger channels.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Dependence of the fracture transmissivity normalized by that of the parallel plate of aperture equal to the rough fracture's aperture, $T/T_‖$, on both the dimensionless external pressure gradient $\overline{\nabla P}/\overline{\nabla P_c}$ and the fracture closure $\sigma_a/\langle a \rangle$, produced for different ratios $L/L_c$: 1, 4, and 16, for panels (a-c), respectively. The plots are otherwise identical to those of Figure 4.</figcaption></figure>

**Behavior of the correlation lengths' ratio:** We discuss in the present subsection how the ratio of the flow's longitudinal correlation length $\ell$ to the aperture field's longitudinal correlation length $\lambda$ depends on the problem's parameters. Note that since $\lambda \propto L_c$, as discussed earlier (see also the numerical simulations by Cintoli et al. (2005)), the ratios $\ell/\lambda$ and $\ell/L_c$ differ by a factor. In particular, Figure 6a shows the influence of the fracture closure and dimensionless pressure gradient on $\ell/\lambda$. On the right-hand side, Figures 6b and 6c display several realizations of the fracture aperture field and the corresponding flow field, for different correlation ratios ($L/L_c \in \{1; 4; 16\}$) and fracture closures: $\sigma_a/\langle a \rangle = 0.8$ (Figure 6b), and $\sigma_a/\langle a \rangle = 0.2$ (Figure 6c). Under a sufficiently small external pressure gradient ($\overline{\nabla P}/\overline{\nabla P_c} \to 0$) the flow correlation tends to the Newtonian behavior depicted in Figure 3. Unlike the results obtained for transmissivity, the flow correlation length in small fractures (i.e., $L \approx L_c$) shows a minor statistical dispersion around the ensemble average, which results in an almost uniform value ($\ell/\lambda \approx 1$) over the parameter space. In cases where $L/L_c \in \{4; 16\}$, the ST effect increases for increasing imposed pressure gradients. Higher dispersion characterizes media of lower fracture closure, where the flow correlation length is larger than the medium's correlation length: $\ell/\lambda \approx \lambda$ ranges from approximately 3 for $\sigma_a/\langle a \rangle \to 1$ up to 6 times $\sigma_a/\langle a \rangle \to 0$, in $L/L_c = 16$. Indeed, when the fracture closure is high, the presence of a contact zone strongly constrains the flow path, reducing the variability of $\ell/\lambda$, and consequently the impact of the fluid's rheology on $\ell/\lambda$. On the other hand, for a low fracture closure, the flow correlation length exhibits a marked effect of the fluid rheology on the shapes of the streamlines, exhibiting a stronger ST effect under an increasing dimensionless pressure gradient.

## 4. Conclusions

This work has characterized quantitatively the impact on fracture scale metrics (i.e., generalized transmissivity and flow correlation length) of (a) the aperture heterogeneity resulting from wall roughness and (b) the fluid's ST rheology. To this aim, we used a 2-D numerical model relying on a lubrication-based theoretical description of ST flow, under a stochastic framework involving ROMs. The ST constitutive law considered, the so-called Ellis model, corresponds to the complex rheology of polymer fluids and most ST fluids used in subsurface industrial applications. It exhibits a quasi-Newtonian behavior at low shear rates and transitions to a ST behavior under increasing shear rates. The dependence of the fractures scale metrics on the problem's parameters, that is, the topographies of these metrics over the parameter space, are obtained through the PCE as a ROM, relying on an

optimized set of MC simulations in the parameter space. This approach allows investigating a three-dimensional parameter space over a wide range of parameter values (imposed pressure gradient over two orders of magnitude, fracture closure between 0 and 1, ratio of the fracture length to the aperture fields' correlation length between 1 and 16) within a fraction (typically 3 orders of magnitude less) of the time that would be required without model reduction.

Results for Newtonian flow (such as that of water) confirm the previous conclusions, that is, that the fracture transmissivity is on average diminished by increasing the closure due to higher viscous energy losses, but with a dispersion over the statistics that is all the larger as the closure is larger; the correlation length hardly impacts the mean behavior, but the dispersion around that mean value is larger for a larger correlation length. The flow correlation length is weakly influenced by the fracture closure, whereas the medium (i.e., aperture field's) correlation length is more markedly impacted by it. Indeed, for small fractures, when $L/L_c \approx 1$, the medium and flow correlation lengths are similar. Conversely, the flow correlation length grows faster than the aperture field's correlation length when larger fractures are considered ($L > L_c$).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Dependence of the fracture transmissivity normalized by that of the parallel plate of aperture equal to the rough fracture's aperture, $T/T_‖$, on both the dimensionless external pressure gradient $\overline{\nabla P}/\overline{\nabla P_c}$ and the fracture closure $\sigma_a/\langle a \rangle$, produced for different ratios $L/L_c$: 1, 4, and 16, for panels (a–c), respectively. The plots are otherwise identical to those of Figure 4.</figcaption></figure>

**Behavior of the correlation lengths' ratio:** We discuss in the present subsection how the ratio of the flow's longitudinal correlation length $\ell$ to the aperture field's longitudinal correlation length $\lambda$ depends on the problem's parameters. Note that since $\lambda \propto L_c$, as discussed earlier (see also the numerical simulations by Cintoli et al. (2005)), the ratios $\ell/\lambda$ and $\ell/L_c$ differ by a factor. In particular, Figure 6a shows the influence of the fracture closure and dimensionless pressure gradient on $\ell/\lambda$. On the right-hand side, Figures 6b and 6c display several realizations of the fracture aperture field and the corresponding flow field, for different correlation ratios ($L/L_c \in \{1; 4; 16\}$) and fracture closures: $\sigma_a/\langle a \rangle = 0.8$ (Figure 6b), and $\sigma_a/\langle a \rangle = 0.2$ (Figure 6c). Under a sufficiently small external pressure gradient ($\overline{\nabla P}/\overline{\nabla P_c} \to 0$) the flow correlation tends to the Newtonian behavior depicted in Figure 3. Unlike the results obtained for transmissivity, the flow correlation length in small fractures (i.e., $L \approx L_c$) shows a minor statistical dispersion around the ensemble average, which results in an almost uniform value ($\ell/\lambda \approx 1$) over the parameter space. In cases where $L/L_c \in \{4; 16\}$, the ST effect increases for increasing imposed pressure gradients. Higher dispersion characterizes media of lower fracture closure, where the flow correlation length is larger than the medium's correlation length: $\ell/\lambda \approx \lambda$ ranges from approximately 3 for $\sigma_a/\langle a \rangle \to 1$ up to 6 times $\sigma_a/\langle a \rangle \to 0$, in $L/L_c = 16$. Indeed, when the fracture closure is high, the presence of a contact zone strongly constrains the flow path, reducing the variability of $\ell/\lambda$, and consequently the impact of the fluid's rheology on $\ell/\lambda$. On the other hand, for a low fracture closure, the flow correlation length exhibits a marked effect of the fluid rheology on the shapes of the streamlines, exhibiting a stronger ST effect under an increasing dimensionless pressure gradient.

## 4. Conclusions

This work has characterized quantitatively the impact on fracture scale metrics (i.e., generalized transmissivity and flow correlation length) of (a) the aperture heterogeneity resulting from wall roughness and (b) the fluid's ST rheology. To this aim, we used a 2-D numerical model relying on a lubrication-based theoretical description of ST flow, under a stochastic framework involving ROMs. The ST constitutive law considered, the so-called Ellis model, corresponds to the complex rheology of polymer fluids and most ST fluids used in subsurface industrial applications. It exhibits a quasi-Newtonian behavior at low shear rates and transitions to a ST behavior under increasing shear rates. The dependence of the fractures scale metrics on the problem's parameters, that is, the topographies of these metrics over the parameter space, are obtained through the PCE as a ROM, relying on an

optimized set of MC simulations in the parameter space. This approach allows investigating a three-dimensional parameter space over a wide range of parameter values (imposed pressure gradient over two orders of magnitude, fracture closure between 0 and 1, ratio of the fracture length to the aperture fields' correlation length between 1 and 16) within a fraction (typically 3 orders of magnitude less) of the time that would be required without model reduction.

Results for Newtonian flow (such as that of water) confirm the previous conclusions, that is, that the fracture transmissivity is on average diminished by increasing the closure due to higher viscous energy losses, but with a dispersion over the statistics that is all the larger as the closure is larger; the correlation length hardly impacts the mean behavior, but the dispersion around that mean value is larger for a larger correlation length. The flow correlation length is weakly influenced by the fracture closure, whereas the medium (i.e., aperture field's) correlation length is more markedly impacted by it. Indeed, for small fractures, when $L/L_c \approx 1$, the medium and flow correlation lengths are similar. Conversely, the flow correlation length grows faster than the aperture field's correlation length when larger fractures are considered ($L > L_c$).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) Dependence of the ratio of the longitudinal correlation length of flow velocities to that of the aperture field, $\ell/\lambda$, on both the dimensionless external pressure gradient $\nabla P/\nabla P_c$ and the fracture closure $\sigma_d/\langle a \rangle$, produced for different ratios $L/L_c$: 1, 4, and 16, illustrated from top to bottom. Panels (b), (c) display fracture aperture realizations on the left-hand side and velocity magnitude maps on the right-hand, for closure $\sigma_d/\langle a \rangle$ equal to 0.8 and 0.2, respectively, and for the three investigated correlation lengths.</figcaption></figure>

When considering the ST rheology exhibited (among other complex fluids) by polymer solutions, the generalized fracture transmissivity exhibits a transition from the Darcian behavior, typical of the linear Newtonian model, to a non-linear power-law regime, when increasing the imposed pressure gradient. Fracture closure and medium's correlation length affect neither the Newtonian nor the power-law trend, but they influence how fast the transition between these two regimes occurs: a higher fracture closure induces a faster transition between the two regimes (i.e., transition to the non-linear regime occurs under higher dimensionless pressure gradients). Besides, a higher dispersion around the mean behavior is observed in small fractures, that is, for a correlation length that is larger as compared to the fracture length ($L/L_c \approx 1$).

Unlike the purely Newtonian rheology, the ST rheology of polymer solutions mitigates the transmissivity reduction observed on average with respect to the corresponding parallel plate (i.e., smooth) geometry; consequently, the average generalized transmissivity for ST flow can be orders of magnitude higher in the rough geological fracture than in the corresponding parallel plate. The Newtonian model thus overestimates fracture transmissivity when used to predict the flow of fluids that exhibit an ST behavior, since the ST rheology is more efficient at minimizing viscous energy losses induced by the aperture variability.

The ratio of the flow correlation length to the medium (i.e., aperture field's) correlation length was also analyzed for the first time in non-Newtonian fracture flow. Fracture closure significantly impacts the flow correlation
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Dependence of the ratio of the longitudinal correlation length of flow velocities to that of the aperture field, $\ell/\lambda$, on both the dimensionless external pressure gradient $\nabla P/\nabla P_c$ and the fracture closure $\sigma_d/\langle a \rangle$, produced for different ratios $L/L_c$: 1, 4, and 16, illustrated from top to bottom. Panels (b), (c) display fracture aperture realizations on the left-hand side and velocity magnitude maps on the right-hand, for closure $\sigma_d/\langle a \rangle$ equal to 0.8 and 0.2, respectively, and for the three investigated correlation lengths.</figcaption></figure>

When considering the ST rheology exhibited (among other complex fluids) by polymer solutions, the generalized fracture transmissivity exhibits a transition from the Darcian behavior, typical of the linear Newtonian model, to a non-linear power-law regime, when increasing the imposed pressure gradient. Fracture closure and medium's correlation length affect neither the Newtonian nor the power-law trend, but they influence how fast the transition between these two regimes occurs: a higher fracture closure induces a faster transition between the two regimes (i.e., transition to the non-linear regime occurs under higher dimensionless pressure gradients). Besides, a higher dispersion around the mean behavior is observed in small fractures, that is, for a correlation length that is larger as compared to the fracture length ($L/L_c \approx 1$).

Unlike the purely Newtonian rheology, the ST rheology of polymer solutions mitigates the transmissivity reduction observed on average with respect to the corresponding parallel plate (i.e., smooth) geometry; consequently, the average generalized transmissivity for ST flow can be orders of magnitude higher in the rough geological fracture than in the corresponding parallel plate. The Newtonian model thus overestimates fracture transmissivity when used to predict the flow of fluids that exhibit an ST behavior, since the ST rheology is more efficient at minimizing viscous energy losses induced by the aperture variability.

The ratio of the flow correlation length to the medium (i.e., aperture field's) correlation length was also analyzed for the first time in non-Newtonian fracture flow. Fracture closure significantly impacts the flow correlation
length for ST flow. This effect is more effective when fracture closure is sufficiently small (approximately smaller than 0.5). Indeed, a higher value of the closure implies larger contact zones, forcing the flow in open channels of constrained pathways. On the other hand, when contact zones are limited, the interplay between aperture variability and fluid rheology may result in a flow correlation length up to 8 times larger than the medium's correlation length. But this ratio is all the larger as the medium's correlation length is smaller. When the a priori structural correlation length $L_c$ is equal to the fracture length, the flow correlation length remains similar to the medium's correlation length at all closures and imposed pressure gradients.

With these results, we have provided a comprehensive panorama of Newtonian and ST flow in geological fractures, as a function of the statistical parameters describing a fracture's geometry, the fluid's rheology, and the imposed pressure gradient. Future works will address solute transport by such flows.

## Appendix A: Relationship Between $\lambda$ and $L_c$

In this work, two characteristic structural lengths, $L_c$ and $\lambda$, are considered. As described in Section 2, the former is an a priori parameter adopted to generate synthetic wall topographies before removing interpenetrations between the two walls when coupled; this parameter also characterizes the resulting aperture field, which exhibits isotropic self-affinity for scales smaller than $L_c$. On the other hand, the latter is an a posteriori estimate of the zero-crossing length of the ACF obtained as the average of the ACFs computed from each longitudinal profile of the aperture field. Note that the ACFs are estimated over aperture fields that may present contact zones, hence a non-uniform discrete Fourier transform-based algorithm is needed to non-uniformly sample each profile. Despite the differences between the two approaches, the Wiener-Khinchin theorem defines the relationship between the ACF and the spectrum, hence we expect a link between the two lengths. Considering that the ACF analyses the signal (i.e., each profile) in the space domain, while the power spectrum considers the data in the frequency domain, differences may also arise due to computational approximation (i.e., discrete representation) and the finite spatial dimensions of the domain. In Figure A1, we analyze the PDFs of the ratio of $\lambda$ to $L_c$, for three different fracture sizes ($L_c$ being constant) and different fracture closures. It turns out that the fracture closure does not affect the PDF of $\lambda/L_c$, but the PDF depends on $L/L_c$. Hence, for a given $L/L_c$, one can consider that $\lambda$ is proportional to $L_c$ with a dispersion of the statistics around the mean ratio (i.e., $\sigma_{L_c/\mu_{L_c}}$) smaller than 15% for $L/L_c = 1$, 13% for $L/L_c = 4$, and 5% for $L/L_c = 16$.

<figure><img src="image_10.png" /><figcaption>**Figure A1.** Probability density functions (PDF) of the ratio $\lambda/L_c$ for different values of $L/L_c$ and different fracture closures: $\sigma_d/(a) = 0.1$ continuous lines, $\sigma_d/(a) = 0.5$ (dashed line), and $\sigma_d/(a) = 1$ (dash-dotted lines). PDFs have been generated from 1,000 aperture realizations with the following statistical parameters: mesh size $N = 2^{10}$, mechanical aperture $\langle a \rangle = 10^{-3}$ m, Hurst exponent $H = 0.8$, and $L_c = 0.1$ m.</figcaption></figure>
length for ST flow. This effect is more effective when fracture closure is sufficiently small (approximately smaller than 0.5). Indeed, a higher value of the closure implies larger contact zones, forcing the flow in open channels of constrained pathways. On the other hand, when contact zones are limited, the interplay between aperture variability and fluid rheology may result in a flow correlation length up to 8 times larger than the medium's correlation length. But this ratio is all the larger as the medium's correlation length is smaller. When the a priori structural correlation length $L_c$ is equal to the fracture length, the flow correlation length remains similar to the medium's correlation length at all closures and imposed pressure gradients.

With these results, we have provided a comprehensive panorama of Newtonian and ST flow in geological fractures, as a function of the statistical parameters describing a fracture's geometry, the fluid's rheology, and the imposed pressure gradient. Future works will address solute transport by such flows.

## Appendix A: Relationship Between $\lambda$ and $L_c$

In this work, two characteristic structural lengths, $L_c$ and $\lambda$, are considered. As described in Section 2, the former is an a priori parameter adopted to generate synthetic wall topographies before removing interpenetrations between the two walls when coupled; this parameter also characterizes the resulting aperture field, which exhibits isotropic self-affinity for scales smaller than $L_c$. On the other hand, the latter is an a posteriori estimate of the zero-crossing length of the ACF obtained as the average of the ACFs computed from each longitudinal profile of the aperture field. Note that the ACFs are estimated over aperture fields that may present contact zones, hence a non-uniform discrete Fourier transform-based algorithm is needed to non-uniformly sample each profile. Despite the differences between the two approaches, the Wiener-Khinchin theorem defines the relationship between the ACF and the spectrum, hence we expect a link between the two lengths. Considering that the ACF analyses the signal (i.e., each profile) in the space domain, while the power spectrum considers the data in the frequency domain, differences may also arise due to computational approximation (i.e., discrete representation) and the finite spatial dimensions of the domain. In Figure A1, we analyze the PDFs of the ratio of $\lambda$ to $L_c$, for three different fracture sizes ($L_c$ being constant) and different fracture closures. It turns out that the fracture closure does not affect the PDF of $\lambda/L_c$, but the PDF depends on $L/L_c$. Hence, for a given $L/L_c$, one can consider that $\lambda$ is proportional to $L_c$ with a dispersion of the statistics around the mean ratio (i.e., $\sigma_{\lambda/L_c} / \mu_{\lambda/L_c}$) smaller than 15% for $L/L_c = 1$, 13% for $L/L_c = 4$, and 5% for $L/L_c = 16$.

<figure><img src="image_9.png" /><figcaption>**Figure A1.** Probability density functions (PDF) of the ratio $\lambda/L_c$ for different values of $L/L_c$ and different fracture closures: $\sigma_d/(a) = 0.1$ continuous lines, $\sigma_d/(a) = 0.5$ (dashed line), and $\sigma_d/(a) = 1$ (dash-dotted lines). PDFs have been generated from 1,000 aperture realizations with the following statistical parameters: mesh size $N = 2^{10}$, mechanical aperture $\langle a \rangle = 10^{-3}$ m, Hurst exponent $H = 0.8$, and $L_c = 0.1$ m.</figcaption></figure>
## Acknowledgments

AL acknowledges the support of both the European Union—NextGenerationEU under the National Recovery and Resilience Plan, Project: RETURN “multi-Risk science for resilient communities under a changing climate.” Code: PE0000005, CUP: J33C22002840002, and the ALMARie Curie SUpER initiative of the Università di Bologna.

# Data Availability Statement

There are no data-sharing issues since all of the numerical information is provided in the figures produced by solving the equations in the paper.

## References

Bodin, J., Delay, F., & de Marsily, G. (2003). Solute transport in a single fracture with negligible matrix permeability: 1. Fundamental mechanisms. *Hydrogeology Journal*, *11*(4), 418–433. https://doi.org/10.1007/s10040-003-0268-2

Boffa, J., Allain, C., Chertcoff, R., Hulin, J., Plouraboué, F., & Roux, S. (1999). Roughness of sandstone fracture surfaces: Profilometry and shadow length investigations. *European Physical Journal B: Condensed Matter and Complex Systems*, *7*(2), 179–182. https://doi.org/10.1007/s100510050602

Boso, F., & Tartakovsky, D. M. (2018). Information-theoretic approach to bidirectional scaling. *Water Resources Research*, *54*(7), 4916–4928. https://doi.org/10.1029/2017wr021993

Bouchaud, E., Lapasset, G., & Planès, J. (1990). Fractal dimension of fractured surfaces: A universal value? *Europhysics Letters*, *13*(1), 73–79. https://doi.org/10.1209/0295-5075/13/1/013

Bour, O., & Davy, P. (1998). On the connectivity of three-dimensional fault networks. *Water Resources Research*, *34*(10), 2611–2622. https://doi.org/10.1029/98wr01861

Brown, S. (1987). Fluid flow through rock joints: The effect of surface roughness. *Journal of Geophysical Research*, *92*(B2), 1337–1347. https://doi.org/10.1029/jb092ib02p01337

Brown, S. (1989). Transport of fluid and electric current through a single fracture. *Journal of Geophysical Research*, *94*(B7), 9429–9438. https://doi.org/10.1029/jb094ib07p09429

Cintoli, S., Neuman, S., & Di Federico, V. (2005). Generating and scaling fractional Brownian motion on finite domains. *Geophysical Research Letters*, *32*(8). https://doi.org/10.1029/2005GL022608

Ciriello, V., & Barros, F. P. J. (2020). Characterizing the influence of multiple uncertainties on predictions of contaminant discharge in groundwater within a Lagrangian stochastic formulation. *Water Resources Research*, *56*(10). https://doi.org/10.1029/2020wr027867

Ciriello, V., Di Federico, V., Riva, M., Cadini, F., De Sanctis, J., Zio, E., & Guadagnini, A. (2012). Polynomial chaos expansion for global sensitivity analysis applied to a model of radionuclide migration in a randomly heterogeneous aquifer. *Stochastic Environmental Research and Risk Assessment*, *27*(4), 945–954. https://doi.org/10.1007/s00477-012-0616-7

Ciriello, V., Edery, Y., Guadagnini, A., & Berkowitz, B. (2015). Multimodel framework for characterization of transport in porous media. *Water Resources Research*, *51*(5), 3384–3402. https://doi.org/10.1002/2015wr017047

Ciriello, V., Guadagnini, A., Di Federico, V., Edery, Y., & Berkowitz, B. (2013). Comparative analysis of formulations for conservative transport in porous media through sensitivity-based parameter calibration. *Water Resources Research*, *49*(9), 5206–5220. https://doi.org/10.1002/wrcr.20395

Ciriello, V., Lauriola, I., Bonvicini, S., Cozzani, V., Di Federico, V., & Tartakovsky, D. M. (2017). Impact of hydrogeological uncertainty on estimation of environmental risks posed by hydrocarbon transportation networks. *Water Resources Research*, *53*(11), 8686–8697. https://doi.org/10.1002/2017wr021368

Ciriello, V., Lauriola, I., & Tartakovsky, D. M. (2019). Distribution-based global sensitivity analysis in hydrology. *Water Resources Research*, *55*(11), 8708–8720. https://doi.org/10.1029/2019WR025844

Comolli, A., Hakoun, V., & Dentz, M. (2019). Mechanisms, upscaling, and prediction of anomalous dispersion in heterogeneous porous media. *Water Resources Research*, *55*(10), 8197–8222. https://doi.org/10.1029/2019wr024919

Cvetkovic, V., Cheng, H., & Wen, X. H. (1996). Analysis of nonlinear effects on tracer migration in heterogeneous aquifers using Lagrangian travel time statistics. *Water Resources Research*, *32*(6), 1671–1680. https://doi.org/10.1029/96wr00278

de Dreuzy, J.-R., Davy, P., & Bour, O. (2002). Hydraulic properties of two-dimensional random fracture networks following power law distributions of length and aperture. *Water Resources Research*, *38*(12), 12–21. https://doi.org/10.1029/2001wr001009

de Dreuzy, J.-R., Méheust, Y., & Pichot, G. (2012). Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN). *Journal of Geophysical Research*, *117*(B11207). https://doi.org/10.1029/2012JB009461

Dentz, M., Kang, P. K., Comolli, A., Borgne, T. L., & Lester, D. R. (2016). Continuous time random walks for the evolution of Lagrangian velocities. *Physical Review Fluids*, *1*(7), 074004. https://doi.org/10.1103/physrevfluids.1.074004

Di Federico, V. (1997). Estimates of equivalent aperture for non-Newtonian flow in a rough-walled fracture. *International Journal of Rock Mechanics and Mining Sciences*, *34*(7), 1133–1137. https://doi.org/10.1016/S1365-1609(97)90205-7

Di Federico, V. (1998). Non-Newtonian flow in a variable aperture fracture. *Transport in Porous Media*, *30*(1), 75–86. https://doi.org/10.1023/a:1006512822518

Di Federico, V. (2001). On non-Newtonian fluid flow in rough fractures. *Water Resources Research*, *37*(9), 2425–2430. https://doi.org/10.1029/2001wr000359

Durham, W., & Bonner, B. (1995). Closure and fluid flow in discrete fractures. In L. Myer, C. Tsang, N. Cook, & R. Goodman (Eds.), *Fractured and jointed rock masses* (pp. 441–446).

Faroughi, S. A., Pruquot, A. J.-C. J., & McAndrew, J. (2018). The rheological behavior of energized fluids and foams with application to hydraulic fracturing: Review. *Journal of Petroleum Science and Engineering*, *163*, 243–263. https://doi.org/10.1016/j petrol.2017.12.051

Focaccia, S., Panini, G., Pedrazzoli, P., & Ciriello, V. (2021). A meta-modeling approach for hydrological forecasting under uncertainty: Application to groundwater nitrate response to climate change. *Journal of Hydrology*, *603*, 127173. https://doi.org/10.1016/j.jhydrol.2021.127173

Ge, S. (1997). A governing equation for fluid flow in rough fractures. *Water Resources Research*, *33*(1), 53–61. https://doi.org/10.1029/96wr02588

Ghanem, R. G., & Spanos, P. D. (1991). *Stochastic finite elements: A spectral approach*. Springer. https://doi.org/10.1007/978-1-4612-3094-6

Hosseini, H., Tsau, J. S., Shafer-Peltier, K., Marshall, C., Ye, Q., & Ghahfarokhi, R. B. (2019). Experimental and mechanistic study of stabilized dry CO₂ foam using polyelectrolyte complex nanoparticles compatible with produced water to improve hydraulic fracturing performance. *Industrial & Engineering Chemistry Research*, *58*(22), 9431–9449. https://doi.org/10.1021/acs.iecr.9b01390

Jiang, J. (2021). Localized nonlinear solution strategies for efficient simulation of unconventional reservoirs. *Applied Mathematical Modelling*, *90*, 977–994. https://doi.org/10.1016/j.apm.2020.10.018

Jung, R. (1989). Hydraulic in situ investigations of an artificial fracture in the falkenberg granite. *International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts*, *26*(3–4), 301–308. https://doi.org/10.1016/0148-9062(89)91978-5
## Acknowledgments

AL acknowledges the support of both the European Union—NextGenerationEU under the National Recovery and Resilience Plan, Project: RETURN “multi-Risk science for resilient communities under a changing climate.” Code: PE0000005, CUP: J33C22002840002, and the ALMARie Curie SUpER initiative of the Università di Bologna.

# Data Availability Statement

There are no data-sharing issues since all of the numerical information is provided in the figures produced by solving the equations in the paper.

## References

Bodin, J., Delay, F., & de Marsily, G. (2003). Solute transport in a single fracture with negligible matrix permeability: 1. Fundamental mechanisms. *Hydrogeology Journal*, 11(4), 418–433. https://doi.org/10.1007/s10040-003-0268-2

Boffa, J., Allain, C., Chertcoff, R., Hulin, J., Plouraboué, F., & Roux, S. (1999). Roughness of sandstone fracture surfaces: Profilometry and shadow length investigations. *European Physical Journal B: Condensed Matter and Complex Systems*, 7(2), 179–182. https://doi.org/10.1007/s100510050602

Boso, F., & Tartakovsky, D. M. (2018). Information-theoretic approach to bidirectional scaling. *Water Resources Research*, 54(7), 4916–4928. https://doi.org/10.1029/2017wr021993

Bouchaud, E., Lapasset, G., & Planès, J. (1990). Fractal dimension of fractured surfaces: A universal value? *Europhysics Letters*, 13(1), 73–79. https://doi.org/10.1209/0295-5075/13/1/013

Bour, O., & Davy, P. (1998). On the connectivity of three-dimensional fault networks. *Water Resources Research*, 34(10), 2611–2622. https://doi.org/10.1029/98wr01861

Brown, S. (1987). Fluid flow through rock joints: The effect of surface roughness. *Journal of Geophysical Research*, 92(B2), 1337–1347. https://doi.org/10.1029/jb092ib02p01337

Brown, S. (1989). Transport of fluid and electric current through a single fracture. *Journal of Geophysical Research*, 94(B7), 9429–9438. https://doi.org/10.1029/jb094ib07p09429

Cintoli, S., Neuman, S., & Di Federico, V. (2005). Generating and scaling fractional Brownian motion on finite domains. *Geophysical Research Letters*, 32(8). https://doi.org/10.1029/2005GL022608

Ciriello, V., & Barros, F. P. J. (2020). Characterizing the influence of multiple uncertainties on predictions of contaminant discharge in groundwater within a Lagrangian stochastic formulation. *Water Resources Research*, 56(10). https://doi.org/10.1029/2020wr027867

Ciriello, V., Di Federico, V., Riva, M., Cadini, F., De Sanctis, J., Zio, E., & Guadagnini, A. (2012). Polynomial chaos expansion for global sensitivity analysis applied to a model of radionuclide migration in a randomly heterogeneous aquifer. *Stochastic Environmental Research and Risk Assessment*, 27(4), 945–954. https://doi.org/10.1007/s00477-012-0616-7

Ciriello, V., Edery, Y., Guadagnini, A., & Berkowitz, B. (2015). Multimodel framework for characterization of transport in porous media. *Water Resources Research*, 51(5), 3384–3402. https://doi.org/10.1002/2015wr017047

Ciriello, V., Guadagnini, A., Di Federico, V., Edery, Y., & Berkowitz, B. (2013). Comparative analysis of formulations for conservative transport in porous media through sensitivity-based parameter calibration. *Water Resources Research*, 49(9), 5206–5220. https://doi.org/10.1002/wrcr.20395

Ciriello, V., Lauriola, I., Bonvicini, S., Cozzani, V., Di Federico, V., & Tartakovsky, D. M. (2017). Impact of hydrogeological uncertainty on estimation of environmental risks posed by hydrocarbon transportation networks. *Water Resources Research*, 53(11), 8686–8697. https://doi.org/10.1002/2017wr021368

Ciriello, V., Lauriola, I., & Tartakovsky, D. M. (2019). Distribution-based global sensitivity analysis in hydrology. *Water Resources Research*, 55(11), 8708–8720. https://doi.org/10.1029/2019WR025844

Comolli, A., Hakoun, V., & Dentz, M. (2019). Mechanisms, upscaling, and prediction of anomalous dispersion in heterogeneous porous media. *Water Resources Research*, 55(10), 8197–8222. https://doi.org/10.1029/2019wr024919

Cvetkovic, V., Cheng, H., & Wen, X. H. (1996). Analysis of nonlinear effects on tracer migration in heterogeneous aquifers using Lagrangian travel time statistics. *Water Resources Research*, 32(6), 1671–1680. https://doi.org/10.1029/96wr00278

de Dreuzy, J.-R., Davy, P., & Bour, O. (2002). Hydraulic properties of two-dimensional random fracture networks following power law distributions of length and aperture. *Water Resources Research*, 38(12), 12–21. https://doi.org/10.1029/2001wr001009

de Dreuzy, J.-R., Méheust, Y., & Pichot, G. (2012). Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN). *Journal of Geophysical Research*, 117(B11207). https://doi.org/10.1029/2012JB009461

Dentz, M., Kang, P. K., Comolli, A., Borgne, T. L., & Lester, D. R. (2016). Continuous time random walks for the evolution of Lagrangian velocities. *Physical Review Fluids*, 1(7), 074004. https://doi.org/10.1103/physrevfluids.1.074004

Di Federico, V. (1997). Estimates of equivalent aperture for non-Newtonian flow in a rough-walled fracture. *International Journal of Rock Mechanics and Mining Sciences*, 34(7), 1133–1137. https://doi.org/10.1016/S1365-1609(97)90205-7

Di Federico, V. (1998). Non-Newtonian flow in a variable aperture fracture. *Transport in Porous Media*, 30(1), 75–86. https://doi.org/10.1023/a:1006512822518

Di Federico, V. (2001). On non-Newtonian fluid flow in rough fractures. *Water Resources Research*, 37(9), 2425–2430. https://doi.org/10.1029/2001wr000359

Durham, W., & Bonner, B. (1995). Closure and fluid flow in discrete fractures. In L. Myer, C. Tsang, N. Cook, & R. Goodman (Eds.), *Fractured and jointed rock masses* (pp. 441–446).

Faroughi, S. A., Pruvot, A. J.-C. J., & McAndrew, J. (2018). The rheological behavior of energized fluids and foams with application to hydraulic fracturing: Review. *Journal of Petroleum Science and Engineering*, 163, 243–263. https://doi.org/10.1016/j petrol.2017.12.051

Focaccia, S., Panini, G., Pedrazzoli, P., & Ciriello, V. (2021). A meta-modeling approach for hydrological forecasting under uncertainty: Application to groundwater nitrate response to climate change. *Journal of Hydrology*, 603, 127173. https://doi.org/10.1016/j.jhydrol.2021.127173

Ge, S. (1997). A governing equation for fluid flow in rough fractures. *Water Resources Research*, 33(1), 53–61. https://doi.org/10.1029/96wr02588

Ghanem, R. G., & Spanos, P. D. (1991). *Stochastic finite elements: A spectral approach*. Springer. https://doi.org/10.1007/978-1-4612-3094-6

Hosseini, H., Tsau, J. S., Shafer-Peltier, K., Marshall, C., Ye, Q., & Ghahfarokhi, R. B. (2019). Experimental and mechanistic study of stabilized dry CO₂ foam using polyelectrolyte complex nanoparticles compatible with produced water to improve hydraulic fracturing performance. *Industrial & Engineering Chemistry Research*, 58(22), 9431–9449. https://doi.org/10.1021/acs.iecr.9b01390

Jiang, J. (2021). Localized nonlinear solution strategies for efficient simulation of unconventional reservoirs. *Applied Mathematical Modelling*, 90, 977–994. https://doi.org/10.1016/j.apm.2020.10.018

Jung, R. (1989). Hydraulic in situ investigations of an artificial fracture in the falkenberg granite. *International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts*, 26(3–4), 301–308. https://doi.org/10.1016/0148-9062(89)91978-5
Kamal, M. S., Sultan, A. S., Al-Mubaiyedh, U. A., & Hussein, I. A. (2015). Review on polymer flooding: Rheology, adsorption, stability, and field applications of various polymer systems. *Polymer Reviews*, 55(3), 491–530. https://doi.org/10.1080/15583724.2014.982821

Kang, P. K., Dentz, M., Le Borgne, T., Lee, S., & Juanes, R. (2017). Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes. *Advances in Water Resources*, 106, 80–94. https://doi.org/10.1016/j.advwatres.2017.03.024

Kim, T., Han, W. S., Piao, J., Kang, P. K., & Shin, J. (2022). Predicting remediation efficiency of LNAPLs using surrogate polynomial chaos expansion model and global sensitivity analysis. *Advances in Water Resources*, 163, 104179. https://doi.org/10.1016/j.advwatres.2022.104179

Lenci, A., Méheust, Y., Putti, M., & Di Federico, V. (2022). Monte Carlo simulations of shear-thinning flow in geological fractures. *Water Resources Research*, 58(9). https://doi.org/10.1029/2022wr032024

Lenci, A., Putti, M., Di Federico, V., & Méheust, Y. (2022). A lubrication-based solver for shear-thinning flow in rough fractures. *Water Resources Research*, 58(8). https://doi.org/10.1029/2021wr031760

Leung, D. Y., Caramanna, G., & Maroto-Valer, M. M. (2014). An overview of current status of carbon dioxide capture and storage technologies. *Renewable & Sustainable Energy Reviews*, 39, 426–443. https://doi.org/10.1016/j.rser.2014.07.093

Méheust, Y., & Schmittbuhl, J. (2000). Flow enhancement of a rough fracture. *Geophysical Research Letters*, 27(18), 2989–2992. https://doi.org/10.1029/1999gl008464

Méheust, Y., & Schmittbuhl, J. (2001). Geometrical heterogeneities and permeability anisotropy of rough fractures. *Journal of Geophysical Research*, 106(B2), 2089–2102. https://doi.org/10.1029/2000jb900306

Méheust, Y., & Schmittbuhl, J. (2003). Scale effects related to flow in rough fractures. *Pure and Applied Geophysics*, 160(5), 1023–1050. https://doi.org/10.1007/pl00012559

Merchán-Rivera, P., Wohlmuth, B., & Chiogna, G. (2021). Identifying stagnation zones and reverse flow caused by river-aquifer interaction: An approach based on polynomial chaos expansions. *Water Resources Research*, 57(12), e2021WR029824. https://doi.org/10.1029/2021WR029824

Morales, V. L., Dentz, M., Willmann, M., & Holzner, M. (2017). Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory. *Geophysical Research Letters*, 44(18), 9361–9371. https://doi.org/10.1002/2017gl074326

Napov, A., & Notay, Y. (2012). An algebraic multigrid method with guaranteed convergence rate. *SIAM Journal on Scientific Computing*, 34(2), A1079–A1109. https://doi.org/10.1137/100818509

Notay, Y. (2012). Aggregation-based algebraic multigrid for convection-diffusion equations. *SIAM Journal on Scientific Computing*, 34(4), A2288–A2316. https://doi.org/10.1137/110835347

Oladyshkin, S., de Barros, F., & Nowak, W. (2012). Global sensitivity analysis: A flexible and efficient framework with an example from stochastic hydrogeology. *Advances in Water Resources*, 37, 10–22. https://doi.org/10.1016/j.advwatres.2011.11.001

Pollock, D. W. (1988). Semianalytical computation of path lines for finite-difference models. *Ground Water*, 26(6), 743–750. https://doi.org/10.1111/j.1745-6584.1988.tb00425.x

Power, W. L., Tullis, T. E., Brown, S. R., Boitnott, G. N., & Scholz, C. H. (1987). Roughness of natural fault surfaces. *Geophysical Research Letters*, 14(1), 29–32. https://doi.org/10.1029/gl014i001p00029

Pruess, K. (2006). Enhanced geothermal systems (EGS) using CO₂ as working fluid—A novel approach for generating renewable energy with simultaneous sequestration of carbon. *Geothermics*, 35(4), 351–367. https://doi.org/10.1016/j.geothermics.2006.08.002

Schmittbuhl, J., Schmitt, F., & Scholz, C. (1995). Scaling invariance of crack surfaces. *Journal of Geophysical Research*, 100(B4), 5953–5973. https://doi.org/10.1029/94jb02885

Thompson, M. E., & Brown, S. R. (1991). The effect of anisotropic surface roughness on flow and transport in fractures. *Journal of Geophysical Research*, 96(B13), 21923–21932. https://doi.org/10.1029/91jb02252

Tosco, T., Marchisio, D. L., Lince, F., & Sethi, R. (2012). Extension of the Darcy–Forchheimer law for shear-thinning fluids and validation via pore-scale flow simulations. *Transport in Porous Media*, 96(1), 1–20. https://doi.org/10.1007/s11242-012-0070-5

Tsang, Y. W., & Tsang, C. F. (1987). Channel model of flow through fractured media. *Water Resources Research*, 23(3), 467–479. https://doi.org/10.1029/wr023i003p00467

Webster, M., Tatang, M. A., & McRae, G. J. (1996). Application of the probabilistic collocation method for an uncertainty analysis of a simple ocean model (Tech. Rep. No. MIT joint program on the science and policy of global change reports series No. 4). MIT.

Wiener, N. (1938). The homogeneous chaos. *American Journal of Mathematics*, 60(4), 897. https://doi.org/10.2307/2371268

Xiu, D., & Karniadakis, G. E. (2002). The Wiener–Askey polynomial chaos for stochastic differential equations. *SIAM Journal on Scientific Computing*, 24(2), 619–644. https://doi.org/10.1137/s1064827501387826

Yang, J., Meng, L., & Guo, L. (2017). In situ remediation of chlorinated solvent-contaminated groundwater using ZVI/organic carbon amendment in China: Field pilot test and full-scale application. *Environmental Science & Pollution Research*, 25(6), 5051–5062. https://doi.org/10.1007/s11356-017-9903-7

Zhang, W.-X. (2003). Nanoscale iron particles for environmental remediation: An overview. *Journal of Nanoparticle Research*, 5(3/4), 323–332. https://doi.org/10.1023/a:1025520116015

Zhong, L., Oostrom, M., Wietsma, T., & Covert, M. (2008). Enhanced remedial amendment delivery through fluid viscosity modifications: Experiments and numerical simulations. *Journal of Contaminant Hydrology*, 101(1–4), 29–41. https://doi.org/10.1016/j.jconhyd.2008.07.007

Zhu, J., Hu, K., Lu, X., Huang, X., Liu, K., & Wu, X. (2015). A review of geothermal energy resources, development, and applications in China: Current status and prospects. *Energy*, 93, 466–483. https://doi.org/10.1016/j.energy.2015.08.098
Kamal, M. S., Sultan, A. S., Al-Mubaiyedh, U. A., & Hussein, I. A. (2015). Review on polymer flooding: Rheology, adsorption, stability, and field applications of various polymer systems. *Polymer Reviews*, 55(3), 491–530. https://doi.org/10.1080/15583724.2014.982821

Kang, P. K., Dentz, M., Le Borgne, T., Lee, S., & Juanes, R. (2017). Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes. *Advances in Water Resources*, 106, 80–94. https://doi.org/10.1016/j.advwatres.2017.03.024

Kim, T., Han, W. S., Piao, J., Kang, P. K., & Shin, J. (2022). Predicting remediation efficiency of LNAPLs using surrogate polynomial chaos expansion model and global sensitivity analysis. *Advances in Water Resources*, 163, 104179. https://doi.org/10.1016/j.advwatres.2022.104179

Lenci, A., Méheust, Y., Putti, M., & Di Federico, V. (2022). Monte Carlo simulations of shear-thinning flow in geological fractures. *Water Resources Research*, 58(9). https://doi.org/10.1029/2022wr032024

Lenci, A., Putti, M., Di Federico, V., & Méheust, Y. (2022). A lubrication-based solver for shear-thinning flow in rough fractures. *Water Resources Research*, 58(8). https://doi.org/10.1029/2021wr031760

Leung, D. Y., Caramanna, G., & Maroto-Valer, M. M. (2014). An overview of current status of carbon dioxide capture and storage technologies. *Renewable & Sustainable Energy Reviews*, 39, 426–443. https://doi.org/10.1016/j.rser.2014.07.093

Méheust, Y., & Schmittbuhl, J. (2000). Flow enhancement of a rough fracture. *Geophysical Research Letters*, 27(18), 2989–2992. https://doi.org/10.1029/1999gl008464

Méheust, Y., & Schmittbuhl, J. (2001). Geometrical heterogeneities and permeability anisotropy of rough fractures. *Journal of Geophysical Research*, 106(B2), 2089–2102. https://doi.org/10.1029/2000jb900306

Méheust, Y., & Schmittbuhl, J. (2003). Scale effects related to flow in rough fractures. *Pure and Applied Geophysics*, 160(5), 1023–1050. https://doi.org/10.1007/pl00012559

Merchán-Rivera, P., Wohlmuth, B., & Chiogna, G. (2021). Identifying stagnation zones and reverse flow caused by river-aquifer interaction: An approach based on polynomial chaos expansions. *Water Resources Research*, 57(12), e2021WR029824. https://doi.org/10.1029/2021WR029824

Morales, V. L., Dentz, M., Willmann, M., & Holzner, M. (2017). Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory. *Geophysical Research Letters*, 44(18), 9361–9371. https://doi.org/10.1002/2017gl074326

Napov, A., & Notay, Y. (2012). An algebraic multigrid method with guaranteed convergence rate. *SIAM Journal on Scientific Computing*, 34(2), A1079–A1109. https://doi.org/10.1137/100818509

Notay, Y. (2012). Aggregation-based algebraic multigrid for convection-diffusion equations. *SIAM Journal on Scientific Computing*, 34(4), A2288–A2316. https://doi.org/10.1137/110835347

Oladyshkin, S., de Barros, F., & Nowak, W. (2012). Global sensitivity analysis: A flexible and efficient framework with an example from stochastic hydrogeology. *Advances in Water Resources*, 37, 10–22. https://doi.org/10.1016/j.advwatres.2011.11.001

Pollock, D. W. (1988). Semianalytical computation of path lines for finite-difference models. *Ground Water*, 26(6), 743–750. https://doi.org/10.1111/j.1745-6584.1988.tb00425.x

Power, W. L., Tullis, T. E., Brown, S. R., Boitnott, G. N., & Scholz, C. H. (1987). Roughness of natural fault surfaces. *Geophysical Research Letters*, 14(1), 29–32. https://doi.org/10.1029/gl014i001p00029

Pruess, K. (2006). Enhanced geothermal systems (EGS) using CO₂ as working fluid—A novel approach for generating renewable energy with simultaneous sequestration of carbon. *Geothermics*, 35(4), 351–367. https://doi.org/10.1016/j.geothermics.2006.08.002

Schmittbuhl, J., Schmitt, F., & Scholz, C. (1995). Scaling invariance of crack surfaces. *Journal of Geophysical Research*, 100(B4), 5953–5973. https://doi.org/10.1029/94jb02885

Thompson, M. E., & Brown, S. R. (1991). The effect of anisotropic surface roughness on flow and transport in fractures. *Journal of Geophysical Research*, 96(B13), 21923–21932. https://doi.org/10.1029/91jb02252

Tosco, T., Marchisio, D. L., Lince, F., & Sethi, R. (2012). Extension of the Darcy–Forchheimer law for shear-thinning fluids and validation via pore-scale flow simulations. *Transport in Porous Media*, 96(1), 1–20. https://doi.org/10.1007/s11242-012-0070-5

Tsang, Y. W., & Tsang, C. F. (1987). Channel model of flow through fractured media. *Water Resources Research*, 23(3), 467–479. https://doi.org/10.1029/wr023i003p00467

Webster, M., Tatang, M. A., & McRae, G. J. (1996). Application of the probabilistic collocation method for an uncertainty analysis of a simple ocean model (Tech. Rep. No. MIT joint program on the science and policy of global change reports series No. 4). MIT.

Wiener, N. (1938). The homogeneous chaos. *American Journal of Mathematics*, 60(4), 897. https://doi.org/10.2307/2371268

Xiu, D., & Karniadakis, G. E. (2002). The Wiener–Askey polynomial chaos for stochastic differential equations. *SIAM Journal on Scientific Computing*, 24(2), 619–644. https://doi.org/10.1137/s1064827501387826

Yang, J., Meng, L., & Guo, L. (2017). In situ remediation of chlorinated solvent-contaminated groundwater using ZVI/organic carbon amendment in China: Field pilot test and full-scale application. *Environmental Science & Pollution Research*, 25(6), 5051–5062. https://doi.org/10.1007/s11356-017-9903-7

Zhang, W.-X. (2003). Nanoscale iron particles for environmental remediation: An overview. *Journal of Nanoparticle Research*, 5(3/4), 323–332. https://doi.org/10.1023/a:1025520116015

Zhong, L., Oostrom, M., Wietsma, T., & Covert, M. (2008). Enhanced remedial amendment delivery through fluid viscosity modifications: Experiments and numerical simulations. *Journal of Contaminant Hydrology*, 101(1–4), 29–41. https://doi.org/10.1016/j.jconhyd.2008.07.007

Zhu, J., Hu, K., Lu, X., Huang, X., Liu, K., & Wu, X. (2015). A review of geothermal energy resources, development, and applications in China: Current status and prospects. *Energy*, 93, 466–483. https://doi.org/10.1016/j.energy.2015.08.098