Depth-integrated Two-dimensional Model for Immiscible Two-phase Flow in Open Rough-walled Fractures with Smoothly Varying Aperture

Two-phase flow in geological fractures holds significant relevance in various applications, including subsurface fluid storage and oil and gas exploitation. The pore-scale modelling of such flows is a challenging task influenced by many factors, such as the complex interplay between the viscous, capillary, gravitational and inertial forces, intricate geometries, as well as molecular scale phenomena such as moving contact lines and thin wetting films. Although various modelling approaches have been used to tackle these challenges, the high computational demands required to accurately capture the three-dimensional fluid-fluid interfacial dynamics often render these models impractical for real-world applications. Consequently, from a practical point of view, the simplification of these 3D models may become imperative to facilitate efficient and reasonably accurate predictions of flow quantities.

Depth-integrated two-dimensional modelling is one such approach which enables saving computational time and effort at the expense of not resolving the third dimension. Here, the governing equations are solved in two dimensions, the influence of the third dimension being incorporated through appropriate additional terms. While such models have been used previously, they have so far been restricted to either permanent single-phase flow in rough fractures or two-phase flow in 2D porous media of homogeneous depth. In a rough fracture, the fluid-fluid interface possesses not only an in-plane curvature but also an out-of-plane curvature, which must be accommodated in the 2D depth-integrated model. Therefore, to address the immiscible flows in rough fractures it is essential to reformulate the 2-D depth-integrated approach from the first principles.

To perform the depth integration, we proceed from the traditional direct numerical simulation (DNS) approach, where the Navier-Stokes equations, coupled with an interface capturing technique, which in our case is the Volume of Fluid (VOF), are solved numerically. We integrate the governing flow equations in the vertical direction while expressing the flow fields in terms of 2D depth-averaged flow quantities. To account for the out-of-plane curvature and the wall shear stress arising from the no-slip conditions on the fracture walls, we assume locally a plane Poiseuille configuration (Hele-Shaw). 

The derived 2D depth-integrated model is implemented in the open-source CFD code OpenFOAM. We validate our model using the Saffman-Taylor instability case, comparing predictions with experiments and full 3D model results. We then extend our study to two numerically generated rough fractures with (a) smoothly and periodically varying aperture and (b) a more realistic aperture field with a larger roughness. We investigate drainage (i.e., the displacement of the wetting fluid by the non-wetting fluid) over a range of Capillary numbers spanning more than three orders of magnitude. We compare our 2D model predictions of both, pore-scale and macroscopic flow variables, with those obtained using 3D simulations. Our 2D model accurately estimates key statistical indicators with a tenfold reduction in computation time, offering an excellent compromise between solution accuracy and computational efficiency. We also discuss the limitations of the depth-averaged model depending on flow ranges.