The capillary pressure vs. saturation curve for a fractured rock mass: fracture and matrix contributions

The fractal topography of fracture surfaces challenges the upscaling of laboratory test results to the field scale, therefore the study of rock masses often requires numerical experimentation. We generate digital fracture analogues and model invasion percolation to investigate the capillarity-saturation Pc-Sw fracture response to changes in boundary conditions. Results show that aperture is Gaussian-distributed and the coefficient of variation is scale-independent. The aperture contraction during normal stress increments causes higher capillary pressures and steeper Pc-Sw curves, while shear displacement results in invasion anisotropy. The three-parameter van Genutchen model adequately fits the fracture capillary response in all cases; the capillary entry value decreases with fracture size, yet the fracture Pc-Sw curve normalized by the entry value is size-independent. Finally, we combine the fracture and matrix response to infer the rock mass response. Fracture spacing, aperture statistics and matrix porosity determine the rock mass capillarity-saturation Pc-Sw curve. Fractures without gouge control the entry pressure whereas the matrix regulates the residual saturation at high capillary pressure Pc.

Normal stress σ y [Pa] Material yield stress φ i [m] Sinusoid phase ω i [m] Sinusoid frequency Mixed fluid conditions are common in fractured rocks and affect both engineering applications and natural processes 1,2 . Examples include energy resource extraction such as hydrocarbons and geothermal 3,4 , the vadose zone and environmental remediation 5,6 , geological storage of CO 2 and nuclear waste 7,8 , and hydrothermal alteration and ore deposition 9 .
In the absence of gouge, fractures provide preferential flow pathways that dominate fluid flow. Then, invasion can readily bypass wetting (and even non-wetting) fluids in the matrix and fast injection studies tend to conclude low displacement efficiency (e.g., CO 2 injection due to viscous fingering, and aggravated by fractures and heterogeneity). On the other hand, long-term analyses must be based on capillarity-saturation curves at thermodynamic equilibrium.
Pore size distribution, pore interconnectivity and spatial variability determine the relationship between the capillary pressure and the degree of saturation at equilibrium P c -S w 10 . These constitutive equations capture the entry pressure and the sensitivity of saturation to changes in capillary pressure. Most numerical codes rely on these expressions to analyze coupled processes 11,12 .
Measurement techniques for the intact rock matrix are well established and rely on porous plates and tensiometers 13 , mercury injection 14 , centrifuge methods 15 , controlled suction 16,17 , dew point hydrometer 18 , and thermal/electrical conductivity measurements on calibrated porous stones 19 . The measured capillary pressure-saturation curves are relatively smooth for intact rocks and can be fitted with simple functions 20,21 .
In general, we can anticipate that the geometric aperture and its spatial correlation control capillary phenomena in fractures 22,23 . However, stress-sensitive aperture [24][25][26] , fractal fracture topography 27 , specimen size limitations and fracture-matrix interaction hinder the experimental determination of P c -S w curves in fractures and fractured rocks [28][29][30] . Furthermore, P c -S w curves would take weeks to months to reach thermodynamic equilibrium, particularly when trapped zones equilibrate through dissolution-diffusion processes 31,32 ; in fact, most published laboratory studies report rapid fluid invasion tests 33 .
In view of spatial and temporal scale limitations, numerical experimentation emerges as a necessary tool to study mixed fluid conditions in fractured rocks under long-term thermodynamic equilibrium. Previous numerical studies have explored the link between the fracture aperture and capillary pressure [34][35][36] , and their associated changes with normal stress 37,38 . Yet, realistic fracture surface generation algorithms, the evolution of aperture with normal stress and shear displacement, and implications on the capillarity of fractured rock masses require further research.
This study explores the relationship between capillary pressure and saturation in fractured rocks, with emphasis on fracture-matrix interaction at equilibrium following the invasion of a non-wetting phase (Note: fast immiscible fluid invasion and transient saturation are not part of this study).

Fracture studies: methodology
First, we investigate the effects of normal stress, shear displacement and fracture surface topography on the fracture capillary response. We generate digital fracture analogues that resemble natural fractures and bring the two rough surfaces together to create the fracture pore space. Then, we use network model simulations to investigate the effects of surface roughness and matedness on the fracture capillarity-saturation curve, and its evolution during normal loading and shear displacement. Fracture generation, response to normal stress and shear displacement, and fluid invasion modeling are described next. where the reference wavelength is ref = 1m . This power equation implies a fractal topography and provides a convenient framework for numerical studies. The wave amplitude X(λ) [m] for a certain wavelength λ is related to the spectral density as G(λ) = C |X(λ)| 2 , where the scaling factor C = N∆x/4 [m] depends on the selected Fourier pair and transform definition (i.e., one-sided vs. two-sided). We use the Inverse Fast Fourier Transform to synthesize the surface roughness profile in space by assuming a uniformly distributed random phase (see similar surface generation models in Refs. 38,41 ). The fractal topography implies lack of a characteristic scale or www.nature.com/scientificreports/ representative equivalent size. However, fractal geometries have limits in natural systems; in our case, the longest wavelength is the fracture size under consideration. Two contiguous rock surfaces define a fracture. We capture matedness by considering the correlation between them. There are two end-members: "perfectly mated fractures" consist of two identical surfaces with opposite orientations, and "unmated fractures" which are made of two uncorrelated surfaces. Figure 1b shows the aperture size distribution of unmated fracture surfaces generated using various α-factors (Eq. (1)). where A f [m 2 ] is the total fracture area. The two fracture surfaces interpenetrate a vertical distance δ v to satisfy the computed contact area A c (Eq. (2)-see Fig. 1c). While mass is not conserved at contacts, its effect on aperture distribution is negligible due to the small contact area A c . The selected rigid-plastic model depends on a single variable -the yield stress σ y -in line with Ockham's principle of parsimony 42 . For comparison, the simplest elastoplastic model with post-peak softening requires ≥ 4 variables 38,43 and additional parameters would be needed to track the evolution of the yield strength with contact deformation 44 . The rigid-plastic model is simple to implement, supports robust analyses and is adequate for capillary studies.
Shear displacement. We impose shear by displacing the two surfaces as rigid bodies without changing their surface topographies; natural processes are more complex and involve asperity shearing and/or overriding. The periodic boundary condition maintains a constant fracture area A f during shear displacement whereby the moving surface re-appears on the opposite side. This boundary condition assumes the shear displacement is a subset of an infinite medium.
Invasion percolation for drainage. After shear displacement and normal loading, we map the resulting N × N aperture field h i,j with spatial resolution Δx = Δy onto a N p × N p square lattice with fourfold connectivity where N p < N (i.e., a "checkerboard" model-see Refs. [45][46][47] ). The height h avg of each cell is the average aperture in the original fracture within the corresponding area Δl × Δl, so that h i,j x 2 ≈ h avg l 2 .
The selected cell size Δl × Δl defines the aperture resolution. While the pressure-dependent saturation of a given pore depends on the pore geometry, the global trends exhibit limited sensitivity to local details: results from a focused numerical study reveal that aperture averaging in Δl × Δl affects only the lower end of the capillarysaturation P c -S w curve, i.e., at high capillary pressures and low degrees of saturation. Furthermore, apertures h are much smaller than wavelengths in natural fractures (h/λ < < 1-Ref. 27 ); therefore, the in-plane radius of curvature is negligible and capillarity-saturation results are unaffected by aperture averaging when Δl < 10Δx. Previous studies show that the in-plane curvature can influence invasion patterns 34,38,48 ; however, our fracture generation model ensures h/λ < < 1 and the capillary pressure is only a function of the aperture. www.nature.com/scientificreports/ The invasion percolation algorithm assumes equilibrium at any given pressure (i.e., neither viscous forces nor time effects), and non-trapping of the wetting phase. Various pore-scale phenomena justify the non-trapping assumption, including (Fig. 2): corner flow along rough surfaces 31 , fluid transport into the matrix 49,50 , water evaporation and vapor pressure equilibration 51 . While these processes have characteristic time scales, the equilibrium assumption disregards any transient trapping of the wetting phase.
We implement the Young-Laplace equation in terms of the aperture-induced curvature, P c = T s /h (Note: this expression assumes cylindrical interfaces given that h < Δl and a perfectly wetting mineral surface θ = 0°). The largest aperture connected to the inlet is invaded first and defines the entry value. All apertures connected to the non-wetting phase throughout the medium are candidates for invasion. Invasion proceeds by monotonically increasing the capillary pressure to define the capillary pressure versus saturation P c -S w curve.

Results
The fracture P c -S w curve. Fractures vary over a wide range of length scales. Figure 3 shows the aperture size distribution for three uncorrelated and unmated fractures of size L x × L y = 1 × 1 m, 0.1 × 0.1 m and 0.01 × 0.01 m subjected to zero normal stress (power spectral density parameters: α = 1 × 10 -6 m 3 , β = 2.9). The selected fracture variables represent natural conditions and highlight size effects. Each curve represents an average of 100 numerical realizations in order to obtain statistically representative results. The mean aperture size distribution  www.nature.com/scientificreports/ is Gaussian-distributed in the three cases, with almost the same coefficient of variation s h /μ h defined in terms of the aperture standard deviation s h and the mean aperture μ h (see Fig. 3). The capillary pressure P c at a given saturation S w increases as the fracture size decreases due to inverse relationship between aperture and fracture size. We fit the mean P c -S w curve using the three-parameter van Genutchen model 52 : where P 0 [Pa] relates to the entry value, and m 1 and m 2 capture the sensitivity of saturation to capillary pressure. This three-parameter model provides an excellent fit compared to the two-parameter van Genutchen and Brooks-Corey models. While the P 0 value is inversely proportional to the fracture size L, the m 1 and m 2 parameters remain constant regardless of the fracture size for surfaces generated with the same power spectral density α and β parameters (Eq. (1)); hence, the normalized (P c /P 0 ) − S w curves are scale-independent.
The effect of normal stress. As the normal stress increases, the aperture size distribution shifts towards smaller values and a cutoff at zero aperture emerges as contact yield results in zero aperture contact points (Fig. 4a). The coefficient of variation of the fitted truncated Gaussian distributions increases with normal stress, and the aperture field exhibits higher variability.
Reduced apertures require higher capillary pressures for the same degree of saturation and the slope of the P c -S w curve increases (Fig. 4b). The evolution of the van Genutchen model parameters P 0 , m 1 and m 2 are shown in the inset (Fig. 4b). Note that higher normal stress and fracture closure increase the degree of saturation of the wetting phase at constant capillary pressure (see arrow in Fig. 4b).
The effect of shear displacement. The shearing of unmated, uncorrelated fracture surfaces results in statistically identical aperture fields and the capillarity-saturation P c -S w trends remain the same (in the absence of asperity shearing and gouge formation).
Shear displacement causes aperture changes only when there is some initial degree of matedness between surfaces 27 . Figure 5 shows the capillarity-saturation P c -S w curves for initially mated fractures: as the shear displacement increases, the mean aperture increases and the capillary pressure decreases for a given degree of saturation.
Shear displacement induces anisotropy in the aperture field of initially mated fractures as aperture ridges emerge transverse to the shear direction. Results in Fig. 5 show that the ensuing anisotropy in aperture www.nature.com/scientificreports/ connectivity produces slightly different P c -S w curves -particularly at low P c values-when the capillary pressure is controlled at a boundary that is normal or parallel to ridges.
The P c -S w curve for the fractured rock mass. The fracture and matrix P c -S w curves combine to define the capillary response of the fractured rock mass. Consider three orthogonal fracture sets with the same spacing d in the three directions so that the repetitive unit is of size d 3 (Fig. 6a). The matrix porosity n, the mean aperture size μ h and the fracture spacing d determine the volume of voids in the matrix V M v ≈ nd 3 and in fractures The volume of voids in the matrix V M v becomes a significant fraction of the total volume of voids V M v +V F v for small aperture to spacing ratios μ h /d and high matrix porosity n: Conversely, the storativity in fractures gains relevance in rock masses with low matrix porosity n where η F = 1 − η M . The fracture and the matrix share the same capillary pressure at equilibrium, therefore, the resultant P c − S w curve is a void-volume average of the saturation contributed by the fracture and the matrix Let's consider the matrix and the fracture P c -S w curves and combine them for various matrix void fractions η M to estimate the rock mass capillary response. We compute the fracture P c -S w curve using the algorithm described above, and a similar algorithm for the matrix where pores are represented as connected tubes with capillary pressure P c = 2T s /r (see model details in Ref. 53 ). Results in Fig. 6 show that fractures -without gouge-control the entry pressure, whereas the matrix dominates the behavior at high capillary pressure. As the η M fraction increases, the rock mass capillary-saturation curves approach the matrix P c -S w curve.

Discussion
The long-term saturation of a fractured rock mass will depend on the capillary pressure and the degree of saturation at equilibrium P c -S w . This will determine the original oil saturation profile in the reservoir and the residual oil after production, the long-term CO 2 and H 2 storage capacity, the distribution of LNAPLs and DNAPLs contaminants and environmental remediation strategies, the residual distribution of hydraulic fracturing fluids and the relative permeability of the resulting fractured rock mass.
The numerical study revealed surprising emergent properties. In particular, why is the aperture coefficient of variation s h /μ h independent of fracture size? (Fig. 3). The fracture surface topography z is the sum of k independent sinusoids with amplitudes a i that follow the power law in Eq. (1) and random phase φ i , z = k i z j = i a i sin(ω i x j + ϕ i ) . The probability density function of the sum of random variables f z is the convolution of the individual density functions, f z = f z 1 * f z 2 … * f z k . In the limit k → ∞, the central limit theorem emerges and f z converges to a Gaussian distribution 54 . While the central limit theorem often requires f z i to be www.nature.com/scientificreports/ independent and identically distributed, f z can converge to a Gaussian distribution for non-identical density functions (Lyapunov's Central limit theorem-Ref. 55 ). Then, the fracture surface topography z satisfies a Gaussian distribution, z ~ N (μ z , s z 2 ). The fracture geometric aperture is the subtraction of the top z t and bottom z b surface topographies, h = z t − z b,. Therefore, the fracture aperture values also exhibit a Gaussian distribution when the two surfaces are uncorrelated h ~ N (μ zt − μ zb , s zt 2 + s zb 2 ). The scale invariant features in our numerically generated fractures result from the fractal nature of the surface roughness 56 .
The aperture Gaussian distribution f h allows us to obtain the lower bound of the P c -S w curve: the pore volume distribution is f v = (h × Δl 2 )f h , then, the cumulative distribution function of f v (from largest to smallest h) normalized by the integral of (h × Δl 2 ) corresponds to the wetting phase saturation S w with capillary pressure P c = T s /h.

Conclusions
The capillarity pressure vs. saturation response of fractured rock masses is needed for long-term analyses. However, the time needed to reach thermodynamic equilibrium and size-dependent fracture topology limit our ability to experimentally gather relevant capillarity-saturation curves.
Numerical experiments show that the capillary pressure versus saturation P c -S w curve for a fractured rock mass is determined by pore-scale characteristics in the fractures and the matrix, including aperture and pore size statistics, spatial variability and connectivity.
The fracture surface roughness can be synthesized as a sum of independent sinusoids. A power law relates the sinusoidal amplitudes to their wave wavelength. The central limit theory emerges and the aperture size follows a Gaussian distribution. The mean aperture increases with fracture-size, yet, the coefficient of variation s h /μ h is scale independent.
While the wetting phase could remain occluded in the matrix during non-wetting fluid invasion, trapping in fractures is limited. Wetting fluid connectivity allows for the wetting fluid to escape through a network of connected corners formed by the surface roughness or via the permeable porous matrix. Long-term equilibrium in liquid-vapor systems further promotes non-trapping conditions in fractures.
The increase in normal stress contracts the aperture and results in higher capillary pressure, increased aperture variability, increased true contact area and steeper P c -S w curves. Capillary invasion reflects a slight anisotropy www.nature.com/scientificreports/ during shear and affects saturation at low capillary pressures. The three-parameter van Genutchen model captures the fracture P c -S w response and its evolution with normal stress and shear displacement. While the capillary entry value is inversely proportional to the fracture size, the normalized capillarity-saturation curves (P c /P 0 ) − S w follow similar trends and are scale independent. At equilibrium, fracture spacing and aperture statistics combine with the matrix porosity to determine the capillary pressure versus saturation curve for a rock mass and its storativity. In the absence of gouge, the fracture P c -S w curve controls the entry pressure, whereas the matrix regulates the rock mass residual saturation at high capillary pressures.

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author.