We propose an innovative numerical model for single-phase flow in carbonate reservoirs containing karstified layers, characterized by partially enlarged fractures and high-permeability conduits formed through superimposed chemical dissolution along fracture intersections. The computational framework is based on a discontinuous Galerkin (dG) formulation applied to a modified system of mixed-dimensional flow equations, which explicitly incorporates permeability enhancement due to localized fracture enlargement near intersections. The methodology is designed with a hierarchical structure, allowing for a consistent and unified multiscale representation of the porous domain, encompassing 3D flow in the rock matrix, 2D flow within fractures, 1D flow along the partially dissolved conduits, and 0D interaction at fracture-conduit junctions. This formulation proves highly effective in capturing the complex, multidimensional flow dynamics induced by karstification, and in quantifying its influence on flow patterns and production curves. We also provide a rigorous mathematical analysis, establishing well-posedness and convergence in broken Sobolev spaces. A suite of computational experiments is presented to demonstrate the robustness, accuracy, and predictive capabilities of the model when applied to realistic karstified reservoir scenarios.

We introduce the Multiscale Embedded Discrete Karst Model (MsEDKM), designed to represent the complex mass-exchange mechanisms between karst conduits and the surrounding porous matrix in carbonate formations. The model consists of a generalization of the karst index-based formulation proposed in Murad et al., (2020) and is developed within a three-scale framework that facilitates the accurate computation of macroscopic exchange transmissibilities, obtained by projecting a nonlocal integral kernel onto a coarse computational grid. At the fine (microscale) level, conduit geometries are reconstructed from static geological models generated using geostatistics. The resulting topology exhibits pronounced irregularities characterized by shell waviness and local variations in cross-sectional area. In the micro-to-meso upscaling stage, this intricate conduit geometry is systematically replaced by an Equivalent Elliptical Cylinder (EEC) using a moment-of-inertia-based geometric homogenization strategy. A subsequent meso-to-macro homogenization step is carried out through a flow-based upscaling procedure, which quantifies the matrix–conduit mass transfer in terms of non-neighboring transmissibilities in the context of the finite volume method. To accurately characterize these exchange processes at the quasi-stationary regime, a comprehensive set of mesoscopic simulations is performed using the finite-element solution of a transient diffusion flow equation. These simulations span a wide range of conduit hydraulic properties and spatial karst-conduit configurations embedded within a coarse cell, generating a dataset used to train a surrogate model, in which geometric and hydraulic descriptors are used as input features, and the upscaled transmissibilities obtained from mesoscopic numerical simulations constitute the target outputs. The trained model demonstrates high predictive performance, reproducing the outcomes of direct simulations at a fraction of the computational cost. Finally, macroscopic simulations incorporating machine learning-predicted transmissibilities within the MsEDKM framework demonstrate the potential of the proposed approach for reservoir-scale applications. This integrated framework provides a robust and scalable tool for improving the prediction and management of flow behavior in karstic reservoirs.

We develop an innovative mixed-dimensional 3D/1D flow model in carbonate rocks containing multiple karst cave conduits with underlying heterogeneity in the petrophysical properties stemming from different geological stages of cave-pipe collapse systems. Such geological structures manifest in distinct heterogeneity patterns inherent to the successive stages of burial, mechanical failure, and collapse, resulting in discrete collapsed passages in the conduit network. In addition, breakdown products appear within the cave system associated with chaotic breccia, suprastrata deformation, and vertical tube-like geo-bodies, herein referred to as breccia pipes, containing faults and fractures around the vertical pipe. The input parameters of the mixed-dimensional flow model show the ability to incorporate the complex multiple heterogeneities associated with the geological objects at different stages of collapse. After populating the geo-bodies with proper petrophysical properties, the mixed-dimensional flow equations are discretized by a locally conservative extended version of the mixed-hybrid finite element method, which incorporates the new nonlinear discrete transmission jump conditions between elements adjacent to the breccias within the conduits. Computational simulations are performed for particular configurations of heterogeneous karst conduit systems with distinct geological time scales, illustrating the influence of the karst and solution breccia-pipe deposits upon the flow regimes, streamline patterns, and well productivity in real-case scenarios of hypogenic cave networks.

We propose a new computational model for solving the Black-Oil flow model, incorporating geomechanical coupling within the framework of the fixed-stress-split scheme. The extended flow equations describing the movement of two slightly compressible liquids and a highly compressible gas are recast in terms of multiphase-multicomponent flow. Here, we construct a nonlinear extension of the fixed-stress split proposed in earlier work (Correa and Murad, J. Comput. Phys. v.373, pp. 493-532, 2018), which also allows for the compressibility of the liquid phases, dissolution of the gas in the oil phase, and gas phase appearance and disappearance. Flow and transport subsystems are rephrased in terms of compositions, and two alternative sequential coupling strategies at different levels are introduced to link the flow/transport and mechanics subsystems. These strategies incorporate a suitable definition of a trusted saturation variable within the flow equations, ensuring the full resolution of a three-equation system of conservation laws for the compositions and thereby enhancing the overall stability and accuracy of the proposed scheme. Flow and mechanics subsystems are discretized by mixed finite element formulations, whereas the transport system is solved by an innovative semi-discrete central-upwind finite volume scheme for hyperbolic conservation laws, capable of capturing spatial and temporal variability in the Lagrangian porosity and also obviating the need to adopt operator-splitting schemes for the storativity in the transport equations. The innovative numerical model clearly demonstrates its ability to capture the intricate interaction between geomechanical effects and phase change in the vicinity of the bubble point. Numerical experiments are performed, including an undrained setting upon cyclic loading and water-flooding problems, illustrating precisely the influence of the bubble point pressure upon the evolution of the poromechanical variables and hydrocarbon production.

We develop an enhanced reduced model for single-phase flow in fractured porous media capable of incorporating more realistic interface conditions at the fracture terminations. In addition to the traditional dimensional model reduction, where the elements of the discrete fracture network are treated as lower dimensional manifolds embedded in the porous matrix, we explore the microscale behavior of the boundary layer flow at the entrances of a fracture bounded by two parallel plates to construct a new set of interface conditions of Robin-type, giving rise to localized pressure jumps at the fracture edges. Within this enriched description, sharper reduced flow and tracer transport mixed-dimensional models are constructed in the asymptotic limit ruled by two small parameters related to the ratio between fracture aperture and entrance developing length and a macroscopic length scale.The discrete flow/transport mixed-dimensional model is discretized by a new discontinuous Galerkin(dG)-based formulation. An adequate version of the Galerkin-Newton method is developed for the numerical treatment of the non-linear Robin interface condition. Considering several fracture arrangements, numerical results illustrate the sharper description of the model proposed herein in predicting flow and tracer transport patterns in fractured media.

We construct herein a three-scale coupled mechanical model for naturally fractured coalbed methane reservoir with the ability of describing the stress balance between the solvation force, arising from the gas adsorption in nanopores, and the restoration stress stemming from the elastic response of the cleats. To determine the cleat porosity, the non-linear hyperbolic Barton-Bandis (BB) law, which captures increase in joint stiffness induced by the cleat closure due to matrix swelling, is postulated for the fracture mechanical response. At the microscale, the theory incorporates the coupling between the effects of the solvation force and the elastic response of the matrix. Such system of governing equations is coupled with the fluid pressure in the discrete cleat system with dependency of aperture with the normal stress dictated by the aforementioned BB-model. A reiterated homogenized procedure is pursued and capable of providing the constitutive response of the homogenized poromechanical parameters on gas pressure. Numerical simulations illustrate the performance of the proposed model.

We construct a new three-scale model for single phase incompressible flow in faulted rocks containing multiple damage zones. At the finer scale (O(1m)), flow is influenced by the high-contrast layered heterogeneities inherent to the core and adjacent damage zones, which are populated by geological anomalies, such as compaction bands, debris, and fine sediments and joints. In the first stage of the reiterated homogenization procedure, we construct a lower-dimensional reduced model, where the discontinuity is envisioned as (n−1)-dimensional manifold (n=2,3), topologically attached to multiple layered structures, with flow patterns characterized by various jumps in pressure and velocity fields. Subsequently, by aligning the fault with the interface between adjacent simulation cells of a coarse grid, the upscaling of the reduced flow model gives rise to transmissibility multipliers, whose constitutive response stems naturally from the local flow patterns, exhibiting improved accuracy compared to the traditional harmonic mean, inherent to the two-point flux approximation. Computational simulations obtained with the finite element method with localized discontinuous spaces illustrate the ability of the three-scale model to provide further insight into the behavior of the transmissibility multipliers, which capture the effects of fault zone texture upon the flow discretization. The methodology proposed herein shows enormous potential for the development of more accurate transmissibility preprocessors at relatively low computational costs, consequently overcoming the shortcomings of a direct application of the local high-fidelity approach.

We construct a new computational model to describe coupled 3D/1D flow in carbonate rocks intertwined by a network of karst cave conduits. The proposed approach shows ability to incorporate pointwise velocity-dependent jumps in the pressure field arising from localized partial obstructions due to the presence of collapse-breccia within the discrete conduit network. At the microscale, we postulate single phase viscous flow governed by the Navier-Stokes equations in the conduit network coupled with Darcian flow in the rock matrix and supplemented by transmission conditions at the common interface. Subsequently, we proceed by constructing a sharper lower-dimensional reduced model wherein, in addition to the usual high geometric aspect ratio between the length and hydraulic diameter of the cave system, we introduce an additional small parameter containing the ratio between the localized width of the perturbed flow region, in the vicinity of each breccia, and characteristic length of the network. The asymptotic behavior gives rise to a coupled mixed-dimensional flow, where 1D sub-manifolds appear embedded in the 3D carbonate matrix with coupling ruled by a mass exchange line-source δ -function, acting synergistically with discrete non-linear pressure jumps of Robin type at the discrete set of breccia locations. The mixed-dimensional flow equations are discretized by a locally conservative extended version of the mixed-hybrid finite element method, showing capability of incorporating the new non-linear discrete transmission jump conditions between elements adjacent to the breccias. Computational simulations are performed for particular configurations of well/karst conduit systems, illustrating the influence of the karst and breccia upon the flow regimes, streamline patterns and well productivity.

We construct a new operator splitting scheme to describe a fluid-driven brittle fracture propagation in a Biot medium based on a time-scale separation assumption. In this context, we propose an alternate hybrid time-stepping scheme, where in the injection step, prior to fracture advance, we explore the framework of the fixed stress split scheme with a hydrodynamic subsystem solved ahead of the geomechanics for a frozen total mean stress. Conversely, after convergence of the fixed stress split iterations, when the pore pressure exceeds a critical value, the coupling between poromechanics and fracture propagation is accomplished by considering a fast time scale, with frozen pore-pressure and Darcy velocity fields. Such a latter step is performed in a separate iteration loop with the elasticity subsystem incorporating the Francfort Marigo variational model for a thin damage region. In this setting, the evolution of the damaged zone is governed by the sensitivity of the associated shape functional with respect to the nucleation of a small damaged zone, which is computed within the framework of the topological derivative method. The resulting approach is algorithmically described in detail. A numerical assessment of the model is constructed by performing a series of benchmark examples, showing different features of the proposed approach, such as characterization of fracture-activation pressure, crack path forecast, and the ability to capture kinking and bifurcations and quantifying the effects of the in situ stress field on the crack path.

We construct an innovative static-dynamic integrated workflow capable of bridging the gap between input geological data, inherent to a lacustrine carbonate outcrop containing karst geobodies, and the description of the flow patterns and quantification of the multi-well productivity index (MWPI) for a particular well configuration in the outcrop. The workflow incorporates additional features stemming from the use of Machine Learning-based methods to mitigate lack of data in the locations away from the sections of input signals, along with the construction of new upscaling methods to assess the MWPI matrix. The ML-enhanced geostatic model hinges upon shallow surface geophysical data collected using Ground-Penetrating Radar (GPR) techniques. Furthermore, by discretizing the flow equations and adopting a flow-based upscaling method, we construct correlations between well flow rates and pressure drawdown in a typical five-spot well configuration. In this setting, we analyze the sensitivity of each well productivity with respect to heterogeneity distribution and correlations in the karst system within the outcrop. Computational simulations illustrate the ability of the integrated workflow proposed herein to improve prediction of hydraulic-connectivity between well pairs, which appear manifested in the entries of the MWPI matrix, whose magnitude aims at quantifying the effects of the karst geobodies upon geofluid production.

We devise a new computational model to describe numerically the coarse-scale response of coupled conduit/matrix flow in karstified carbonate rocks. The methodology is based on the combination of a two-scale coupled 1D/3D flow model, rigorously derived by formal homogenization, and an extension of the recursive Mixed Multiscale Method (MuMM), seated on domain decomposition and redesigned herein to more complex bi-modal systems. Within this double-coarsening framework, our target macroscopic scenario consists of a limestone matrix intertwined by a network of karst-conduits of high aspect ratio between longitudinal length and hydraulic diameter. In order to upscale the high-fidelity flow equations to a mesoscale model characterized by the presence of 1D sub-manifolds, associated with the conduit network, embedded in the 3D carbonate matrix, we proceed within a model reduction seated on matched asymptotic expansions. Such an upscaling leads to the appearance of an exchange function described by a line-source $$\delta$$ -function localized in the coordinates intercepted by the conduit network. Subsequently, the extended multiscale method is constructed to numerically capture the macroscopic response of the karstified medium by properly designing multiscale basis functions for the mixed formulation of the reduced problem, adapted to the scenario with presence of the local conduit/matrix exchange functions, numerically approximated by a local regularization of the $$\delta$$ -function. Within the framework of the recursive multiscale domain decomposition methodology, the coarse-grid interface problem, which assembles the local solutions, is replaced by the family of small localized interface linear systems obtained by the clustering of the multiscale basis functions associated with the nearest neighbor subdomains. Parallel implementation of the recursive procedure proposed herein exhibits good efficiency and scalability properties, showing great potential for the resolution of giant karstified reservoirs containing up to billions of cells. Computational simulations are performed considering a particular karst arrangement of branchwork-type in the sense of the pioneering work of Palmer (Geol Soc Am Bull 103: 1–21, 1991). To assess the quality of the coarser-scale multiscale solutions, numerical results of pressure and velocity fields are computed and compared with a reference solution computed by the mixed-hybrid finite element method adopting fine meshes.

We construct a new symmetric interior penalty discontinuous Galerkin (SIPDG)-method for incompressible flow in fractured porous media. The method is developed for computing accurate approximations of the reduced flow model, where fractures are treated as lower dimensional objects. Unlike previous formulations, the methodology proposed herein shows ability to capture the asymptotic limits of highly permeable and fracture seals, also covering a wider range of values of the quadrature integration parameter which appears in the pressure jumps across the fractures. As a first novel contribution we should mention that the method was developed for specific interface condition in the reduced problem for which known in literature DG method are not applicable. As a second novelty we propose a new penalty technique for stabilization of the method and obtain explicit estimate for the penalty parameters associated with flow in matrix and fractures in order to achieve stability. And finally we derive new high order hp type a priori error estimates for the numerical solution in the energy norm. Numerical results illustrate the performance of the proposed SIPDG-method in simulating discrete fracture models.

We propose a new computational model for two-phase immiscible flow in a poroelastic medium overlain by a saline formation displaying creep behavior with viscous strain ruled by a nonlinear dependence of power-law type on the deviatoric stress. Within the framework of the fixed-stress split algorithm, hinged upon freezing the total mean stress in the flow equations, hydrodynamics and geomechanics subsystems are solved by mixed finite element methods whereas the hyperbolic equation for the water saturation is formulated in terms of the Lagrangian porosity and solved adopting an operator splitting fractional-step method combined with a higher-order non-oscillatory finite volume central scheme. Within this typical multi-physics setting, the sequential formulation allows greater flexibility in the choice of the meshes for the subsystems, particularly for the geomechanics module, which can be solved in the extended domain including the adjacent impervious formations (over-, under- and side-burdens). Considering the upper interface of the rock salt dictated by the profile of saline dome geobodies, whose impermeability property precludes hydraulic communication between adjacent formations, the initial in situ undrained state is also modeled by invoking Skempton’s coefficient along with the characterization of the undrained bulk and shear modulus. By constructing a robust evolving algorithm between the discrete counterparts of the three subsystems in the reservoir, in conjunction with an iterative scheme underlying the mixed method for the nonlinear creep problem in the rock salt, numerical simulations of a water-flooding problem in secondary oil recovery are performed for different realizations of the input random fields. Numerical results illustrate the influence of viscoelastic effects in the cap rock upon subsidence, reservoir compaction, finger grow and breakthrough curves. Comparisons with the performance of traditional one-way formulation are also presented.

We develop a new three-scale (micro/meso/macro) computational model based on a reiterated homogenization procedure to describe flow in carbonate rocks containing complex geological structures, such as fractures and solution-collapse breccias in the sense of Loucks (AAPG Bull. 83(11), 1795–1834, 1999). In this setting, we construct a hierarchical karst-fracture model wherein the larger geological objects are incorporated explicitly whereas the higher density microscopic structures are homogenized and replaced by equivalent continua with properties computed from self-consistent homogenization schemes. In the upscaling method, we subdivide the different clastic arrangements in the breccia into crackle, mosaic, and chaotic substructures where equivalent permeability and elastic constants are assigned to each layer within the breccia. After reconstructing the mesoscopic coefficients, we adopt a flow-based upscaling to the macroscale, where the characteristic length is associated with a typical coarse grid cell of a reservoir simulator. The mesoscopic flow equations are constructed based on the discrete fracture model (DFM) and discretized by a robust computationally scheme with the ability to handle strong heterogeneity induced by the collapse breccia and pressure jumps across flow barriers. In addition to the scenario wherein the collapse breccia network is composed of disconnected objects (isolated chambers), we also develop a reduced model for the case of connected karst facies playing the role of a network of enlarged fractures. Numerical results, with input data extracted from outcrops drone images, are presented illustrating the influence of different settings on flow patterns and their effect upon the magnitude of macroscopic properties.

We develop a new multiscale model to compute effective properties such as relative permeability, contact angle and partition coefficients in low salinity enhanced oil recovery processes for two-phase flow in sandstones containing reactive surfaces of kaolinite clay. In this setting, we construct a three-scale approach which entails the local nanoscale description ruled by the electro-chemistry of a confined electrolyte solution containing Na + , Ca^2+ , H^+ , Cl^- and OH^- ions residing between bounded crude-oil droplets at residual saturation and clay substrate. Our analysis focuses on the case of surface complexation geochemical reactions between the ionic species of the invading water and the electrically charged kaolinite and oil–water interfaces. In this scenario, we construct a local electric double layer problem for the electric potential based on a non-symmetric Poisson–Boltzmann equation supplemented by nonlinear boundary conditions with the magnitude of the surface charge strongly dictated by the geochemical reactions. By invoking the local mechanical equilibrium of the electrolyte solution and solving numerically the nonlinear problem using the finite element method, we compute the local ionic profiles and reconstruct numerically the disjoining pressure and adsorption isotherms for each ionic species for a wide range of brine compositions and pH of the water phase. Furthermore, combining the disjoining pressure results with the Frumkin/Derjaguin wetting theory allows to compute the dependence of the contact angle on wettability, pH and salinity. Subsequently, the formal homogenization procedure is adopted to upscale the pore-scale flow and ion transport to the macroscale giving rise to a new Darcy scale coupled flow/transport model. The hyperbolic part of the nonlinear homogenized model is solved analytically in an 1D example of enhanced oil recovery.

We develop a new computational model to describe hydro-mechanical coupling in fractured rocks composed of a linear poroelastic Biot medium and nonlinear elastic joints with constitutive response governed by the Barton–Bandis (BB) law. The model aims at capturing increase in stiffness induced by fracture closure during fluid withdrawal. The nonlinear hydro-mechanical formulation is constructed within the framework of the Discrete Fracture Model, with flow and geomechanical sub-systems coupled through a sequential iterative algorithm. The internal contact constraint arising from non-overlapping between opposite fracture faces is enforced through the weak fulfillment of the BB-law. Such a constraint is captured within the framework of the Augmented Lagrangian formulation, where the non-linear mechanical interaction is enforced adopting successive approximations of the Lagrange multiplier, interpreted as the contact pressure in the joint, supplemented by a penalty component associated with the rock stiffness. Furthermore, adopting a traditional flow based upscaling method, macroscopic permeabilities are numerically reconstructed with magnitude strongly dependent on the local stress state. Such a mechanical dependence of the homogenized properties is represented in a discrete manner through pseudo-coupling tables, with enormous potential to be explored within a preprocessing stage in reservoir simulators to compute multipliers relative to a chosen pre-stressed reference state, where input data is available. Numerical simulations are performed for some fracture arrangements illustrating the potential of the formulation proposed herein in bridging hydromechanical coupling at different scales in jointed rocks.

A micro-structure model of dual porosity type for flow and deformation in deformable porous media that incorporates intra-grain delayed how during secondary consolidation is proposed. The model is derived within the framework of the homogenization technique and is resultant from scaling up the constitutive relations on the finer scales. In particular, from upscaling the poroelastic equations for the grains (matrix blocks) coupled with the Stokesian bulk-water movement in the system of larger fissures. In the homogenized dual porosity model with microstructure the macroscopic porous medium is represented as two porous distinct structures coexisting at each macroscopic point: one representing the matrix blocks (or grains) and the other the bulk water. In this picture, the macroscopic flow is that of the bulk water and matrix blocks act as sources/sinks of mass and momentum to the macroscale bulk phase. In the case of deformable porous media under consolidation processes, this source/sink transfer function governs the creep constitutive equations of the macroscale effective stress tensor as it accounts for secondary deformation of the skeleton arising from the delayed drainage of the fluid within the blocks. The governing equations of dual porosity type are discretized by the finite element method. Numerical simulations of creep during secondary consolidation are presented to illustrate the performance of the proposed approach.