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

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

particular flow localization in areas of higher conductance, where the fracture aperture is larger.The resulting hydraulic behavior of the (rough) geological fracture of interest is usually compared to the parallel plate model where both walls are assumed to be perfectly smooth planes, parallel to each other and separated by a distance that is identical to the mechanical aperture (distance between the mean planes of the two rough walls) of the rough fracture.Flow channeling induces either a flow-enhancing or flow-inhibiting hydraulic behavior with respect to the corresponding parallel plate model (Méheust & Schmittbuhl, 2000).However, over a large population of fractures with identical statistical properties, flow-inhibition dominates (Méheust & Schmittbuhl, 2001), but the statistical dispersion (i.e., variability) over the fractures' hydraulic behavior depends on their correlation length (Méheust & Schmittbuhl, 2003).
This study focuses both on Newtonian and shear-thinning (ST) fluids.Subsurface industrial activities widely adopt complex fluids, such as foams (Faroughi et al., 2018;Hosseini et al., 2019), clayey muds, polymer solutions, and aqueous and non-aqueous suspensions, which are ST, that is, whose viscosity tends to decrease under increasing shear rates.CO 2 -based ST fluids (Jiang, 2021) have also been recently proposed as a water-saving solution to reduce the use of freshwater in subsurface activities.ST fluids have been demonstrated to enhance the cleanup efficiency (Hosseini et al., 2019) and carrying capability of remedial amendments (Tosco et al., 2012;Zhong et al., 2008).
The effects of ST rheology on the transmissivity of variable aperture fractures were first investigated in the past using simplified deterministic (Di Federico, 1997, 2001) or stochastic (Di Federico, 1998) conceptualizations.More recently, the fracture's hydraulic behavior has been shown to be augmented by the interplay between pore-scale heterogeneity and ST rheology.This potentially leads to fracture transmissivities orders of magnitude larger than the Newtonian counterpart (Lenci, Méheust, et al., 2022).The transmissivity is computed from the ratio of the volumetric flow rate to the applied macroscopic pressure gradient.In numerical simulations, the flow rate is obtained from simulating the fluid velocity field.The Reynolds equation, based on the lubrication approximation (Brown, 1987), is often used to model Newtonian flow in geological fractures.This depth-averaged approach relies on the local validity of the Hagen-Poiseuille law, which holds as long as the conditions for Stokes flow are satisfied and spatial variations of the aperture fields are slow.
Recently, an efficient algorithm has been proposed by Lenci, Putti, et al. (2022) to solve a generalized non-linear Reynolds equation that extends the model to ST fluids whose rheology includes a transition from a Newtonian plateau at low shear rates to a power law ST behavior at large shear rates.This numerical strategy has been adopted to produce a Monte Carlo (MC) study of the ST flow through rough fractures (Lenci, Méheust, et al., 2022), overcoming the computation burden that is typical of stochastic analysis of non-linear problems in spatially varying domains and thus allowing to evidence the impact of the fluid's rheology as discussed above.However, exploring the space of controlling parameters for ST flow in geological fractures over a wide range of values of these parameters, using fracture-scale simulations, would inevitably require a vast number of MCs, and thus a tremendous amount of computational time, due to the dimensionality of the parameter space (consisting of at least three dimensions: fracture closure, fracture correlation length, imposed pressure gradient), the nonlinearity of the model, and the variability of the aperture field's topology over different realizations.
In this study, reduced-order models (ROMs) are thus employed to drastically decrease the computational cost of this statistical framework and to obtain the dependence of given quantities of interest (QoIs) describing a fracture's hydraulic behavior on the controlling parameters in a simple mathematical form, rather than using a pure MC approach as previously done by Lenci, Méheust, et al. (2022).The QoIs, whose response surfaces are reconstructed in this study, are the ensemble average and difference between the third and first quartiles (computed over the number of MCs of the aperture field) of a set of MFMs, which are related to the effective transmissivity and flow correlation length.
The polynomial chaos expansion (PCE) technique is used for model reduction since it has demonstrated its effectiveness in different hydrogeological settings.For example, PCE has been used to accelerate variance-based global sensitivity analysis and uncertainty quantification (otherwise prohibitive from a computational point of view) in transport problems in subsurface hydrology (Ciriello & Barros, 2020;Ciriello et al., 2012Ciriello et al., , 2017;;Kim et al., 2022;Oladyshkin et al., 2012) and for a sensitivity-based transport model calibration (Ciriello et al., 2013(Ciriello et al., , 2015)).PCE has also been used to compute distribution-based global sensitivity analysis (Ciriello et al., 2019), and to predict the future behavior of water resources under climate change (Focaccia et al., 2021).Merchán-Rivera et al. (2021) used PCE to quantify the uncertainty in surface water-groundwater interaction estimates.The key advantage of 10.1029/2023WR036026 3 of 14 using PCE is that this technique provides an approximation of the response surface of the QoIs in the form of a polynomial series in the space of variability of the selected governing parameters (Ghanem & Spanos, 1991;Wiener, 1938).To be determined, the PCE requires a negligible number of model runs, thus being advantageous over onerous MCs for mapping out the governing parameter space (Ciriello & Barros, 2020).
Hence, we greatly extend our previous findings on the hydraulic behavior of Newtonian and ST fracture flow, by obtaining a prediction of the generic hydraulic behavior of a geological fracture, under Newtonian and ST flow, which is virtually continuous over the parameter space of practical interest.In addition, we characterize the relationship between the structural correlation length and that of the velocity field; to our knowledge, this has never been investigated for fracture flow, let alone for ST fracture flow.Furthermore, our approach provides both the average behavior over the statistics and the statistical dispersion as a function of the parameters, both for the hydraulic behavior and for the flow correlation length.These findings are only achievable through model reduction, which, to our knowledge, has never been used to address fracture-scale hydraulics.
The paper is organized as follows.In Section 2, the geometry of geological fractures is described and the flow equations, numerical model, and stochastic framework are described; the principle of the ROM is also explained.In Section 3, the results of the analysis are presented and interpreted.We draw conclusions in Section 4.

Problem Statement
Synthetic geological fractures: Geological planar fractures consist of two facing rough walls which are self-affine (Power et al., 1987) at all scales (Schmittbuhl et al., 1995).Assuming statistical isotropy of the topographies, this scale invariance manifests in the Fourier transform of the topographies of the upper and the lower wall, z u and z l respectively, which scale as (Méheust & Schmittbuhl, 2001) where is the wave number and H is the Hurst exponent that controls the scale invariance.The measured values of H fall in the range [0.5; 0.8] (Boffa et al., 1999;Bouchaud et al., 1990), so that the topographies are planar at large scales, the corresponding asymptotic planes being also the mean planes of the topographies.
Assuming that these two mean planes are parallel to each other, the local fracture aperture is the distance between the two fracture walls, measured along the direction normal to the mean planes.Further assuming that  u|l have been defined with a zero mean, the aperture field is defined as where the mechanical aperture a m (Brown, 1987) is the distance between the mean planes of the two walls and is equal to the mean aperture 〈a〉 if the walls are not in contact.When a m is set sufficiently low, contact zones result from eventual interpenetrations between the rough walls; in such cases the negative values of a from Equation 2are set to zero and the mean aperture 〈a〉 becomes larger than a m .To simplify notations, we shall, in the following, use the notation 〈a〉 instead of a m ; 〈a〉 is then to be understood as the distance between the mean planes of the fracture walls, and is only strictly equal to the aperture field's mean value when the walls are not in contact.
Fracture walls are generally not entirely uncorrelated with respect to each other: their Fourier modes are identical above a given characteristic length scale which we denote correlation length L c and assume to be a constitutive parameter of the rock formation, resulting from its history of tectonic activity and weathering.Consequently, the Fourier transform of the aperture field a only possesses the same self-affinity as the walls, expressed by Equation 1, at scales smaller than L c .Depending on the fracture length L, this will include the entire relevant scale range (if L ≤ L c ), or not.
A geological fracture's geometry can thus be described from the following statistical parameters: its length L, the correlation length L c , its mechanical aperture a m , its standard deviation σ a , the Hurst exponent H, and a parameter ω controlling the stochastic nature of the aperture field (e.g., the seed of the random number generator used to generate it numerically).An example of synthetic fracture is represented in Figure 1.The synthetic aperture wall topographies have been generated following the pseudocode described in Lenci, Putti, et al. (2022), which implements the above theoretical geometrical properties.The same pseudocode has been used to generate all the synthetic aperture fields which are considered in this study.

Fluid rheology:
The ST fluid rheology is modeled with the three-parameter Ellis constitutive law, which relates the fluid's dynamics viscosity to the shear stress τ; μ 0 defines the low-shear rate quasi-Newtonian viscosity plateau, n is the ST exponent that defines the power-law trend at large shear rate, and the characteristic shear stress τ 1/2 shapes the transition between those two asymptotic behaviors.These rheological parameters can be derived experimentally from rheometric measurements by least-squares regression, or inferred from other constitutive laws as proposed for the Carreau-Yasuda model by Lenci, Putti, et al. (2022).

Lubrication flow model:
The steady-state isothermal viscous flow of an incompressible Ellis fluid through a geological fracture is modeled using a lubrication-based approach.The local pressure distribution P across the mean fracture plane is governed by the generalized Reynolds equation (Lenci, Putti, et al., 2022): which holds for Stokes flow conditions (Re ≪ 1) as long as the aperture field's gradient remains small (‖∇a‖ ≪ 1).For Newtonian fluids, the classical Reynolds equation can be recovered for n = 1 or τ 1/2 → + ∞.Solving Equation 4 yields the spatial distribution of the pressure (P), and the velocity field is then obtained as Transmissivity: The volumetric flow rate Q can then be computed from the velocity field to estimate the apparent transmissivity, Even though for Newtonian flow the transmissivity is an intrinsic property of the fracture's geometry, for ST fluids it depends on both that geometry, the magnitude of the external pressure gradient, and the rheological parameters n and τ 1/2 .
Hereinafter, the subscript "0" is adopted to denote quantities related to the Newtonian model, while the subscript ‖ stands for quantities estimated for the parallel plate geometry, that is, assuming a uniform aperture field a = 〈a〉.
In this work two different dimensionless formulations of the transmissivity are thus used: T/T 0 to study the deviation, due to the ST rheology, from the hydraulic behavior of the fracture under Newtonian flow, and T/T ‖ to analyze the compound effect of fracture heterogeneity and fluid rheology.
In the literature, the transition from the quasi-Newtonian (T/T 0 ≈ 1) to the ST hydraulic behavior (T/T 0 > 1) has been observed to occur when the external pressure gradient reaches the crossover value where τ c is a crossover shear stress obtained from ]−1 (Lenci, Putti, et al., 2022).

Correlation length of the velocity field:
The finite volume scheme (Lenci, Putti, et al., 2022) provides edge velocities, which can be used to define the cell-centered flow velocity by bilinear interpolation scheme, as proposed by Pollock (1988).Considering a regular partition of the domain, the longitudinal and transverse components of the flow velocity located at the center of each mesh cell are obtained by arithmetic averaging of the edge velocities at opposite cell faces.The correlation length of the velocity magnitude ℓ is then analyzed in relation to the correlation length of the aperture field λ.Both correlation lengths are defined as the first zero-crossing lag of the autocorrelation function (ACF), estimated by averaging the ACFs of all longitudinal profiles (respectively of the aperture and velocity fields) of the fracture, and discarding the contact zones.Figure 2 provides examples of such average ACFs of the velocity magnitude for the flows of a Newtonian fluid and a ST fluid in two fracture realizations, one with L/L c = 1 and one with L/L c = 16.
Note that λ is related to L c , which is the upper limit scale for the self-affinity of the aperture field and is imposed when generating the geometry, but the former is proportional to the latter.Indeed, when a random fractal is truncated beyond a given length scale, say L c , the correlation length of the resulting field is proportional to L c , as demonstrated numerically by Cintoli et al. (2005).The proportionality factor between λ and L c depends on L c but not on the fracture closure.More details are provided in Appendix A.

Numerical Method
Stochastic apertures fields defined by the set of parameters (L, L c , a m , σ a , H, ω) are generated using a random generator and shaping the Fourier transform into the product of a power law cone (see Equation 1) and of the Fourier transform of a random noise (see Lenci, Putti, et al., 2022;Méheust & Schmittbuhl, 2003).
We adopt a depth-averaged approach to solve Equation 4 over a 2D square domain of the mean fracture plane.Permanent flow occurs under the effect of a longitudinal pressure drop ΔP imposed through Dirichlet boundary conditions on the eastern and western boundaries and "no-flow-through-the-boundary" conditions imposed on the northern and southern boundaries.The latter are Neumann boundary conditions that stipulate a zero pressure gradient across the boundary; they are well suited to simulating fractures with a finite lateral extent and are particularly relevant to simulating geological fracture flow.The other type of boundary conditions that could be applied on the lateral boundary are periodic boundary conditions, which are much less relevant to real fractures.
Equation 4 is solved on a 2 10 × 2 10 grid with the finite volume scheme proposed by Lenci, Putti, et al. (2022), whose implementation relies on an inexact preconditioned conjugate gradient (PCG) Newton-Krylov method combined with a parameter continuation strategy (Lenci, Putti, et al., 2022).The mesh size is selected based on previous studies.In particular, the nonlinearity of the constitutive law results in a higher computational expense as compared to Newtonian (and hence, Darcian) flow, so we choose a convenient value for the grid size (2 10 = 1,024) to keep an adequate level of detail while allowing to perform the required number of simulations in a reasonable time.Note that larger modes of the fracture aperture field dominate smaller modes in controlling the flow heterogeneity; as such, further grid refinement is unnecessary because it wouldn't impact the predicted hydraulic behavior.
Here, we have replaced the PCG method with the aggregation-based algebraic multigrid method (Napov & Notay, 2012;Notay, 2012), which performs better in achieving global convergence.As already assessed by Lenci, Méheust, et al. (2022), this highly efficient scheme allows a stochastic analysis of the phenomenon via MCs, resulting in a reasonable amount of computational time.

Stochastic Framework
The MFMs which we evaluate with the model are: (a) the fracture transmissivity (T) and (b) the flow correlation length (ℓ).In particular, the former is expressed by means of the ratios T 0 /T 0‖ and T/T ‖ ; these are used to measure the effect of the aperture field's heterogeneity on the transmissivity with respect to Newtonian and ST fluids, respectively.The parallel plate model, characterized by parallel planar walls and thus a uniform aperture a = 〈a〉, is taken as the homogeneous reference case for the transmissivity.As for the flow's characteristic length scale, it is relevant for upscaling transport in heterogenous media (Comolli et al., 2019), where the medium's correlation length is generally used as a surrogate to relate macroscopic solute transport to flow conditions (Dentz et al., 2016).
In particular, the flow's correlation length has been adopted as the only parameter present in Markov models for predicting the evolution of the spatially sampled Lagrangian velocities in Bernoulli and Ornstein-Uhlenbeck processes (Dentz et al., 2016;Kang et al., 2017;Morales et al., 2017).However, the incomplete characterization of the flow's correlation length in both porous and fractured media, caused it to be replaced with the medium's correlation length, which is supposed to be of the same order of magnitude.Here, we focus on the ratios between the flow's and medium's correlation lengths for Newtonian and ST fluids, ℓ 0 /λ and ℓ/λ respectively, to quantify how they differ under different geometrical, rheological, and external forcing conditions.Indeed, if for heterogeneous porous media the two characteristic lengths diverge significantly (Cvetkovic et al., 1996), little has been investigated for fractured media.Finally, to explore the transition from the quasi-Darcian to the non-linear ST regime over the entire parameter space, the ratio between the transmissivities of the two fluids, T/T 0 , is also analyzed for different flow conditions and degrees of heterogeneity.This allows us to extend the results obtained by Lenci, Méheust, et al. (2022), which were constrained to a limited set of scenarios.
Once defined the set of MFMs in this study, that is,  {0∕0‖,  ∕‖, 0∕, ∕,  ∕0} , the QoIs are key statistics of these selected metrics, computed over 100 MCs of the aperture field.We focus, in particular, on the ensemble average (hereinafter this will be indicated directly with the name (⋅) of each MFM for the sake of brevity) and the difference between the third and first quartiles (hereinafter Δq (⋅) ).Among these QoIs, T 0 /T 0‖ and ℓ 0 /λ (and the corresponding quartiles) are functions only of the set of parameters dealing with the heterogeneity: While we span continuously the QoIs over a proper range of variation of the remaining parameters by modeling them as uniformly distributed random variables.Specifically, we consider: • The fracture closure  ∕⟨⟩ ∈ [0; 1] , with 〈a〉 = 0.001 m, • The dimensionless pressure gradient  ∇ ∕∇c ∼∈ [0.1; 10] to consider an interval ranging from naturally occurring groundwater head gradients to industrial operating conditions (Jung, 1989).
In order to be able to explore such a vast parameter space, instead of relying solely on the MC simulations of the numerical flow model, we apply model reduction as reported in the next section.

Reduced-Order Models
In order to reconstruct the dependence of the QoIs on the selected parameters, that is, the topography of these QoIs in the selected random parameter space, continuously and at a negligible computational cost, ROMs are used according to the PCE theory.
Let Ω denote the generic QoI of this study.Let f be the high-fidelity model (HFM) described in Sections 2.1 and 2.2, which provides Ω as a function of the aforementioned set of parameters (see Section 2.3), M of them being random and collected in vector p: Ω = f(p), with dim(p) = M.If the variance of Ω is finite, the approximation of the HFM provided by the PCE is admitted (Ghanem & Spanos, 1991;Wiener, 1938;Xiu & Karniadakis, 2002) as In Equation 8, multi-indices   = {1, . . .,  } ∈ ℕ  are associated with multivariate polynomials Ψ a of degree  || = ∑  =1  , which constitute an orthonormal basis with respect to the joint probability density function (PDF) of p (Xiu & Karniadakis, 2002); while the PCE coefficients s a are deterministic coordinates of the spectral decomposition (Ghanem & Spanos, 1991).In practice equation 8 is truncated as where q is the maximum degree of the expansion, that is,  || ≤ q for all   ∈ ℕ  .The PCE coefficients are computed through a regression-based approach consisting in the minimization of the variance of the residual,  = | Ω − Ω| , with respect to the s i ; the set of regression points in the random parameter space is provided by the probabilistic collocation method (PCM) (Webster et al., 1996).
In our study, with M = 2 and q = 4, thus resulting in a number of HFM runs equal to P = 15 according to the PCM (in place of a full MC series in the parameter space).Since the parameters in p are distributed uniformly, the Legendre polynomials are used as the basis in the ROMs given by Equation 9 (Xiu & Karniadakis, 2002).

Newtonian Flow
We first focus the analysis on the QoIs for Newtonian fluids, T 0 /T 0‖ and ℓ 0 /λ, depicted in the left and right panels of Figure 3, respectively.Figure 3 shows, for different correlation lengths  (∕c ∈ {1; 4; 16}) distinguishable by color variations, the MC ensembles provided by the PCM and adopted to generate the predictions of QoIs (continuous lines) as a function of σ a /〈a〉 in order to calibrate the PCE.In particular, each MC ensemble is provided with the first and third quartile (bars), 100 single realizations (dots with colored edges), and the relative average (circle symbols with black edges).In Figure 3a, single realizations exhibit two opposite hydraulic behaviors referring to the parallel plate model transmissivity: flow-inhibiting if T 0 /T 0‖ < 1, or flow-enhancing if T 0 /T 0‖ > 1.When the fluid rheology is Newtonian, the transmissivity is an implicit feature of the fracture geometry, hence these behaviors are solely affected by the morphology of the aperture field.The flow-inhibiting behavior is dominant in the statistics and promoted by higher fracture closure, as it induces major viscous energy losses due to roughness.Obstructions, such as contact areas, hinder fluid motion and increase flow path tortuosity.On the other hand, the presence of areas with large aperture aligned parallel to the external pressure gradient is associated with a flow-enhancing behavior, although this geometrical feature is not favored statistically (Lenci, Putti, et al., 2022;Méheust & Schmittbuhl, 2001).The effect of the fracture closure on the transmissivity is similar among different values of L/L c .However, small fractures (i.e., L ≈ L c ), present a higher statistical dispersion around the ensemble average with several outliers displaying prominent flow-enhancing behavior.
Figure 3b shows the same plots as Figure 3a but for the ratio ℓ 0 /λ of the flow correlation length to the geometrical correlation length.To our knowledge, ℓ (including ℓ 0 for Newtonian flow) has hardly been investigated in the literature.Neither the ensemble average nor the width of the confidence interval are observed to be significantly sensitive to fracture closure, hence L/L c is the most significant parameter that substantially affects ℓ 0 /λ.In particular, the higher the ratio L/L c , the higher the ensemble average and the statistical dispersion.

Shear-Thinning Flow
Comparison to Newtonian transmissivity: We now extend the analysis to complex fluids that may develop a ST behavior.Figure 4 shows the dependency on fracture closure and dimensionless pressure gradient of the PCE-based predictions of QoIs: the ensemble average of the ratio of the transmissivity to the corresponding Newtonian transmissivity, T/T 0 , and the corresponding amplitude of the confidence interval,  Δ   0 . The former is adopted to quantify the departure of the ST hydraulic behavior from the Newtonian transmissivity, while Under low external pressure gradients, the Ellis fluid exhibits a quasi-Darcian behavior , as the transmissivity is comparable with the Newtonian transmissivity (T ≈ T 0 ).This results from low shear rates across the flow domain remaining under the crossover value below which the Ellis fluid does not develop a relevant ST behavior (μ ≈ μ 0 , low-shear rate plateau).Consequently, the hydraulic behavior is controlled solely by medium heterogeneity, as for Newtonian fluids (Figure 3).Increasing external pressure gradients triggers a transition to a non-linear regime, which is present in all panels of Figure 4 and exhibits a power law trend of slope  1  − 1 when the dimensionless pressure gradient is sufficiently higher than 1.The power law trend of T/T 0 is weakly affected by the fracture closure and the correlation length, although heterogeneity contributes to the regulation of the transition between the two regimes.In particular, the establishment of the non-linear regime due to the fluid's ST behavior, as well as the statistical dispersion around the ensemble average, are promoted by higher fracture closures and lower values of L/L c , while  ∇ ∕∇c is the control parameter of the transition.Indeed, for a given  ∇ ∕∇c , higher closures mean a higher degree of flow channeling, which results in a wider range of hydraulic behaviors (from flow-inhibiting to flow-enhancing). .The figure's layout is otherwise identical to that of Figure 4. Again, for low dimensionless pressure gradients the hydraulic behavior tends to the Newtonian transmissivity (see Figure 3a), as T → T 0 and T ‖ → T 0‖ .For sufficiently high dimensionless pressure gradient, indicatively  ∇ ∕∇c > 3 (a transition roughly occurs between 0.1quasi-Newtonian behavior-and 4), almost horizontal contour lines are seen in the T/T ‖ topography, indicating that the fracture closure's impact on the transmissivity scales equivalently in the rough-walled fracture and the corresponding parallel plate geometry.This is all the truer as the correlation length L c is larger.Moreover, as depicted in the maps on the right-hand side, the fracture closure does not significantly impact the statistical dispersion of T/T ‖ when the flow is strongly ST.Furthermore, when the ST hydraulic behavior is developed, the observed relation between the transmissivity and the fracture closure is non-monotonic, for any value of L/L c .The ST nature of the fluids mitigates the viscous energy losses induced by the higher heterogeneity due to fracture closure, increasing the statistical likelihood of flow-enhancing hydraulic behavior.As the aperture field becomes more heterogeneous, both the transmissivity and statistical dispersion increase, with many realizations characterized by a transmissivity that is orders of magnitude higher than that of the corresponding parallel plate case.This can be explained by the formation of channels where the fluid shear-thins due to high shear rates, so that the resulting apparent viscosity decreases and velocities increase, resulting in preferential flow paths of very high conductance due to minor viscous energy losses.The higher statistical dispersion observed when fracture closure increases and L/L c = 1 depends on the more chaotic nature of the medium, which then affects the flow by generating larger channels.Behavior of the correlation lengths' ratio: We discuss in the present subsection how the ratio of the flow's longitudinal correlation length ℓ to the aperture field's longitudinal correlation length λ depends on the problem's parameters.Note that since λ ∝ L c , as discussed earlier (see also the numerical simulations by Cintoli et al. (2005)), the ratios ℓ/λ and ℓ/L c differ by a factor.In particular, Figure 6a shows the influence of the fracture closure and dimensionless pressure gradient on ℓ/λ.On the right-hand side, Figures 6b  and 6c display several realizations of the fracture aperture field and the corresponding flow field, for different correlation ratios  (∕c ∈ {1; 4; 16}) and fracture closures: σ a /〈a〉 = 0.8 (Figure 6b), and σ a /〈a〉 = 0.2 (Figure 6c).

Comparison to parallel plate transmissivity:
Under a sufficiently small external pressure gradient  ( ∇ ∕∇ c → 0 ) the flow correlation tends to the Newtonian behavior depicted in Figure 3. Unlike the results obtained for transmissivity, the flow correlation length in small fractures (i.e., L ≈ L c ) shows a minor statistical dispersion around the ensemble average, which results in an almost uniform value (ℓ/λ ≈ 1) over the parameter space.In cases where  ∕c ∈ {4; 16} , the ST effect increases for increasing imposed pressure gradients.Higher dispersion characterizes media of lower fracture closure, where the flow correlation length is larger than the medium's correlation length: ℓ/λ ≈ λ ranges from approximately 3 for σ a /〈a〉 → 1 up to 6 times σ a /〈a〉 → 0, in L/L c = 16.Indeed, when the fracture closure is high, the presence of a contact zone strongly constrains the flow path, reducing the variability of ℓ/λ, and consequently the impact of the fluid's rheology on ℓ/λ.On the other hand, for a low fracture closure, the flow correlation length exhibits a marked effect of the fluid rheology on the shapes of the streamlines, exhibiting a stronger ST effect under an increasing dimensionless pressure gradient.

Conclusions
This work has characterized quantitatively the impact on fracture scale metrics (i.e., generalized transmissivity and flow correlation length) of (a) the aperture heterogeneity resulting from wall roughness and (b) the fluid's ST rheology.To this aim, we used a 2-D numerical model relying on a lubrication-based theoretical description of ST flow, under a stochastic framework involving ROMs.The ST constitutive law considered, the so-called Ellis model, corresponds to the complex rheology of polymer fluids and most ST fluids used in subsurface industrial applications.It exhibits a quasi-Newtonian behavior at low shear rates and transitions to a ST behavior under increasing shear rates.The dependence of the fractures scale metrics on the problem's parameters, that is, the topographies of these metrics over the parameter space, are obtained through the PCE as a ROM, relying on an optimized set of MC simulations in the parameter space.This approach allows investigating a three-dimensional parameter space over a wide range of parameter values (imposed pressure gradient over two orders of magnitude, fracture closure between 0 and 1, ratio of the fracture length to the aperture fields' correlation length between 1 and 16) within a fraction (typically 3 orders of magnitude less) of the time that would be required without model reduction.
Results for Newtonian flow (such as that of water) confirm the previous conclusions, that is, that the fracture transmissivity is on average diminished by increasing the closure due to higher viscous energy losses, but with a dispersion over the statistics that is all the larger as the closure is larger; the correlation length hardly impacts the mean behavior, but the dispersion around that mean value is larger for a larger correlation length.The flow correlation length is weakly influenced by the fracture closure, whereas the medium (i.e., aperture field's) correlation length is more markedly impacted by it.Indeed, for small fractures, when L/L c ≈ 1, the medium and flow corre lation lengths are similar.Conversely, the flow correlation length grows faster than the aperture field's correlation length when larger fractures are considered (L > L c ).When considering the ST rheology exhibited (among other complex fluids) by polymer solutions, the generalized fracture transmissivity exhibits a transition from the Darcian behavior, typical of the linear Newtonian model, to a non-linear power-law regime, when increasing the imposed pressure gradient.Fracture closure and medium's correlation length affect neither the Newtonian nor the power-law trend, but they influence how fast the transition between these two regimes occurs: a higher fracture closure induces a faster transition between the two regimes (i.e., transition to the non-linear regime occurs under higher dimensionless pressure gradients).Besides, a higher dispersion around the mean behavior is observed in small fractures, that is, for a correlation length that is larger as compared to the fracture length (L/L c ≈ 1).
Unlike the purely Newtonian rheology, the ST rheology of polymer solutions mitigates the transmissivity reduction observed on average with respect to the corresponding parallel plate (i.e., smooth) geometry; consequently, the average generalized transmissivity for ST flow can be orders of magnitude higher in the rough geological fracture than in the corresponding parallel plate.The Newtonian model thus overestimates fracture transmissivity when used to predict the flow of fluids that exhibit an ST behavior, since the ST rheology is more efficient at minimizing viscous energy losses induced by the aperture variability.
The ratio of the flow correlation length to the medium (i.e., aperture field's) correlation length was also analyzed for the first time in non-Newtonian fracture flow.Fracture closure significantly impacts the flow correlation

Figure 1 .
Figure 1.Representation of a geological fracture, with definitions of the various relevant geometric quantities and flow boundary conditions.Wall topographies are both generated over a 1,024 × 1,024 nodes mesh with the following parameters: L/L c = 8, σ a /〈a〉 = 0.2, L c = 0.1 m, 〈a〉 = 1 mm, and H = 0.8.

Figure 2 .
Figure 2. Average longitudinal autocorrelation functions of the velocity magnitude for the flow of a Newtonian and a shear-thinning fluid in two fracture realizations, one with L/L c = 1 (a) and one with L/L c = 16 (b).The parameters adopted otherwise for the simulations are: σ a /〈a〉 = 0.5, 〈a〉 = 10 −3 m, H = 0.8, and L c = 0.1 m.

Figure 3 .
Figure 3. Impact of aperture heterogeneity on the ratio T 0 /T 0‖ and normalized flow correlation length along the direction of the external pressure gradient, ℓ 0 /λ-panels (a, b) respectively.Parameters that describe the aperture heterogeneity are the fracture closure σ a /〈a〉 (x-axes) and the ratio L/L c .Different colors refer to different values of L/L c : yellow is 1, orange is 4, and 16 for purple.Solid lines are the polynomial chaos expansion estimations, circles denote data from the MC simulations (among which black edges denote ensemble averages), and bars indicate the first and third quantiles.
Figure 5 depicts the PCE-predicted behavior of the ratio of the fracture transmissivity T to that of the corresponding parallel plate geometry, T ‖ , and of the amplitude of the relative confidence interval  Δ   ‖

Figure 4 .
Figure 4. Dependence of the fracture transmissivity normalized by the corresponding transmissivity for Newtonian flow, T/T 0 , on both the dimensionless external pressure gradient ∇P/∇P c and the fracture closure σ a /〈a〉, produced for three different L/L c ratios: 1, 4, and 16, illustrated respectively in panels (a-c).For each of these panels, the inset within the left side figure shows the topography of T/T 0 in the (  ∇ ∕∇c , σ a /〈a〉) parameter space, while the main plots are 1D profiles extracted from this topography at given closures, the black dashed line denoting the parallel plate model transmissivity T/T 0 and the solid black line the asymptotic power law trend; the right hand side smaller figure is a plot similar to the aforementioned inset, but showing the amplitude of the confidence interval.

Figure 5 .
Figure 5. Dependence of the fracture transmissivity normalized by that of the parallel plate of aperture equal to the rough fracture's aperture, T/T ‖ , on both the dimensionless external pressure gradient  ∇ ∕∇c and the fracture closure σ a /〈a〉, produced for different ratios L/L c : 1, 4, and 16, for panels (a-c), respectively.The plots are otherwise identical to those of Figure 4.

Figure 6 .
Figure 6.(a) Dependence of the ratio of the longitudinal correlation length of flow velocities to that of the aperture field, ℓ/λ, on both the dimensionless external pressure gradient  ∇ ∕∇c and the fracture closure σ a /〈a〉, produced for different ratios L/L c : 1, 4, and 16, illustrated from top to bottom.Panels (b, c) display fracture aperture realizations on the left-hand side and velocity magnitude maps on the right-hand, for closure σ a /〈a〉 equal to 0.8 and 0.2, respectively, and for the three investigated correlation lengths.