Direct numerical simulations of immiscible two-phase flow in rough fractures: Impact of wetting film resolution

The immiscible displacement of a wetting fluid by a non-wetting fluid (i.e., drainage) in rough-walled, open geological fractures holds significant relevance for various applications, encompassing subsurface fluid storage (Kim , 2011; Muller , 2019), subsurface contamination remediation (Lee , 2010), and oil and gas exploitation (Bikkina , 2016; Zhang , 2017). In hydrogeological frameworks, the subsurface environment is often conceptualized as a rigid, porous medium. However, many of the above applications involve a fractured or fissured subsurface, such as rock formations (Birkholzer and Tsang, 2000; Xue , 2018). Accurate predictions of two-phase flow phenomena in these applications necessitate resolving hydrodynamic-scale physics, in which interfacial interactions profoundly impact macroscopic two-phase flow dynamics (Meakin and Tartakovsky, 2009). Although significant research has been dedicated to studying the rich phenomenology of two-phase flows in porous media over the past 20 years, the modeling of immiscible displacement in rough fractures is a comparatively emerging field.

Rough fractures consist of the space between two solid walls that result from the solid material having been fractured. In the subsurface, in particular, fracturing occurs due to constraints imposed by tectonic activity or the quick cooling of solid rock. Irrespective of the solid material and the fracturing mode, the resulting fracture surfaces have a universal geometrical property, which is a self-affine scale invariance characterized by a roughness exponent (or Hurst exponent) (Bouchaud, 1997; Schmittbuhl , 1995). This property is valid at scales ranging from micrometers in metals and ceramics to thousands of kilometers in tectonic faults. The resulting geometry is asymptotically planar at large scales but gets all the rougher as one examines smaller scales. In a freshly fractured solid material, if there is no loose material, the two facing walls are perfectly matched, i.e., they have the same geometry. However, any shift in the walls' positions along their mean planes renders the distance between the two facing walls, denoted local fracture aperture, spatially dependent. In geological fractures, chemical weathering over geological times, as well as tectonic constraints and fracturing around contact points, also results in apertures that vary spatially, with fracture walls that have maintained their self-affine geometry.

The complex fracture geometry leads to complex flow, both for single-phase Newtonian (Brown, 1987; Zimmerman and Bodvarsson, 1996; and Méheust and Schmittbuhl, 2001) and non-Newtonian (Lavrov, 2014; Lenci , 2022a, 2022b; and Zhang , 2023) fluids, with in-plane channeling, which impacts the fracture transmissivity. This structural heterogeneity also impacts immiscible two-phase flow (Glass , 2001). For two-phase flows, the displacement patterns have been studied extensively for capillary-dominated flows. Depending on the roughness and correlation structure, the patterns change from invasion percolation structures to more compact patterns, as the in-plane component of the capillary pressure smoothes out small-scale fluid clusters (Glass , 2001; Neuweiler , 2004; and Yang , 2012). The patterns were investigated in detail by Yang (2016), who showed that closing a fracture more leads to less compact patterns due to an increase in the out-of-plane component of the capillary pressure, relative to the in-plane component, in agreement with the previous findings. They also showed that the latter in-plane component suppresses the trapping of fluid clusters below the correlation length of the fractures and changes the characteristics of the size distribution of trapped clusters. Yang (2019) studied the fluid distribution for a wider range of capillary numbers, for which the patterns transitioned from capillary fingers to viscous fingers. The in-plane curvature influenced patterns for cases with low capillary numbers and large viscosity ratios, leading to not only more compact clusters but also a change in trapping behavior; the patterns' fractal dimension varied between 1.55 and 1.85 depending on the capillary number and viscosity ratio between the two fluids. Chen (2017) characterized the transition between capillary fingering and viscous fingering experimentally in detail from various observables, in a transparent replica of a fracture in granite rock; in particular, they observed different dependencies of the tip of the most advanced finger on the capillary number, depending on the displacement regime. The same group later extended these experimental findings by examining the corresponding flow velocity fields (Chen , 2018), using direct numerical simulations in the same geometry as that of Chen (2017), validated by comparisons with their experiments.

Numerical approaches for modeling two-phase flows in rough fractures have developed gradually in the past 20 years. The first models addressed very low capillary (quasi-stationary) displacements with invasion percolation (IP) schemes (Glass , 1998; Detwiler , 2009) inspired by similar approaches on two-dimensional (2D) porous media, with a constant in-plane component of the capillary pressure due to an in-plane curvature radius that was considered a characteristic length scale of the geometry. Yang (2016) also used an IP approach, but with an in-plane curvature that was computed in time and space as a purely geometrical quantity based on the geometry of the interface. Yang (2019) proposed a mesoscopic model treating each fluid phase at the continuum scale with proper boundary conditions at the fluid–fluid interface, including an in-plane curvature defined from the interface geometry. Pore network models (Fatt, 1956) have also been used for simulating fluid flow in permeable media, including rough fractures (Ferer , 2011) or fracture networks (Hughes and Blunt, 2001). Although computationally efficient, these models employ simplified representations of pore geometry and, especially in multiphase flow, simplified physics (Meakin and Tartakovsky, 2009). Recent advancements in computational capabilities have enabled the resolution of the Navier–Stokes (NS) equations, either directly (direct numerical simulations) or using methods that are equivalent to solving the NS equations, such as smoothed particle hydrodynamics (SPH) (Tartakovsky and Meakin, 2005) and lattice Boltzmann methods (LB) (Dou , 2013; Guiltinan , 2021). However, these methods face challenges in relating interaction forces to fluid properties. For example, LB methods require case-specific model calibration to match density ratio and surface tension based on LB parameters, with complex parameter interdependencies (Porter , 2012; Huang , 2011).

The direct numerical simulation (DNS) of the NS equations has superior numerical efficiency and is capable of handling large viscosity and density ratios (Meakin and Tartakovsky, 2009), hence its use in recent studies (Ferrari , 2015; Chen , 2018). Hydrodynamic scale modeling of immiscible flows in rough fractures is a complex and challenging task influenced by a multitude of factors, such as the complex interplay between the viscous, capillary, gravitational, and inertial forces (Detwiler , 2009; Glass , 1998); characteristic wall roughness and wetting properties (Auradou, 2009); intricate geometries and aperture variability (Glass , 2003); molecular scale phenomena such as moving contact lines and thin wetting films (Shi , 2018; Qiu , 2023); and interfacial fragmentation and coalescence (Amundsen , 1999; Santiago , 2016). Among these challenges, a crucial yet largely unexplored one is capturing the thin film deposited on the fracture walls by the displaced wetting fluid. In subsurface applications such as enhanced oil recovery, the wetting films have a significant impact on the macroscopic transport properties (Blunt , 1992; Hui and Blunt, 2000). For instance, films can speed up invading fluid's flow through constrictions, resulting in a lubricating effect that, depending on the viscosity ratio, can yield relative permeabilities greater than 1 (Blunt , 1992). Resolving thin films in hydrodynamics-scale models is computationally challenging, primarily due to the large aspect ratio of the different length scales (Zhao , 2019). For context, fracture dimensions along the mean fracture plane typically span centimeters to kilometers in the subsurface, with aperture thicknesses up to a few centimeters, while in laboratory replicas, the dimensions span up to tens of centimeters and the aperture up to a few millimeters. For comparison, film thicknesses range from a micrometer to ten micrometers. While some numerical studies have reported satisfying comparisons between the experimental measurements and numerical predictions using numerical models that do not resolve wetting films (Chen , 2018), to the best of our knowledge, the consequences of ignoring the film on the predicted flow and spatial distribution of fluid phases for drainage in open rough fractures have not been studied systematically. The main focus of the present study centers around this critical question.

To resolve the fluid–fluid interface, we have chosen to couple a NS solver with the volume of fluid (VOF) method (Hirt and Nichols, 1981). It is one of several interface-capturing methods, which also include the Level-Set (LS) (Sethian, 1996) and phase-field (PF) (Sun and Beckermann, 2007) methods. These Eulerian-based methods are well-suited for flows with complex interface motion and can handle large interface deformations such as breakup or coalescence (Gopala and van Wachem, 2008). The VOF method, in particular, has been effectively applied to model immiscible flow in both fractured (Chen , 2018) and porous media (Ferrari and Lunati, 2013; Ferrari , 2015). Experimental validation from these studies has substantiated the VOF model's capacity to reasonably predict invasion structures and faithfully reproduce phenomena such as viscous meniscus deformation, snap-off events, and coalescence. Valuable insights on the issue of resolving wetting films have emerged from porous media studies, where the issue of wetting films was successfully tackled using the VOF approach. For instance, Horgue (2013) conducted experimental and numerical investigations of liquid jets spreading across in-line arrays of cylinders confined within a Hele–Shaw cell. Additionally, Ferrari (2015) employed the VOF model to replicate primary drainage experiments in a 2D porous medium consisting of a Hele–Shaw cell filled with cylindrical grains. Inspired by prior studies on Taylor flow in microchannels (Gupta , 2009; Horgue , 2012), these two studies showed that to numerically capture the thin wetting film, the liquid film region must be covered by at least two grid cells.

Transferring the film resolution approach from these studies to the case of rough fractures is not straightforward. Complex fracture geometries present significant challenges for numerical grid generation procedures. As the location of the fluid–fluid interface cannot be known a priori, selective region refinements are not feasible. Abrupt changes in local aperture distribution and the presence of zero-aperture regions where fracture walls touch each other further add to the difficulties. Previous studies on bubble motion in tubes with wetting fluid film have revealed that film thickness scales as $Ca2/3$⁠, where Ca is the capillary number (Hodges , 2004; Bretherton, 1961). As Ca decreases and gets into ranges $<10−5$⁠, film thickness nears the nanoscale, requiring ultra-fine grid resolutions. To accurately describe multiphase flow in permeable media, another crucial consideration is the behavior of moving contact lines resulting from fluid–fluid interface interactions with solid surfaces (Meakin and Tartakovsky, 2009). Modeling fluid behavior around a moving contact line is complicated, as the no-slip boundary condition at the solid–fluid interface leads to a force singularity at the contact line (Pasandideh-Fard , 1996). Therefore, it is important to include “slip” when dealing with moving contact lines (Dussan , 1991). The magnitude of this slip is quantified by the so-called slip length, which is usually in the molecular range and is much smaller than any numerically achievable grid sizes (see Renardy , 2001 for details). In practice, numerical calculations implicitly introduce a slip length, approximately $Δxwall/2$⁠, with $Δxwall$ representing the typical cell dimension near solid walls (Renardy , 2001). Consequently, grid resolution near the walls significantly influences the contact line motion, making complete grid independence for a VOF code challenging in such scenarios (Rosengarten , 2006). Hence, conducting a systematic validation of the VOF model becomes paramount to ensure the accuracy and reliability of investigations of immiscible displacements in rough fractures across varying grid resolutions.

The main goal of this study is to comprehensively assess how resolving thin wetting films impacts model predictions of drainage in open rough fractures compared to not resolving them. Our study is restricted to the range of flow cases where the wetting film thicknesses do not approach the nanoscale and can be well resolved with realistic grid resolutions. We aim to identify scenarios where wetting films significantly affect results, necessitating their inclusion in the model. Our approach involves initial model validation considering a Hele–Shaw geometry based on the classical experiments of Saffman and Taylor (1958). Subsequently, we perform numerical drainage experiments with two distinct rough fracture geometries. We then compare outcomes from two configurations: one resolving the film and the other not. We analyze various hydrodynamic and macroscopic scale properties, including invasion patterns and average pressure drops, while discussing the computational requirements. Our numerical simulations utilize interFoam, a finite-volume, VOF-based code developed by OpenFOAM (OpenFOAM, 2012), which we have self-modified to address the specific numerical challenges.

The paper is structured as follows: In Sec. II, the numerical model is described in detail. The model is then validated on a Hele–Shaw geometry in Sec. III, while the results and discussions of the numerical experiments in rough fractures are presented in Sec. IV. Finally, Sec. V presents the conclusions and future outlooks.

In this section, for the sake of completeness, we briefly outline the numerical model as implemented in OpenFOAM (2012) and employed to perform the simulations presented hereafter. We describe the fundamental governing equations (Sec. II A) as well as an enhancement to the standard volume-of-fluid (VOF) OpenFOAM solver (Sec. II B) that we have implemented, and briefly present the solution procedure (Sec. II C). Readers interested in the details of the model implementation, numerical discretizations, and solution procedure are encouraged to refer to the provided supplementary material and references therein.

Note that this smoothing procedure is not included in the original OpenFOAM VOF solver, which thus had to be modified in this respect.

The numerical model presented above has been implemented in the open-source CFD software OpenFOAM (OpenFOAM, 2012). For all the flow domains under consideration, the 3D computational mesh consisting of purely hexahedral volume elements was generated using the blockMesh utility provided in OpenFOAM. The governing equations [Eqs. (1), (2), and (6)] were discretized spatially using the cell-centered-based finite volume technique (Berberović , 2009) with a second-order accuracy in space. The implicit Euler scheme was used to perform the time integration, resulting in first-order accuracy in time.

The solution procedure starts by solving the phase fraction advection equation [Eq. (6)] with the compression velocity term [Eq. (7)] to restrict numerical diffusion. The resulting phase fraction is then used to update the bulk density and viscosity [using Eq. (4)] and to calculate the surface tension forces using the CSF model as described in Sec. II A. To limit oscillations in the surface tension force calculation and to suppress parasitic velocities, we have modified the CSF model of OpenFOAM to implement the smoothing of the phase fraction using Eq. (9) (Sec. II B). Finally, to solve for the pressure and the velocity field from the momentum equation [Eq. (2)] and the continuity equation (1), the Pressure Implicit with Splitting Operators (PISO) algorithm by Issa (1986) is used. For further details on the implementation of the PISO algorithm and the finite volume discretization of the governing equations in OpenFOAM, see Deshpande (2012), Berberović (2009), and Rusche (2003) and the supplementary material. For all simulations, the initial time step is set to a very low value (⁠ $∼10−8$ s), with a run-time adjustable time step scheme, which automatically adjusts the time steps based on the Courant–Friedrichs–Lewy (CFL) criterion (Rusche, 2003). The maximum CFL number for all the cases was set to a value of 0.5. All numerical simulations were conducted in parallel on the cluster system at Leibniz University of Hannover, Germany.

To validate the two-phase VOF model described in Sec. II, we consider the immiscible displacement of a wetting fluid, which initially saturates the flow cell, by a non-wetting fluid (i.e., primary drainage), in a long Hele–Shaw cell. This cell is thus a rectangular cuboid flow cell with one of its dimensions transverse to the mean flow that is very small compared to the other two dimensions. The Hele–Shaw cell (or channel) can be considered the limit model for a perfectly smooth fracture. While smooth-walled fractures are unrealistic, using Hele–Shaw cells offers valuable insights into processes that could occur within more complex fracture geometries (Meakin and Tartakovsky, 2009). In this primary drainage configuration, the displacing fluid has a lower viscosity than the displaced fluid, so the process controlling the dynamics of the fluid–fluid interface is viscous instability, which leads to the formation of a fingerlike pattern, first studied experimentally and theoretically in the seminal work by Saffman and Taylor (1958). The most important observation in their experiments was the emergence of a single dominant finger in the channel, whose width's ratio to the channel width, λ, is close to 0.5 unless the flow is very slow. At very low mean flow velocities, the measured values of λ diverge from 0.5, going to the limiting case of 1.0 as the mean velocity tends to zero.

The numerical setup for our validation case (shown in Fig. 1) is inspired by the original experimental configuration used by Saffman and Taylor (1958), with slight modifications on account of the computational effort and time, e.g., the width to thickness ratio for our case is $Ly/b=31.25$⁠, which is very close to the original experimental setup (31.75). The channel length is 10 cm, which is 8 times its width. We impose an initial perturbation on the fluid interface, as shown in Fig. 1, in order to suppress the growth of smaller fingers and accelerate the propagation of one single dominant finger in the channel. We designate the two fluids as fluid 1 and fluid 2, where fluid 1 is the invading, less viscous non-wetting fluid, while fluid 2 is the defending, more viscous fluid, wetting the channel walls. The material properties for the pair of fluids are listed in Table I.

Depending on the prescribed injection velocities of the invading fluid 1 at the inlet, $u(x=0,y,z)=(U,0,0)$⁠, and using the definitions of the Reynolds number as $Re=ρ1Ub/μ1$⁠, which quantifies the nonlinearity of the flow, and capillary number as $Ca=μ1U/σ$⁠, which is an estimate of the typical ratio of the viscous forces' magnitude to that of the capillary forces, we obtain a range of different flow configurations as listed in Table I. The viscosity ratio $M=μ1/μ2$ is equal to 1/30, which is the value used in Saffman and Taylor's experiments where water penetrates an oil of higher viscosity. The boundary conditions that are imposed on the different domain boundaries (Fig. 1) are listed in Table II.

The numerically induced slip presents a challenge for achieving a truly grid-independent solution with the VOF model. Therefore, in our grid convergence study, our aim is to reasonably predict the two-phase flow behavior while considering computational constraints. We conduct two sets of numerical simulations, which differ in how they handle the film of the defending fluid attached to the top and bottom walls of the Hele–Shaw channel. The first set of simulations, which we denote “resolved-film” (RF), resolves that film, while the other set, denoted “unresolved-film” (UF), does not. The corresponding computational meshes used for these two configurations are shown in Fig. 2. The RF mesh [Fig. 2(a)] has a much higher number of cells in the z direction, with a uniform grading near the walls, while the UF mesh [Fig. 2(b)] has cells of uniform thickness in the z direction, equal to their horizontal size. Previous studies (Horgue , 2012, 2013) indicate that to capture the wetting fluid film accurately, it should be covered by at least two grid cells, and that having more than two cells in the liquid film does not significantly affect solution accuracy. Our simulations, conducted with varying vertical resolutions in a reduced domain, confirm this. However, it is essential to note that the vertical resolution at the walls impacts the numerical slip of the contact line, which can influence other flow behaviors, such as the finger growth velocity.

The horizontal grid resolution is defined by the uniform horizontal mesh size $Δx=Δy$⁠, which is chosen as a fraction of the constant aperture b: $Δx=b/4, b/8$⁠, or $b/16$⁠, depending on the simulation runs. For both RF and UF cases, these three levels of horizontal discretization were investigated while maintaining a constant vertical resolution. Figure 3 illustrates the impact of the horizontal discretization on finger evolution for the UF case with Re = 0.06 and $log (Ca)=−2.60$⁠. For the coarse grid (⁠ $b/Δx=4$⁠), corresponding to the dashed curves in Fig. 3, insufficient grid resolution leads to numerical inaccuracies, preventing the formation of a stable finger. As the grid is refined to $b/Δx=8$ and $b/Δx=12$⁠, a stable dominant finger evolution is observed (solid curves in Fig. 3). Similar findings regarding the impact of the horizontal grid resolution were obtained for the RF simulations (not shown here). Additionally, pressure drops $ΔPio=Pin−Pout$ across medium's length were recorded at fixed times. The pressure drops for $b/Δx=8$ and $b/Δx=12$ differed by less than 2%. Similar trends were observed for the RF simulations. Consequently, we have selected a horizontal discretization level of $b/Δx=8$ for all our simulations.

We shall now compare the results of numerical simulations conducted for the immiscible two-phase displacement in the long Hele–Shaw channel. We specifically focus on two aspects: (i) comparing the numerical results of both RF and UF cases with the experimental observations and theoretical findings of Saffman and Taylor (1958) and (ii) examining the differences between the results of the RF and UF cases, particularly in predicting the macroscopic flow properties.

We first study the differences in invasion dynamics between the RF and UF configurations. Naturally, it is anticipated that the RF simulations are expected to better describe the physics at play and thus better predict the experimental findings, but they are more costly computationally. Hence, we want to test the comparative predictive efficiency of UF simulations.

Figure 4 depicts the initial stage of the fluid–fluid interface's displacement and highlights a crucial difference between the RF and UF predictions concerning the dynamics of the contact line at the top and bottom walls. In the RF case (solid curves in Fig. 4), where refined grid cells are employed at the walls, the contact line hardly moves, resulting in the formation of a distinct liquid film of the displaced wetting fluid on the top and bottom walls. This is what is observed experimentally. Conversely, in the UF case (dashed curves in Fig. 4), the much coarser vertical grid resolution in the vicinity of the walls causes the contact line to slip, preventing film development due to the fact that the vertical resolution is too coarse for the film to be resolved by the numerical simulation.

Figure 5 displays the fluid–fluid interface in vertical cross sections of the channel, both longitudinal [Fig. 5(a)] and transverse [Fig. 5(b)] to the mean flow direction, for three Ca values: high, low, and intermediate (Table I). The defending fluid's film thickness decreases with lower capillary numbers (Ca), which is consistent with prior numerical and experimental studies (Bretherton, 1961; Horgue , 2012; and Rosengarten , 2006). Notably, when examining the fluid–fluid interface in the transverse vertical plane [Fig. 5(b)], upward bulges become apparent at the four corners, alongside a comparatively flatter interface at the center. Previous experimental and theoretical studies on viscous fingering in conventional Hele–Shaw cells (Tabeling and Libchaber, 1986; Jensen , 1987; Tabeling , 1987; Reinelt, 1987; and Tanveer, 2000) have well acknowledged the presence of such a thin wetting film and its distinct transverse shape. It is important to note that this detail is specifically captured in the RF configuration, which further strengthens our model validation efforts. However, a detailed analysis of the transverse shape and variable film thickness falls outside the scope of our current study. For readers interested in a more in-depth examination, the work of De Lozar (2008) may provide valuable insights.

The presence of a film of defending fluid on the horizontal walls reduces the effective aperture available to the invading fluid finger for flowing through the channel. Therefore, for a given injection flow rate, the mean velocity of the displacing fluid is necessarily larger for a RF simulation than for the corresponding UF simulation. This is reflected in the breakthrough time (⁠ $t*$⁠), i.e., the time at which the displacing fluid reaches the outlet. In Fig. 6, we show the relative difference in breakthrough times between UF and RF simulations across the entire range of the investigated Ca values. RF simulations consistently predict lower breakthrough times compared to UF simulations. Importantly, as Ca decreases, the relative difference in breakthrough times diminishes, dropping to less than 5% for the lowest Ca values.

Figure 7 illustrates the time dynamics of the mid-plane cross section of the fluid–fluid interface from t = 0 until $∼0.9t*$⁠, when the Hele–Shaw cell is initially filled by the wetting fluid, for both RF [Fig. 7(a)] and UF [Fig. 7(b)] simulations, with Re = 0.016 and $log (Ca)=−3.18$⁠. Both simulations exhibit a behavior consistent with experimental observations, with a singular stable finger of displacing fluid that exhibits a uniform width as it advances at constant velocity along the channel. Expectedly, the different color scales in Figs. 7(a) and 7(b) show that the finger in the RF simulation reaches the outlet earlier than that in the UF simulation. Furthermore, the differing shapes of the finger constrictions near the inlet illustrate the influence of the film and of the corresponding contact line motion. For the cases from Table I that are not shown here, similar observations are reported.

In order to further assess the reliability of these results, we compare the width and tip shape of the simulated fingers with Saffman and Taylor's (1958) experimental observations of the finger width and theoretical prediction of the finger shapes for fingers of water displacing two different oils (Diala and Talpa). To determine the widths of our simulated fingers, we first obtain the intersection curve between the horizontal mid-plane of the Hele–Shaw cell and the fluid–fluid interface at time $0.9t*$⁠. The finger width is subsequently obtained as the mean width of the segment of the curve located in the region $x/Lx∈[0.6,0.8]$⁠.

We characterize the finger width by λ, which is its ratio to the Hele–Shaw channel width. In Fig. 8, we plot λ as a function of $μ2Uf/σ$⁠, as presented by Saffman and Taylor (1958), where $Uf$ is the finger-tip velocity of a stable finger. The black data points represent the historical experimental measurements involving direct and photographic assessments of finger positions and profiles, while the colored ones represent our simulation results. The colored data points in Fig. 8(b), correspond to our simulations, conducted at different levels of vertical discretization while maintaining the same horizontal grid refinement level of $(b/Δx=8)$⁠. As mentioned earlier in Sec. III B, vertical refinement of the RF simulations' numerical mesh has allowed us to reach grid convergence. All filled symbols indicate simulations for which we consider that grid convergence has been reached. At some values of Ca, unfilled symbols are shown, closely clustered together and with the filled symbol; they correspond to successive steps of the grid refinement procedure. Our investigations show that for the lowest Ca case, where the film thickness is approximately 3 μm, a grid distribution consisting of 40 cells in the vertical direction, a wall-adjacent cell thickness of 0.33 μm, and a cell-to-cell expansion ratio of 1.2 provide a sufficient resolution of the thin film. In contrast, the UF cases require a uniform vertical cell thickness, using the three levels of vertical discretization: $b/Δz=4$⁠, 8, and 16.

The RF numerical results closely match the experimental data across the entire range of investigated Ca, which points to an excellent validation of our OpenFOAM implementation of the VOF-based two-phase flow model. In comparison, the UF simulations tend to overpredict the finger width (λ) and accordingly yield lower finger velocity magnitudes (⁠ $Uf$⁠) compared to both experimental data and the RF simulations [Fig. 8(a)]. These discrepancies are more noticeable when examining Fig. 8(b), where λ values are plotted against Ca, defined using the inlet velocity magnitude U. As expected, these differences diminish as Ca decreases. At higher Ca values, where a relatively thicker defending fluid film exists, the limited vertical resolution in the UF simulations fails to capture film effects. In the RF simulations, the contact line is pinned to the walls from the start of the injection [Fig. 4(a)], while the UF configuration shows contact line slippage [Fig. 4(b)], as mentioned in Sec. III C 1. Although minor in comparison to the longitudinal direction, contact line slippage also occurs in the transverse direction in the UF simulations. Consequently, with a given imposed flow rate, lateral contact line motion in UF simulation leads to an overestimation of λ, further resulting in a reduced finger velocity. As the wetting film thickness decreases with decreasing Ca, the film's influence diminishes, enhancing the agreement between UF and RF simulations.

The VOF model described in Sec. II and validated in Sec. III is now used to conduct numerical experiments on rough fracture geometries. Our objective is to compare the results from RF simulations to those from the corresponding UF simulations, aiming to elucidate how film resolution influences predictions of various flow-related quantities during drainage in a geological fracture. It allows us to quantify the prediction errors arising from the absence of resolved films. We aim to identify situations where such resolution is essential and those where, despite the lack of film resolution, predictions largely remain reliable. Additionally, this study highlights the trade-off between accuracy and computational expenses, guiding our decisions in selecting film-resolution approaches and strategies for drainage simulations within rough fractures.

Figures 11 and 12 show the schematics of the two numerically generated rough fractures that were used for the numerical simulations of drainage. The first fracture (Fig. 11) has an aperture field that varies sinusoidally (and is thus periodic); in the following, we shall refer to it as the sinusoidal fracture. The gradient of its aperture field is moderate (⁠ $||∇a(x,y)||≪1$⁠). It is important to note that no actual fracture, regardless of the fracture material, will exhibit a geometry similar to this one. This fracture geometry serves as a transitional case of increased complexity from the parallel plate (Hele–Shaw) setup used as our validation (Sec. III) case to the realistic geometry shown in Fig. 12. This second fracture exhibits wall topographies that seem to be distributed at random (in particular, they are not periodic), but the fracture has statistical geometrical properties that are consistent with those of naturally occurring geological fractures. Its aperture field is also rougher at small scales (⁠ $max(||∇a(x,y)||) ⪆ 1$⁠) than that of the sinusoidal fracture. Hereafter, we shall refer to that fracture as the geological fracture.

Both fractures consist of a top and a bottom wall whose horizontal mean planes are parallel and separated by a distance $a¯$⁠, which is thus the mean aperture of the fracture (the walls are never in contact). In the case of the sinusoidal fracture, the bottom surface coincides with the plane z = 0, while the top surface is represented by a two-dimensional sinusoidal function (see Fig. 11). On the other hand, the top and bottom walls of the geological fracture are self-affine isotropic topographies, which exhibit long-range correlations (Bouchaud, 1997; Plouraboué , 1995) characterized by the Hurst exponent H (here H = 0.8). The two walls have identical large-scale variations above a characteristic length $Lc$⁠, which we denote correlation length. In addition to (i) $a¯$⁠, (ii) H, and (iii) $Lc$⁠, the aperture field of such geological fractures is defined primarily by (iv) its standard deviation σ_a and (iv) the horizontal dimensions L_x and L_y of the fracture (Méheust and Schmittbuhl, 2003; Dewangan , 2022). For details on the geometrical properties of such geological fractures, see also Lenci (2022b) and Dewangan (2022). The fracture wall topographies of Fig. 12 are generated using a spectral method initially proposed by Méheust and Schmittbuhl (2003). The generation of the correlated random fields corresponding to the wall topographies is done by first generating two independent self-affine walls of zero mean and identical standard deviation, $h+$ and $h−$⁠; then attributing the same large Fourier modes (for wave number larger than $2π/Lc$⁠), taken from $h+$⁠, to $h−$ (or conversely); then imposing to both a new standard deviation such that the $h+−h−$ has the desired standard deviation σ_a; and finally adding $a¯/2$ to $h+$ and subtracting $a¯/2$ from $h−, a¯$ being the desired mean aperture.

For both fracture geometries, we simulate the same drainage configuration, described in Sec. III, where a less viscous fluid 1 displaces a fluid 2 of higher viscosity (see Table III). To accurately capture the onset of instability of the fluid–fluid interface, we initialize the interface (flat) up to a small distance of 0.0025 m  $=Lx/20$ from the inlet, meaning at t = 0 fluid 1 occupies a small volume of the fracture domain. Here we use the same nomenclature for the different domain boundaries (see Figs. 11 and 12) as used in Sec. III A for the Hele–Shaw channel, and thus the boundary conditions for our simulation setup are the same as prescribed in Table II. We investigate five different cases of Ca which span over three orders of magnitude: $log10Ca=−3.0,−3.5,−4.0,−4.5,−5.0$⁠. Table III lists the inlet velocities $U=σCa/μ1$ obtained for these five cases and the resulting $Re=ρ1Uina¯in/μ1$ (we have used inlet mean aperture $a¯in$ to calculate Re) for both the sinusoidal and geological fractures.

We again conduct two sets of simulations, namely with resolved-film (RF) and unresolved-film (UF) for both the sinusoidal and geological fracture geometries. The computational mesh employed for our numerical calculations is shown in Fig. 13, where, similarly to our validation on Hele–Shaw geometry (Sec. III), the RF simulations have a substantially higher degree of resolution in the vertical direction than the UF simulations. To ensure grid independence in our numerical simulations, we used our findings from the Hele–Shaw validation case (Sec. III B) and chose three horizontal discretization levels: $Δx=100$ μm, $Δx=50$ μm, and $Δx=25$ μm. We then performed grid sensitivity analyses on a smaller domain of the fracture geometries while maintaining a fixed vertical resolution. For both geometries and the corresponding RF and UF configurations, we found that the difference between the results corresponding to the discretization levels of $Δx=50$ μm and $Δx=25$ μm (e.g., average pressure drops $ΔPio$⁠) were less than 1%–2%. Consequently, for all the simulations, we have chosen a horizontal discretization corresponding to $Δx=50$ μm, which results in 1000 × 1000 cells in the horizontal plane of the fracture geometries.

As previously indicated in Secs. III B and III C, ensuring convergence, particularly in the vertical resolution, poses challenges. Furthermore, with rough fractures, whose apertures vary spatially, excessively large grid cell numbers in the z direction are impractical due to the emergence of tiny volume cells in small aperture regions. For RF simulations, numerical experiments on reduced domains with varying numbers of cells in the z direction suggested that for both geometries, 32 uniformly graded cells (1.33 cell-to-cell expansion ratio) sufficed to resolve wall-adjacent films, even at the smallest capillary numbers. This arrangement results in mean cell thicknesses of 0.33 and 0.39 μm for the sinusoidal and geological fractures, respectively. Conversely, considering the dependence of wetting-film thickness on the Ca, we conducted analogous UF numerical experiments (i.e., without vertical grading in the z direction). For higher Ca values [ $log(Ca)≥−3.5$], we used eight cells; for lower Ca cases [ $log(Ca)<−3.5$], 16 cells were considered.

The fluid displacement morphology, commonly referred to as invasion patterns, determines the macroscopic residual saturations and thus, for example, the sweep efficiency. It is controlled by the interplay between viscous and capillary forces, and possibly gravity, as well as the local aperture conditions. When the interface is viscously unstable (M < 1), this interplay results in fingerlike patterns. In porous media, if capillary forces are dominant (low Ca), the so-called capillary fingering occurs, which features patterns resulting from growth occurring at various places along the fluid–fluid interfaces, depending on where the capillary forces are more favorable to the interfaces' displacement (Lenormand and Zarcone, 1989). Conversely, at high Ca values, viscous fingering prevails, characterized by more branched fingers (Måløy 1985). A similar behavior is observed in rough fractures (Chen , 2017), but the phase diagrams of the flow patterns as a function of Ca and M are not necessarily identical to those known for porous media (Lenormand, 1990) and remain an open question. Figures 14 and 15 depict the displacement morphologies for the sinusoidal and geological fractures for different capillary numbers. The images display superimposed invasion patterns, using distinct colors to show agreement (yellow) or discrepancy (blue and red) between patterns from RF and UF simulations. In line with what was observed for the Hele–Shaw configuration (Sec. III), finger advancement time differs between film and without-film cases. To emphasize the patterns, the snapshots for RF and UF simulations are taken at different time steps to maintain approximately equidistant positions from the outlet.

For the sinusoidal fracture case (Fig. 14), an immediate observation is the presence of longitudinal preferential paths of advancement in the invasion front's progression at all capillary numbers. These pathways form because the initial invasion is easier along these paths as they feature a large aperture zone in contact with the fracture inlet. Once these high aperture regions have been invaded, parallel fingers have been initiated, which will then retain an advantage and thus prevent the fluid from invading other areas in the vicinity of the inlet. Hence, only these fingers grow. Due to the viscous instability of the interface, a finger that gets ahead of another neighboring one at a given time will keep this advantage at later times, so that the differences in finger propagation increase in time. For the geological fracture (Fig. 15), where aperture variations are no longer periodic but disordered, with a correlation length of the apertures that is a quarter of the fracture size, the same phenomenon will lead to the formation of a smaller number of disordered fingers exhibiting irregular patterns resembling those seen in natural fractures or their laboratory replicas, as reported in previous studies (e.g., Chen , 2017, 2018).

For both fracture geometries, a noticeable broadening of the fingers in the transverse direction is seen as the Ca decreases. Additionally, when the Ca is decreased, the capillary forces start impacting the invasion process more strongly, so that increased local surface tension effects lead to stronger shielding of the least advanced fingers' growth by the most advanced fingers, thus reducing the number of fingers that proceed effectively toward the outlet. These features are reproduced with both RF and UF simulations. The enhanced shielding mechanism at low Ca is more prominent in the sinusoidal fracture, due to the fact that the various fingers encounter the same geometries as they travel forward, so that no stochastic geometry heterogeneity competes with the shielding mechanism to control the differential propagation of fingers. In this geometry, and at the two lower Ca values of $log Ca=−4.5$ and –5.0, a considerable disagreement is thus observed between UF and RF invasion patterns, as UF simulations lack the ability to capture the full 3D effects of dominant surface tension forces. This effect is much less present in the geological fracture, due to the stochastic heterogeneity of the aperture field, as discussed above.

At higher Ca values [ $log (Ca)>−4.0$], discrepancies between the invasion patterns predominantly appear in the figures as extended blue regions (RF) aligned with the flow direction, for both geometries. These differences tend to become more pronounced in the later stages of the invasion process [e.g., Figs. 14(a″) and 15(a″)]. The relatively thick wetting film at high Ca limits the effective aperture available for the invading fluid, accelerating the movement of the advancing fluid fingertips. The widening of the fingers in the transverse direction, though present in both RF and UF simulation results across all Ca cases, is predominately more marked in the UF cases. This is consistent with our observations in the case of a single finger propagating in the Hele–Shaw channel, where UF simulations yielded finger widths larger than those of the RF simulations for the entire range of Ca. Analogously, the lateral slip of contact lines in the UF simulations translates to slightly broader finger widths compared to the corresponding RF simulations. This is also depicted in Fig. 16, illustrating the fluid–fluid interface in small cross sections of the two fracture geometries.

In the sinusoidal fracture at the two lower Ca values [ $log (Ca)=−4.5$ in Fig. 14(d-d″) and $log (Ca)=−5.0$ in Fig. 14(e-e″)], we also observe cutting-off of the invading fluid, both in the RF and UF simulations, by pinch-off in low aperture areas. As shown in Fig. 17 for $log (Ca)=−4.5$⁠, this occurs in fingers that are left behind and is more frequent and quicker in the RF simulations. In the latter configurations, such pinch-off events are associated with later re-merging of these isolated clusters with the base of the short finger from which they had separated; this is not observed in UF simulations. Further details are provided in the  Appendix, with additional pictures of the flow patterns for each of the flow configurations, and links to online videos of the corresponding experiments (Multimedia available online).

The overall relative deterministic discrepancy between the UF and RF simulations can be quantified from the ratio $δA/A$ of the joint area of the blue and red regions in Figs. 14 and 15 (for the sinusoidal and geological fractures, respectively) to the area $A$ of the region of the fracture plane that is occupied by the displacing fluid in the RF simulation (i.e., the joint area of the blue and yellow regions in Figs. 14 and 15). This ratio is plotted at breakthrough time in Fig. 17(d), as a function of the capillary number and for the two types of fracture geometries. Its dependence on Ca is well explained by the phenomenology discussed above. In particular, at the lowest capillary number (⁠ $log Ca=−5$⁠), the aforementioned shielding of less advanced fingers, coupled with the marked effect of capillary forces, leads to fingers partially developing at different places in the fracture plane in the UF simulations, as compared to the RF simulations, and thus to $δA/A$ values between 20% and 30%. At the largest capillary number (⁠ $log Ca=−5$⁠), on the contrary, the interface's enhanced instability leads to a similar discrepancy and $δA/A$ values as high as 50% for the geological fracture. At intermediate Ca values, for which none of these two effects is strong, the $δA/A$ values are typically close to 10%. Note that deterministic discrepancies between the two displacement patterns that are larger than 10%, essentially result from finger branches developing at different places in the UF and RF simulations, which is expected due to the unstable nature of the fluid–fluid interface and the relatively ordered geometries. We shall see that comparing quantities that correspond to more statistical measures of the flow behavior (e.g., breakthrough times, fluid–fluid interface length, longitudinal pressure profiles, longitudinal saturation profiles) will yield much smaller discrepancies between the invasion patterns predicted by the UF and RF simulations. This was also observed, for similar reasons, by Ferrari (2015) when comparing numerical predictions to experimental measurements of two-phase flow in random two-dimensional porous media.

For the two fracture geometries, an essential observation emerges from the RF simulations when examining the respective film thicknesses. We estimate the film thickness at a given position along the fracture plane in the following manner. We first depth-average the indicator function $γ(x,y,z)$ to obtain a two-dimensional saturation field $s1(x,y)$ for the invading fluid 1, and thence the corresponding depth-averaged saturation of the defending fluid film, $s2(x,y)=1−s1(x,y)$⁠. This 2D saturation s₂ is either 1, in the portion of the fracture plane where no invading fluid is present, or significantly smaller than 1 in a region that we denote Ω. The value of s₂ in Ω results from the presence of the two films of displaced fluid (on each of the fracture walls), and the average film thickness can be estimated as the average of $(s2×a)/2$ over Ω, which we write as $⟨s2 a⟩Ω$⁠.

Note that this estimate of the mean film thickness tends to slightly decrease in time. Its average over the three different times at which the invasion patterns of Figs. 14 and 15 are shown is plotted in Fig. 18 as a function of the capillary number for the two fracture geometries. For both geometries, it exhibits a power law behavior close to $Ca1/3$⁠, with values ranging from 2.5 to 25 μm. In the figure, we also indicate the Bretherton law for wetting film thickness measured when a gas displaces a liquid (Bretherton, 1961); it is not expected to hold here due to the much lower viscosity ratio between the two fluids. Across all Ca, the sinusoidal fracture consistently exhibits larger film thickness than the geological fracture, but the thickness ratio between the two geometries is overall all the smaller as Ca is larger, reaching 1 at $Ca=10−3$⁠. At $Ca=10−5$⁠, it is larger than 2. A detailed study of the films' temporal and spatial fluctuation is out of the scope of the present study.

Two key observations emerge when examining the relative differences between breakthrough times $t*$ in the UF simulations and RF simulations, $Δt*=tUF*−tRF*$⁠. They are plotted against $log (Ca)$ values in Fig. 19 for both fracture geometries. First, $Δt*$ is always positive, which means that the UF simulations always overpredict the breakthrough times, i.e., underpredict the finger-tip propagation. Second, decreasing Ca leads to diminishing the relative differences in $t*$ for both fracture geometries, which aligns with our Hele–Shaw cell validation findings (Sec. III). Third, the relative differences in $t*$ for the sinusoidal fracture geometry, while comparable to those obtained for the Hele–Shaw cell case (see Fig. 6), are considerably larger (nearly one order of magnitude) than those for the geological fracture. This cannot be only explained by the differences in wetting film thickness presented in Sec. IV D 2. For the sinusoidal fracture in particular, the widening of the invading fluid pattern plays an important role in slowing down finger-tip propagation. In addition, not resolving the wetting films leads to the predominant propagation of a different finger in Fig. 17(a), which can result in a significant discrepancy in the breakthrough times, whereas for the geological fracture, the finger that propagates faster is the same in the UF and RF simulations.

For the evolution of fluid pressure, we examine the variation of the area-weighted average pressure drop $ΔPio$ across the inlet and outlet of the fracture geometries with the invading fluid saturation S₁ (Fig. 20). Higher capillary numbers result in elevated and relatively steady pressure drops due to increased viscous forces. At lower Ca values, where capillary forces become dominant, considerable fluctuations in pressure drops are seen. Qualitatively, these fluctuations also mark the transition from a viscous-dominated flow to a more capillary-dominated invasion process. Similar behavior of the pressure drops has been reported in previous studies on drainage in both porous media (Ferrari and Lunati, 2013) and rough fractures (Chen , 2018; Yang , 2019).

For the pressure drops, a generally good agreement is observed between those obtained from RF and UF simulations, particularly in the case of the sinusoidal fracture. This improved match is primarily attributed to a thinner film of defending fluid in the sinusoidal fracture geometry. At $log (Ca)=−3.0$⁠, when the film thickness is the largest, the pressure drops obtained in both fractures from the UF simulations are considerably lower than those obtained from the RF simulations. Another noteworthy observation is the larger $ΔPio$ values in the sinusoidal fracture compared to the geological fracture. In addition to the mean aperture of the geological fracture being 1.3 times larger than that of the sinusoidal fracture, the larger quantity of wetting film in the latter case further reduces the effective aperture available to the invading fluid. Consequently, higher pressure drops are obtained in the sinusoidal fracture under the same inlet flow conditions.

We conduct another comparison of the invasion morphologies by computing the length L of the fluid–fluid interface within the mean aperture plane (xy) of the fracture geometries. For this, the fluid distributions are mapped onto a constant-thickness z grid, and the interface length is obtained as the intersection length between the interface and the mid-cell plane of this grid. Starting with a flat interface $(L(t=0)=Ly)$ and some initial saturation of the invading fluid $S1(t=0)$⁠, we are interested in the evolution of the interfacial length as a function of the saturation S₁. Figure 21 depicts this evolution for the two fracture geometries across three representative Ca values spanning from highest to lowest. A relatively linear increase in L with saturation is observed across all cases. Higher Ca values correspond to a pattern close to viscous fingering, leading to a larger L. These trends in L are captured well by the RF simulations and are consistent with prior studies (Chen , 2018; Ferrari , 2015). Interface lengths predicted by the UF simulations match this behavior only at low capillary numbers when the influence of the wetting film's thickness diminishes. Notably, in the geological fracture at the highest Ca [ $log (Ca)=−3.0$], a rapid growth in L is seen [Fig. 21(b)]. At large Ca, breakup events dominate the invasion morphology, leading to an increase in L (Pak , 2015; Chen , 2018). At low Ca, cutoff and trapping of the fluid cluster also impact the evolution of L. This is seen in the sinusoidal fracture at $log (Ca)=−5.0$ [Fig. 21(a)], where trapping occurs early in the invasion process [see Fig. 14(e-e″)], temporarily increasing the interfacial length. For further details, see the  Appendix, which provides additional pictures of the flow patterns for each of the flow configurations and links to online videos of the corresponding simulations (Multimedia available online).

Figure 22 shows the interface tip velocities $Utip$ plotted as a function of the invading fluid saturation S₁ for the two geometries. The tip velocities are the velocities of the most advanced interface point. With decreasing capillary numbers (Ca) and injection velocities (U), the average tip velocities also decrease. In the sinusoidal fracture, $Utip$ exhibits a quasi-periodic behavior, while in the geological fracture, unsurprisingly, its spatial variations are irregular. Additionally, higher fluctuations in $Utip$ can be noticed in the case of the sinusoidal fracture. At low Ca values [ $log (Ca)=−4.5,−5.0$], strong fluctuations indicative of increased capillary fingering appear in the tip velocities. In the non-smooth fracture, a clear crossover is observed between $log(Ca)=−4.0$ and –4.5, signifying a comparable influence of capillary and viscous forces (Chen , 2017). This behavior is well captured by the UF simulations, although not quite as well as by the RF simulations.

The aforementioned parameters, such as interfacial lengths and average pressure drops, while revealing insightful flow behaviors, do not encompass the entire information regarding the overall fluid distributions. For instance, statistically different fracture geometries might exhibit contrasting invasion patterns yet exhibit the same temporal evolution of their interface lengths (Ferrari , 2015). To complement the analysis, we thus examine longitudinal saturation profiles of the invading fluid, obtained as averages of the indicator function across the transverse (y) direction. In Fig. 23, these profiles are presented as functions of the normalized longitudinal coordinate, $x/Lx$⁠, for the two fracture geometries. It is noteworthy that the longitudinal saturation profile also shares a close connection with the Darcy scale depiction of flow through permeable media, which is based on the method of volume-averaging (see, e.g., Bear, 2013). Regarding the observed distinctions in saturation profiles between the RF and UF simulations, these dissimilarities can be attributed to one of the following factors: (i) the presence (or absence) of the defending fluid film, which could lead to a reduction (or increase) in the saturation of the invading fluid; (ii) the lateral extent occupied by the advancing fingers; or (iii) a mismatch between the resultant invasion patterns in the two scenarios, as observed in the sinusoidal fracture at low Ca. At high Ca [ $log (Ca)=−3.0$], when the wetting film thickness is substantial, for the sinusoidal fracture case [Fig. 23(a)], the saturation profiles from the UF simulations exhibit higher saturation levels at all locations (x), but with larger differences at longitudinal positions where the average aperture is larger. Conversely, at the same Ca value for the geological fracture [Fig. 23(d)], the dissimilarities in the saturation profiles primarily arise due to the mismatch in the branched invasion patterns [Fig. 15(a″)] and are smaller. For both fracture geometries, at intermediate Ca values [ $log (Ca)=−4.0$⁠, Figs. 23(b) and 23(e)], when the influence of the wetting film is moderate and the observed alignment between the corresponding invasion patterns from the RF and UF simulations is most prominent [Figs. 14(c″) and 15(c″)], the saturation profiles also exhibit the highest level of agreement. At the lowest Ca value [ $log (Ca)=−5.0$⁠, Fig. 23(c)], the disparities in saturation profiles for the sinusoidal fracture scenario are attributed to the pronounced mismatch between the invasion morphologies resulting from the two film resolution configurations [refer to Fig. 14(e″)]. On the contrary, for the geological fracture, although at $log (Ca)=−5.0$ the patterns exhibit considerable resemblance at large scales, these disagreements predominantly arise due to small-scale discrepancies, where the advancing fluid front becomes more compact and occupies greater aperture space in the transverse direction [see Fig. 15(e″)].

We have investigated to which extent resolving the film left behind by the defending fluid on the fracture walls or not influences the immiscible fluid–fluid displacements in rough fracture geometries, and the impact of not resolving the film on the predicted flow behavior. The grid resolution employed in these simulations plays a significant role in capturing the intricate dynamics of fluid flow within the fractures. Thus, it is important to recognize the distinct computational demands imposed by each of the two approaches, namely, the RF and the UF. As stated in Secs. III B and IV C, the RF simulations demand a more refined grid structure compared to their UF counterparts. Specifically, the RF simulations involve nearly twice the number of cells in the vertical z direction. Moreover, the grid resolution is vertically graded, resulting in higher cell density near the top and bottom walls.

The computational time required for these simulations is proportional to the number of cells in the grid. Consequently, RF simulations demand significantly more computational time due to their higher cell count. Additionally, the computational effort also scales with the number of time steps and the time step size, which in turn is constrained by the CFL stability criterion. For such flow cases where surface tension forces are relevant and involve complex fluid interfaces, to maintain numerical stability and mitigate the occurrence of parasitic currents, a CFL value of 0.5 is typically imposed (Deshpande , 2012). The CFL criterion necessitates smaller time steps, which implies that as grid resolution increases, requiring smaller cell sizes, the time step size decreases accordingly. Thus, the increased cell count and the reduction in time step sizes compound the total computational effort required for the RF simulations. For instance, in the case of the sinusoidal fracture geometry at the lowest capillary number [ $log (Ca)=−5.0$], the RF simulation employed a total cell count of $3.2×107$⁠, twice that of the UF simulation. Given a CFL criterion of 0.5, this led to time steps approximately an order of magnitude smaller, around $∼2×10−5$ s, compared to those in the UF counterpart. While the UF simulation utilized 120 processors and demanded 927.62 CPU-hours, the RF simulation ran on 250 processors, requiring 14 461.04 CPU-hours, i.e., 15 times the CPU-hours needed for the UF simulation. This stark difference highlights the substantial increase in computational demand associated with resolving the defending fluid films. It underscores the trade-off between computational cost and film resolution accuracy, a crucial consideration when choosing the simulation strategy in rough fracture drainage studies.

A challenging yet essential consideration in hydrodynamic-scale modeling of drainage within open rough fractures is the accurate resolution of the film of defending fluid, which adheres to the fracture walls. This study employs a VOF-based CFD approach, using a custom-modified version of the open-source code OpenFOAM, to evaluate the impact of film resolution on hydrodynamic-scale drainage simulations.

We have presented a detailed validation of the numerical model utilizing the classical case of viscous instability within a conventional Hele–Shaw geometry. Depending on whether the computation mesh resolves the defending fluid film or not, the results from two sets of simulations; “resolved-film” (RF) or “unresolved-film” (UF), were compared with the experimental and analytical results of Saffman and Taylor (1958). This comparison revealed that the RF simulations demonstrated excellent agreement with experimental results across the entire range of capillary numbers (Ca), while UF simulations at higher Ca resulted in slower finger propagation velocities and increased finger widths. With smaller Ca values (⁠ $Ca≤10−4$⁠) and smaller film thickness, UF results aligned more closely with the experimental and RF simulation findings.

The investigation was then extended to two numerically generated rough fracture geometries, where the first has a sinusoidally varying aperture and the other mimics a naturally occurring geological fracture with a stochastic self-affine roughness of the fracture walls and correlations between the walls' topographies at large scales. We conducted drainage simulations on these fractures, employing both RF and UF simulations and covering a broad range of Ca which spanned over three orders of magnitude (⁠ $10−5≤Ca≤10−3$⁠). To assess the effect of film resolution, we performed an extensive comparison of hydrodynamic-scale and macroscopic quantities of interest between the UF and RF simulations. Consequently, our analysis leads us to the following findings:

• In the sinusoidal fracture case, throughout all capillary numbers, the observed invasion patterns had distinct preferential flow paths along channels with overall larger apertures, leading to flow-aligned fingers. Unsurprisingly, the geological fracture, on the other hand, displayed more irregular invasion patterns, consistent with the stochastic geometry. In this case, the displacement patterns transition from more to less flow-aligned as the Ca decreases.

• At higher capillary numbers (⁠ $Ca>10−4$⁠) with a considerable film of the defending fluid present on the walls, overall invasion patterns are similar, but the UF simulations notably underestimate other hydrodynamic-scale predictions such as interfacial lengths and macroscopic pressure drops, but overpredict invading fluid saturations and the breakthrough times.

• At very low capillary numbers (⁠ $Ca=10−5$⁠), while most hydrodynamic-scale predictions from the UF simulations align well with those from the RF simulations, the detailed invasion patterns, especially at smaller scales, can vary.

• Regarding macroscopic flow variables such as the average pressure drop and saturation, differences do exist between the RF and UF simulations, especially at higher Ca. However, not resolving the film only results in serious mispredictions at large Ca (⁠ $Ca≥10−3$⁠). At low Ca (⁠ $Ca<10−4$⁠), these differences can be considered negligible.

• The more abrupt spatial variations in wall topographies in the geological fracture lead to a reduction in wetting film thickness as compared to the sinusoidal fracture, resulting in less pronounced film effects and thus a better match between the UF simulations' predictions and those from the RF simulations.

• During the invasion process, the thin film on the fracture walls reduces the effective aperture perceived by the invading fluid. This film thickness is not uniform, even in the ideal parallel plate geometry. One potential approach to partially accounting for the film's influence without resolving it is to adjust macroscopic parameters such as the effective apertures available to the invading fluid. However, it is important to note that this cannot fully rectify all discrepancies. This limitation arises from the inherent variability in film thicknesses, which are not known a priori, and the broadening of the fingers in the absence of film resolution.

In summary, the choice of film resolution in modeling two-phase flow in rough fractures involves a trade-off between solution accuracy and computational resources. Our study demonstrates that accurately resolving the film can require up to tenfold more computational time compared to simulations where the wetting film is not resolved. Resolved-film simulations differ from unresolved-film simulations; however, these differences become smaller at small capillary numbers. In any case, taking these considerations into account, our study demonstrates that the employed VOF model is reliable and has the potential to serve as a valuable complement to experimental observations. Moreover, it enables exploration across a wide range of parameters and scenarios that may be difficult to emulate with experiments and field observations.

Prospects to this study are numerous and include the detailed study of the wetting film's heterogeneity in space and time and of its dependence on Ca, M, and the fracture geometry, as well as the characterization of flow regimes in the (Ca, M, Bo) parameter space (Bo being the Bond number).

See the supplementary material which presents a detailed description of the VOF model employed in this study. Readers interested in detail of the governing equations, discretization, and solution procedure implemented in OpenFOAM are encouraged to see the references included therein.

We gratefully acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Project No. NE 824/20-1 and the Agence Nationale de la Recherche (ANR, French National Research Agency) under Project No. ANR-20-CE92-0026 for the DFG-ANR project “2PhlowFrac.”

The authors have no conflicts to disclose.

R. Krishna: Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Y. Méheust: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (supporting); Methodology (equal); Project administration (supporting); Resources (supporting); Software (supporting); Supervision (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal). I. Neuweiler: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (supporting); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Validation (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The following figures show the snapshots of the fluid–fluid invasion patterns in the RF (invading fluid shown in blue) and UF (invading fluid shown in red) simulations. The light blue color indicates the defending fluid. The corresponding multimedia video files are available online.

• The video linked to Fig. 24 (Multimedia available online) shows the momentary pinch-off and merging of the invading phase seen at high Ca [ $log (Ca)=−3.0$] in the geological fracture when the film is resolved (RF simulation). These phenomena are not observed in the corresponding UF simulation, as shown in Fig. 25 (Multimedia available online).

• The video linked to Fig. 26 (Multimedia available online) shows the progressive pinch-off (isolated blobs of invading fluid) of the invading phase in the RF simulation of the displacement in the sinusoidal fracture for $log (Ca)=−4.5$⁠. However, as shown in Fig. 27 (Multimedia available online), this phenomenon is completely absent in the UF configuration. A similar observation is made for the sinusoidal fracture at $log (Ca)=−5.0$ where the RF simulation as shown in Fig. 28 (Multimedia available online), in addition to pinch-off, also leads to merging of two fingers of the invading phase. In the corresponding UF simulation shown in Fig. 29 (Multimedia available online), this pinch-off and the resulting isolation of the invading phase are present, but to a much lesser extent.