Influence of Local Aperture Heterogeneity on Invading Fluid Connectivity During Rough Fracture Drainage

Determining the (in)efficiency of wetting phase displacement by an invading non-wetting phase (drainage) in a single fracture is key to modelling upscaled properties such as relative permeability and capillary pressure. These constitutive relationships are fundamental to quantifying the contribution, or lack thereof, of conductive fracture systems to long-term leakage rates. Single-fracture-scale modelling and experimental studies have investigated this process, however, a lack of visualization of drainage in a truly representative sample at sufficient spatial and temporal resolution limits their predictive insights. Here, we used fast synchrotron X-ray tomography to image drainage in a natural geological fracture by capturing consecutive 2.75 μm voxel images with a 1 s scan time. Drainage was conducted under capillary-dominated conditions, where percolation-type patterns are expected. We observe this continuously connected invasion (capillary fingering) only to be valid in local regions with relative roughness, λ b ≤ 0.56. Fractal dimension analysis of these invasion patterns strongly aligns with capillary fingering patterns previously reported in low λ b fractures and porous media. Connected invasion is prevented from being the dominant invasion mechanism globally due to high aperture heterogeneity, where we observe disconnected invasion (snap-off, fragmented clusters) to be pervasive in local regions where λ b ≥ 0.67. Our results indicate that relative roughness has significant control on flow as it influences fluid conductivity and thus provides an important metric to predict invasion dynamics during slow drainage.


Introduction
Subsurface fluid flow primarily transpires in porous rocks, however, in low-permeability formations such as mudrocks, interconnected rock fractures can govern flow.Here we focus on drainage of fractures, i.e. the immiscible displacement of a wetting phase, wp, by a non-wetting phase, nwp.This process is fundamental to many subsurface engineering applications, from contamination remediation of non-aqueous phase liquids in groundwater (Detwiler et al. 2009), enhanced oil recovery (Nguyen et al. 2018), and the displacement of resident brine (wp) during CO 2 (nwp) migration from a storage reservoir through a naturally fractured sealing formation (Phillips et al. 2020;Rizzo et al. 2024;Snippe et al. 2022).
Modelling field-scale fracture drainage requires an understanding of the macroscopic effective properties that purport to describe single-fracture two-phase flow, namely capillary pressure, P c , and relative permeability, k r .These constitutive relationships are commonly modelled to depend on wetting and non-wetting phase saturations, S w and S nw , but are also influenced by fluid properties (Glass et al. 2003), surface wetting (Dou et al. 2013;Lee et al. 2010), and gravitational, viscous and capillary forces (Glass et al. 1995(Glass et al. , 1998;;Loggia et al. 2009).All these properties are also central to determining porous media twophase flow properties (Méheust et al. 2002;Toussaint et al. 2012).Unlike porous media, natural fractures comprise opposing surfaces with variable roughness (Brown 1987(Brown , 1995)), which are highly sensitive to stress alterations (from geological confinement and natural/engineered fluid pressure alterations) that produce discrete contact points, and consequently, yield heterogeneous apertures that complicate flow (Kang et al. 2016;Phillips et al. 2021; Pyrak-Nolte and Morris 2000; Pyrak-Nolte and Nolte 2016; Wang and Cardenas 2016;Zimmerman et al. 1992).It has therefore been debated whether employing porous-media models to simulate two-phase fracture flow is valid as flow predominantly transpires in 2D rather than the 3D pore space of a rock matrix (Huo and Benson 2016), meaning increased phase interference (disruption in continuous flow) in fractures (Fourar and Bories 1995;Persoff and Pruess 1995;Watanabe et al. 2015), which also increases with higher roughness.Numerous fracture k r models reported in the literature (cf.Gong et al. 2021) are commonly used in commercial simulators.However, the absence of a unifying relationship describing drainage in a single rough fracture, coupled with the scarcity of successful matching between experimental and numerically simulated invasion patterns (Chen et al. 2018b;Yang et al. 2019), suggests that flow behaviour not captured by currently available k r correlations exists.
For two-phase fracture flow, internal fracture geometry heavily influences fluid invasion geometries (Babadagli et al. 2015;Hatami and Walsh 2022;Karpyn et al. 2007;Neuweiler et al. 2004;Nicholl et al. 2000;Yang et al. 2013).Central to robust predictions is identifying the link between fracture roughness and nwp invasion (in)efficiency.To do so, statistical aperture characterisation is required (Yang et al. 2012), as this captures the defining features of internal fracture geometry that impact flow (wall roughness, contact area, R c , aperture distribution).Determining the mean aperture, < b > , and its standard deviation, σ b , enables the calculation of relative roughness, λ b = σ b / < b > , a metric which has been shown to describe fracture geometry (Glass et al. 2003), and to influence nwp front stability and displacement efficiency (Hu et al. 2019;Wang and Bayani Cardenas 2017;Ye et al. 2015).Previous two-phase flow studies which have investigated the link between λ b and invasion geometry (Chen et al. 2017(Chen et al. , 2018b;;Hu et al. 2019;Yang et al. 2016) have provided insights into invasion characteristics.However, they are limited in the sense that: (i) They investigate drainage in laboratory-induced, synthetically generated or replica fractures.This raises concerns towards the applicability of the results for natural fractures generated via geological processes; (ii) Due to the lack of confinement, and the nature of the materials studied, fractures examined have large apertures on the scale of hundreds of µm to mm.This is likely not representative of fracture apertures found in fine-grained, low permeability mudrocks, at reservoir depths (and therefore stress conditions) of >1000 m, relevant for geoenergy applications.Smaller fracture apertures would result in a higher impact of fracture roughness on fluid flow and fluid distribution in the fracture; (iii) In contrast to steady-state investigations, experimental time-resolved imaging of fracture drainage has only been reported using camera-based imaging systems, which are resolution-limited (> 40 μm), meaning smaller-scale fluid invasion events that could have a significant impact on k r have not yet been visualised; (iv) λ b values range between ~ 0.05 and 0.4, while numerical studies have investigated a broader λ b range (~ 0-2) (Brush and Thomson 2003;Wang and Bayani Cardenas 2017).
Given the vast ranges in mineralogy, grain sizes and fracture modes found in nature, broader λ b ranges may be more representative of natural fractures as they generate immobile flow zones (note that published data on natural fractures are sparce).This increases the likelihood that especially during Haines jumps, inertial forces can influence the flow field within a fracture and yield more complex invasion patterns.For example, high aperture variability and roughness increase the propensity for complex invasion processes (snap-off, trapping) (Gong et al. 2021;Tokan-Lawal et al. 2015), which requires understanding due to its influence on k r and P c ; (v) Recent porous media drainage studies (Andrew et al. 2015;Herring et al. 2018;Mascini et al. 2021) demonstrate the propensity for complex invasion events under capillarydominated flow.However, for fracture drainage, such processes have yet to be visualised dynamically and thus visco-capillary balance during slow drainage and its ultimate influence on multiphase flow in fractures is still to be understood.
In this study, we investigate multiphase flow in a natural fracture under confinement and flow rates that are deemed relevant for flow through low-permeability siltstone.We utilise, for the first time, fast 4D X-ray tomography to image drainage over ~150 min in a natural geological fracture by capturing consecutive 2.75 μm voxel images at a 1 s scan time.We investigated a conceptual drainage model.This describes capillary-dominated unsteady-state two-phase fluid displacement as a connected-phase invasion percolation process for continuous connectivity of the nwp and possible controls on it.Global aperture heterogeneity (λ b = 1.1) prevents the invading nwp from propagating through the fracture as an interconnected fluid body, demonstrating that this traditional conceptualisation of drainage for this natural fracture might be a simplification only.Through visualisation of local fluid displacement events, we identify a correlation between connected or disconnected flow and local aperture λ b .The reduced fluid connectivity and therefore conductivity slows down the overall invasion process and thus influences the visco-capillary balance.

Sample Material
The experiment was performed on a cylindrical siltstone rock (see Table 1) from the Carmel Formation, a regional caprock sequence overlying a naturally leaking CO 2 -charged reservoir (Kampman et al. 2014).The sample was obtained from a cored interval (depth = 198 m) recovered from the CO2W55 well in Green River, Utah (USA) (Kampman et al. 2013).Several natural fractures were identified in the parent material, which were cored to satisfy size requirements for our flow apparatus.The optimum sample was selected by imaging each fracture using a laboratory μ-CT scanner at Ghent University's Centre for X-ray computed tomography (UGCT; (Bultreys et al. 2016)) (see Text S1 and Figure S1 in Supplementary Material (SI)).

Fluid Flow Apparatus
Drainage was performed by first fully saturating the fracture with the wp, brine (S w = 1), followed by injection of the nwp, decane (see Table 1 and Text S2 for fluid properties).Decane was injected from bottom to top at a constant flow rate (Q = 200 nL/min), which yielded a capillary-dominated regime analogous to most geological flows (Piri and Karpyn 2007), quantified by a macroscopic capillary number (describing the ratio of viscous to capillary forces) of Ca = 1.65 * 10 -7 (cf.Table 1 for formulation).Flow was continuous throughout the experiment, meaning that the pressure gradient and visco-capillary balance was maintained during imaging.See Text S3 and Figure S2 for the full experimental procedure and flow apparatus.

Synchrotron Imaging and Image Processing
4D X-ray imaging was performed using the X02DA TOMCAT beamline at the Swiss Light Source, Paul Scherrer Institute (Villigen, Switzerland).We imaged drainage over ~ 150 min at a 2.75-μm voxel size and 1 s scan time (plus an additional 2 s to rotate the sample back to its initial position before the following scan could proceed).Time-resolved (unsteadystate flow) imaging was supplemented with steady-state imaging before the experiment, which acted as high-quality reference images (Table S3).See Text S4 for full image acquisition details.
Image processing was conducted using ImageJ (https:// imagej.nih.gov/ ij/) and Per-Geos 2020.2 (Thermo Fischer Scientific).Reconstructed tomograms were cropped to only include the fracture and immediately bordering rock matrix (see Figure S3) and then binned (2 × 2 × 2) to yield tomograms measuring 936 × 228 × 740 voxels with a 5.5-μm voxel size which reduced noise and computational load.A 2D noise-reducing median filter was applied to all reconstructed volumes.The fracture (pore) and matrix (solid) were  a The discrepancy between the actual total scans (i.e. 8 sequences × 401 scans per sequence = 3,208) and the quoted number of scans here is the result of discarding redundant scans from the final analysis.b, c Value obtained from PubChem, open chemistry database.d Value obtained from Zeppieri et al. (2001).e Estimated as Ca = Qμ/γA c (Ferer et al. 2011), where A c is the average cross-sectional area measured from the 2D image slices.f Estimated as Bo = ((ρ w -ρ nw )g < b > 2 /γ) (Loggia et al. 2009).g Estimated as Bo * = Bo-Ca (Méheust et al. 2002).h Estimated using the method of Andrew et al. (2014).See Text S8 and Figure S8.i Estimated as = < b > 2 /12 (Nicholl et al. 2000).j Denote fracture properties measured directly from the field of view (cf. Figure 1a).k Estimated using the method of Barton et al. (1985).See Text S12.l Calculated as λ b = σ b / < b > (Glass et al. 2003 , λ b 1.1 segmented from the high-quality reference image (Figure S4).From this, aperture characteristics were measured (Fig. 1c, Table 1, Text S9).All time-resolved images were subtracted from the high-quality brine-saturated images, which yielded differential images with only decane remaining.2D median filtering was again performed on the subtracted images to further increase the signal-to-noise ratio.The segmented matrix was applied as a mask on all time-resolved drainage images, followed by intensity-based thresholding to yield the final segmented decane.See Text S5 for full image processing details.Sensitivity analysis was performed to investigate the impact of image binning on nwp segmentation.Global nwp saturation and connectivity trends were similar for the 5.5 μm and 2.75 μm voxel size images (see Figure S5 and Text S7), but local connectivity was higher at certain time intervals in the 2.75 μm voxel size images.Consequently, the global connectivity (Fig. 2) and contact angle analysis (Figure S8) was performed on the 2.75 μm voxel size images, while all other analysis was performed on the binned (5.5 μm) images.

Experimental Considerations
During prolonged imaging, high-energy synchrotron X-rays can cause gas bubble formation in the brine due to the presence of potassium iodide (KI).To ensure that gas bubbles did not impact the nwp (decane) analysis, dark phase greyscale values (which could be attributed to decane or gas bubbles) were examined throughout the sequences.At t = 152 min, gas bubbles began to form, which appeared different in morphology and greyscale (darker) compared to decane (see Text S10 and Figure S6).Consequently, all images after t = 152 min were discarded, verifying that all analyses before t = 150 min focused solely on the brine and decane dynamics.
Capillary pressure P c was measured continuously throughout the experiment (see Text S3 and Figure S7), with a mean P c = 4.1 kPa, which aligns with P c estimation calculated via the incorporation of contact angle measurements (Hu et al. 2018) (see Figure S8), giving P c = 2cosϴ/ < b > = ~ 4 kPa (depending on < b >).We do not assign P c fluctuations to specific fluid redistribution events as the entire sample was not imaged, meaning dynamics occurring outside the field of view could contribute to certain P c fluctuations (Berg et al. 2013).
Nwp flow was concentrated through the fracture as (i) the fracture permeability k f was orders of magnitude greater than the matrix permeability k m (Table 1); (ii) the matrix pore sizes were significantly smaller than the mean aperture < b > , meaning nwp invasion into the matrix was hindered by a capillary entry pressure barrier, and (iii) for drainage configurations (i.e.water-wet system where the wp predominantly coats fracture walls and contact angles < 90° (Figure S8)), nwp invasion into the matrix is unlikely.Visual and greyscale μ-CT image inspection showed no nwp (or wp) in the matrix.
Here, fluid injection proceeded from fracture bottom to top.An experimental oil-water study showed a range of two-phase flow regimes depending on combinations of Bo and Ca (Loggia et al. 2009).Bo for our study is < < 1 (Table 1), suggesting capillarity dominates and buoyancy has a negligible impact on interface geometry (Yang et al. 2016).

Global Overview of Fracture Drainage
Our main observation is the dynamic nwp evolution during slow fracture drainage.This, to our knowledge, has not yet been visualized with high spatial and temporal resolution using synchrotron imaging in the literature.Drainage progresses as several disconnected nwp bodies that adopt variable invasion mechanisms depending on location within the fracture (Fig. 1b), not as an interconnected body of emergent capillary fingers.
We observe rapid (millisecond) interfacial Haines jumps to be a common invasion mechanism throughout the fracture, evidenced by rapid (1 s) P c fluctuations (Figure S7).Haines jumps under capillary-dominated conditions have been observed in a 3D Navier-Stokes simulation of rough fracture drainage (Chen et al. 2018a) and are commonly reported in porous media (Andrew et al. 2015;Armstrong and Berg 2013;Berg et al. 2013;Biswas et al. 2018;Haines 1930;Måløy et al. 1992;Moebius and Or 2014).Alongside Haines jumps, other distinctive invasion mechanisms identified include, but are not limited to, capillary fingering, Roof snap-off and fragmented clusters.These processes can unfortunately not be visualized because they occur at smaller time scales (i.e., < seconds) than the time it takes to scan the fractured samples (i.e., seconds) with resulting phase distributions presented in Fig. 2a-f.Furthermore, distal snap-off events have been observed for porous media (Andrew et al. 2015) which can be used to explain the absence of a connected nonwetting phase in the section that is shown in Fig. 3d, h.We discuss the propensity for more connected invasion (e.g.capillary fingering) or more disconnected invasion (e.g.snap-off) in Sect.3.2.We attribute this wide array of invasion mechanisms to the high global aperture variability (a lognormal distribution with λ b = 1.1) (Fig. 1c), which is supported by previous numerical investigations that suggest more compact invading geometries with lower λ b (Yang et al. 2019).However, we also acknowledge that local surface wettability variations could contribute.The adherence to/violation of capillary dominance in such proximity suggests that simpler aperture geometries that do not represent the roughness of geological fractures (e.g.λ b < 0.5) provide an incomplete (and therefore potentially unphysical) account of slow drainage in geological fractures.Consequently, nwp flow connectivity may be overestimated, resulting in a possible overestimation of CO 2 or H 2 leakage rates from subsurface storage reservoirs.
Considering the classical conceptualization of capillary fingering and invasion percolation (IP), originally proposed for porous media (Chandler et al. 1982;Wilkinson and Willemsen 1983), we would expect an initial increase in the number of nwp clusters at low S nw .As S nw increases, individual clusters start merging, whereby the number is decreasing due to the establishment of more connections.Here, we find a near linear increase (R 2 = 0.9) in the number of disconnected clusters with increasing nwp complexity (i.e.pervasive fragmentation, variable cluster orientation & size; Figs.2g and h) with S nw .This finding, which deviates from the expectation of drainage under IP, has also been reported in a recent study of slow drainage in porous media (Herring et al. 2018).Here, increasing nwp disconnections have been observed experimentally as drainage proceeds (for S nw < 50%) and have been verified via Lattice Boltzmann simulations.However, to our knowledge, this behaviour has not yet been reported for rough fracture drainage.We postulate that understanding drainage at low S nw is important for characterising the onset of fracture leakage along sealing layers in subsurface storage settings.We here provide evidence that the nwp has propagated through the field of view without exceeding a S nw of ~ 5%, implying the establishment of an intermittently connected pathway (that is continuously breaking and reforming) and that can transport fluid in its current configuration.This may also be the case for capillarydominated fracture leakage at low S nw , which further imaging studies may be required to confirm.
Our observations support previous time-resolved imaging studies in porous media (Andrew et al. 2015;Singh et al. 2017), which report persistent disconnection-reconnection events, and highlight the similarities between slow drainage in porous media and rough fractures.This further supports the notion that modelling rough fracture drainage as an invasion percolation process is incomplete for low S nw .While the shortcomings of an IP approach have been discussed previously (thus leading to its modification) for other drainage aspects such as in-plane curvature, aperture geometry and capillary-to-viscous forces transition (Chen et al. 2017(Chen et al. , 2018a;;Glass et al. 1998;Yang et al. 2019), our study provides evidence for phase disconnections of the nwp (particularly at low S nw ), even at low capillary numbers.As a result, modelling efforts that do not account for disconnected phase invasion mechanisms (see above), may yield inaccurate fluid distributions.It also challenges P c estimations, as boundary P c measurements cannot account for disconnected clusters.Here, clusters comprise 77% of the total nwp volume (0.053 mm 3 ) at the end of the experiment.
We quantify invasion patterns using the fractal dimension, D f (-), which is a classical measure of the degree to which a pattern fills a 2D space.Increasing complexity with time is evidenced by an increasing D f (1.58 to 1.7) (Fig. 2g & Figure S9a -9f).Notably, the mean D f (1.68) aligns with previously reported measurements (D f = 1.66-1.68)for viscous fingering in rough fractures (Chen et al. 2017;Hu et al. 2019), rather than capillary fingering.It should be noted that these studies measured D f for higher S nw (~ 20%), compared here to S nw ~ 5% (Fig. 1g).This further illustrates the challenges of predicting nwp invasion morphologies using a capillary number that is based on mean quantities such as the average cross-sectional area, even in rough fractures with high capillarity.However, we emphasise that this does not mean that viscous forces dominate, but likely have an influence.It rather seems to be the case that the dynamics are governed by the heterogeneity and not the mean value of the aperture.The Bond number and generalized Bond number obtained for this study (Table 1) are identical and small, indicating that surface tension is dominant, and gravity can be neglected.Consequently, stable fluid displacement is likely.The authors further describe this stable displacement as a process where gravity limits fast-growing fingers from developing much faster than slower ones.While it may not be conclusive in the context of fluid connectivity, it reasonably suggests important dynamics beyond the dominant unstable displacement (fingering).Within this complex flow environment, we think that the fracture geometry controls the strength and speed of the fingers (for instance, the rougher a fracture the slower the flow) irrespective of acting forces.In fact, Mascini et al. (2021) found that fluid conductivity (in mixed wet and highly heterogeneous porous media) can locally dominate the flow dynamics over capillary forces, even at a low capillary number.This further informs on the impact of roughness which limits connection of the nwp.

Roughness Control of Invasion Path Connectivity
Despite a globally low Ca number, local aperture variations promote diverse nwp invasion morphologies that vary in terms of connectivity.While we studied the wider aperture region above, we identify four separate aperture regions within the global aperture field below.These represent observations across the entire aperture field and facilitate either connected flow (i.e.capillary fingering) or disconnected flow (snap-off) and characterise their respective aperture characteristics (Fig. 3).Here, regions for these analyses were selected based on visual inspection, where we observed changes in saturations or fluid occupancy over time.

Locally Connected Invasion
For the case study presented here, we observe connected pathway flow (capillary fingering; CF1 & CF2) to occur in rather low λ b regions (0.56 and 0.53) (Figs. 3b,.Compared to the global aperture field, these regions exhibit narrow aperture ranges (~ 0-85 μm), low σ b and < b > , and a smaller skewness in the aperture distribution (Fig. 3b).Both CF1 and CF2 emerge initially as rapid nwp Haines jumps into the large apertures (Fig. 3d, h), which occur below the temporal resolution.Note that Haines jumps are recorded through a small (< 1 kPa) decrease in P c (Figure S7).Following this, invasion evolves radially as interconnected capillary fingers that are unconstrained by the macroscopic flow direction (i.e.forward and reverse front propagation), which aligns with previously reported behaviour for rough fractures under high capillarity (Chen et al. 2017(Chen et al. , 2018a;;Hu et al. 2019).These fingers preferentially invade the centre of the largest apertures, while the wp primarily resides on the fracture walls and in corners and crevices (Figure S11).This behaviour adheres to typical drainage invasion percolation patterns in fractures and porous media (Arshadi et al. 2018;Mascini et al. 2020), where regions with larger apertures are preferentially invaded as they require a lower P c to invade.
D f analysis of nwp invasion patterns in these local regions vary from the global D f values (Fig. 2g).The D f of CF1 and CF2 measure 1.82 and 1.83 (Figure S9g,h), respectively.This is consistent with experimentally observed capillary fingering in low (< 0.3) λ b fractures (D f = 1.78-± 0.06) (Hu et al. 2019), with theoretical values in porous media based on percolation theory (D f = 1.82) (Wilkinson and Willemsen 1983), and with numerical and experimental porous media values (D f = 1.82-1.83)(Holtzman and Juanes 2010;Løvoll et al. 2004;Toussaint et al. 2005;Zhao et al. 2016).Based on the match of the fractal dimensions we find that for the case study presented here, invasion percolation is able to describe nwp invasion in aperture regions with λ b ≤ 0.56 at low Ca numbers.It is also comparable to porous media, despite the 2D pore network geometry of a rough fracture.While for our case example, invasion percolation is observed for lower λ b (≤ 0.56) and disconnected flow for larger λ b (≥ 0.67), it is difficult to generalise this finding and further samples need to be tested for confirmation.
Throughout their respective evolutions, both capillary fingering bodies are continuously supplied by non-local nwp Haines jumps which are indicated by small changes in capillary pressure towards lower values (SI, Figure S7).This implies that even in low λ b aperture regions, which are conducive to classical percolation-type capillary fingering (i.e.CF1 & CF2), the high degree of aperture heterogeneity surrounding these regions limits the ability for capillary fingering to be the pervasive invasion mechanism globally.Notably, both instances of capillary fingering develop early in the experiment (see temporal color-coded Figs.2a-f), and the connected pathways they establish remain unaltered for the remainder of the experiment, barring non-local periodic interfacial recessions in response to P c changes (Figure S11).This suggests that lower λ b regions are the most prone for invasion in terms of P c threshold and are less prone to fluid redistribution events that alter the global S nw .

Locally Disconnected Invasion
Disconnected invasion (snap-off; SO1 & SO2) occurs in higher λ b regions (0.67 and 0.71 in these examples) (Figs.3b, k-n, o-3) that comprise larger aperture ranges (~ 0-230 μm), higher σ b and < b > , and a higher positive skewness in the aperture distribution (Fig. 3b).Notably, higher contact area R c does not yield more disconnected flow, with both instances of capillary fingering exhibiting higher R c (8.6 and 5%) than for snap-off (3.3 and 3.1%), suggesting that the configuration of apertures open to flow (non-zero apertures), rather than the percentage of no-flow regions (zero apertures), exert a more significant influence on nwp connectivity.
Snap-off arises from the loss of connectivity when the nwp flows through aperture bottlenecks, that is, narrow constrictions that form a bridge between two larger aperture regions (Figs.3k, o).This aligns with previous experimental studies which attribute the trapping of nwp clusters in larger aperture pockets to aperture regions with significant heterogeneity (Karpyn et al. 2007).Here, similar to capillary fingering, snap-off events are not confined to the forward-most nwp interface or the macroscopic flow direction and can hence be described as local (Roof snap-off) and non-local (distal snap-off) events (cf.Andrew et al. 2015).The presence of distal snap-off is especially important to consider as these events are associated with more persistent fluid configurations and can alter displacement pathways (i.e.connected nwp can lose access to pore space), potentially decreasing displacement efficiency (Zacharoudiou et al. 2018).One of the few studies to incorporate snap-off in an invasion percolation model successfully replicated experimental fluid distributions (e.g.Amundsen et al. 1999) but was limited by an inability to predict snap-off based upon aperture heterogeneity.We propose that knowledge of λ b can aid in determining when invasion percolation models should incorporate disconnected invasion mechanisms.

Conclusions
We experimentally investigated slow (capillary-dominated) drainage in a natural rough fracture by using 4D X-ray imaging.This is of relevance for predicting the risk of leakage from subsurface gas (CO 2 , H 2 ) storage reservoirs.We show that modelling slow rough fracture drainage as a continuous-phase invasion percolation is likely incomplete for low non-wetting phase saturations (~ 5%).Consequently, such model approaches may yield physically unrealistic fluid distributions and therefore unrealistic capillary pressure and relative permeability estimations.Further evidence of this deviation from an IP process is demonstrated by an increase in the number of disconnected non-wetting phase clusters with increasing non-wetting phase saturation.Adherence to percolation-based models is observed in aperture regions with a lower relative roughness, where connected capillary fingers propagate in a comparable manner (fractal dimension) to that in low λ b fractures and porous media.
Conversely, aperture regions with a larger λ b promote disconnected invasion.Practically, the prevalence of disconnected invasion during low capillary number drainage has implications for subsurface engineering applications, such as CO 2 or hydrogen (H 2 ) storage, where snap-off (capillary trapping) of the non-wetting phase (e.g.CO 2 , H 2 ) may reduce leakage rates through fracture networks, which is advantageous for secure subsurface storage.While further research is required on natural fractures in low-permeability rocks, our results indicate that the aperture heterogeneity (and thus relative roughness) encountered influences of fluid conductivity and overall invasion dynamics during drainage, altering the visco-capillary balance which dictates the capillary number.This makes λ b an important parameter to consider beyond capillary number for characterizing flow.As a result, a

Fig. 1
Fig.1Global fracture overview.a 3D aperture rendering from the field of view used in the analysis.The surrounding rock volume is rendered transparent apart from the bottom section which is shown for reference.The white space denotes contact points (zero-aperture regions).See Table1for fracture properties.b 3D rendering of two-phase (brine & decane) flow through the fracture at the end of the experiment (normalized time = 1, see Fig.2g).See Movie S1 (https:// doi.org/ 10. 6084/ m9.figsh are.24050 940.v1) for the full drainage sequence.c Aperture size distribution and inset of cumulative size distribution

Fig. 2
Fig. 2 Global two-phase flow experiment overview.a-f 2D visualizations of nwp progression over the experiment.Temporal color-coded visualizations of the corresponding saturation maps, showing only the nwp.These visualization windows are cropped to only include the nwp invasion path and thus appear smaller than the full field of view shown in.Fig.1g Nwp complexity over the experiment, quantified by fractal dimension (Figure S9) and the number of individual nwp clusters measured every 50th scan (see Text S7).h Inset of nwp saturation and cumulative volume over the experiment measured every 10th scan

Fig. 3
Fig.3Local aperture regions for Capillary Fingering (CF1 and CF2) and. a 2D view of the global aperture field, with the locations of CF1, CF2, SO1 and SO2 highlighted.See FigureS10for 2D views of the local apertures.b Corresponding aperture size distributions and characteristics for CF1, CF2, SO1 and SO2.c-f Capillary Fingering example 1 (CF1).The local aperture region is shown in c and the evolution of capillary fingering is shown in d-f.The wp is shown in blue, nwp in red and the surrounding rock is rendered transparent.See Movie S2 (https:// doi.org/ 10. 6084/ m9.figsh are.24050 964.v1) for the full image sequence.g-j Capillary Fingering example 2 (CF2).The local aperture region is shown in g and the evolution of capillary fingering is shown in h-j.Movie S3 (https:// doi.org/ 10. 6084/ m9.figsh are.24050 955.v1) for the full image sequence.k-n Snap-Off example 1 (SO1).The local aperture region is shown in k and the evolution of snap-off is shown in (l-n).See Movie S4 (https:// doi.org/ 10. 6084/ m9.figsh are.24050 958.v1) for the full image sequence.o-r Snap-Off example 2 (SO2).The local aperture region is shown in o and the evolution of snap-off is shown in p-r.See Movie S5 (https:// doi.org/ 10. 6084/ m9.figsh are.24050 961.v1) for the full image sequence ▸

Table 1
Experimental Parameter ;Wang and Cardenas 2018) l Relative roughness