Geomodeling Insights for Numerical Simulations of Naturally Fractured Reservoirs

Abstract

The continuous evolution of technologies in subsurface petroleum exploration and characterization makes imperative the creation of multidisciplinary up-to-date technical workflows aimed to represent complex reservoirs with robust reliability. The objective of this chapter is to promote an innovative and integrated collection of the latest technologies available for the study of naturally fractured hydrocarbon reservoirs. We assembled some of the present-day advanced technologies ranging from high-resolution digital core analysis, borehole fracture characterization, numerical deformation and geomechanics modeling, stochastic fracture modeling and the numerical simulation of dual-porosity, dual-permeability fractured reservoirs. The efficient application of robust fracture models in field development is also discussed. This chapter is intended to serve as a practical and structured handbook for professionals working on subsurface reservoir modeling, with emphasis on fractured reservoirs in the petroleum industry.

1. Introduction

Natural fractures distribution is the result of independent natural processes affecting rock masses at different times during the geological evolution of a region. Whereas tectonic deformation and the structural style influences fracture orientation and the spatial distribution of the fracture network, the distribution of intrinsic rock mechanical properties controls mainly vertical extent and intensity in layered rocks. Additionally, in-situ stress regime affects the hydraulic conductivity of the fractured media. This chapter describes an advanced reservoir fracture modeling workflow [1] that incorporates main geological processes and integrates conventional and modern technologies for the study of relevant aspects of the natural fracture systems:

• Core and wellbore fracture characterization

• 3D numerical deformation modeling

• Seismic discontinuity attributes

• 3D Mechanical Earth Modeling

• Dual-Porosity/Dual-Permeability numerical simulation

Commonly, subsurface naturally fractured reservoirs are modeled explicitly in a 3D pillar-grid. Discrete fracture planes of which position, spacing, distribution, orientation and shape as well as permeability and aperture, are modeled with advanced geostatistical workflows. The modeling process described herewith combines continuous and discrete variables to construct a stochastically distributed discrete fracture network (DFN).

Existing 3D deformation modeling solutions use paleo-structural reconstruction methods aimed to identify paleo tectonic events that originated natural fractures. The 3D deformation algorithms are based on boundary element and finite element simulation techniques, such as Kine3D [2], which solve for strains and stresses related mainly to tectonic deformation. The resultant tectonic drivers are used as an input to predict spatial distribution and intensity variability of the fracture network. Solid understanding of the structural setting and tectonic evolution is a key factor when deciding on the 3D paleo-structural reconstruction method to use. A 3D mechanical earth model includes intrinsic rock mechanical properties and the in-situ stress regime, represented as continuous variables. Elastic and strength properties as well as stress magnitudes are the input to compute 3D discrete geo-mechanical facies, which represent mechanical units within the reservoir with differential response to the stress/strain deformation. The geo-mechanical facies can be validated with direct observation or measurement of fracture density from borehole image log tools. Additionally, the in-situ stress regime is utilized to identify fluid-flow pathways by applying critical stress state methodology [1].

Combining these major elements into a stochastic fracture modeling workflow requires sequential and assertive selection of simulation methods and geostatistical techniques. Furthermore, optimization solutions, such as genetic algorithms with machine learning and probabilistic methods to generate multiple scenarios, facilitate the evaluation of uncertainty associated with parametrization to assist and accelerate dual porosity-dual permeability (DPDP) numerical reservoir simulation and field development.

2. Dual porosity dual permeability geomechanics model

The main parameters that describe subsurface natural fracture systems are orientation, density, distribution, aspect ratio, permeability and aperture. Advanced techniques, such as high-resolution X-ray micro computed tomography scan, stress dependent permeability tests (SDK) and Azimuthal shear anisotropy analysis (ASAA) among others, allow detailed characterization of those key elements ([3] and references therein). These tools, combined with the conventional borehole image log analysis, represent an integrated way to characterize the natural fractures at wellbore scale.

The geo-mechanical workflow consists of several steps to characterize the natural fractures and mechanical properties:

• 1D Mechanical Earth modeling and fracture characterization is a process that integrates a full description of natural fractures from borehole images including fracture types, fracture intensity, orientation and aperture. The process deploys an extensive rock mechanical test program to obtain the most reliable information about the mechanical properties. The program includes ultrasonic test, multi-stage tri-axial test, pore volume compressibility (PVC), stress dependent permeability and fracture porosity measurements. Additionally, total volume strain is recorded from the tri-axial test to create in the geo-mechanical facies classification. Where available, walkaway vertical seismic profiles (VSP) are used to extract acoustic anisotropy properties, showing the fast shear direction with the natural fractures strikes similar orientation.

• 3D Seismic discontinuity analysis combines several attributes (e.g., similarity or coherence and curvature) to determine the possible discontinuity trends for the entire field. The cluster analysis is used to divide the field into sectors, and within these sectors, conduct a comparison in terms of fracture and seismic discontinuity orientation.

• The 3D Mechanical Earth model: Elastic rock mechanical properties and rock strength properties are populated using the porosity model from the 3D grid geological model to maintain consistency between petrophysical and mechanical properties. The in-situ stress tensor and strain tensor are calculated by using finite element model (FEM) geomechanics simulator capturing features such as possible stress rotation due to faults configurations. The geo-mechanical facies are modeled using principal component analysis (PCA) and neural network-based classification techniques: these facies are contrasted against the density of natural fractures to establish concordance between the brittle facies and high-density fracture.

• Paleo-deformation model is used to evaluate different tectonic stages related to faults creation to find the best fit values that can represent most of the natural fractures from the set observed in the borehole image interpretation.

• Fracture modeling: fractures are populated using the geo-mechanical facies and the driver indicator generated by the stress inversion method. Additionally, the fractures are constrained utilizing in-situ stress tensor to evaluate the state of criticality [1] of natural fractures. Permeability calculations are based on the cubic law function.

2.1 3D structural deformation model

The seismic discontinuity clustering analysis is a technique that is regularly deployed for the spatial analysis of variations in the lineament’s trends, giving valuable clues on the complexity of the fracture network at seismic resolution scale. Advanced 3D seismic discontinuity attributes such as curvature and curvature gradients have long been used in the detection of subsurface structural lineaments [4].

Combined with deformation modeling techniques, seismic discontinuities attributes can be a supportive tool to explore the impact of faulting or folding deformation on fracture spatial distribution [5]. Furthermore, special attributes such as tensor/SO-semblance/Dip can reveal lineament patterns that show congruence with known faults and fracture orientations. This can be valuable for adequate identification and modeling of different tectonic environments.

2.2 Intrinsic mechanical properties model

The modeling of mechanical properties applies different algorithms such as sequential Gaussian simulation (SGS) to represent the elastic and rock strength properties defined in the wellbore scale analysis. The 3D stress analysis involves calculating stresses that represent the pre-production conditions throughout the reservoir and its surroundings. Due to the complex variation of structure and properties within the model, the stress equilibrium must be solved numerically. Geo-mechanical simulation software uses a finite element method to determine the required solution, producing a 3D stress tensor whose magnitude and orientation vary laterally and vertically. The model uses the 3D structure and the rock mechanical properties together with the loads that govern stresses (gravitational, pore pressures and boundary conditions) to simulate the initial stress state of the field.

3D elastic and strength properties as well as stress magnitudes are inputs to compute 3D discrete geo-mechanical facies, which represent mechanical units within the reservoir with differential response to the stress/strain. Neural network classification is applied to distribute mechanical rock types in 3D space.

The key aspect in mechanical Earth model construction, for example while building the 3D grid geo-cellular model, is to preserve the intrinsic relation between mechanical properties and petrophysical properties such as porosity and Young Modulus. This preservation is important because it will honor the distribution of geo-mechanical facies within the layering model and can be achieved by using the co-kriging interpolation algorithm [6].

3. Stochastic discrete fracture model

This section explores several techniques developed to embed a discrete DFN model into a 3D grid geo-cellular model, with emphasis on stochastic methods and parametrization and ranking of main variables such as fracture length, aperture, density or length/height ratio and flow capacity calibration for (closed loop) dynamic reservoir simulation model inversion. The integrated fracture modeling workflow used for 3D stochastic fracture model parametrization is shown in Figure 1. The workflow incorporates four main components:

1. The first component as shown in Figure 1 encompasses the fracture modeling process where main geo-mechanical drivers are constructed by using mechanical stratigraphy, paleo-deformation model, and 3D stress modeling to model the natural fractures.

2. The second component establishes the initial fracture distribution parameters, such as fracture length and geometry, fracture intensity as well as aperture distribution and intrinsic fracture permeability.

3. The third component leverages techniques of global, multi-objective optimization, such as genetic algorithm to minimize the differences between the dynamic properties such as flow capacity as measured data and the predicted data from the model [7].

4. The fourth component is a process to identify the natural fracture parameters that render consistent realizations, used for simultaneous closed loop inversion in dual porosity dual permeability numerical simulation.

Figure 1.

Workflow for generating constrained DFN model realizations and integration into the reservoir simulation model.

Fracture networks usually serve as the major pathways for fluid transport in subsurface rocks, especially if the matrix is relatively impermeable compared to the fractures. The partitioning of fluid flow within a fracture population relies on the spatial connectivity of fracture geometries and the transmissivity of individual fractures, both of which can be affected by the geo-mechanical conditions. The flow paths and fracture conductivity based on critical stress are in general defined following the Mohr Coulomb criterion [8]:

$$\text{Critical Fractures}=\left(\tau -\sigma n^{\ast }\tan \phi \right)\left(\phi \right)\ge 0$$E1

where τ, σ and φ correspond to shear stress on fracture plane, normal stress on fracture plane and the angle of internal friction, respectively. Under this concept, the orientation of the in-situ stress affects the normal and shear stresses on the fracture plane. Natural fractures optimally oriented with respect to in situ stress orientation will be close to shear failure, and therefore, more hydraulically conductive in contrast with those not optimally oriented. The comparison between full and critical fracture sets generated by the fracture modeling workflow is shown in Figure 2.

Figure 2.

The comparison between full and critical fracture sets generated with fracture modeling workflow.

To represent the continuous subsurface medium, equivalent to the DFN system, the fracture modeling workflow performs the property scale-up which allows the calculation by grid block of the fracture porosity, effective fracture permeability tensor (fracture KX, KY and KZ), and fracture block dimensions (fracture LX, LY and LZ) that are used to calculate the fracture-matrix mass transfer function (sigma, σ; [9]):

$$\sigma =4\left(\frac{KX}{L{X}^2}+\frac{KY}{L{Y}^2}+\frac{KZ}{L{Z}^2}\right)$$E2

where KX, KY and KZ represent the components of the matrix permeability tensor.

The generated fracture realizations represent the input for the dual porosity, dual permeability (DPDP) reservoir simulation model. An example of a DFN realization superimposed on fracture KX distribution is given in Figure 3, while Figure 4 depicts the vertical cross-section of the fracture model indicating the contrast in fracture intensity, populated in ductile anhydrite layers versus brittle hydrocarbon-rich layers.

Figure 3.

An example of a DFN realization superimposed on fracture KX distribution, as generated with the fracture modeling workflow.

Figure 4.

Vertical cross-section of the fracture model indicating the contrast in fracture intensity/density as populated in ductile anhydrite layers vs. brittle hydrocarbon bearing layers.

3.1 Calibration of fracture model flow capacity

The permeability-stress relationship to constrain fracture network permeability is a key component in this iterative process. The goal is to minimize the error between the computed fracture network flow capacity and measured dynamic data such as flow capacity from PTA and flowmeter production log. During the iterative process, fracture network upscaling is conducted for each model realization.

We apply a hierarchy methodology to calculate the equivalent permeability for the three media usually found in carbonate reservoirs, where the natural fractures have the major impact for the fluid-flow movement, followed by high permeability streaks (HPS), characterized by the core interpretation, and the matrix permeability. Using this hierarchy, the flow capacity for each rock media component is computed and optimized according to the reservoir dynamic response. The reservoir dynamic data used to validate quantitatively the fracture model realizations are pressure transient analysis (PTA), and production log test (PLT). Additional qualitative validation can be performed using tracers, water encroachment, cumulative production, etc. that indicate possible fluid-flow behavior due to natural fractures.

For quantitative validation or calibration, a workflow for multiple fracture realizations can be used to optimize fracture properties by reducing the difference between the flow capacity from the pressure transient analysis (KH_PTA) and the total flow capacity produced by each component of the geological model: fractures (KH_{Fract-simulate}), matrix (KH_Facim) and high permeability streaks (KH_HPS) as shown in Figure 5. The HPS and matrix flow capacity calibration steps are depicted as grouped because both components effectively represent the same rock medium.

Figure 5.

Flow capacity for each element contributing to the fluid flow.

During the calibration of fracture network, several variables are parametrized based on uncertainty analysis. These variables can be conceptually defined as fracture size, fracture shape, intrinsic fracture properties and criticality of these fractures. The variability that can be selected for these variables will depend on the geo-mechanical modeling as well as deformation model, determining for example, the fracture size distribution within a given reservoir. The differences between each fracture realization are reflected in an upscaled grid fracture property which is iteratively compared with the observed KH_PTA.

Figure 6 shows multiple realizations and the calculated error distribution. The error minimization using a genetic algorithm renders a first population used to explore the range of variability within the parameters established. After the population is created, a minimization process by iterating within the range of variables is conducted to achieve a stationary solution where, within the error tolerance of parameters, insignificant variability is expected.

Figure 6.

Multiple realization algorithm to minimize the flow capacity error.

The range of parameters can be defined by exploring the area of “best set of parameters” bracketed in Figure 6. The ranges of variables defined within the minimized solution can be limited by exploring the extreme values in that area. Extreme values are defined and utilized for the simultaneous closed-loop inversion of the reservoir simulation model, entailing the uncertainty quantification and the dynamic model calibration with computer-assisted history matching (AHM), as demonstrated in Section 4.

4. Dynamic calibration of DPDP model

The dynamic calibration of the DFN model is performed with synthetic simulation model of a field with four main stacked naturally fractured reservoirs under natural aquifer drive, denoted as Reservoir 1, 2, 3 and 4 in ascending order from top to bottom, as illustrated in Figure 7. The field is a salt-related anticlinal structure and has several predominantly normal faults with fault patterns consistent with the structural dome. The reservoirs are separated by anhydritic zones as a vertical flow barrier. Reservoirs 1 and 2, however, indicate vertical pressure communication through the faults and are, in terms of dynamic pressure response, combined into one flow unit, denoted as Reservoir 1–2.

Figure 7.

Schematic illustration of the modeled reservoir structure.

The simulation model is a 3-phase flow, black oil, dual-porosity dual-permeability (DPDP) model and has an approximate grid size of 10 million cells. In computational terms, the grid size if effectively double since in DPDP systems, the matrix and fracture media are each represented with individual grids, hydrodynamically connected via fracture-matrix mass transfer function (see Eq. (2)). To avail acceptable simulation runtimes for models of comparable size and complexity it is recommended to utilize parallel and highly scalable simulators [10], deployed on high-performance computing (HPC) platforms.

4.1 Uncertainty quantification, sensitivity and dynamic variability

In dynamic model reconciliation studies, one of the crucial tasks is to understand the process of generating multiple (fracture) representative model realizations under uncertainty (i.e., the definition of sampling distributions and uncertainty intervals) that render sufficient dynamic response variability as well as the integration of fracture models into the workflow for dynamic model updating and assisted history matching (AHM).

The fracture modeling workflow stochastically samples from probability density distribution functions (pdf’s) quantified for all principal components of the DFN model, such as fracture length, L/H ratio, orientation, fracture aperture with intrinsic permeability as well as fracture density and intensity. Following the data analysis and uncertainty quantification, the sampled distributions have assigned the types/shapes and uncertainty boundaries. Examples of pdf’s, corresponding to selected DFN parameters as modeled and sampled herewith, are given in Figure 8.

Figure 8.

Examples of probability density functions (pdf’s) corresponding to selected uncertain DFN parameters as sampled in the process of uncertainty quantification. Sampled frequency corresponds to 2-level DoE (Plackett-Burman) and 3-level DoE (Latin hypercube).

The convergence effectiveness of iterative and computationally intensive algorithms for AHM under uncertainty depends largely on the complexity and volume of the parameter sampling domain. Particularly with the analyses of real life-like reservoirs, represented by multi-million, high resolution simulation grids, the screening, optimization and dimensionality reduction of the prior sample space, based on a thorough engineering understanding and sensitivity analysis is not only relevant but also mandatory.

The local, also referred to as univariate, sensitivity analysis uses design of experiments (DoE) technique One-Variable-at-the-Time (OVAT) [11]. It is based on calculating the effect on the model output of small perturbations around baseline parameter. OVAT scenarios start at baseline values and then each parameter is successively varied over its range with other parameters held constant. The OVAT design with N parameters and 3 levels (min/base/max) renders 3*2(N + 1) scenarios. The local sensitivity analysis has the advantage of being highly intuitive in its interpretation. Tornado charts are routinely used to rank the parameters and quantitatively visualize their importance.

A tornado chart corresponding to reservoir pressure, as a targeted dynamic variable is shown in Figure 9. In this example, the chart confirms that the upper 50% of parameter domain can be captured by the 12 most influential and impactful parameters. This significantly simplifies the stochastic sampling process, improves the convergence of AHM process and enhances the understanding of reservoir behavior and its dynamics:

• The aquifer strength, aquifer response time (modeled by using aquifer porosity or pore volume and aquifer permeability multipliers, respectively), and fracture permeability of reservoirs 1 and 2 in communication represent the top three most influential parameters for reservoir connectivity and the corresponding pressure behavior. This confirms that the field is under strong natural aquifer drive and that the oil legs of reservoirs 1 and 2 are well connected with the aquifer through the fracture network. Varying these parameters manifests a positive correlation with reservoir pressure behavior.

• Three parameterized variables render negative correlation (or inverse proportion) with reservoir pressure behavior: the gas-oil contact, fracture permeability of reservoir 3 and correlation length in vertical variogram direction of the DFN model. This can be interpreted as the connectivity of reservoir 3 with the aquifer through the fractures is less (e.g. due to non-conductive, cemented fractures) and the vertical pressure communication through fractures is less pronounced. Increasing the value for any of these parameters will result in reduction of reservoir pressure and vice versa.

• Seven out of the twelve most influential parameters for reservoir pressure behavior correspond to fracture properties. This confirms the spatial, geometrical and petrophysical properties of the DFN model as the highly uncertain set of reservoir properties and the main focus in dynamic model calibration.

Figure 9.

Tornado charts corresponding to reservoir pressure: Upper 50% of the parameter domain can be captured by the 12 most impactful parameters.

To reduce the risk of anchoring the baseline bias in local (OVAT) analysis one may alternatively resort to the global or generalized sensitivity analysis (GSA) [12, 13] that explores the input parameters space across its range of variation and then quantify the input parameter importance by characterizing the model response surface. It is worth mentioning though that the interpretation of GSA results may be more challenging and less intuitive.

The sensitivity analyses are followed by dynamic model variability assessment using multi-level DoE techniques. The Plackett-Burman (P-B) design, a standard two-level screening method [11] and one of the most used fractional factorial designs to sample the envelope of parameters’ uncertainty domain. This design can be used for studying up to k = (N-1)/(L-1) factors, where L is the number of levels and N (a multiple of four) is the number of simulations. In P-B design, the importance of all main effects appears at the same precision.

The Latin Hypercube (LHC), a standard three-level screening design is used to refine the sampling of uncertainty domain within the envelope of parameters’ tolerances. The LHC of n runs for m parameters is an n x m matrix, where each column consists of n equally spaced intervals, sampled by the Monte Carlo method [11]. The overall objective of the dynamic variability analyses is to design and validate the model uncertainty sampling scheme that brackets the spread (variance) of the observed/measured data with an uncertainty envelope and uniformly samples the uncertainty space within the envelope.

Figure 10 illustrates the outcomes of dynamic variability sampling: the observed/measured field pressure data are shown as green dots, the simulated uncertainty envelope is depicted with red and blue lines, while the cyan-colored lines correspond to simulated refined sampling within the uncertainty envelope. The two green dots shown at the end of the simulated timeline represent pressure data, measured after a prolonged period of the field shut-in. Since the reservoir is under aquifer drive, the field-re-pressurization is strongly attributed to aquifer strength, which validates the high rank of corresponding aquifer attributes on the tornado chart (Figure 9). Moreover, the narrowing of uncertainty bound toward the end of simulated timeline, reflected in dynamic response of all OVAT scenarios (P-B and LHC), further validates the modeled spatial parameterization of aquifer strength and response time.

Figure 10.

Dynamic model variability of the field pressure: Top - results of two-level P-B screening design; bottom - results of three-level LHC screening design.

The designed uncertainty workflow is successfully capturing multi-variate complexity of the modeled system at the field level. Arguably, quantified variable uncertainty intervals represent a valid starting point for rigorous AHM dynamic model calibration.

4.2 Closed loop DPDP model inversion

The workflow for simultaneous closed loop inversion of DFN simulation models [14, 15] is schematically depicted in Figure 11. It is initialized with rigorous multi-variate uncertainty quantification (UQ), where probability distribution functions with associated tolerances of sampling intervals for all relevant modeling parameters and variables are defined in close collaboration between the geo-modeler, fracture modeler and simulation engineer.

Figure 11.

Schematic flowchart rendering of the workflow for simultaneous closed loop inversion of fracture simulation models.

The inversion workflow automatically interfaces between applications for geo-model and fracture model update, the full-physics, fluid-flow reservoir simulator and the global, multi-objective optimizer, which communicates with model update applications via embedded geo- and fracture modeling workflows. The simulation deck combines the full set of observed (historical) well data, ranging from pressure data (e.g., flowing bottom-hole, static), production rates (e.g., oil, water, gas), injection rates (e.g., water, gas) as well as well pressure or fluid-flow profiles (RFT/PLT). The observed data are used to define the multi-objective global misfit objective function (OF) (see Eqs. (3)-(5)), per specific well and per specific simulation time-step. In the model at hand, the global OF incorporates misfit terms for the well static pressure.

The optimization algorithm iteratively evaluates the misfit objective function and closes the inversion loop by generating updates of geo- and fracture models and the variables parameterized implicitly in the simulation model. The closed loop inversion is completed when enough optimization iterations are achieved to minimize the misfit OF, following an automated convergence criterion (e.g., maximum entropy criterion, [16]) or based on engineering judgment.

An example of a multi-objective stochastic optimization shown herewith is the technique of Evolutionary Strategies with Covariance Matrix Adaptation (ES-CMA), which belongs to the class of population-based, direct search methods [17]. ES-CMA is a gradient-free optimization method and requires the objective function (OF) value to determine a new search step. In the framework of Bayesian inference, the OF is defined as a least-square term:

$$OF=\left(1-\mu \right)O{F}^{Likelihood}+\mu O{F}^{Prior}$$E3

where the prior term corresponds to the prior geo-model realization and the likelihood term is evaluated with forward reservoir simulation. The prior term of the OF is usually mapped on the geo-cellular grid of the 3D geo-model. An additional weight factor μ is introduced to prioritize either contribution to the global OF. In multi-Gaussian representation, the likelihood term of the OF is generically defined as:

$$O{F}^{Likelihood}=\sum O{F}_i=\frac{{\omega }_i}{N_i}\sum _j^{N_i}{\omega }_i^j\frac{{\left(si{m}_i^j-ob{s}_i^j\right)}^2}{{\sigma }_i^2}$$E4

where terms “sim” and “obs” correspond to simulated and observed/measured data, respectively. Index i runs over the number of wells (N_i) accounted for in the misfit calculation, while index j runs over the number of time-steps used to discretize the timescale of dynamic parameter. The prior term of the OF is defined as:

$$O{F}^{Prior}=\sum _m{\left(x_m-{\overline{x}}_m\right)}^T{\left[C_x^{Prior}\right]}^{-1}\left(x_m-{\overline{x}}_m\right)$$E5

where m represents an arbitrary prior model parameter, with ${\overline{x}}_m$ as mean or expected value. The prior covariance matrix is calculated for each grid-cell of the 3D geo-cellular grid.

It is herewith worth mentioning that many other population-based optimization techniques exist that are used in state-of-the-art AHM workflows, e.g., multi-objective genetic algorithm (MOGA) [18], particle swarm optimization (PSO) [19], Markov chain Monte Carlo (McMC)-based proxy modeling [16] and ensemble-based multiple data assimilation (ES-MDA) [20].

In the given example, the joint global misfit OF function is following the definition in Eq. (4) and combines well response vectors for static well pressure, with more than 70 misfit terms. Figure 12 outlines static pressure variability for nine producing wells, represented by 75 DoE runs. Adequate variability and bracketing of historic/observed pressure data points (red dots) by simulated response profiles (blue lines) indicate a rigorous definition and sampling of the underlying parameter uncertainty space.

Figure 12.

An example of well pressure variability (prior to history matching) for nine producing wells, represented by 75 DoE (25 P-B and 50 LHC) runs. Historic data are depicted with red dots, simulation runs with blue lines.

Five variability realizations are identified as the starting point for global stochastic optimization and model inversion. Figure 13 depicts the monotonic decrease of the global misfit OF for a total number of six iterations performed, while Figure 14 indicates a significant improvement in the precision of well pressure profiles, while retaining the accuracy, relative to observed data points.

Figure 13.

A depiction of the monotonic decrease of the global misfit OF with the increasing number of iterations. The red circles represent the mean value of the misfit function per iteration, the blue lines represent the associated variance.

Figure 14.

An example of well pressure history match for nine producing wells, represented by 30 ES-CMA runs. Historic data are depicted with red dots, simulation runs with blue lines. The quality of history match is quantified through significant improvement in pressure profile precision, while retaining accuracy, relative to observed data points.

As derived from tornado chart in Figure 9, fracture properties represent highly impactful parameters for the field pressure dynamics in DPDP models. Consequently, the geologically consistent reconciliation of fracture permeability, with historically observed dynamic data, plays a vital role in calibrating spatial reservoir connectivity. Figure 15 shows maps of spatial fracture X-permeability distribution, before (left) and after (right) AHM model inversion. The locations of nine wells used in model pressure calibration are depicted as well. As indicated, the permeability updates are moderate and well-conditioned, resulting in desired fine-tuning and high predictability of the model.

Figure 15.

An example of average spatial fracture X-permeability distribution in the vicinity of nine wells depicted in Figures 12 and 14, before (prior) and after (posterior) AHM model inversion.

5. Discussion

This chapter describes a practical and structured approach to subsurface reservoir modeling in the petroleum industry, with emphasis on fractured reservoirs. It outlines the relevant aspects of multiple static and dynamic data integration and modeling processes for natural fracture systems, covering core and wellbore fracture characterization, 1D/3D mechanical Earth modeling and stochastic natural fracture prediction. The loop is closed with dynamic calibration and reconciliation of a dual-porosity dual-permeability (DPDP) reservoir simulation model of a fractured and faulted reservoir developed under natural aquifer drive with embedded discrete fracture network (DFN).

The chapter further demonstrates the application of state-of-the-art probabilistic algorithm for simulation model inversion with robust multi-variate quantification of subsurface uncertainty. The outlined approach emphasizes the value of statistical inference and parameterization using Design of Experiments to assess the boundaries of the uncertainty domain and to preserve the unbiased sampling of parameters’ probability distributions. Particularly when modeling real-life reservoirs, represented by multi-million, high resolution simulation grids, the implementation of rigorous sensitivity analyses is necessary to simplify the model by eliminating input parameters deemed insignificant. The intuitive quantitative visualization of sensitivity analysis results (e.g., by using Tornado charts) also enhances the interpretation and understanding of the impact that most influential model parameters may have on its dynamic response.

But, most importantly, the described integrated fracture modeling and dynamic calibration approach generates geologically consistent subsurface model updates and renders measurable improvements in history match quality. This improves the overall reliability and predictability of the model for field development planning and reserves management tasks. The presented workflows and practices have been further generalized and benchmarked in reservoir simulation studies of complex, large-scale reservoirs under water-flooding strategies [21].