Naturally fractured reservoirs (NFR), such as the large carbonate reservoirs in the Middle East, contain a major part of the world's remaining conventional oil reserves, but recovering these is especially challenging as the fractures only constitute fluid conduits while the oil is trapped in a low-permeability rock matrix. Recovery factors are therefore difficult to estimate, permeability anisotropy is high, size and shape of drainage areas are difficult to constrain, early water breakthrough is likely to be associated with a high and irreversible water cut, and secondary recovery behaviour is unusual. Outcrop-analogue model-based discrete fracture and matrix (DFM) simulations have emerged recently, helping us to disentangle and rationalize this erratic production behaviour. They allow us to understand the emergent flow behaviour and resulting saturation patterns in NFRs. Thus, classical simulation approaches, such as dual-continua conceptualizations, can be critically evaluated and improved where they fail to capture the flow behaviour of interest. This paper discusses recent advances in DFM simulation of single- and multi-phase flow processes in geologically realistic outcrop-analogue models, and solved with finite-element (FE) and finite-volume (FV) methods. It also reviews key results from recent DFM simulation studies, in particular how new measures such as the fracture-matrix flux ratio and velocity spectra can provide new means to analyse flow behaviour in heterogeneous domains or how results from outcrop-based simulations can be used to test the suitability of conventional upscaling approaches for NFR and guide the development of new ones. We close by enlisting outstanding challenges in outcrop-based flow simulations such as the need to capture the fracture-matrix transfer processes due to capillary, gravity and viscous forces accurately, which often implies detailed grid refinement at the fracture-matrix interface and small time-steps to resolve the physical processes adequately. Thus, we explore how outcrop-based flow simulation could be applied more routinely in NFR reservoir characterization and simulation workflows.@en

A joint AAPG-Society of Petroleum Engineers-Society of Exploration Geophysicists Hedberg Research Conference was held in Saint-Cyr sur Mer, France, on July 8 to 13, 2012, to review current research and explore future research directions related to improved production from carbonate reservoirs. Eighty-seven scientists from academia and industry (split roughly equally) attended for five days. A primary objective for the conference was to explore novel connections among different disciplines (primarily within geoscience and reservoir engineering) as a way to define new research opportunities. Research areas represented included carbonate sedimentology and stratigraphy, structural geology, geomechanics, hydrology, reactive transport modeling, seismic imaging (including four-dimensional seismic, tomography, and seismic forward modeling), geologic modeling and forward modeling of geologic processes, petrophysics, statistical methods, numerical methods for simulation, reservoir engineering, pore-scale processes, in-situ flow experiments (e.g., x-ray computed tomography), visualization, and methods for data interaction.@en

Hydrocarbon reservoirs commonly contain an array of fine-scale structures that are below the resolution of seismic images. These features may impact flow behavior and recovery, but their specific impacts may be obscured by the upscaling process for sector and field-scale reservoir simulations. It is therefore important to identify those situations in which subseismic structures can introduce significant departures from full-field flow predictions. Using exposures of Jurassic carbonate outcrops near the village of Amellago in the High Atlas Mountains of Morocco, we have developed a series of flow simulations to explore the interactions of a hierarchical fracture network with the rock matrix of carbonate ramp strata. Model geometries were constructed in CAD software using field interpretations and LiDAR1 data of an outcrop area that is 350 m long by 100 m high. The impact of water injection on oil recovery between an injector and producer pair was investigated. Simulations were performed by a single medium reservoir simulator using a single mesh to represent fracture planes as well as rock-matrix volumes. The effects of changing scenarios for rock permeability and porosity as well as facture permeability distributions were investigated. First-order results show that the best recovery was achieved by a model with a high permeability, homogeneous matrix combined with a heterogeneous fracture network. The worst recovery scenario was given by a model with low, homogeneous permeability and high fracture permeabilities. The results highlight the importance of the permeability contrasts between the matrix and the fractures for overall recovery and the very significant impact that fractures can have on recovery by creating shadow zones and providing critical connections between permeable layers. The presence of the hierarchical fracture network developed strong fingering even in homogeneous matrix cases and evolving velocity patterns reveal competing fluid pathways among matrix and fracture routes. Insights from these models can help to develop production strategies to improve recovery from fractured carbonate reservoirs and provide an initial platform from which to extend further evaluations of different populations of conductive and baffling structures, spatial variations in wettability and capillary pressures and well positions.@en

Discrete-fracture and rock matrix (DFM) modelling necessitates a physically realistic discretisation of the large aspect ratio fractures and the dissected material domains. Using unstructured spatially adaptively refined finite-element meshes, we find that the fastest flow often occurs in the smallest elements. Flow velocity and element size vary over many orders of magnitude, disqualifying global Courant number (CFL)-dependent transport schemes because too many time steps would be necessary to investigate displacements of interest. Here, we present a higher-order accurate implicit pressure-(semi)-implicit transport scheme for the advection-diffusion equation that overcomes this CFL limitation for DFM models. Using operator splitting, we solve the pressure and the transport equations on finite-element, node-centred finite-volume meshes, respectively, using algebraic multigrid methods. We apply this approach to field data-based DFM models where the fracture flow velocity and mesh refinement is 2-4 orders of magnitude greater than that of the matrix. For a global CFL of ≤10,000, this implies sub-CFL, second-order accurate behaviour in the matrix, and super-CFL, at least first-order accurate, transports in fast-flowing fractures. Their greater refinement, however, largely offsets this numerical dispersion, promoting a highly accurate overall solution. Numerical and fracture-related mechanical dispersions are compared in the realistic DFM models using second-order accurate runs as reference cases. With a CFL histogram, we establish target error criteria for CFL overstepping. This analysis indicates that for extreme fracture heterogeneity, only a few transport steps can be sufficient to analyse macro-dispersion. This makes our implicit method attractive for quick analysis of transport properties on multiple realisations of DFM models.

@en

Discrete-fracture modeling and simulation of two-phase flow in realistic representations of fractured reservoirs can now be used for the design of improved-oil-recovery (IOR) and enhanced-oilrecovery (EOR) strategies. Thus far, however, discrete-fracture simulators usually do not include a third compressible gaseous phase. This hinders the investigation of the performance of gas gravity drainage, water alternating gas injection, and blowdown in fractured reservoirs. We present a new numerical method that expands the capabilities of existing black-oil models for three-component, three-phase flow in three ways: (a) It uses a finite-element/finite-volume discretization generalized to unstructured hybrid element meshes, (b) It employs higher-order accurate representations of the flux terms, (c) Flash calculations are carried out with an improved equation of state allowing for a more realistic treatment of phase behavior. We illustrate the robustness of this numerical method in several applications. First, quasi-lD simulations are used to demonstrate grid convergence. Then, 2D discrete-fracture models are used to illustrate the effect of mesh quality on predicted production rates in discrete-fracture models. Finally, the proposed method is used to simulate three-component, three-phase flow in a realistic 2D model of fractured limestone mapped in the Bristol Channel, UK, and create a 3D stochastically generated discretefracture model.@en

Simulation grid blocks of naturally fractured reservoirs contain thousands of fractures with variable flow properties, dimensions, and orientations. This complexity precludes direct incorporation into field-scale models. Macroscopic laws capturing their integral effects on multiphase flow are required. Numerical discrete fracture and matrix simulations show that ensemble relative permeability as a function of water saturation (k_ri[S_w]), water breakthrough, and cut depend on the fraction of the cross-sectional flux that occurs through the fractures. This fracture-matrix flux ratio [q _fqm] can be quantified by steady-state computation. Here we present a new semianalytical model that uses q_fqm and the fracture-related porosity fjf) to predict k_ri(S_w) capturing that, shortly after the first oil is recovered, the oil relative permeability (kro) becomes less that that of water (k_rw), and k _rwk_ro approaches q_f/q_m as soon as the most conductive fractures become water saturated. To include a capillary-driven fracture-matrix transfer into our model, we introduce the nonconventional parameter A_fw/_Sw), the fraction of the fracture-matrix interface area in contact with the injected water for any grid-block average saturation. The A_fw/_Sw) is used to scale the capillary transfer modeled with conventional transfer functions and expressed in terms of a rate- and capillary-pressure- dependent k_ro. All predicted parameters can be entered into conventional reservoir simulators. We explain how this is accomplished in both, single- and dual-continua formulations. The predicted grid-block-scale fractional flow (f _i[S_w]) is convex with a near-infinite slope at the initial saturation. The upscaled flow equation therefore does not contain an S _w shock but a long leading edge, capturing the progressively widening saturation fronts observed in numerical experiments published previously.

@en

We have been able to solve a reservoir simulation problem which was previously thought of as intractable: We simulated multiphase displacement, including viscous, capillary, and gravitational forces, for highly resolved and geologically realistic models of naturally fractured reservoirs (NFR) at the sector, i.e. kilometre, scale with very reasonable runtime. This has been possible because we used massive parallelisation and hierarchical solvers in conjunction with a new discrete fracture and matrix modelling (DFM) technique that is based on mixed-dimensional unstructured hybrid-element discretisations. High-resolution DFM simulations are important to resolve the non-linear coupling of small scale capillary - viscous and large scale gravitational - viscous processes adequately for sector scale NFR. Cross-scale process coupling in NFR controls oil recovery and NFR often exhibit power-law fracture length distributions, i.e. they do not possess an REV, and highly permeable fractures can extend over the full hydrocarbon column height. As a consequence, emergent displacement patterns have been observed which are difficult to quantify using traditional means of upscaling. However, such patterns could now be used as benchmarks to reach a better consensus on the correctness of promising new upscaling techniques. The parallel DFM technologies presented here allow us to to obtain these results much more efficiently and hence explore the parameter space in greater detail. We observed a linear scaling behaviour for up to 64 processes and a significant decrease in runtime when applying our parallel DFM approach to three highly refined NFR simulations. These contain thousands of fractures, up to 5 million elements, and have local grid-refinements below 1 m for model dimensions between I and 10 kilometres. We achieved this excellent speedup because we reduced inter-processor communication by minimising the overlap between individual domains and decreased idle time of individual processors by distributing the number of unknowns equally among the processors.@en

Field data-based simulations of geologic systems require much computational time because of their mathematical complexity and the often desired large scales in space and time. To conduct accurate simulations in an acceptable time period, methods to reduce runtime are required. A parallelization approach is attractive because fast multi-processor clusters are nowadays readily available. Here we report on our recent efforts to parallelize our multiphysics code CSMP ++ (Complex System Modelling Platform). In particular, we describe a parallel finite element-finite volume method for multi-phase fluid flow in heterogeneous porous media. We take a domain partitioning approach where the finite element mesh is partitioned into sub-domains, assigning each of them to a single processor. For each sub-domain a local finite volume mesh is constructed. We can now solve advection-dispersion type equations taking an operator splitting approach: Pressure diffusion is calculated with an implicit finite element method and advection with an implicit or explicit finite volume scheme. We have tested the accuracy, robustness and computational speedup of our new parallel scheme on a Linux cluster by means of three geologic applications. All tests give excellent computational speedup with increasing number of up to 32 processors. These results broaden the range of possible simulations in terms of spatial and temporal scale and resolution as well as numerical accuracy up to two orders of magnitude.@en

Realistic simulation of structurally complex reservoirs (SCR) is challenging in at least three ways: (1) geological structures must be represented and discretized accurately on vastly different length scales; (2) extreme ranges and discontinuous variations of material properties have to be associated with the discretized structures and accounted for in the computations; and (3) episodic, highly transient and often localized events such as well shut-in have to be resolved adequately within the overall production history, necessitating a highly adaptive resolution of time. To facilitate numerical experiments that elucidate the emergent properties, typical states and state transitions of SCRs, an application programmer interface (API) called complex systems modelling platform (CSMP ++) has been engineered in ANSI/ISO C ++. It implements a geometry and process-based SCR decomposition in space and time, and uses an algebraic multigrid solver (SAMG) for the spatio-temporal integration of the governing partial differential equations. This paper describes a new SCR simulation workflow including a two-phase fluid flow model that is compared with ECLIPSE in a single-fracture flow simulation. Geologically realistic application examples are presented for incompressible 2-phase flow, compressible 3-phase flow, and pressure-diffusion in a sector-scale model of a structurally complex reservoir.

@en

Discrete-fracture modeling and simulation of two-phase flow in realistic representations of fractured reservoirs can now be used for the design of IOR and EOR strategies. Thus far, however, discrete fracture simulators fail to include a third compressible gaseous phase. This hinders the investigation of the performance of gas-gravity drainage, water alternating gas injection, and blow-down in fractured reservoirs. Here we present a new numerical method that expands the capabilities of existing Black-Oil models for threecomponent - three-phase flow in three ways: i) It utilizes a finite element - finite volume discretization generalized to unstructured hybrid element meshes, ii) It employs higher-order accurate representations of the flux terms, iii) Flash calculations are carried out with an improved equation of state allowing for a more realistic treatment of phase behavior. We illustrate the robustness of this numerical method in .several applications. First, quasi-ID simulations are used to demonstrate grid convergence. Then, 2D discrete fracture models are employed to illustrate the impact of mesh quality on predicted production rates in discrete fracture models. Finally, the proposed method is used to simulate three-component - three-phase flow in a realistic 2D model of fractured limestone mapped in the Bristol Channel, U.K. and a 3D stochastically generated discrete fracture model.@en

We present the benchmarking of a new finite element - finite volume (FEFV) solution technique capable of modeling transient multiphase thermohaline convection for geological realistic p-T-X conditions. The algorithm embeds a new and accurate equation of state for the NaCl-H2O system. Benchmarks are carried out to compare the numerical results for the various component-processes of multiphase thermohaline convection. They include simulations of (i) convection driven by temperature and/or concentration gradients in a single-phase fluid (i.e., the Elder problem, thermal convection at different Rayleigh numbers, and a free thermohaline convection example), (ii) multiphase flow (i.e., the Buckley-Leverett problem), and (iii) energy transport in a pure H2O fluid at liquid, vapor, supercritical, and two-phase conditions (i.e., comparison to the U.S. Geological Survey Code HYDROTHERM). The results produced with the new FEFV technique are in good agreement with the reference solutions. We further present the application of the FEFV technique to the simulation of thermohaline convection of a 400°C hot and 10 wt.% saline fluid rising from 4 km depth. During the buoyant rise, the fluid boils and separates into a high-density, high-salinity liquid phase and a low-density, low-salinity vapor phase.@en

We present a new finite element - finite volume (FEFV) method combined with a realistic equation of state for NaCl-H₂O to model fluid convection driven by temperature and salinity gradients. This method can deal with the nonlinear variations in fluid properties, separation of a saline fluid into a high-density, high-salinity brine phase and low-density, low-salinity vapor phase well above the critical point of pure H₂O, and geometrically complex geological structures. Similar to the well-known implicit pressure explicit saturation formulation, this approach decouples the governing equations. We formulate a fluid pressure equation that is solved using an implicit finite element method. We derive the fluid velocities from the updated pressure field and employ them in a higher-order, mass conserving finite volume formulation to solve hyperbolic parts of the conservation laws. The parabolic parts are solved by finite element methods. This FEFV method provides for geometric flexibility and numerical efficiency. The equation of state for NaCl-H₂O is valid from 0 to 750°C, 0 to 4000bar, and 0-100 wt.% NaCl. This allows the simulation of thermohaline convection in high-temperature and high-pressure environments, such as continental or oceanic hydrothermal systems where phase separation is common.

@en

We present new, accurate numerical simulations of 2D models resembling hydrothermal systems active in the high-permeability axial plane of mid-ocean ridges and show that fluid flow patterns are much more irregular and convection much more unstable than reported in previous simulation studies. First, we observe the splitting of hot, rising plumes. This phenomenon is caused by the viscous instability at the interface between hot, low-viscosity fluid and cold, high-viscosity fluid. This process, known as Taylor-Saffman fingering could potentially explain the sudden extinguishing of black smokers. Second, our simulations show that for relatively moderate permeabilities, convection is unsteady resulting in transiently varying vent temperatures. The amplitude of these fluctuations typically is 40 °C with a period of decades or less, depending on the permeability. Although externally imposed events such as dike injections are possible mechanisms, they are not required to explain temperature variations observed in natural systems. Our results also offer a simple explanation of how seismic events cause fluctuating temperatures: Earthquake-induced permeability-increase shifts the hydrothermal system to the unsteady regime with accompanying fluctuating vent temperatures. We demonstrate that realistic modelling of these high-Rayleigh number convection systems does not only require the use of real fluid properties, but also the use of higher order numerical methods capable of handling high-resolution meshes. Less accurate numerical solutions smear out sharp advection fronts and thereby artificially stabilize the system.@en

Water flooding of fractured reservoirs is risky because water breakthrough can occur early, leading to a prohibitively high water cut. In mixed or oil-wet carbonates, capillary drive is negligible or absent. For this scenario, we investigate fluid-pressure-driven displacement of oil by water in two-phase flow numerical models based on naturally fractured limestone beds mapped along the British Channel coast. These reservoir analogs are represented by unstructured finite-element grids with discrete representations of intersecting fractures. We solve the governing equations for slightly compressible two-phase flow with our original control-volume finite-element method. This permits the direct examination of displacement patterns in fractures and rock matrix. We find that the irreducible saturation in the fractured carbonate is much higher than the value prescribed to the rock matrix. The shape of water invasion fronts is highly sensitive to the viscosity ratio of oil and water. When tho Brooks-Corey relative permeability model is applied to the rock matrix at a viscosity ratio of 1, the total mobility, λ_t, is low at intermediate saturations. This stabilizes displacement fronts where a girdle of reduced λ _t develops, but this effect disappears as the viscosity ratio increases. For an idealized model with a water-wet matrix, we have also evaluated the effect of countercurrent capillary-pressure-driven flow across fracture-matrix interfaces. The rate of this countercurrent imbibition scales with the specific fracture surface area and decays exponentially as intermediate saturation zones develop adjacent to the fractures. The resulting reduced λ_t feeds back into the fluid-pressure-driven displacement process. Copyright

@en

The permeability of the Earth's crust commonly varies over many orders of magnitude. Flow velocity can range over several orders of magnitude in structures of interest that vary in scale from centimeters to kilometers. To accurately and efficiently model multiphase flow in geologic media, we introduce a fully conservative node-centered finite volume method coupled with a Galerkin finite element method on an unstructured triangular grid with a complementary finite volume subgrid. The effectiveness of this approach is demonstrated by comparison with traditional solution methods and by multiphase flow simulations for heterogeneous permeability fields including complex geometries that produce transport parameters and lengths scales varying over four orders of magnitude.

@en