# A Lubrication-Based Solver for Shear-Thinning Flow in Rough Fractures

Alessandro Lenci, M. Putti, V. Di Federico, Yves Méheust

► To cite this version:

Alessandro Lenci, M. Putti, V. Di Federico, Yves Méheust. A Lubrication-Based Solver for Shear-Thinning Flow in Rough Fractures. Water Resources Research, 2022, 58 (8), pp.e2021WR031760. 10.1029/2021WR031760 . insu-03760218

HAL Id: insu-03760218

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

Submitted on 25 Aug 2022

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/2021WR031760

**Key Points:**

* A computationally efficient code is implemented to solve the flow of shear-thinning (ST) fluids in geological fractures

* Increasing the ST property promotes flow heterogeneity and stronger localization patterns

* Transmissivity attenuation due to flow channeling is mitigated by the ST rheology

**Correspondence to:**

A. Lenci,
alessandro.lenci@unibo.it

**Citation:**

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, e2021WR031760. https://doi.org/10.1029/2021WR031760

Received 17 JAN 2022

Accepted 28 JUN 2022

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

## A Lubrication-Based Solver for Shear-Thinning Flow in Rough Fractures

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

¹Department of Civil, Chemical, Environmental and Materials Engineering, Alma Mater Studiorum Università di Bologna, Bologna, Italy, ²University Rennes, CNRS, Géosciences Rennes – UMR 6118, Rennes, France, ³Department of Mathematics, Università di Padova, Padua, Italy

**Abstract** The depth-averaged (2D) lubrication theory is often adopted to simulate Newtonian flow in rough fractures. This approach, which is computationally much less expensive than using 3D CFD solvers, allows addressing large ensembles of stochastic fracture realizations. For creeping flow, the degree of approximation introduced is limited as long as the apertures vary relatively smoothly. We propose the first generalization of this approach addressing the flow of fluids whose rheology, described by the Ellis model, is shear-thinning (ST) above a crossover shear stress and Newtonian (of viscosity $\mu_0$) below. The resulting nonlinear Reynolds equation for pressures is solved for a vast range of realistic rheological parameter values using a novel and specifically designed finite volume-based numerical model. The spatial discretization takes inspiration from the graph p-Laplacian to yield a symmetric Newton Jacobian, allowing for a highly efficient inexact implementation of the preconditioned conjugate gradient-based Newton-Krylov method. This is combined with a parameter continuation strategy to increase code robustness and ensure global convergence for flow indices as low as 0.1 with an excellent efficiency. This original solver is used to investigate realistic synthetic rough fracture geometries, which exhibits both self-affinity and a correlation length. The results show that the ST rheology mitigates the effects of aperture heterogeneities, increasing fracture transmissivity by several orders of magnitudes as compared to the Newtonian flow of viscosity $\mu_0$ if the imposed macroscopic gradient is sufficiently large, and even rendering the rough fracture up to 10 times more permeable than a smooth fracture of identical mean aperture.

### 1. Introduction

The comprehension of the hydraulic, mechanical, and chemical behaviors of fluids in geological formations is fundamental to the success of subsurface techniques aimed at a variety of applications: resources recovery (Ciriello et al., 2021; Suleimanov et al., 2011), geothermal exploitation (Sanner et al., 2003), carbon sequestration (Leung et al., 2014), and soil reclamation (Mulligan et al., 2001). The decreasing prospects of new oil reservoir discoveries are stimulating companies to invest in unconventional reservoir exploitation (Sheng et al., 2019) and enhanced oil recovery (EOR) to maximize the recovery factor of mature oilfields (Muggeridge et al., 2014). Reservoir stimulation via hydraulic fracturing is a well-established approach to increase formation permeability, which allows extending existing reserves. In unconventional reservoirs, where the pore space is poorly connected, induced stimulation permits the production of oil and gas from formations of low permeability (e.g., shale) (Curtis, 2002), reactivating natural fractures (Gale et al., 2014), and generating new ones (Cipolla et al., 2010). Induced stimulation is also utilized in enhanced geothermal systems (EGS), which constitutes an innovative power system (Lu, 2018): they involve the injection of fluids in artificially fractured hot rocks to exploit an abundant renewable heat source. CO₂-based EGS (Aminu et al., 2017) or CO₂-EOR (Dowell et al., 2017) in fractured geological formations has also been proposed to offset the costs of the subsurface storage of CO₂, which is currently considered a viable approach to reduce anthropogenic CO₂ emissions, responsible for two thirds of the increased greenhouse effects (Leung et al., 2014).

These applications have led in the last 10 years to renewed scientific interest for flow in subsurface porous media, in particular, fractured porous media. In igneous rocks in particular, and more generally in low-permeability formations, fractures provide preferential pathways of high conductance with respect to the almost impervious rock matrix. The fractures are organized in connected networks (Bonnet et al., 2001), and the overall hydraulic behavior of the medium subjected to Newtonian flow is mainly governed by their connectivity (Bour & Davy, 1998) and by the distribution of fracture permeabilities throughout the network (de Dreuzy et al., 2002).
In discrete fracture networks (DFN), the parallel plate conceptualization has traditionally been used to model Newtonian flow in a single fracture using the simplified geometry of two planar and parallel fracture walls, leading to a fracture permeability proportional to the square of its aperture. However, fracture wall topographies are in fact rough as a result of the fracturing process and exhibit self-affine and long-range-correlated topographies (Brown, 1987; Brown & Scholz, 1985; Schmittbuhl, Schmitt, & Scholz, 1995; Schmittbuhl, Vilotte, & Roux, 1995). This geometric property is transferred to the aperture field at scales smaller than the correlation length between the two fracture walls (Brown, 1995). The resulting spatial variability of the aperture field controls the heterogeneity of Newtonian flow, all the more as the fracture is more closed (i.e., as aperture fluctuations relative to the mean aperture are larger) (Brown, 1987; Glover et al., 1998; Méheust & Schmittbuhl, 2000, 2001), while the ratio of the correlation length to the fracture size controls the impact of the flow heterogeneity on the fracture's hydraulic behavior (Méheust & Schmittbuhl, 2003). Moreover, flow heterogeneity below the fracture scale can in some cases modify the (Newtonian) flow connectivity at the network scale, thus strongly impacting the permeability of the fractured formation (de Dreuzy et al., 2012).

Consequently, Newtonian flow in a single rough fracture has been the topic of many past studies, some of which are cited above. And though simple deterministic wall geometries (Di Federico, 1997; Elsworth & Goodman, 1986; Zimmerman & Bodvarsson, 1996) or an aperture probability density function without spatial correlations (Felisa et al., 2018; Lenci & Di Federico, 2020) allow for analytical or semi-analytical investigation of the flow, the study of the aforementioned realistic geometries must rely on numerical modeling. A vast number of studies, including the seminal work of Brown (1987), have relied on solving the Reynolds equation. The latter work states that the local flux, defined as the integral of the fluid velocity over the local fracture aperture, is conservative and can be expressed in terms of the local in-plane pressure gradient according to a local cubic law, that is, a Darcy law where the local fracture transmissivity is proportional to the cube of the local aperture (Brown, 1987; Zimmerman & Bodvarsson, 1996). Other studies, following the recent increases in computational power, have simulated Newtonian flow in the 3D space between the fracture walls (among the first such studies, see Mourzenko et al. (1995); Brush and Thomson (2003)), by means of commercial or open source software, either able to numerically solve the Navier-Stokes equation in 3D flow domains or adopting the Lattice Boltzmann method (LBM) (Wang et al., 2016) to obtain the velocity field. They have rightfully pointed to the moderate underestimation, by the Reynolds equation, of the fracture transmissivity's deviation from that of a smooth fracture of identical mean aperture; this underestimation results from the inability to model out-of-plane flow tortuosity. However, such 3D simulations are computationally expensive, limiting applications to the study of a few fracture realizations, and making a full stochastic analysis impossible. The lubrication theory, which can be assumed to be valid as long as the variation of the aperture field $w$ is sufficiently smooth ($\nabla w \ll 1$) and the fluid can be considered in creeping motion ($Re \ll 1$), reduces the flow problem to two dimensions and leads, for Newtonian flow, to the Reynolds equation. This approximation thus allows for very efficient numerical solvers. Since the variability in the hydraulic behaviors of rough fractures described by the same statistical geometrical parameters (ratio of mean aperture to length, ratio of apertures' standard deviation to mean aperture, and ratio of correlation length to length) can be very large (in particular when the correlation length is not much smaller than the fracture size; Méheust & Schmittbuhl, 2001, 2003), determining a generic hydraulic behavior for a given set of such geometric parameters requires a stochastic Monte-Carlo analysis, in which the use of an efficient lubrication-based solver comes very handy.

In this paper, the fluids of interest are those used in the aforementioned applications; they are water-based but contain surfactants or macrom polymers (e.g., xanthan gum) that boost fluid viscosity and lower surface tension. Various additives can be introduced depending on the application and the geological formation: (a) crosslinkers to increase molecular weight and proppant-carrying capacity; (b) friction reducers to increase pump efficiency and keep proppant particles evenly distributed in the suspension; (c) breakers to provide rheoreversibility, favoring flowback or fluid-disposal process; (d) biocides to kill bacteria; and (e) gellants to adjust the viscosity and form a gel (Barati & Liang, 2014). A vast range of fluids have been adopted in the various steps of the hydraulic fracturing process: natural polysaccharides, synthetic polymer solutions, organic or chemical gels, emulsions, muds, micellar surfactant solutions, and aqueous physical (Barbati et al., 2016). The complex microstructure of these fluids inevitably induces a non-Newtonian mechanical behavior. Solid particles such as proppants will typically lead to jamming effects and thus to shear thickening of the flow, while most additives will provide the fluid with a soft matter-type microstructure (colloids, polymers, droplets, bubbles, and micelles), resulting in a shear-thinning behavior (hereinafter ST). This behavior is by far prevalent in fluids engineered for subsurface applications
(Barbati et al., 2016; Economides, 2000). Recently, nontoxic inexpensive CO₂-based alternatives to conventional fluids have been considered because of their high viscosity and rheoreversibility; these fluids show the desirable ST behavior without incorporating toxic chemical agents (Chen et al., 2016; Jung et al., 2015; Li et al., 2017).

Hence, better understanding the flow of ST fluids in rough fractures is crucial to predicting the outcome of various subsurface operations. The coupling of ST rheology and domain heterogeneity is known to produce higher flow localization and channeling in porous media (Sullivan et al., 2006; Zami-Pierre et al., 2016) as compared to Newtonian flow in the same geometry. However, the topic of ST fluid flow in rough fracture geometries has been little addressed in the literature. Morris et al. (2015) proposed a lubrication-based model for the flow of Herschel-Bulkley fluids in simple deterministic geometries, while Lavrov (2013a, 2013b) used a similar model to study the flow of power law (PL) fluids in relatively small (33 × 33) size realizations of self-affine rough fractures (without a correlation length) and provided considerations on the impact of the PL exponent on the flow geometry. Zhang et al. (2019) performed 3D CFD simulations of ST laminar flow, focusing on the impact of the fluid rheology on inertial effects and the transition to a Forchheimer-like regime. To the best of our knowledge, the flow of ST fluids featuring a rheological transition from a Newtonian plateau at low-shear rates to a PL ST behavior at high-shear rates in large-scale and realistic fracture geometries (featuring a correlation length) has never been addressed in the literature, and no code implementation that allows easily performing Monte Carlo simulations over a large number of aperture fields has been proposed.

In this study, we propose a model simulating the flow of an ST fluid in a variable aperture fracture based on the depth-averaged lubrication theory and considering a three-parameter Ellis constitutive model for the fluids' rheology. The latter accounts for both low-shear rate quasi-Newtonian and high-shear rate ST behaviors. This approach yields a quasilinear partial differential equation (PDE) that is the sum of a linear heterogeneous diffusion term and a nonlinear generalized (heterogeneous) *p*-Laplace operator. The latter is a mathematically well-studied prototype of nonlinear PDE and is at the heart of many models of nonlinear diffusion (Philip, 1961), including non-Newtonian flows, turbulent filtration and reactive-diffusion (Diaz & De Thelin, 1994), and nonlinear Darcy flows (Firdaouss et al., 1997). The reader is referred to Benedikt et al. (2018) for a historical perspective. Similar mathematical models have been adopted in groundwater hydraulics (e.g., (Cainelli et al., 2012) and pore network modeling (Perrin et al., 2006)). Recently, the interest on *p*-Laplacian operators has increased, resulting in significant developments of efficient numerical solvers (Facca et al., 2021; Loisel, 2020). Taking into consideration the aim of running Monte Carlo simulations, we propose an efficient Finite Volume (FV)-based discretization for the nonlinear PDEs written as a generalized graph *p*-Laplacians in a uniform square mesh. The proposed scheme turns out to be a slight variant of the lowest order method of the family of schemes analyzed in Andreianov et al. (2004). The developed approach has the advantage that it leads to a symmetric Jacobian, thus enabling the implementation of a very efficient Newton-Krylov iteration for the solution of the ensuing nonlinear system, employing the Preconditioned Conjugate Gradient (PCG) linear solver. A parameter continuation method is added to increase robustness for wider ranges of *p*. This allows addressing a vast range of geometrical and rheological parameters of practical interest, in particular, small rheological indices (i.e., ST exponents) of order 0.1, without resorting to more sophisticated and computationally expensive methods such as those proposed in Facca et al. (2021). The ensuing solver achieves a favorable compromise between accuracy and computational cost that enables us to address the study of the flow behavior in the fracture within a Monte-Carlo framework as discussed above. The impact of the ST rheology on the spatial distribution of the velocity field in the fracture plane and on the resulting fracture scale hydraulic behavior are then investigated. The solver is two to three orders of magnitude faster than the few alternative lubrication-based algorithms presented in the literature to solve the flow of ST fluids in fractures. In fact, the Newton algorithm shows a superlinear rate of convergence, while only a linear convergence can be achieved by the Picard-based alternatives adopted by Lavrov (2013a) and Morris et al. (2015) to simulate the flow of PL fluids. The numerical scheme permits overcoming the actual computational limits so as to allow for a stochastic analysis in a reasonable amount of time, considering the large meshes (1024 × 1024 squares) required by the fracture statistical description, and accounting for the entire ranges of the parameters of practical interest.

Based on these premises, the organization of the paper follows naturally: Section 2 describes the geometric properties of fracture apertures and how realistic synthetic aperture fields can be generated. Section 3 provides the derivation of the generalized Reynolds equation for an Ellis fluid in a variable aperture fracture; Section 4 presents the implementation of this theoretical model into a novel numerical solver, while Section 5 reports on
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Representation of a fracture's geometry with the top wall positioned above the bottom wall. $z_u(x)$ and $z_1(x)$ described the topographies of the rough walls with respect to their mean planes. The separation between those planes is the mechanical aperture $\langle w \rangle$. Large-scale spatial fluctuations of the two walls' topographies are identical, while their small-scale roughnesses are different. (a) Three-dimensional view. (b) Two-dimensional view in the plane indicated by the red lines in panel (a).</figcaption></figure>

numerical results regarding the impact of the fluid rheology on fracture flow and the fracture's hydraulic aperture; Section 6 lists our conclusions and prospects for future work.

## 2. Synthetic Aperture Fields

As depicted in Figure 1, the aperture field $w(x)$ of a fracture is typically estimated as the distance between the two rough fracture walls, whose mean planes are parallel. Denoting the mean aperture of the fracture, and $z_u(x)$ and $z_1(x)$, respectively, the upper and lower topographies of the walls with respect to their mean planes $(x = (x_1, x_2)^T$ being the position vector in the fracture plane), $w(x)$ is simply

$$w(x) = z_u(x) - z_1(x) + \langle w \rangle. \qquad (1)$$

Typically, the walls of a fracture surface can be described as isotropic self-affine topographies, exhibiting long-range correlations (Bouchaud, 1997; Schmittbuhl, Vilotte, & Roux, 1995) characterized by an exponent $H$ (the so-called Hurst exponent), up to their linear size. For a fresh fracture, the two wall topographies are identical, so the fracture aperture depends on the shear slip history of the walls. Shear slip induces an anisotropic decorrelation of the two topographies along the shear direction and at scales smaller than the slip distance (Plouraboué et al., 1995). In geological fractures, however, the aperture field is usually measured to be isotropic, but also with a decorrelation between the wall topographies at a scale smaller than a correlation length $L_c$ (Brown, 1995). Above that scale, the walls can be considered to be matched, but below it, the aperture field, being a linear combination of two independent self-affine topographies, is also self-affine. In the Fourier space, this translates into a power spectral density in the form of a PL of exponent $-2(H + 1)$ at scales smaller than $L_c$ and flat above that length (see Figure 2a). Note that the Hurst exponent $H$ has been measured to a quasi-universal value

of around 0.8 over a wide range of scales and materials, including ceramics, metals, and rocks, such as granite and basalt (Bouchaud et al., 1990), with a few exceptions, such as sandstones ($H \sim 0.45$, see Boffa et al. (2000)). This property of the fracturing process in sandstones has been attributed to its intergranular nature.

Recursive algorithms can be implemented to generate rough surfaces with spatial correlations (by successive random addition (Liu et al., 2004)), but the generation of large fractures with these algorithms is inefficient. Alternatively, we adopt a method proposed by Méheust and Schmittbuhl (2003) and relying on the Fourier spectrum's properties (see also Barnsley et al. (1988)). Generating a random Fourier transform of the aperture field with the spatial correlation properties described above allows us to obtain an appropriate synthetic aperture field without the need to generate the two wall topographies (Méheust & Schmittbuhl, 2003). The procedure is very fast as it relies on the fast Fourier transform. We start with a spatially uncorrelated white noise and multiply the modulus of its Fourier transform by the modulus of the wave number, $|k| = \sqrt{k_{x_1}^2 + k_{x_2}^2}$, raised to the power $-(H + 1)$ (Figure 2c). The resulting aperture field can then be scaled and translated so as to tune its mean value and standard deviation, with consequent negative values set to 0, corresponding to an assumption of perfect plastic closure. Zero values are then further regularized to a minimum cutoff $w_0 = 10^{-8}$ m to avoid ill-posedness of the flow model. The aperture field generation is implemented according to Algorithm 1 in Appendix A. Figure 2b shows two representations of the aperture field of the same fracture at two different degrees of closure, or coefficients of variations, $\sigma_w/\langle w \rangle$.

## 3. Generalized Reynolds Equation

Engineered fluids used in subsurface operations typically display an overall shear-thinning (ST) behavior, wherein the apparent viscosity $\mu$ of the fluid is a decreasing function of shear rate ($\dot{\gamma}$), that is, the fluid thins with increasing deformation. More specifically, the typical flow curve $\mu = \mu(\dot{\gamma})$ for an ST fluid exhibits two
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Mean radial profile (obtained by angular integration) of the 2D Fourier transform of a fracture's aperture field, displaying a low-frequency cutoff for wave numbers smaller than $k_c = 2\pi/L_c$, and a power law (PL) trend of negative exponent $-2(1 + H)$ for higher frequencies. (b) Aperture field representations for different values of the closure: 0.5 (upper panel) and 1 (bottom panel). (c) Representation of the wave number to the power $-(1 + H)$. The set of parameter adopted for the generations is $L = 80$ cm, $L_c = 10$ cm, $\langle w \rangle = 1$ mm, $H = 0.8$.</figcaption></figure>

Newtonian plateaus, one of high viscosity at low-shear rates, the other one of low viscosity at high-shear rates. These plateaus are separated by a ST regime starting for shear rates higher than a critical threshold value $\dot{\gamma}_c$. The Carreau-Yasuda (CY) five-parameter model (Yasuda et al., 1981) and its original four-parameter version, the Carreau model (see Appendix B for both), exhibit both plateaus and are the most common constitutive laws able to fully capture such a rheology. On the other hand, the PL two-parameter model represents the simplest non-Newtonian constitutive equation, but is able to reproduce only the ST behavior typical of the intermediate range of shear rates in a CY or Carreau fluid (Escudier et al., 2001). The CY model is commonly used in CFD applications, but is difficult to handle analytically. Conversely, the more tractable PL model oversimplifies the constitutive law because of its unrealistic apparent viscosity, which is unbounded as the shear rate approaches zero and tends to zero at high-shear rates. When dealing with creeping flow in rough fractures, the velocities are usually sufficiently low for high-shear rates, corresponding to the low viscosity plateau of the rheological curve, to be rare. In other words, the apparent viscosity attains the low viscosity plateau value in a limited number of spatial locations, so that a rheological model only accounting for the low-shear rate plateau and the PL decrease at higher shear rates is an acceptable simplification in practical cases. Such a model is perfectly described by the three-parameter Ellis constitutive law, first introduced by Bird et al. (1987), where, at variance with most generalized Newtonian fluids, the apparent viscosity depends on the local shear stress $\tau$ according to

$$ \mu = \frac{\mu_0}{1 + \left(\frac{\tau}{\tau_{1/2}}\right)^{\frac{1}{n}-1}} \quad (2) $$

As shown in Figure 3, $\mu_0$ is the plateau viscosity (having the same dimensions as the dynamic viscosity of a Newtonian fluid), $n$ defines the PL exponent ($n-1$) of the ST regime, and the characteristic $\tau_{1/2}$, defined by $\mu(\tau_{1/2}) = \mu_0/2$,

shapes the transition between the two tendencies. The Ellis rheologic parameters depicted in Figure 3 and listed in Table 1 are those of four real ST fluids taken from the literature that will be used in the following to demonstrate the code.

To derive the flow rate of an Ellis fluid between parallel plates, one of the main building blocks of the code, we recall the equations governing isothermal flow of an incompressible fluid of density $\rho$, that is, the momentum and continuity equations, reading:

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = \rho \mathbf{g} + \nabla \cdot \mathbf{T}, \quad (3a) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (3b) $$

where $\mathbf{u}$ is the velocity, $\mathbf{g}$ is the body force vector, and $\mathbf{T}$ is the stress tensor. For a generalized Newtonian fluid, the latter can be written as $\mathbf{T} = -p'\mathbf{I} + 2\mu\mathbf{D}$, where $p'$ is the pressure, $\mathbf{I}$ is the identity tensor, $\mathbf{D} = \frac{1}{2}\nabla\mathbf{u} + \frac{1}{2}(\nabla\mathbf{u})^T$ is the strain rate tensor, and $\mu$ is the apparent viscosity, which depends on the scalar shear rate $\dot{\gamma}$, defined from the tensor $\mathbf{D}$ according to $\dot{\gamma} = \sqrt{2\mathbf{D}} : \mathbf{D}$, where : stands for the tensor double-dot product.

In the following, we consider Stokes flow, that is, we assume that the Reynolds number is much smaller than 1 and thus that the nonlinear term in the right-hand side of Equation 3a can be neglected in comparison to viscous forces. For permanent Stokes flow, Equation 3a becomes

$$ [\nabla \cdot (\mu \nabla)] \mathbf{u} = \nabla P, \quad (4) $$

where the reduced pressure is defined as $P = p' + \rho gz$, $z$ being the coordinate along the direction perpendicular to the fracture plane (i.e., the vectorial plane parallel to the mean planes of the fracture's walls).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Dependence of the dimensionless apparent viscosity $\mu/\mu_0$ on the shear rate $\dot{\gamma}$, for the four Ellis fluids listed in Table 1: F1 (yellow solid line), F2 (orange solid line), F3 (red solid line), and F4 (purple solid line). The parameters of these fluids are deduced from data originally interpreted with the Carreau-Yasuda model (see Appendix B); the original Carreau model and its low-shear rate approximation are reported as dashed and dotted lines, respectively. The high-shear rate power law exponents for all solid curves are $n-1$.</figcaption></figure>

We now assume the lubrication approximation, that is, that the gradient of the aperture field is very small ($\nabla w \ll 1$). This means that the typical horizontal scale of variation of the velocity field is much larger than its typical vertical scale of variation. It follows that due to flow incompressibility (Equation 3b), the vertical velocity is negligible with respect to the horizontal velocity and that the projection of Equation 4 on the $z$ direction reduces to $\partial P/\partial z = 0$, so that $P$ can be considered to not depend on $z$. The projection of Equation 4 on the fracture plane is then

$$[\nabla \cdot (\mu \nabla)] \mathbf{u} = \nabla_{\parallel} P, \qquad (5)$$

where $\mathbf{u}$ is considered parallel to the fracture plane and $\nabla_{\parallel}$ is the projection of the gradient operator on that plane. From the aforementioned considerations on the horizontal and vertical typical scales of variations of the velocity field, it follows that the term involving the derivative with respect to $z$ in the left-hand term of Equation 5 dominates the others, so that this equation can be approximated as

$$\frac{\partial}{\partial z} \left( \mu \frac{\partial \mathbf{u}}{\partial z} \right) = \nabla_{\parallel} P, \qquad (6)$$

which holds at any position $\mathbf{x}$ along the fracture plane. Integrating this equation once over $z$ between the vertical mid-position $\bar{z}(\mathbf{x})$ (where the $\partial \mathbf{u}/\partial z$ must be $\mathbf{0}$ by symmetry) and $z$ yields a relation between the shear stress $\tau$ and the horizontal pressure gradient:

$$\tau = \mu \frac{\partial \mathbf{u}}{\partial z} = (z - \bar{z}) \nabla_{\parallel} P. \qquad (7)$$

Replacing $\tau$ by $(z - \bar{z}) \|\nabla_{\parallel} P\|$ in Equation 2 and inserting the resulting expression for the viscosity into Equation 7, we then obtain

$$\frac{\partial \mathbf{u}}{\partial z} = \frac{1}{\mu_0} \left[ 1 + \left( \frac{(z - \bar{z}) \|\nabla_{\parallel} P\|}{\tau_{1/2}} \right)^{\frac{1}{n}-1} \right] (z - \bar{z}) \nabla_{\parallel} P. \qquad (8)$$

The dependence between the vertical velocity profile at position $\mathbf{x}$ in the plane and the local pressure gradient $\nabla_{\parallel} P$ can be obtained by integrating Equation 8 between the value of $z$ at the lower wall and $z$. Introducing a new centered vertical variable $\bar{z} = z - \bar{z}$, we obtain

$$\mathbf{u}(\bar{z}) = \int_{-w/2}^{z} \left[ 1 + \left( \frac{s \|\nabla P\|}{\tau_{1/2}} \right)^{\frac{1}{n}-1} \right] \frac{s \nabla P}{\mu_0} ds \\ = - \left[ \frac{1}{8\mu_0} (w^2 - 4\bar{z}^2) + \frac{n}{n+1} \left( \frac{1}{2^{n+1}\mu_0^n \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}} (w^{\frac{1}{n}+1} - |2\bar{z}|^{\frac{1}{n}+1}) \|\nabla P\|^{\frac{1}{n}-1} \right] \nabla P. \qquad (9)$$

We define the local flux $\mathbf{q}$ as the integral of the fluid velocity along the vertical coordinate and over the local fracture aperture:

$$\mathbf{q} = - \left[ \frac{w(\mathbf{x})^3}{12\mu_0} + \frac{n}{2n+1} \left( \frac{1}{2^{n+1}\mu_0^n \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}} w(\mathbf{x})^{\frac{2n+1}{n}} \|\nabla P\|^{\frac{1}{n}-1} \right] \nabla P. \qquad (10)$$

It turns out that $\mathbf{q}$ is conservative (Zimmerman & Bodvarsson, 1996): $\nabla \cdot \mathbf{q} = 0$, which leads to the following nonlinear lubrication equation:
<table><caption>**Table 1**
Rheologic Parameters $\mu_0$, $n$, and $\tau_{1/2}$ Related to the Four Shear-Thinning Fluids Interpreted With the Ellis Models</caption><thead><tr><th rowspan="2">Fluid ID</th><th rowspan="2">Solution</th><th>μ<sub>0</sub></th><th>τ<sub>1/2</sub></th><th>n</th></tr><tr><th>(Pa · s)</th><th>(Pa)</th><th>(-)</th></tr></thead><tbody><tr><td>F1</td><td>CMC at 0.3 wt%</td><td>0.0510</td><td>4.07</td><td>0.72</td></tr><tr><td>F2</td><td>CMC at 0.5 wt%</td><td>0.2203</td><td>2.50</td><td>0.51</td></tr><tr><td>F3</td><td>CMC at 1.0 wt%</td><td>2.9899</td><td>5.14</td><td>0.40</td></tr><tr><td>F4</td><td>VES</td><td>49</td><td>1.07</td><td>0.10</td></tr></tbody></table>

Note. The fluids considered are: three different Carboxymethylcellulose (CMC) water-based solutions from Sousa et al. (2005) and a viscoelastic surfactant (VES) from Moukhtari and Lecampion (2018). Experimental data have been originally fitted with the CY model (parameters listed in Table B1) and interpreted with the Ellis model as described in Appendix B.

$$-\nabla \cdot \left[ \left( \frac{w(\mathbf{x})^3}{12\mu_0} + \frac{n}{2n+1} \left( \frac{1}{2^{n+1}\mu_0^n \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}} w(\mathbf{x})^{\frac{2n+1}{n}} \|\nabla P\|^{\frac{1}{n}-1} \right) \nabla P \right] = 0. \quad (11)$$

This equation reduces to the classical, linear, Reynolds equation for a Newtonian fluid ($n = 1$), the solution depending only on the fracture aperture in that case. For $n < 1$, conversely, the non-Newtonian pressure field also depends on the fluid's rheology, which renders the problem nonlinear.

Note that $\mathbf{q}$ acts as the velocity of an incompressible 2D flow occurring in the fracture plane. The mean 3D velocity at position $(x)$ can be computed as

$$\bar{\mathbf{u}} = \frac{\mathbf{q}}{w} = \frac{1}{w} \int_{-\frac{w}{2}}^{\frac{w}{2}} \mathbf{u}(z) \ dz = - \left[ \frac{w^2}{12\mu_0} + \frac{n}{2n+1} \left( \frac{1}{2^{n+1}\mu_0^n \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}} w^{\frac{1+n}{n}} \|\nabla P\|^{\frac{1}{n}-1} \right] \nabla P . \quad (12)$$

The fraction of the fracture plane that is concerned by the nonlinear rheology, and to which extent it is, can be visualized from maps of the local depth-averaged apparent viscosity, which must be estimated a posteriori (i.e., once the flow velocity field has been computed):

$$\bar{\mu}(\mathbf{x}) = \frac{\mu_0}{w} \int_{-w/2}^{w/2} \left[ 1 + \left( \frac{\|\nabla P\|z}{\tau_{1/2}} \right)^{\frac{1}{n}-1} \right]^{-1} dz. \quad (13)$$

This 2D viscosity field allows introducing the Reynolds number for a generalized Newtonian fluid ($Re_G$), as

<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Fracture plane representation with boundary conditions. (b) Finite volume method: five-point stencil. Pressures are defined at the center of each cell, while the fracture's local apertures are defined on the edges.</figcaption></figure>
$$Re_G = \frac{\rho \langle \bar{u} \rangle \langle w \rangle}{2 \langle \bar{\mu} \rangle} . \qquad (14)$$

where $\rho$ is the fluid density (in the following, it will be assumed equal to the water density $10^3$ kg/m$^3$) and $\langle \bar{\mu} \rangle$ is the average over the fracture plane of the depth-averaged apparent viscosity; similarly, $\langle \bar{u} \rangle$ is the average over the fracture plane of the depth-averaged velocity (or, equivalently, the average of the 3D velocity over the entire fracture volume). That definition of the Reynolds number allows us to verify a posteriori that the Stokes flow assumption is still verified in the computed ST flow, a condition that is requested for the lubrication approach to be valid.

Note that an estimate of the critical macroscopic pressure gradient magnitude $\overline{\nabla P} = \|\overline{\nabla P}\|$ above which non-Newtonian effects begin to appear in the flow, which we denote $\overline{\nabla P}_c$, can be obtained by considering the parallel plate configuration whose aperture is equal to the rough fracture's mean aperture. From the rheological law (2), it is seen that the viscosity begins to deviate from the Newtonian plateau viscosity $\mu_0$ when the largest shear stress locally occurring in the fracture, $\tau_{max}$, is on the same order as the crossover shear stress $\tau_c$ in the CY rheological law relating the apparent viscosity to the shear rate (see Equation B1). The crossover shear stress $\tau_c$ can be estimated numerically from the transition shear stress $\tau_{1/2}$ and plateau viscosity $\mu_0$ of the Ellis model (Equation 2) (see Appendix B).

From Equation (7), it is possible to see that at any position $x$ in the fracture plane, the largest shear stress (in absolute value) occurs at the walls for any fluid. An estimate of it can be obtained from the parallel plate configuration of aperture equal to the mean fracture aperture $\langle w \rangle$ and is equal to $\tau_{max} = \langle w \rangle \overline{\nabla P}/2$. Hence, the condition $\tau_{max} \approx \tau_c$ provides the following crossover value for the applied macroscopic pressure gradient:

$$\overline{\nabla P}_c = \frac{2\tau_c}{\langle w \rangle} . \qquad (15)$$

If the imposed $\overline{\nabla P}$ is much smaller than $\overline{\nabla P}_c$, the flow is expected to be Newtonian, whereas if it is much larger than the crossover pressure gradient, the flow can be expected to be strongly ST.

# 4. Numerical Modeling

## 4.1. Finite Volume Scheme Formulation

We consider an ST fluid flowing through a fracture whose projection on its mean plane is square ($\Omega = (0, L) \times (0, L)$) as depicted in Figure 1. As depicted in Figure 4a, the flow in the fracture is driven from left to right by an externally applied pressure drop $\Delta P$ corresponding to an average pressure gradient magnitude $\overline{\nabla P} = \Delta P/L$. Without loss of generality, we assume the following Dirichlet boundary conditions over $\partial \Omega_D$ (see Figure 4a): for $0 < x_2 < L$, the pressure is taken to be zero at the downstream end ($P(L, x_2) = 0$), while upstream $P(0, x_2) = L \overline{\nabla P}$ is assumed. The upper and lower sides of the domain ($\partial \Omega_N$) are both assumed impervious, imposing Neumann boundary conditions such that $\partial_2 P(x_1, 0) = \partial_2 P(x_1, L) = 0$, for $0 < x_1 < L$.

The generalized Reynolds Equation 11 results in a quasilinear PDE problem of the second order written as

$$-\nabla \cdot \left[ \left( a(\mathbf{x}) + b(\mathbf{x}) \|\nabla P\|^2 \right) \nabla P \right] = 0 \quad \mathbf{x} \in \Omega \qquad (16a)$$

$$P = P_0 \quad \mathbf{x} \in \partial \Omega_D \qquad (16b)$$

$$\nabla P \cdot \mathbf{v} = 0 \quad \mathbf{x} \in \partial \Omega_N \qquad (16c)$$

where $\mathbf{x}$ is the position vector, $\mathbf{v}$ is the outward unit normal vector defined on the boundary, and the coefficients $a(\mathbf{x})$ and $b(\mathbf{x})$ are defined as

$$a(\mathbf{x}) = \frac{w(\mathbf{x})^3}{12\mu_0} \qquad (17)$$

$$\text{and} \quad b(\mathbf{x}) = \frac{n}{2n+1} \left( \frac{1}{2^{1+n} \mu_0^n \tau_{1/2}^{1-n}} \right)^{\frac{1}{n}} w(\mathbf{x})^{\frac{2n+1}{n}} . \qquad (18)$$
Equation 16a is the sum of two terms, the first one involving a linear heterogeneous Laplacian differential operator and the second one involving a nonlinear heterogeneous p-Laplacian operator with $p = 1/n + 1$. To ensure well-posedness of the problem (Equations 16a, 16b, and 16c), the two coefficients $a$ and $b$ must be strictly positive, and the exponent $p$ must fall within the range $1 < p < \infty$. The latter condition implies that the ST index $n$ remains in the range $0 < n < +\infty$; this well-posedness region contains the physical range of $n$ ($0 < n \le 1$). Hence, the only condition that needs to be imposed is that $a$ and $b$ be strictly positive and thus, that the aperture be strictly positive everywhere. This is the reason for the introduction of a nonzero aperture threshold $w_0$ enforced in the contact zones between the two fracture walls in Algorithm 1 of Section 3.

Next, we describe the FV-based discretization of Equation 16a. The design of the scheme is grounded on the need to make the overall solver computationally efficient so that relatively large ensembles of flow data can be built with a relatively large spatial resolution of the stochastic aperture realizations. The development of the algorithm proceeds as follows. In the first step, the fracture domain $\Omega$ is partitioned into a set $Q$ of $4^M$ nonoverlapping square control volumes $\mathcal{K}$ (i.e., $\forall j \mathcal{K}_j \subset \Omega$ and $\cup_j \mathcal{K}_j = \Omega$), where $M$ indicates the mesh level. Let $\mathbf{x}_j$ be the center of the $\mathcal{K}_j$ control volume and let $\sigma(j) = \{E, W, N, S\}$ be the set of neighboring control volumes. The edge between two cells can be identified by

$$ \mathcal{K}_i \cap \mathcal{K}_j = \begin{cases} e_{ij} & \text{if } i \in \sigma(j) \\ \emptyset & \text{otherwise.} \end{cases} \qquad (19) $$

Note that $|e_{ij}| = h$ is the length of a side of the square control volume and is constant across the mesh. Equation 16a is integrated over each cell $\mathcal{K}_j$ to obtain up to an irrelevant plus or minus sign:

$$ \forall j \quad \int_{\mathcal{K}_j} \nabla \cdot \left[ \left( a(\mathbf{x}) + b(\mathbf{x}) \|\nabla P\|^{\frac{1}{n}-1} \right) \nabla P \right] d\mathbf{x} = 0 \qquad (20) $$

Applying the divergence theorem in each finite volume $\mathcal{K}_j$ (Figure 4b), leads to

$$ \int_{\partial \mathcal{K}_j} \left( a(\mathbf{x}) + b(\mathbf{x}) \|\nabla P\|^{\frac{1}{n}-1} \right) \nabla P \cdot \mathbf{v} \, ds = \sum_{i=1}^{4} \int_{e_{ij}} \mathbf{q}_i^{(j)} \cdot \mathbf{v}_i \, ds = 0, \quad j = 1, \dots, 4^M, \qquad (21) $$

where $\mathbf{v}_i$ is the outward unit normal vector on the cell edge $e_{ij}$. Numerical integration by means of the mid-point quadrature rule over $e_{ij}$ leads to a second-order approximation of the flux through each finite volume edge:

$$ \sum_{i=1}^{4} \int_{e_{ij}} \mathbf{q}_{i}^{(j)} \cdot \mathbf{v}_{i} \, ds \cong +h \sum_{i=1}^{4} \mathbf{q}_{i}^{(j)} \cdot \mathbf{v}_{i} = 0 \quad j = 1, \dots, 4^{M}. \qquad (22) $$

where the edge fluxes are defined as $\mathbf{q}_i^{(j)} = (a_i^{(j)} + b_i^{(j)} \|\nabla P_i^{(j)}\|^{\frac{1}{n}-1}) \nabla P_i^{(j)}$. We then consider a numerical approximation of the pressure gradient's normal component on the edge, defined as

$$ |\nabla P_i^{(j)} \cdot \mathbf{v}_i^{(j)}| = \left| \frac{P_i - P_j}{h} \right| . \qquad (23) $$

To maximize computational efficiency, we use this first-order approximation as well to evaluate the Euclidean norm of the gradient evaluated at edge $e_{ij}$, that is, $\|\nabla P_i^{(j)}\| \approx |P_i - P_j|/h$. This scheme results in a consistent approximation of the p-Laplacian in Equations 16a as long as we maintain our discretization on a square mesh. Indeed, the approximation (23) is exact for affine functions on $\mathcal{K}$ (Andreianov et al., 2004) and the gradient approximation based on graph-Laplacians converges to the continuous counterpart for both the heterogeneous linear Laplace equation (Singer, 2006) and the nonlinear p-Laplace equation (Calder, 2018). Note that this FV approximation, being typical of graph-p-Laplacians, has already been used in general graph-based applications, and, in the context of hydrological sciences, it has been used for pore-scale modeling in subsurface hydrology (see e.g., Balhoff and Thompson (2006); Sochi and Blunt (2008)). To the best of our knowledge, this is the first time that such a scheme is used to model a continuous problem based on a p-Laplacian operator. Since the approximation of the norm of the gradient causes the convergence rate to be in any case only first order, we chose to estimate the fracture aperture at the cell interface in the expressions for the coefficients $a_i^{(j)}$ and $b_i^{(j)}$ using the arithmetic mean:
$$w_i^{(j)} = \frac{w_i + w_j}{2}. \qquad (24)$$

Note that the use of the harmonic mean would preserve the energy of the scheme. In spite of this, we opted for the arithmetic mean because it is less affected by the ill-conditioning induced by the aperture field's variability, which may span several orders of magnitude (Mazzia et al., 2011).

## 4.2. Model Implementation

In this section, we describe the implementation details of our FV formulation of the nonlinear lubrication Equation 16a. For each cell $K_j$, we have:

$$\left( a_E^{(j)} + b_E^{(j)} \left| \frac{P_E - P_j}{h} \right|^{\frac{1}{n}-1} \right) (P_E - P_j) + \left( a_W^{(j)} + b_W^{(j)} \left| \frac{P_W - P_j}{h} \right|^{\frac{1}{n}-1} \right) (P_W - P_j) + \\ + \left( a_N^{(j)} + b_N^{(j)} \left| \frac{P_N - P_j}{h} \right|^{\frac{1}{n}-1} \right) (P_N - P_j) + \left( a_S^{(j)} + b_S^{(j)} \left| \frac{P_S - P_j}{h} \right|^{\frac{1}{n}-1} \right) (P_S - P_j) = 0, \qquad (25)$$

valid for $j = 1, ..., 4^M$. After all the terms have been put together, the following system of nonlinear equations is obtained:

$$\mathbf{F}(\mathbf{p}) = \mathbf{A}(\mathbf{p})\mathbf{p} - \mathbf{f} = 0, \qquad (26)$$

where $\mathbf{p}$ is the unknown pressure vector, and, for any given $\mathbf{p}$, the matrix $\mathbf{A}(\mathbf{p})$ is symmetric positive definite and pentadiagonal with coefficients $A_{ij}$ given by

$$A_{ij} = \begin{cases} -\sum_{k \in \sigma(j)} \left( a_k^{(j)} + b_k^{(j)} \left| \frac{P_k - P_j}{h} \right|^{\frac{1}{n}-1} \right) & \text{if } i = j; \\ a_i^{(j)} + b_i^{(j)} \left| \frac{P_i - P_j}{h} \right|^{\frac{1}{n}-1} & \text{if } i \in \sigma(j); \\ 0 & \text{otherwise,} \end{cases}. \qquad (27)$$

The components $f_j$ of the right-hand-side vector $\mathbf{f}$ implementing the nonzero Dirichlet boundary conditions are given by

$$f_j = \begin{cases} - \left( a_W^{(j)} + b_W^{(j)} \left| \frac{P_W - \overline{\nabla P L}}{h} \right|^{\frac{1}{n}-1} \right) & \text{if } j = k(M-2) \text{ for } k = 1, \dots, t-2, \\ 0 & \text{otherwise.} \end{cases} \qquad (28)$$

The system of nonlinear algebraic Equation 26 is solved by the Newton method, starting with an initial guess $\mathbf{p}_0$ that is the pressure field corresponding to Newtonian flow (exponent $n = 1$). Denoting with $k$ the Newton iteration number and with $J(\mathbf{p})$ the Jacobian matrix, Newton's scheme takes on the form:

$$\left\{ \begin{array}{l} J(\mathbf{p}_k) \mathbf{s}_k = -\mathbf{F}_k(\mathbf{p}_k); \\ \mathbf{p}_{k+1} = \mathbf{p}_k + \mathbf{s}_k. \end{array} \right. \qquad (29)$$

where for the $k$th Newton iteration the Jacobian matrix is formally defined as

$$J(\mathbf{p}) = \mathbf{F}'(\mathbf{p}) = \mathbf{A}(\mathbf{p}) + \mathbf{A}'(\mathbf{p})\mathbf{p}, \qquad (30)$$

where the $\mathbf{A}'$ symbol denotes the derivative of matrix $\mathbf{A}$ with respect to the pressures at each mesh node.

The calculation of the elements of the Jacobian matrix is provided in Appendix C with its generic $ij$-th element reported here:
<table><caption>**Table 2**

Parameters for Numerical Integration</caption><tr><td>Maximum number of PCG iterations (k<sub>max</sub>)</td><td>10<sup>3</sup></td></tr><tr><td>Drop tolerance (ε)</td><td>10<sup>-4</sup></td></tr><tr><td>Diagonal shift compensation coefficient (δ)</td><td>10<sup>-3</sup></td></tr><tr><td>PCG tolerance upper limit (η<sub>max</sub><sup>lin</sup>)</td><td>0.99</td></tr><tr><td>PCG tolerance scaling factor (η<sub>0</sub>)</td><td>0.90</td></tr></table>

$$J_{ij}(\mathbf{p}) = \begin{cases} \sum_{s \in \sigma(i)} \left( a_s^{(i)} + \frac{1}{n} b_s^{(i)} \left| \frac{P_s - P_i}{h} \right|^{\frac{1}{n}-1} \right) & \text{if } i = j; \\ -a_i^{(j)} - \frac{1}{n} b_i^{(j)} \left| \frac{P_i - P_j}{h} \right|^{\frac{1}{n}-1} & \text{if } i \in \sigma(j); \\ 0 & \text{otherwise.} \end{cases} \quad (31)$$

The resulting $J(\mathbf{p})$, at any given $\mathbf{p}$, is a symmetric, pentadiagonal, and positive definite matrix. We would like to remark here that the symmetry of $J(\mathbf{p})$ is a consequence of our specific choice of discretization method. Indeed, the

symmetry of $J(\mathbf{p})$ arises from the symmetry of the particular operation of differentiation of the absolute value function. This symmetry allows the use of the Preconditioned Conjugate Gradient (PCG) method for the solution of the Newton linear system (29), with obvious improvements in computational performance. Row equilibration of the Jacobian matrix is achieved via diagonal scaling to improve the problem conditioning. This means that the left- and right-hand-side terms of Equation 29 are both left-multiplied by $D^{-1/2}(\mathbf{p})$, $D = \text{diag}(J)$ being a matrix consisting solely of the Jacobian's main diagonal. Applying the scaling transformation results in a more accurate solution and reduces computational time. We obtain

$$[D(\mathbf{p})]^{-1/2} J(\mathbf{p}) [D(\mathbf{p})]^{-1/2} [D(\mathbf{p})]^{1/2} \mathbf{s} = -[D(\mathbf{p})]^{-1/2} F(\mathbf{p}), \quad (32)$$

The preconditioner of choice is a variable-fill-in incomplete Cholesky factorization. This preconditioning constitutes a problem transformation, which results in a smaller condition number, and consequently, improves the rate of convergence of the PCG. The entire procedure detailed above is encapsulated within an inexact Newton-Krylov (INK) framework (Kelley, 1987). To avoid oversolving, the PCG exit tolerance is decreased as iterations progress using the following rule:

$$\eta_k^{\text{lin}} = \begin{cases} \eta_{\text{max}}^{\text{lin}} & \text{if } k = 0; \\ \min \left\{ \eta_{\text{max}}^{\text{lin}}, \eta_0 \times \left( \frac{\|\mathbf{F}(\mathbf{p}_k)\|_2}{\|\mathbf{F}(\mathbf{p}_{k-1})\|_2} \right)^2 \right\} & \text{otherwise,} \end{cases} \quad (33)$$

where the choice of $\eta_{\text{max}}^{\text{lin}} = 0.99$ and $\eta_0 = 0.9$ allows extracting most of the information from the inner iteration and guarantees an upper limit to the sequence $\{\eta_k^{\text{lin}}\}$. Table 2 summarizes these choices with the related values of the parameters.

When the exponent of the Ellis rheological model ($n$) is small, and/or the applied macroscopic pressure gradient ($\overline{\nabla P}$) is large, the pressure field of the solution for the Newtonian rheology ($n = 1$) may not fall within the Newton basin of attraction, so that some sort of globalization techniques must be employed. A parameter continuation strategy has thus been implemented for such strongly nonlinear cases. A sequence $\{n_d\}$ converging to the desired fluid ST index/exponent $n_D$ can be constructed to approach the basin of attraction. The solution for the case $n_{d-1}$ is used as initial guess for solving the case $n_d$. The sequence $\{n_d\}$ is obtained by subdividing the interval from 1 to $n_D$, according to the following rule:

$$n_d = \begin{cases} 1 & \text{if } d = 0; \\ n_1 \left( \frac{n_D}{n_1} \right)^{\frac{d}{D}} & \text{for } d = 1, \dots, D, \end{cases} \quad (34)$$

where $d$ indicates the parameter continuation iteration, while $n_1$ and $n_D$ are the second and final ST indices of the sequence. In general, the index $n_1 = 1$ (Newtonian case) represents a valid candidate, but for strongly nonlinear cases, indicatively when $n < 0.5$, the second ST index can be imposed to 0.5 to reduce the computational time needed to achieve overall convergence. The overall algorithm is summarized in Algorithm 2 in Appendix A.
The parameter continuation strategy estimates $\mathbf{p}^{(d)}$ following the sequence of ST indices $\{n_D\}$ and uses $\eta^{\text{lin}}$ and $\eta^{\text{nl}}$ as tolerances for the linear and nonlinear iterations, respectively. The tolerance adopted for the continuation strategy is also reduced while approaching the final ST index $n_D$, such that

$$ \eta_d^{\text{nl}} = \eta_{\text{max}}^{\text{nl}} \left( \frac{\eta_D^{\text{nl}}}{\eta_{\text{max}}^{\text{nl}}} \right)^{\frac{d}{D}} \quad \text{for } d = 1, \dots, D, \qquad (35) $$

with $\eta_{\text{max}}^{\text{nl}} = 10^{-3}$ and $\eta_D^{\text{nl}} = 10^{-8}$. This strategy allows achieving convergence in a limited number of steps even for strongly nonlinear cases (e.g., $n \approx 0.1$), covering essentially all fluids of practical interest.

# 5. Results

## 5.1. Experimental Convergence of the Proposed Method

We first examine the experimental convergence of the proposed solver in practical applications. The test considers the domain described at the beginning of Section 4.1. A 2 × 2 aperture field (mesh level 1) is generated and kept constant as the mesh is refined, resulting in an aperture field of only four different values. A sequence of 8 mesh refinements is then obtained by uniform subdivision while maintaining the original four-value aperture field for all the refined meshes originating from the same level-1 mesh. Since, as mentioned before, the scheme is consistent, the solution at mesh level 9 is considered as a “proxy” analytical solution against which we can calculate the error at different mesh levels. At each level $M$, the $L^2$ error norm is evaluated as

$$ \|e_M\| = \sqrt{\Delta x^2 \sum_j (P_{M,j} - P_{9,j})^2} \approx \sqrt{\int_{\Omega} (P_M - P_9)^2} \qquad (36) $$

Convergence of this sequence of errors together with the scheme consistency implies convergence to the true solution of the overall solver. We repeat these simulations for all fluid types listed in Table 1. The results are shown in Figure 5, which show the log-log plot of the error versus the mesh level as identified by the value of $M$.

We now investigate the nonlinear convergence of the solver. Recall that to minimize nonlinear iterations, the general strategy is to reduce the pre-asymptotic phase of the Newton method as much as possible, thus engaging efficient quadratic convergence as quickly as possible. For fluids characterized by $n$ values in the range $0.5 \le n \le 1$, Newton quadratic convergence is always observed (see fluids F1–F3 in Table 3). Conversely, in the most difficult situations, that is, for strongly ST fluids characterized approximately by the $n$-range $0.1 \le n < 0.5$, we resort to parameter continuation since Newton fails to converge when the initial guess ($n_0$) is not in the Newton basin of attraction. As depicted in Figure 6 and reported in Table 3 for the fluid F4, the convergence speed of this latter strategy controls the computational efficiency for these difficult cases. The convergence is shown in Figure 6, where the $L^2$ norms of the difference between two consecutive iterations ($s_p^{(d)}$) and of the nonlinear residual ($F^{(d)}$) are plotted against the parameter continuation step. These norms are scaled by a factor $1/2^M$ to remove the dependency on the mesh size. The results show an initial linear convergence often followed by superlinear convergence. This behavior can be intuitively expected as a consequence of the adopted scaled tolerance proposed in Equation 33. In Table 3, the performance of the numerical scheme is reported in terms of the nonlinear residual and difference between the two last iterations versus normalized computational time. Different parameter continuation strategies are adopted for fluid F4, showing higher accuracy in the solution for a larger number of continuation iterations and the possibility of reducing computational time by considering a different value of $n_1$ in rule (34) for strongly ST fluids.

## 5.2. Impact of the Rheology on the Velocity Field, Apparent Viscosity, and Fracture Transmissivity

The simulated flow field is on average cocurrent along the direction of $\overline{\nabla P}$, as an obvious consequence of the imposed boundary conditions, which forces the flow from the left-hand side to the right-hand side of $\Omega$. However, aperture field heterogeneity results in channeling of the flow in the fracture plane along paths of less resistance.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Logarithmic plot of $||e_M||$, the $L^2$ norm of the solution error calculated with respect to the finest mesh solution $P_9$ (mesh level 9), plotted as a function of the mesh level $M$ (i.e., the 2-logarithmic mesh size). The different mesh levels are uniform refinements of a $2 \times 2$ initial aperture field. Different curves are related to the different fluids listed in Table 1: F1 (yellow line), F2 (orange line), F3 (red line), and F4 (purple line). The superimposed lines show the first-order accuracy of the proposed solver.</figcaption></figure>

This phenomenon, and how it is impacted by fracture closure, has been studied for Newtonian flow for decades (Brown, 1987; Méheust & Schmittbuhl, 2001). The fracture closure $\sigma$ and the ratio $L_c/L$ both control the aperture heterogeneity. The former quantifies aperture fluctuations with respect to the mean fracture aperture, as well as the amount of contact (if the closure is sufficiently large), while the correlation length defines the size of regions of correlated large apertures (or, equivalently, of correlated small apertures), which essentially controls the spatial

<caption>**Table 3**</caption>

Performance of the Parameter Continuation Strategy for Different Numbers of Iterations and Initial Shear-Thinning Index $n_0$

<table><thead><tr><th>Simulation ID</th><th>Fluid ID</th><th>||F<sup>(D)</sup>||/2<sup>M</sup> (m<sup>2</sup>/s)</th><th>||s<sup>(D)</sup>||/2<sup>M</sup> (Pa/m)</th><th>t/t<sub>0</sub> (-)</th><th>Parameter continuation N. Iterations (initial ST index)</th></tr></thead><tbody><tr><td>Sim 1</td><td>F4</td><td>1.80 × 10<sup>-8</sup></td><td>9.10 × 10<sup>-9</sup></td><td>1021</td><td>5 (n<sub>1</sub> = 1)</td></tr><tr><td>Sim 2</td><td>F4</td><td>3.44 × 10<sup>-9</sup></td><td>2.58 × 10<sup>-9</sup></td><td>1662</td><td>10 (n<sub>1</sub> = 1)</td></tr><tr><td>Sim 3</td><td>F4</td><td>7.75 × 10<sup>-9</sup></td><td>4.08 × 10<sup>-9</sup></td><td>1095</td><td>10 (n<sub>1</sub> = 0.5)</td></tr><tr><td>Sim 4</td><td>F4</td><td>1.15 × 10<sup>-10</sup></td><td>9.56 × 10<sup>-11</sup></td><td>2041</td><td>20 (n<sub>1</sub> = 0.5)</td></tr><tr><td>Sim 5</td><td>F4</td><td>2.00 × 10<sup>-12</sup></td><td>6.11 × 10<sup>-9</sup></td><td>5960</td><td>50 (n<sub>1</sub> = 0.5)</td></tr><tr><td>Sim 6</td><td>F1</td><td>4.97 × 10<sup>-16</sup></td><td>2.82 × 10<sup>-10</sup></td><td>617</td><td>—</td></tr><tr><td>Sim 7</td><td>F2</td><td>1.03 × 10<sup>-14</sup></td><td>4.01 × 10<sup>-11</sup></td><td>678</td><td>—</td></tr><tr><td>Sim 8</td><td>F3</td><td>1.94 × 10<sup>-14</sup></td><td>1.45 × 10<sup>-9</sup></td><td>527</td><td>—</td></tr></tbody></table>

*Note.* The algorithm performance is measured in terms of residual and absolute error scaled norms at the final continuation step. The normalized computational time is obtained by dividing the computational time ($t$) by the one of the Newtonian solver ($t_0 = 0.13$ s). The aperture field is generated with $L = 0.4$ m, $\langle w \rangle = 10^{-3}$ m, $H = 0.8$, $\sigma_w/\langle w \rangle = 1$, and $L/L_c = 8$, considering a dimensionless pressure gradient of $\overline{\nabla P}/\overline{\nabla P_c} = 4.81$.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Convergence of the parameter continuation strategy toward the solution, for the most strongly shear-thinning fluid, F4 ($n = 0.1$). Panels (a and b) show the Newton error and the residual scaled norm at each continuation parameter step, respectively (see Table 3). The aperture field (mesh size $2^M = 2^{10}$) adopted for this simulation has been generated with $L_c = 0.1$ m, $\langle w \rangle = 10^{-3}$ m, $H = 0.8$, $\sigma_u/\langle w \rangle = 1$, $\overline{\nabla P}/\overline{\nabla P_c} = 4.81$, and $L/L_c = 8$. Results refer to simulations reported in Table 3: blue (Sim 1), green (Sim 2), red (Sim 3), yellow (Sim 4), and purple (Sim 5). The dotted line in panel (a) represents the tolerance $\eta_{\max}^{nl}$.</figcaption></figure>

patterns of flow channeling (Méheust & Schmittbuhl, 2003). In particular, correlated large aperture channels can form if $L \sim L_c$; if such a channel is aligned with the imposed macroscopic pressure gradient, the resulting permeability is larger than that of the smooth fracture of identical mean aperture (Méheust & Schmittbuhl, 2000). Note however that such configurations are not favored statistically (Méheust & Schmittbuhl, 2001), so, on average, heterogeneity induces higher energy losses across the fracture, which inevitably results in a lower transmissivity as compared to the parallel plate model. If the closure is sufficiently high, contact areas tend to exacerbate the behavior observed at closures that do not allow for fracture wall contact. Fractures where $L_c \sim L$ present few large contact zones, while for $L_c \ll L$, contacts are sparser across $\Omega$ and the fracture resembles a quasi-2D porous medium.

When a non-Newtonian, ST rheology is involved, all these effects are relevant, but additionally they are impacted by the fluid's rheology. Figure 7 illustrates the compound effect on flow localization of the ST fluid nature and of fracture heterogeneity. For each fluid listed in Table 1, the left-hand column shows maps of the ratio of the velocity magnitude $\bar{u} = \|\bar{u}(x)\|$ to the average velocity $\langle \bar{u} \rangle$, while the right-hand column shows the ratio of the depth-averaged apparent viscosity $\bar{\mu}(x)$ to the crossover viscosity $\mu_c = \bar{\mu}(\tau_c)$ (see definition of $\tau_c$ in Appendix B).

Figures 7a, 7c, 7e, and 7g show that the flow localization increases for decreasing values of $n$ (i.e., increasing fluid ST property). High velocities concentrate in areas of higher conductance, and low velocities in the proximity of contact zones, which typically exhibit higher resistance to flow due to aperture heterogeneity. This is similar to what is observed with Newtonian flow, but in addition, Figures 7b, 7e, 7f, and 7h show that the high velocity regions coincide with low-apparent viscosity values as expected due to the ST nature of the fluid. Hence, the flow of the ST fluid tends to be even more localized than that of the corresponding Newtonian flow. Note that the numerical solution suffers from mesh-alignment issues at small scales. This behavior develops preferentially at larger $n$ values when concentration effects dominate (Facca et al., 2021) and is expected because of the low-order reconstruction of the velocity field. However, while the local (i.e., small scale) behavior of the velocity magnitude seems to be impacted by this effect, the consistency and stability (i.e., convergence) of the FV global solution guarantees that the global behavior (i.e., at scales significantly larger than the mesh cell size) is indeed a solution of the proposed model.

Let us first discuss weakly ST cases. In these cases, it is seen that both velocity and apparent viscosity maps show a relatively slight dispersion around their ensemble average value. Small velocities are located in a limited portion of the fracture near the closed regions, where the apparent viscosity is close to $\mu_0$, implying a quasi-Newtonian behavior. On the other hand, low values of apparent viscosity are visible in small spots between contact zones, where the fluid is forced to flow under high-shear rate conditions. Strongly ST cases exhibit a different behavior with a higher dispersion around ensemble averages for both velocity and apparent viscosity. The areas with lower apparent viscosity cover a larger percentage of the fracture plane and appear more elongated with a more channel-like shape as the exponent $n$ is closer to its lowest investigated value, 0.1.

Next, we consider the fracture's global transmissivity defined as

$$T = \frac{Q\mu_0}{\overline{\nabla P}}. \qquad (37)$$

Figure 8 depicts the transition from a quasi-Darcian flow ($T \propto \overline{\nabla P}$) to a nonlinear regime ($T \propto \overline{\nabla P}^{1/n}$) for increasing dimensionless pressure gradients $\overline{\nabla P}/\overline{\nabla P_c}$. As expected, the crossover occurs at the characteristic pressure gradient defined in Equation 15. The Newtonian transmissivity $T_0$ is obtained from solving the flow for an ST index $n = 1$ and dynamic viscosity $\mu_0$. The flow regime tends to the quasi-Newtonian behavior at
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Maps of velocity magnitude $\bar{u}/\langle\bar{u}\rangle$ (left column) and depth-averaged apparent viscosity $\bar{\mu}/\mu_c$ (right column) for different fluids. Each row is associated with a fluid listed in Table 1 with the fluid rheological index *n* decreasing from top to bottom. Streamlines (white continuous lines) are shown superimposed to the velocity maps in the left column. In the viscosity maps (right column), the quasi-Newtonian areas ($\bar{\mu} \ge \mu_c$) are colored in yellow. The aperture field is generated adopting the following parameters: $M = 10$, $H = 0.8$, $L/L_c = 4$, $\langle w \rangle = 1$ mm, $\sigma_u/\langle w \rangle = 1$, $L = 0.4$ m, and $\nabla P/\nabla P_c = 10$ (fluids F1–F3) or $\nabla P/\nabla P_c = 3$ (fluid F4).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

**Figure 8.** Transition from the quasi-Newtonian to the shear-thinning behavior: ratio of the actual transmissivity $T$ to the transmissivity $T_0$ for a Newtonian fluid of viscosity $\mu_0$ as a function of the dimensionless pressure gradient $\overline{\nabla P}/\overline{\nabla P_c}$, evaluated for fluids F1–F4 in Table 1. Simulations have been conducted considering an aperture field generated with parameters: $M = 10$, $H = 0.8$, $\langle w \rangle = \sigma_w = 10^{-3}$ m, and $L/L_c = 4$. The color code for the fluids is the following: yellow line for F1, orange line for F2, red line for F3, and purple line for F4.

low-imposed macroscopic pressure gradients $\overline{\nabla P}$ and diverges from it as $T \propto \overline{\nabla P}^n$ when a pressure gradient higher than $\overline{\nabla P_c}$ is imposed.

Figure 9 illustrates how the fracture transmissivity $T$, normalized by the transmissivity $T_{pp}$ of the equivalent parallel plate fracture (i.e., of aperture equal to the rough fracture's mechanical aperture), evolves as a function of the fracture closure, again for the fluids F1–F4 listed in Table 1. The behavior of the Newtonian case corresponding to a constant $\mu_0$ viscosity, $T_0/T_{0,pp}$, is also shown for comparison. Note that, for a fracture subjected to ST fluid flow, whether the walls are rough or not, $T$ is not an intrinsic property of the fracture; it also depends on the imposed macroscopic pressure gradient $\overline{\nabla P}$ due to the nonlinearity of the hydraulic response; on the contrary, $T_0$ is independent of $\overline{\nabla P}$. Here, the results for fluids F1–F3 refer to a dimensionless pressure gradient $\overline{\nabla P}/\overline{\nabla P_c} = 10$, while for F4, it is $\overline{\nabla P}/\overline{\nabla P_c} = 3$ as setting it to 10 would have brought the Reynolds number $Re_G$ above the 1 upper limit for creeping flow. The results are shown for 200 fracture realizations generated with the same set of parameters (see caption of Figure 9). The curves go to 1 at small fractures closures as expected. The ST property of the ST fluids increases continuously from (a) to (d) and indeed, the deviation of the ST fluid's behavior from that of Newtonian flow increases all the more as the fracture closure is larger. Furthermore, this deviation is always positive and can be so large at sufficiently large closures and for a sufficiently strong ST property (see, e.g., Figures 9c and 9d), that the median behavior of a rough fracture becomes much more permeable than that of the equivalent parallel plate by a factor which reaches an order of magnitude. Indeed, as the flow becomes more localized in channels of low-apparent viscosity and high velocity, conveying most of the volumetric flow rates of the fracture, the viscous dissipation within the fluid becomes less than what it would be in the homogeneous flow of the equivalent parallel plate configuration. In other words, the ST rheology contrasts the median tendency of transmissivity to decrease due to increasing aperture heterogeneity and even reverts it in a spectacular manner.

In addition, the dispersion over the statistics of the ratio $T/T_{pp}$ also increases much more dramatically with fracture closure for ST flow than for Newtonian flow. This reflects the impact of the spatial arrangement of contact zones and of large permeability regions on the variability of the flow among the individual fracture realizations.

Note also that $T_0/T_{0,pp}$, the ratio of the parallel plate transmissivity for the ST fluid to that for Newtonian flow, is also strongly dependent on the exponent $n$ of the ST rheological law. Its values for the four fluids F1–4 are shown in Table 4. They vary from less than 3 to nearly 1000 as the ST exponent $n$ varies from 0.72 down to 0.1. Finally, note that for all the fracture flow configurations addressed here, the generalized Reynolds number $Re_G$ remains smaller than 1.

## 6. Discussion and Conclusions

This paper has presented a study of the combined effects of a nonlinear ST rheology and aperture variability on low Reynolds number flow in a single rough fracture. A Fourier transform based method has been used to generate realistic synthetic aperture fields of geological fractures. A novel flow solver has been derived to provide the first lubrication-based numerical model for the flow of an ST fluid whose rheology features a transition in the apparent viscosity between a low-shear rate Newtonian plateau and an ST behavior of index $n$ at larger shear rates. This rheological behavior was modeled with an Ellis constitutive law.

A numerical model has been set up based on a novel nonlinear generalized Reynolds equation describing the flow of an Ellis fluid in a fracture geometry under lubrication assumptions. The computational cost often associated with the solution of the nonlinear system of equations has been optimized by the developed finite volume-based solver, which yields a symmetric Jacobian that allows for the implementation of a PCG-based INK algorithm,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Dependence on the closure $\sigma_w/\langle w \rangle$ of the ratio of the actual transmissivity $T$ to the Ellis parallel-plate model transmissivity $T_{pp}$, plotted in colors for fluids F1–F3 with $\overline{\nabla P}/\overline{\nabla P}_c = 10$ and for fluid F4 with $\overline{\nabla P}/\overline{\nabla P}_c = 3$. Confidence intervals are estimated over 200 simulations. Fracture realizations are generated with: $M = 10$, $H = 0.8$, $\langle w \rangle = 10^{-3}$ m, $L/L_c = 4$ and $L = 0.4$ m. Panels (a–d) correspond to the different fluids reported in Table 1: (a) F1 in yellow, (b) F2 in orange, (c) F3 in red, and (d) F4 in purple. The data plotted in gray shows the Newtonian behavior (which is identical for the four panels). The disks linked by a continuous line represent the median behavior, while the dashed lines show the confidence interval (25th and 75th percentiles).</figcaption></figure>

to avoid oversolving. Furthermore, a continuation parameter strategy has been used to handle strongly nonlinear cases with low $n$ values. The resulting code is robust for a wide range of ST index values (from 0.1 to 1), even in strongly heterogeneous cases implying a significant percentage of closed areas (e.g., $\sigma_w/\langle w \rangle = 1$), and for pressure gradients typical of forced flow in subsurface industrial applications (e.g., $\overline{\nabla P} = 10^5$ Pa/m), much higher than typical groundwater natural gradients. Computational efficiency is maximized by ensuring that convergence is always achieved. Linear convergence is typically displayed only during the initial Newton iterations and continuation parameters steps, while quadratic asymptotic convergence is always observed. The solver is two to three orders of magnitude faster than the few alternative lubrication-based algorithms presented in the literature so far to solve the flow of ST fluids in fractures. In comparison to simulating the (Navier-)Stokes equation in the 3D space of the fracture, it is likely four to five orders of magnitude faster.

Concerning the physical effects jointly controlled by fluid rheology and heterogeneity, the ST behavior and aperture variability both contribute to flow localization. Elongated zones of high velocity and low apparent viscosity tend to form in a continuous pattern from the inlet to the outlet of the fracture for low ST index $n$ and high relative closure of the aperture field. Streamlines are affected by an increasing ST behavior of the fluid with the appearance of stronger localization patterns and zones where the flow is almost stagnant.

<table><caption>**Table 4**
Parallel-Plate Transmissivity Ratios Between the Ellis Fluids Listed in Table 1 ($T_{pp}$) and the Corresponding Newtonian Fluid ($T_{o,pp}$), of Dynamic Viscosity $\mu_\sigma$</caption><thead><tr><th>Fluid ID</th><th>F1</th><th>F2</th><th>F3</th><th>F4</th></tr></thead><tbody><tr><td><span style="text-decoration: overline;">&nabla;P</span>/<span style="text-decoration: overline;">&nabla;P<sub>c</sub></span></td><td>10</td><td>10</td><td>10</td><td>3</td></tr><tr><td><span style="text-decoration: overline;">T<sub>pp</sub></span>/<span style="text-decoration: overline;">T<sub>o,pp</sub></span></td><td>2.72</td><td>5.34</td><td>12.15</td><td>972.81</td></tr></tbody></table>

Note. A constant separation $\langle w \rangle = 10^{-3}$ m is considered between the plates.
At the fracture scale, a sufficiently high externally imposed macroscopic pressure gradient induces a non-Darcian flow regime leading to a non-Newtonian transmissivity that is higher than its Newtonian counterpart, possibly by several orders of magnitude. Such a strong nonlinear tendency is mainly due to the bulk of the flow being conveyed in marked preferential channels, typically along the path of lowest resistance, where the velocity is high and the apparent viscosity is low. The transition from a (quasi-)Darcian to a nonlinear regime occurs when the magnitude of the average shear stress is larger than a characteristic value, which is itself a function of fluid and fracture properties. Equivalently, this transition occurs close to a characteristic value of the applied macroscopic pressure gradient, which is related analytically to the aforementioned characteristic shear stress. Both quantities can be defined analytically from the Ellis constitutive law.

When multiple realizations are examined, an interesting observation becomes evident. The median fracture transmissivity decreases with increasing aperture heterogeneity (and thus closure) as expected (and well-known for Newtonian flow), but the fluid ST behavior tends to mitigate this effect as a result of the aforementioned enhanced flow channeling. Thus, an increase in non-Newtonian behavior contrasts the tendency of transmissivity to decrease due to increasing aperture heterogeneity. So much so that for sufficiently large fracture closures and ST properties, the transmissivity is seen to increase again with fracture closure, and to exceed the Newtonian transmissivity by up to one order of magnitude.

The performance achieved with this numerical scheme allows overcoming the limits of the current numerical alternatives. Future prospects include a systematic stochastic analysis of ST flow in geological fractures, considering numerous realizations for the same set of statistical geometric parameters, as well as the impact of the combined effect of fracture heterogeneity and complex rheology on anomalous transport.

## Appendix A: Numerical Algorithms

In this Appendix, the algorithms for the aperture field generator and flow solver are summarized, respectively, in the following two pseudocodes: Algorithm 1 and Algorithm 2.

<caption>**Algorithm 1.** Pseudo-code for the fracture aperture field generator</caption>

```
[k_{x1}, k_{x2}] = Grid(M);
k = \sqrt{k_{x1}^2 + k_{x2}^2};
k(k < k_c) = k_c;
z = invFFT(FFT(Rand(M)) \times k^{-(1+H)});
z = (z - Mean(z))/Std(z);
w = z \times \sigma_w + \langle w \rangle;
w(w < w_0) = w_0;
```

In Algorithm 1, the parameters are $\langle w \rangle$ and $\sigma_w$, which are, respectively, the desired distance between the walls' mean planes and the standard deviation of the aperture field prior to implementing the perfect plastic closure; $H$, which is the Hurst exponent controlling the self-affinity; and $k_c$, which is the crossover wave number related to the correlation length through $k_c = 2\pi/L_c$. The following functions are adopted:

1. Grid returns a $2^M \times 2^M$ square grid of coordinates $k_{x1}$ and $k_{x2}$;

2. Rand generates a $2^M \times 2^M$ square matrix of random numbers extracted from a uniform distribution;

3. FFT computes the 2D Fast Fourier Transform and rearranges by shifting the zero-frequency components to the center of the array;

4. invFFT computes the inverse 2D Fast Fourier Transform.

5. Mean computes the average value.

6. Std computes the standard deviation.
In the pseudo-code, the parameter $w_0$ is a nonzero lower cutoff imposed to the aperture field in fractures with closed regions. This regularization parameter is introduced to ensure well-posedness of the flow PDE and guarantee solvability in the numerical solver also in the presence of contact points (zero aperture). Typically, $w_0$ is chosen to be sufficiently small so as not to affect the accuracy of the solution. In our experiments, we employed $w_0 = 10^{-8}$ m.

In Algorithm 2, the INK algorithm and the continuation strategy can be summarized as follows.

The functions in Algorithm 2 have the following meaning:

<caption>**Algorithm 2.** Pseudo-code for the flow solver</caption>

p^{(0)} = A\f;
for d = to D do
    for k = 1 to k_{max} do
        J ← DiagonalScaling;
        J ← Reordering;
        C = Michol(J,ε,δ);
        s = PCG(J, F, C, η_k^{lin});
        if ||s|| ≤ η_d^{nl} then
            p^{(d+1)} = p^{(d)} + s;
            break

1. DiagonalScaling implements Equation 32;

2. Reordering applies the reverse Cuthill-McKee ordering;

3. Michol generates the PCG preconditioner C via incomplete Cholesky factorization;

4. PCG solves the linear system of Equation 30 via PCG, adopting C as a preconditioner and $\eta_d^{\text{lin}}$ as a tolerance stopping criterion.

## Appendix B: Estimation of the Parameters of the Ellis Model From the Carreau-Yasuda Model

This Appendix illustrates how the parameters of the fluids used in the main body of the paper were derived. In general, simple non-Newtonian rheological models, such as the Ellis constitutive law (Equation 2), are empirical nonlinear relationship between stress and strain rate that can be adopted to model ST rheology.

<caption>**Table B1**</caption>

Rheologic Parameters Related to the Four Fluids: Experimental Data Are Fitted With the Carreau-Yasuda Models

<table><thead><tr><th rowspan="2">Fluid ID</th><th rowspan="2">Solution</th><th>μ'<sub>0</sub></th><th>μ'<sub>∞</sub></th><th>n'</th><th>γ'<sub>c</sub></th><th>a'</th></tr><tr><th>(Pa · s)</th><th>(Pa · s)</th><th>(-)</th><th>(s<sup>-1</sup>)</th><th>(-)</th></tr></thead><tbody><tr><td>F1</td><td>CMC at 0.3 wt%</td><td>0.0510</td><td>0.001</td><td>0.72</td><td>17.67</td><td>0.71</td></tr><tr><td>F2</td><td>CMC at 0.5 wt%</td><td>0.2203</td><td>0.001</td><td>0.51</td><td>15.85</td><td>0.565</td></tr><tr><td>F3</td><td>CMC at 1.0 wt%</td><td>2.9899</td><td>0.001</td><td>0.40</td><td>2.74</td><td>0.668</td></tr><tr><td>F4</td><td>VES</td><td>49</td><td>0.0003</td><td>0.10</td><td>0.10</td><td>2.00</td></tr></tbody></table>

Note. Parameters for the Carboxymethylcellulose (CMC) solutions are provided in Sousa et al. (2005), while those for the viscoelastic surfactant (VES) are taken from Moukhtari and Lecampion (2018). Corresponding Ellis rheological parameters are listed in Table 1.
When two different models are fitted to experimental data, rheological parameters are estimated by best fitting and no given correspondence between the two sets of parameters exists. Here, due to the scarcity of Ellis model parameters fitted on rheological data in the literature, and in the interest of comparison and simplification, we inferred the Ellis parameters from existing data hitherto interpreted with the Carreau-Yasuda (CY) model. The latter is a five-parameter model, frequently adopted because it typically reproduces the rheology of ST fluids well; its apparent viscosity is expressed as a function of the shear rate $\dot{\gamma}$ as

$$ \mu = \mu'_{\infty} + \frac{\mu'_{0} - \mu'_{\infty}}{\left[1 + (\dot{\gamma}/\dot{\gamma}_{c})^{a'}\right]^{\frac{1-n'}{a'}}} \quad (B1) $$

where $\mu'_{0}$ and $\mu'_{\infty}$ are, respectively, the low- and high-shear rate apparent viscosity plateaus, $\dot{\gamma}_{c}$ is the characteristic shear rate separating the low-shear, pseudo-Newtonian regime from the intermediate PL regime, $a'$ is an index that influences the shape of the transition between the intermediate PL behavior and the high-shear rate plateau $\mu'_{\infty}$, and $n'$ is a ST index. Fixing the value of $a'$ to 2 results in a four-parameter model, termed the Carreau model.

In low Reynolds number applications, the high-shear rate plateau can be neglected ($\mu'_{\infty} \approx 0$) and the low-shear rate approximation of Equation B1 is obtained as

$$ \mu = \frac{\mu_0}{\left[1 + (\dot{\gamma}/\dot{\gamma}_c)^{a'}\right]^{\frac{1-n'}{a'}}} \quad (B2) $$

We consider a given CY flow curve and proceed to analytically find the Ellis model (Equation 3) that best fits the corresponding flow curve $\mu(\dot{\gamma})$. This implies that the low-shear rate viscosity plateaus be identical, and the PL regimes as well. The former constraint immediately implies that $\mu_0 = \mu'_{0}$, while the latter implies that the CY PL trend, which for $\dot{\gamma} \gg \dot{\gamma}_c$ can be simplified to

$$ \mu = \mu'_{0} \left( \frac{\dot{\gamma}_{c}}{\dot{\gamma}} \right)^{1-n'} \quad (B3) $$

be identical to that of the Ellis model, which for $\tau \gg \tau_{1/2}$ can be approximated to

$$ \mu = \mu_0 \left( \frac{\tau_{1/2}}{\tau} \right)^{\frac{1-n}{n}} \quad (B4) $$

Substituting the constitutive equation $\tau = \mu\dot{\gamma}$ in Equation B4 and recalling that $\mu_0 = \mu'_{0}$ yields

$$ \mu = \mu'_{0}^{n} \left( \frac{\tau_{1/2}}{\dot{\gamma}} \right)^{1-n} \quad (B5) $$

which can now be identified to Equation B3. We thus immediately obtain $n = n'$ and

$$ \tau_{1/2} = \mu'_{0} \dot{\gamma}_{c} \quad (B6) $$

Conversely, when determining the crossover pressure gradient $\overline{\nabla P_c}$ from Equation 15, the crossover shear stress $\tau_c$ to be considered is that corresponding to the crossover shear rate $\dot{\gamma}_c$ of the CY model. It can thus be estimated from the Ellis model parameters by solving the following implicit equation numerically:

$$ \frac{\tau_c}{\tau_{1/2}} = \left[ 1 + \left( \frac{\tau_c}{\tau_{1/2}} \right)^{\frac{1}{n} - 1} \right]^{-1} \quad (B7) $$
## Appendix C: Analytical Formulation of the Jacobian

The present Appendix provides details on the analytical evaluation of the Jacobian matrix. The relationship between the Jacobian matrix J(p) and the matrix A is given in Equation 30. In terms of components, we can write the Jacobian matrix as follows:

$$J_{ij} = A_{ij} + \sum_{u=1}^{4M} \frac{\partial A_{iu}}{\partial P_j} P_u = \begin{cases} A_{ij} + \frac{\partial A_{ij}}{\partial P_j} P_j + \frac{\partial A_{ii}}{\partial P_j} P_i & \text{if } i \in \sigma(j) \\ A_{ii} + \frac{\partial A_{ii}}{\partial P_i} P_i + \sum_{k \in \sigma(i)} \frac{\partial A_{ik}}{\partial P_k} P_k & \text{if } i = j \\ 0 & \text{otherwise.} \end{cases} \quad (C1)$$

The Jacobian has the same sparsity pattern as A, resulting in a pentadiagonal matrix. The nonzero components in Equation C1 can be derived from those of A as follows:

$$\frac{\partial A_{iu}}{\partial P_j} = \left(\frac{1}{n} - 1\right) \begin{cases} -\sum_{k \in \sigma(u)} b_k^{(u)} \left|\frac{P_k - P_u}{h}\right|^{\frac{1}{n}-1} \frac{1}{P_k - P_u} & \text{if } u = i = j; \\ +b_i^{(j)} \left|\frac{P_i - P_j}{h}\right|^{\frac{1}{n}-1} \frac{1}{P_i - P_j} & \text{if } u = j \text{ and } i \in \sigma(j); \\ -b_i^{(j)} \left|\frac{P_i - P_j}{h}\right|^{\frac{1}{n}-1} \frac{1}{P_i - P_j} & \text{if } i = u \neq j; \\ +b_u^{(j)} \left|\frac{P_u - P_j}{h}\right|^{\frac{1}{n}-1} \frac{1}{P_u - P_j} & \text{if } i = j \neq u; \\ 0 & \text{otherwise,} \end{cases} \quad (C2)$$

which defines the matrix A' mentioned in Section 4.2. Substituting the expression of the components of Equation C2 in Equation C1 leads to the following formulation for the Jacobian

$$J_{ij} = \begin{cases} -a_i^{(j)} - \frac{1}{n} b_i^{(j)} \left| \frac{P_i - P_j}{h} \right|^{\frac{1}{n}-1} & \text{if } i \neq j \\ \sum_{k \in \sigma(j)} \left( a_k^{(j)} + \frac{1}{n} b_k^{(j)} \left| \frac{P_k - P_j}{h} \right|^{\frac{1}{n}-1} \right) & \text{if } i = j \\ 0 & \text{otherwise,} \end{cases} \quad (C3)$$

where it can be noted that the Jacobian is symmetric (i.e., $J_{ij} = J_{ji}$) since $a_i^{(j)} = a_j^{(i)}$ and $b_i^{(j)} = b_j^{(i)}$.

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

Aminu, M. D., Nabavi, S. A., Rochelle, C. A., & Manovic, V. (2017). A review of developments in carbon dioxide storage. *Applied Energy*, 208, 1389–1419. https://doi.org/10.1016/j.apenergy.2017.09.015

Andreianov, B., Boyer, F., & Hubert, F. (2004). Finite volume schemes for the p-Laplacian on Cartesian meshes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 38(6), 931–959. https://doi.org/10.1051/m2an:2004045

Balhoff, M. T., & Thompson, K. E. (2006). A macroscopic model for shear-thinning flow in packed beds based on network modeling. *Chemical Engineering Science*, 61(2), 698–719. https://doi.org/10.1016/j.ces.2005.04.030

Barati, R., & Liang, J.-T. (2014). A review of fracturing fluid systems used for hydraulic fracturing of oil and gas wells. *Journal of Applied Polymer Science*, 131(16), 40735. https://doi.org/10.1002/app.40735

Barbati, A. C., Desroches, J., Robisson, A., & McKinley, G. H. (2016). Complex fluids and hydraulic fracturing. *Annual Review of Cell Biology*, 7(1), 415–453. https://doi.org/10.1146/annurev-chembioeng-080615-033630

Barnsley, M. F., Devaney, R. L., Mandelbrot, B. B., Peitgen, H.-O., Saupe, D., & Voss, R. F. (1988). Algorithms for random fractals. In H.-O. Peitgen & D. Saupe (Eds.), *The science of fractal images*. Springer. https://doi.org/10.1007/978-1-4612-3784-6

## Acknowledgments

A. Lenci acknowledges support from University of Bologna Marco Polo Grant 2019. V. Di Federico acknowledges support from University of Bologna ALMAIDEA Grant 2017. Y. Méheust acknowledges support from University of Bologna through the Visiting Professor Fellowship Excellence Department Project. Open Access Funding provided by Universita degli Studi di Bologna within the CRUI-CARE Agreement.
Benedikt, J., Girg, P., Kotrla, L., & Takac, P. (2018). Origin of the *p*-Laplacian and A. Missbach. *The Electronic Journal of Differential Equations*, **2018**(16), 1–17.

Bird, R. B., Armstrong, R. C., & Hassager, O. (1987). *Dynamics of polymeric liquids* (Vol. 1). John Wiley & Sons. Retrieved from https://www.ebook.de/de/product/3635666/bird_armstrong_hassager_liquids_2e_v1.html

Boffa, J., Allain, C., & Hulin, J. (2000). Experimental analysis of self-affine fractured rock surfaces through shadow length measurements. *Physica A*, **278**(1–2), 65–86. https://doi.org/10.1016/s0378-4371(99)00598-1

Bonnet, E., Bour, O., Odling, N. E., Davy, P., Main, I., Cowie, P., & Berkowitz, B. (2001). Scaling of fracture systems in geological media. *Reviews of Geophysics*, **39**(3), 347–383. https://doi.org/10.1029/1999rg000074

Bouchaud, E. (1997). Scaling properties of cracks. *Journal of Physics: Condensed Matter*, **9**(21), 4319–4344. https://doi.org/10.1088/0953-8984/9/21/002

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. https://doi.org/10.1029/jb092ib02p01337

Brown, S. (1995). Simple mathematical model of a rough fracture. *Journal of Geophysical Research*, **100**(B4), 5941–5952. https://doi.org/10.1029/94jb03262

Brown, S., & Scholz, C. (1985). Broad bandwidth study of the topography of natural rock surfaces. *Geophysical Research Letters*, **9**(B14), 12575–12582. https://doi.org/10.1029/jb090ib14p12575

Brush, D. J., & Thomson, N. R. (2003). Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations. *Water Resources Research*, **39**(4), 1085. https://doi.org/10.1029/2002wr001346

Cainelli, O., Bellin, A., & Putti, M. (2012). On the accuracy of classic numerical schemes for modeling flow in saturated heterogeneous formations. *Advances in Water Resources*, **47**, 43–55. https://doi.org/10.1016/j.advwatres.2012.06.016

Calder, J. (2018). The game theoretic *p*-Laplacian and semi-supervised learning with few labels. *Nonlinearity*, **32**(1), 301–330. https://doi.org/10.1088/1361-6544/aae949

Chen, Y., Elhag, A. S., Reddy, P. P., Chen, H., Cui, L., Worthen, A. J., et al. (2016). Phase behavior and interfacial properties of a switchable ethoxylated amine surfactant at high temperature and effects on CO₂-in-water foams. *Journal of Colloid and Interface Science*, **470**, 80–91. https://doi.org/10.1016/j.jcis.2016.02.028

Cipolla, C. L., Lolon, E. P., Erdle, J. C., & Rubin, B. (2010). Reservoir modeling in shale-gas reservoirs. *SPE Reservoir Evaluation and Engineering*, **13**(4), 638–653. https://doi.org/10.2118/125530-pa

Ciriello, V., Lenci, A., Longo, S., & Di Federico, V. (2021). Relaxation-induced flow in a smooth fracture for Ellis rheology. *Advances in Water Resources*, **152**, 103914. https://doi.org/10.1016/j.advwatres.2021.103914

Curtis, J. B. (2002). Fractured shale-gas systems. *AAPG Bulletin*, **86**(11), 1921–1938. https://doi.org/10.1306/61eeddbe-173e-11d7-8645000102e1865d

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), 1276. 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

Diaz, J. I., & De Thelin, F. (1994). On a nonlinear parabolic problem arising in some models related to turbulent flows. *SIAM Journal on Mathematical Analysis*, **25**(4), 1085–1111. https://doi.org/10.1137/s0036141091217731

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

Dowell, N. M., Fennell, P. S., Shah, N., & Maitland, G. C. (2017). The role of CO₂ capture and utilization in mitigating climate change. *Nature Climate Change*, **7**(4), 243–249. https://doi.org/10.1038/nclimate3231

Economides, M. (2000). *Reservoir stimulation*. Wiley.

Elsworth, D., & Goodman, R. (1986). Characterization of rock fissure hydraulic conductivity using idealized wall roughness profiles. *Int. J. Rock Mech. Min. Sci. Geomech. Abs.*, **23**(3), 233–243. https://doi.org/10.1016/0148-9062(86)90969-1

Escudier, M., Gouldson, I., Pereira, A., Pinho, F., & Poole, R. (2001). On the reproducibility of the rheology of shear-thinning liquids. *Journal of Non-Newtonian Fluid Mechanics*, **97**(2–3), 99–124. https://doi.org/10.1016/s0377-0257(00)00178-6

Facca, E., Cardin, F., & Putti, M. (2021). Branching structures emerging from a continuous optimal transport model. *Journal of Computational Physics*, **447**, 110700. https://doi.org/10.1016/j.jcp.2021.110700

Felisa, G., Lenci, A., Lauriola, I., Longo, S., & Di Federico, V. (2018). Flow of truncated power-law fluid in fracture channels of variable aperture. *Advances in Water Resources*, **122**, 317–327. https://doi.org/10.1016/j.advwatres.2018.10.024

Firdaouss, M., Guermond, J.-L., & Le Quéré, P. (1997). Nonlinear corrections to Darcy's law at low Reynolds numbers. *Journal of Fluid Mechanics*, **343**, 331–350. https://doi.org/10.1017/s0022112097005843

Gale, J. F., Laubach, S. E., Olson, J. E., Eichhuble, P., & Fall, A. (2014). Natural fractures in shale: A review and new observations. *AAPG Bulletin*, **98**(11), 2165–2216. https://doi.org/10.1306/08121413151

Glover, P. W. J., Matsuki, K., Hikima, R., & Hayashi, K. (1998). Fluid flow in synthetic rough fractures and application to the Hachimantai geothermal hot rock test site. *Journal of Geophysical Research: Solid Earth*, **103**(B5), 9621–9635. https://doi.org/10.1029/97jb01613

Jung, H. B., Carroll, K. C., Kabilan, S., Heldebrant, D. J., Hoyt, D., Zhong, L., et al. (2015). Stimuli-responsive/rheoreversible hydraulic fracturing fluids as a greener alternative to support geothermal and fossil energy production. *Green Chemistry*, **17**(5), 2799–2812. https://doi.org/10.1039/c4gc01917b

Kelley, C. T. (1987). Iterative methods for linear and nonlinear equations. Retrieved from https://www.ebook.de/de/product/5993482/c_t_kelley_iterative_methods_for

Lavrov, A. (2013a). Numerical modeling of steady-state flow of a non-Newtonian power-law fluid in a rough-walled fracture. *Computers and Geotechnics*, **50**, 101–109. https://doi.org/10.1016/j.compgeo.2013.01.004

Lavrov, A. (2013b). Redirection and channelization of power-law fluid flow in a rough-walled fracture. *Chemical Engineering Science*, **99**, 81–88. https://doi.org/10.1016/j.ces.2013.05.045

Lenci, A., & Di Federico, V. (2020). A channel model for bi-viscous fluid flow in fractures. *Transport in Porous Media*, **134**(1), 97–116. https://doi.org/10.1007/s11242-020-01438-5
Leung, D. Y., Caramanna, G., & Maroto-Valer, M. M. (2014). An overview of current status of carbon dioxide capture and storage technologies. *Renewable and Sustainable Energy Reviews*, **39**, 426–443. https://doi.org/10.1016/j.rser.2014.07.093

Li, C., Huang, Y., Sun, X., Gao, R., Zeng, F., Tontiwachwuthikul, P., & Liang, Z. (2017). Rheological properties study of foam fracturing fluid using CO₂ and surfactant. *Chemical Engineering Science*, **170**, 720–730. https://doi.org/10.1016/j.ces.2017.03.022

Liu, H.-H., Bodvarsson, G. S., Lu, S., & Molz, F. J. (2004). A corrected and generalized successive random additions algorithm for simulating fractional Levy motions. *Mathematical Geology*, **36**(3), 361–378. https://doi.org/10.1023/b:matg.0000028442.71929.26

Loisel, S. (2020). Efficient algorithms for solving the *p*-Laplacian in polynomial time. *Numerische Mathematik*, **146**(2), 369–400. https://doi.org/10.1007/s00211-020-01141-z

Lu, S.-M. (2018). A global review of enhanced geothermal system (EGS). *Renewable and Sustainable Energy Reviews*, **81**, 2902–2921. https://doi.org/10.1016/j.rser.2017.06.097

Mazzia, A., Manzini, G., & Putti, M. (2011). Bad behavior of Godunov mixed methods for strongly anisotropic advection–dispersion equations. *Journal of Computational Physics*, **230**(23), 8410–8426. https://doi.org/10.1016/j.jcp.2011.07.021

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

Morris, J. P., Chochua, G. G., & Bogdan, A. V. (2015). An efficient non-Newtonian fluid-flow simulator for variable aperture fractures. *Journal of Chemical Engineering*, **93**(11), 1902–1915. https://doi.org/10.1002/cjce.22314

Moukhtari, F.-E., & Lecampion, B. (2018). A semi-infinite hydraulic fracture driven by a shear-thinning fluid. *Journal of Fluid Mechanics*, **838**, 573–605. https://doi.org/10.1017/jfm.2017.900

Mourzenko, V. V., Thovert, J. F., & Adler, P. (1995). Permeability of a single fracture; validity of the Reynolds equation. *Journal de Physique II*, **5**(3), 465–482. https://doi.org/10.1051/jp2:1995133

Muggeridge, A., Cockin, A., Webb, K., Frampton, H., Collins, I., Moulds, T., & Salino, P. (2014). Recovery rates, enhanced oil recovery and technological limits. *Philos. Trans. R. Soc. A: Mathematical, Physical and Engineering Sciences*, **372**(2006), 20120320. https://doi.org/10.1098/rsta.2012.0320

Mulligan, C., Yong, R., & Gibbs, B. (2001). Surfactant-enhanced remediation of contaminated soil: A review. *Engineering geology*, **60**(1–4), 371–380. https://doi.org/10.1016/s0013-7952(00)00117-4

Perrin, C. L., Tardy, P. M., Sorbie, K. S., & Crawshaw, J. C. (2006). Experimental and modeling study of Newtonian and non-Newtonian fluid flow in pore network micromodels. *Journal of Colloid and Interface Science*, **295**(2), 542–550. https://doi.org/10.1016/j.jcis.2005.09.012

Philip, J. R. (1961). N-diffusion. *Australian Journal of Physics*, **14**(1), 1–13. https://doi.org/10.1071/ph610001

Plouraboué, F., Kurowski, P., Hulin, J.-P., Roux, S., & Schmittbuhl, J. (1995). Aperture of rough cracks. *Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics*, **51**(3), 1675–1685. https://doi.org/10.1103/physreve.51.1675

Sanner, B., Karytsas, C., Mendrinos, D., & Rybach, L. (2003). Current status of ground source heat pumps and underground thermal energy storage in Europe. *Geothermics*, **32**(4–6), 579–588. https://doi.org/10.1016/s0375-6505(03)00060-9

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

Schmittbuhl, J., Vilotte, J.-P., & Roux, S. (1995). Reliability of self-affine measurements. *Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics*, **51**(1), 131–147. https://doi.org/10.1103/physreve.51.131

Sheng, G., Su, Y., & Wang, W. (2019). A new fractal approach for describing induced-fracture porosity/permeability/compressibility in stimulated unconventional reservoirs. *Journal of Petroleum Science and Engineering*, **179**, 855–866. https://doi.org/10.1016/j.petrol.2019.04.104

Singer, A. (2006). From graph to manifold Laplacian: The convergence rate. *Applied and Computational Harmonic Analysis*, **21**(1), 128–134. https://doi.org/10.1016/j.acha.2006.03.004

Sochi, T., & Blunt, M. J. (2008). Pore-scale network modeling of Ellis and Herschel-Bulkley fluids. *Journal of Petroleum Science and Engineering*, **60**(2), 105–124. https://doi.org/10.1016/j petrol.2007.05.009

Sousa, R., Riethmuller, M., Pinto, A., & Campos, J. (2005). Flow around individual Taylor bubbles rising in stagnant CMC solutions: PIV measurements. *Chemical Engineering Science*, **60**(7), 1859–1873. https://doi.org/10.1016/j.ces.2004.11.035

Suleimanov, B., Ismailov, F., & Veliyev, E. (2011). Nanofluid for enhanced oil recovery. *Journal of Petroleum Science and Engineering*, **78**(2), 431–437. https://doi.org/10.1016/j.petrol.2011.06.014

Sullivan, S., Gladden, L., & Johns, M. (2006). Simulation of power-law fluid flow through porous media using lattice Boltzmann techniques. *Journal of Non-Newtonian Fluid Mechanics*, **133**(2–3), 91–98. https://doi.org/10.1016/j.jnnfm.2005.11.003

Wang, M., Chen, Y.-F., Ma, G.-W., Zhou, J.-Q., & Zhou, C.-B. (2016). Influence of surface roughness on nonlinear flow behaviors in 3d self-affine rough fractures: Lattice simulation of power-law fluid flow through porous media using Lattice Boltzmann simulations. *Advances in Water Resources*, **96**, 373–388. https://doi.org/10.1016/j.advwatres.2016.08.006

Yasuda, K., Armstrong, R. C., & Cohen, R. E. (1981). Shear flow properties of concentrated solutions of linear and star branched polystyrenes. *Rheologica Acta*, **20**(2), 163–178. https://doi.org/10.1007/bf01513059

Zami-Pierre, F., de Loubens, R., Quintard, M., & Davit, Y. (2016). Transition in the flow of power-law fluids through isotropic porous media. *Physical Review Letters*, **117**(7), 074502. https://doi.org/10.1103/physrevlett.117.074502

Zhang, M., Prodanović, M., Mirabolghasemi, M., & Zhao, J. (2019). 3D microscale flow simulation of shear-thinning fluids in a rough fracture. *Transport in Porous Media*, **128**(1), 243–269. https://doi.org/10.1007/s11242-019-01243-9

Zimmerman, R., & Bodvarsson, G. (1996). Hydraulic conductivity of rock fractures. *Transp. Porous Media*, **23**(1), 1–30. https://doi.org/10.1007/bf00145263