COMPREHENSIVE WORKFLOW TO MODEL NATURALLY FRACTURED RESERVOIRS

Abstract

A method includes modeling a reservoir using a lab scale set of models and a field scale set of models. The reservoir is modeled using the lab scale set of models by scanning a sample of the reservoir into the lab scale set of models to create modeled fractures, estimating hydraulic properties of the modeled fractures, estimating multi-phase dynamic properties of the modeled fractures, and determining characteristics of a flow regime of a fluid flowing through the modeled fractures. The reservoir is modeled using the field scale set of models by modeling a discrete fracture network of the reservoir, upscaling the discrete fracture network, and calibrating the discrete fracture network. An enhanced oil recovery operation is designed and performed on the reservoir using the calibrated discrete fracture network.

Description

BACKGROUND

Hydrocarbons are located in porous reservoirs beneath the Earth's surface. Wells are drilled into the reservoir to access and produce the hydrocarbons. Many reservoirs are naturally fractured, such as carbonate reservoirs. Natural fractures in reservoirs are beneficial as they can enhance the effective permeability, stimulate secondary porosity, and improve reservoir connectivity.

However, high conductive streaks created by the natural fractures may lead to well conformance issues, premature water breakthrough, and eventually result in degraded volumetric sweep efficiency and recovery factor. The presence of natural fractures and their dynamic contribution can have a critical effect on the reservoir management especially when considering tertiary recovery alternatives such as carbon dioxide or water based enhanced oil recovery schemes.

SUMMARY

This summary is provided to introduce a selection of concepts that are further described below in the detailed description. This summary is not intended to identify key or essential features of the claimed subject matter, nor is it intended to be used as an aid in limiting the scope of the claimed subject matter.

This disclosure presents, in accordance with one or more embodiments methods and systems for modeling natural fractures in a naturally fractured reservoir. The method includes modeling a reservoir having natural fractures using a lab scale set of models and a field scale set of models. The reservoir is modeled using the lab scale set of models by scanning a sample of the reservoir into the lab scale set of models to create modeled fractures that represent the natural fractures, estimating hydraulic properties of the modeled fractures using analytical models or artificial intelligence based models, estimating multi-phase dynamic properties of the modeled fractures based on geometric properties of the modeled fractures and the estimated hydraulic properties, and determining characteristics of a flow regime of a fluid flowing through the modeled fractures using the estimated hydraulic properties and the estimated multi-phase dynamic properties. The reservoir is modeled using the field scale set of models by modeling a discrete fracture network of the reservoir using the modeled fractures, the estimated hydraulic properties, the estimated multi-phase dynamic properties, and the characteristics of the flow regime using artificial intelligence and stochastic methods to capture an anisotropy and heterogeneity of each individual modeled fracture, upscaling the discrete fracture network by determining a shape factor using a physics-informed neural network, and calibrating the discrete fracture network using machine learning algorithms to create a calibrated discrete fracture network. An enhanced oil recovery operation is designed and performed on the reservoir using the calibrated discrete fracture network.

The system includes an enhanced oil recovery system having an injection well and a production well drilled into a reservoir having natural fractures, at least one sample of the reservoir; and a computer system configured to model the reservoir using a lab scale set of models and a field scale set of models. The reservoir is modeled using the lab scale set of models by scanning a sample of the reservoir into the lab scale set of models to create modeled fractures that represent the natural fractures, estimating hydraulic properties of the modeled fractures using analytical models or artificial intelligence based models, estimating multi-phase dynamic properties of the modeled fractures based on geometric properties of the modeled fractures and the estimated hydraulic properties, and determining characteristics of a flow regime of a fluid flowing through the modeled fractures using the estimated hydraulic properties and the estimated multi-phase dynamic properties. The reservoir is modeled using the field scale set of models by modeling a discrete fracture network of the reservoir using the modeled fractures, the estimated hydraulic properties, the estimated multi-phase dynamic properties, and the characteristics of the flow regime using artificial intelligence and stochastic methods to capture an anisotropy and heterogeneity of each individual modeled fracture, upscaling the discrete fracture network by determining a shape factor using a physics-informed neural network, and calibrating the discrete fracture network using machine learning algorithms to create a calibrated discrete fracture network.

Other aspects and advantages of the claimed subject matter will be apparent from the following description and the appended claims.

BRIEF DESCRIPTION OF DRAWINGS

Specific embodiments of the disclosed technology will now be described in detail with reference to the accompanying figures. Like elements in the various figures are denoted by like reference numerals for consistency. The sizes and relative positions of elements in the drawings are not necessarily drawn to scale. For example, the shapes of various elements and angles are not necessarily drawn to scale, and some of these elements may be arbitrarily enlarged and positioned to improve drawing legibility. Further, the particular shapes of the elements as drawn are not necessarily intended to convey any information regarding the actual shape of the particular elements and have been solely selected for ease of recognition in the drawing.

FIG. 1 shows a flowchart in accordance with one or more embodiments.

FIG. 2 shows an example of a natural fracture technology suite that may be used to model the naturally fractured reservoir in accordance with one or more embodiments.

FIG. 3 shows the modeled fractures created using the FracGeo module in accordance with one or more embodiments.

FIG. 4 shows a convolutional neural network being used to estimate hydraulic properties of the modeled fractures in accordance with one or more embodiments.

FIG. 5 shows a physics-guided neural network being used to estimate hydraulic properties of the modeled fractures in accordance with one or more embodiments.

FIG. 6 shows a high-fidelity Navier-Stokes simulation being used to estimate hydraulic properties of the modeled fractures in accordance with one or more embodiments.

FIG. 7 shows a stress-dependent permeability workflow under coupled flow-normal-shear condition in accordance with one or more embodiments.

FIG. 8 shows a multi-phase simulator based on the lattice Boltzmann method in accordance with one or more embodiments.

FIG. 9 shows a multi-phase flow simulator in accordance with one or more embodiments.

FIG. 10 shows a fracture geometric property simulation upscaled into a fracture permeability map in accordance with one or more embodiments.

FIG. 11 shows a streamlined diagnostic of a DFN, used to assess the fracture connectivity and the volumetric flow velocity, in accordance with one or more embodiments.

FIG. 12 shows the coupled PINN approach used to upscale the DFN in accordance with one or more embodiments.

FIG. 13 shows a hybrid numerical model used for DFM modeling in accordance with one or more embodiments.

FIG. 14 shows a workflow for calibrating the DFN in accordance with one or more embodiments.

FIG. 15 shows an example EOR operation being performed on a reservoir in accordance with one or more embodiments.

FIG. 16 shows a computer system in accordance with one or more embodiments.

DETAILED DESCRIPTION

In the following detailed description of embodiments of the disclosure, numerous specific details are set forth in order to provide a more thorough understanding of the disclosure. However, it will be apparent to one of ordinary skill in the art that the disclosure may be practiced without these specific details. In other instances, well-known features have not been described in detail to avoid unnecessarily complicating the description.

Throughout the application, ordinal numbers (e.g., first, second, third, etc.) may be used as an adjective for an element (i.e., any noun in the application). The use of ordinal numbers is not to imply or create any particular ordering of the elements nor to limit any element to being only a single element unless expressly disclosed, such as using the terms “before”, “after”, “single”, and other such terminology. Rather, the use of ordinal numbers is to distinguish between the elements. By way of an example, a first element is distinct from a second element, and the first element may encompass more than one element and succeed (or precede) the second element in an ordering of elements.

Due to the effect of natural fractures when producing hydrocarbons from a well, it is beneficial to be able to model the natural fractures as accurately as possible so that effective enhanced oil recovery schemes can be planned and executed. However, modeling naturally fractured reservoirs is a complex task that involves fracture recognition and characterization, static and dynamic data acquisition, dynamic fracture impact on well productivity and recovery performance, building naturally fractured reservoir geological models using stochastic or deterministic approaches, upscaling, multiphase effects on matrix/fracture interaction, reservoir simulation with single-, multiple-continua, or discrete-fracture approaches, model calibration with field observations, and ultimately decision making.

Understanding the fundamental mechanisms of multiphase fluid flow in fractured reservoirs that include matrix-fracture interactions from capillarity, gravity, diffusivity, and viscous forces is crucial for enhanced oil recovery (EOR)/improved oil recovery (IOR) optimization. Bridging the gaps in scales from pore, to plug, to near wellbore, to field scales is a persisting challenge in modeling naturally fractured reservoirs. As such, the present disclosure outlines methods and systems that use a set of lab scale models and a set of field scale models that integrate with one another to model the naturally fractured reservoir.

FIG. 1 shows a flowchart in accordance with one or more embodiments. The flowchart outlines a method for modeling a naturally fractured reservoir and designing/performing an enhanced oil recovery operation on the reservoir. While the various blocks in FIG. 1 are presented and described sequentially, one of ordinary skill in the art will appreciate that some or all of the blocks may be executed in different orders, may be combined or omitted, and some or all of the blocks may be executed in parallel. Furthermore, the blocks may be performed actively or passively.

In S100, a reservoir having natural fractures is modeled using a lab scale set of models and a field scale set of models. FIG. 2 shows an example of a natural fracture technology (NAFT) suite (200) that may be used to model the naturally fractured reservoir. The NAFT suite (200) is divided into a NAFT lab (202) set of models and a NAFT field (204) set of models. The NAFT lab (202) electronically models physical samples taken from the naturally fractured reservoir and the NAFT field (204) applies the sampled portions of the reservoir to an overall model of the reservoir. S100 of FIG. 1 is described using the NAFT lab (202) set of models.

FIG. 2 shows that NAFT lab (202) includes a fracture geometry (FracGeo) module (206), a fracture static properties (FracSPro) module (208), a fracture dynamic properties (FracDPro) module (210), a fracture mechanistic modeling (FracMech) module (212), and a fracture pressure transient analysis (FracPTA) module (214).

FIG. 2 further shows that NAFT field (204) includes a fracture discrete network (FracDNet) module (216), a fracture diagnostic (FracDiag) module (218), a fracture scale-up (FracScale) module (220), a fracture pilot-scale (FracPilot) module (222), and a fracture uncertainty quantification (FracUQ) module (224).

In accordance with one or more embodiments, the NAFT lab (202) is a lab-scale mechanistic modeling tool focused on local scales of fractured formation analysis and modeling. In accordance with one or more embodiments, the objective of the NAFT lab (202) is to provide a workflow with innovative solutions to enable engineers and researchers to streamline the study of fundamental recovery mechanisms of multiphase flow in naturally fractured formation at scales ranging from plug-scale (in centimeters) to near-wellbore scale (in meters).

In accordance with one or more embodiments, NAFT lab (202) may be developed in Python and may run as a standalone tool. NAFT lab (202) may be customized and hooked up with other tools, such as an in-house or commercial reservoir simulator. NAFT lab (202) consists of the modules outlined above that can be customized and designed to streamline the workflow.

In accordance with one or more embodiments, the FracGeo module (206) creates 2D and 3D modeled fractures that mimic the natural fractures in terms of geometric properties, including mean aperture, surface roughness, tortuosity, aperture distribution, contact areas, minimum & maximum aperture, standard deviation of aperture fields, geometric flow directions, etc.

Turning back to FIG. 1, in S102, a sample of the reservoir is scanned into the lab scale set of models to create modeled fractures that represent the natural fractures. In accordance with one or more embodiments, a CT-scan or profilometer may be used to scan the sample of the reservoir and upload the scanned images into the FracGeo module (206).

The sample of the reservoir may come from a core sample taken from the reservoir, a sample taken from an outcrop of the reservoir, and/or borehole images taken while drilling a well into the reservoir. Once the images are scanned and compiled into the FracGeo module (206), the FracGeo module (206) may analyze fracture density, fracture porosity, fracture spacing, fracture orientation (characterized by rose or stereo net diagrams) using imbedded mathematical equations and measured sizes of the modeled fractures.

FIG. 3 shows the modeled fractures created using the FracGeo module (206) in accordance with one or more embodiments. Specifically, FIG. 3 shows a roughness scale (300) that shows the magnitude of roughness of the modeled fractures. Two different fractures are shown: a first modeled fracture (302) and a second modeled fracture (304). The first modeled fracture (302) is shown in 3D in the left image and in 2D in the right image at the top of FIG. 3. The second modeled fracture (304) is shown in 3D in the left image and in 2D in the right image at the bottom of FIG. 3.

The 2D versions of the modeled fractures (302, 304) are cross sections from the 3D versions of the modeled fractures (302, 304) showing different levels of roughness and tortuosity. The roughness of the first modeled fracture (302) and the second modeled fracture (304) may be measured using a profilometer, not pictured. While FIG. 3 only shows two samples of modeled fractures (302, 304), a person skilled in the art will appreciate that there may be more than two samples needed to create an accurate model of the entire relevant portion of the naturally fractured reservoir.

Turning back to FIG. 1, in S104, hydraulic properties of the modeled fractures (302, 304) are estimated using analytical models or artificial intelligence based models and the lab scale set of models. In accordance with one or more embodiments, the FracSPro module (208) is used to estimate the hydraulic properties of the modeled fractures (302, 304).

The fracture permeability is a complex function of various parameters including the fracture mechanical aperture, surface roughness, and contact areas, which are subject to alterations by the stress field acting on the fracture walls. In accordance with one or more embodiments, the FracSPro module (208) is used to estimate the hydraulic properties of the modeled fractures (302, 304), including the fracture hydraulic aperture and the effective permeability.

Several options are embedded in the FracSPro module (208) to estimate the hydraulic properties, including analytical models (e.g., in-house developed corrected cubic law), AI-based models (e.g., support vector machine: SVM, polynomial chaos expansion: PCE, artificial neural network: ANN, convolutional neural network: CNN, physics-informed neural network: PINN, physics-guided neural network: PGNN), and high-fidelity Navier-Stokes simulations. The FracSPro module (208) can also be used to estimate the stress-dependent permeability (SDK) under effective normal stress and coupled flow-normal-shear conditions, respectively.

FIG. 4 shows a convolutional neural network (CNN) which is an image-based deep learning module used to estimate hydraulic properties of the modeled fractures (302, 304) in accordance with one or more embodiments. FIG. 5 shows a physics-guided neural network (PGNN) which is a vector-based machine learning module used to estimate hydraulic properties of the modeled fractures (302, 304) in accordance with one or more embodiments. FIG. 6 shows a high-fidelity Navier-Stokes Equation (NSE) simulation being used to estimate hydraulic properties of the modeled fractures (302, 304) in accordance with one or more embodiments. The NSE simulation in FIG. 6 is used to produce datasets used to train the CNN in FIG. 4 or the PGNN in FIG. 5.

Turning to the CNN module, FIG. 4 shows three steps of the CNN being used to determine the hydraulic properties. In S400, a data set is generated from the scanned modeled fractures (302, 304). In S402, the CNN model undergoes deep training. In S404, the CNN model is evaluated.

In accordance with one or more embodiments, the CNN is trained using high-resolution fracture images, such as the scanned modeled fractures (302, 304), as inputs. Specifically, the input high-fidelity discrete rock fracture image is pre-processed to generate the input data for the deep CNN model. The input high-resolution rock fracture is formed by two lines representing the roughness and tortuousness of the fractures (see the graphs on the right of FIG. 3 representing the modeled fractures (302, 304) in 2D).

In this case, there exists a substantial imbalance between the fracture-occupied and non-fracture-occupied pixels, which could lead to difficult feature extraction during the training process. As such, a filling-in process is conducted to mark the zone bounded by the fracture lines as shown in S400 of FIG. 4. An image-coarsening process is followed, where a suitable resolution (such as 50×1000 pixel-resolution) is applied. A person skilled in the art will appreciate that the size of the coarse resolution may be changed depending on the scenario.

In accordance with one or more embodiments, the coarse resolution is found to be a compromise between accuracy and efficiency, which is sufficient for characterizing the roughed fractures without significant loss of accuracy. On the other hand, it dramatically accelerates the feature extraction and training-validating processes. Finally, the region of interest is captured with a binary image using Otsu's method, in which the fracture-occupied pixels and non-fracture-occupied pixels are represented by 1 and 0, respectively.

Further in S400, training and validation datasets each consisting of a coarse-resolution binary fracture image created through pre-processing (outlined above) are generated. The data sets include the coarse-resolution binary fracture image and their corresponding effective hydraulic aperture. In accordance with one or more embodiments, the full-physics NSE using the mixed finite element (FE) formulation within the FEniCS platform is used to compute the effective hydraulic aperture. Considering the steady-state of incompressible, Newtonian laminar flow with no gravity effects, the full-physics NSE can be given in Equations (1) and (2) below.

$\begin{array}{cc}\rho \left(\stackrel{\rightharpoonup }{u}·\nabla \stackrel{\rightharpoonup }{u}\right)=-\nabla p+\mu {\nabla }^2\stackrel{\rightharpoonup }{u}& \mathrm{Equation}\left(1\right)\end{array}$ $\begin{array}{cc}\nabla ·\stackrel{\rightharpoonup }{u}=0& \mathrm{Equation}\left(2\right)\end{array}$

In the above equations, is velocity vector, ρ is density, p is pressure, and μ is viscosity.

In S402, the training datasets generated in S400, are used to train the deep CNN model. A person skilled in the art will appreciate that any network architecture may be used to train the datasets, however, the suitable network architecture dramatically affects the quality of the CNN model predictions. As such, a deep CNN with a total of 22 layers, including 20 hidden layers, may be used herein to obtain the highest quality CNN model predictions. The training process seeks the minimization of the loss corresponding to the error between the ground truth and CNN predictions using an optimizer function. The optimizer updates the weights and biases of each layer during the training process to improve the prediction accuracy using a back-propagation algorithm.

Attention should be given to analyze the convergence behavior of validation samples. It often occurs that the trained model performs well for the training samples but produces poor predictions for the validation samples. In this work, a trial-error analysis of the coupling training-validating process is applied to seek optimum weights and biases until the trained model achieves an accuracy of 90% or above for both training and validation datasets. In S404, the effective hydraulic apertures for raw cases is estimated using the trained deep CNN model from S402.

Turning to the PGNN model, FIG. 5 shows input layers (500), hidden layers (502), and output layers (504) including the fracture hydraulic properties. The PGNN adds a physical inconsistency term (506) to the loss function in order to guide the neural network to predict physically consistent fracture hydraulic aperture. In accordance with one or more embodiments the input layers (500) may include mean aperture, relative roughness, tortuosity, ratio of minimum to mean apertures, Reynolds number, etc. In accordance with one or more embodiments the output layers (504) may include the fracture hydraulic aperture.

Turning to the high-fidelity Navier-Stokes simulation, FIG. 6 shows the magnitude of velocity profiles (600) of flow streamlines. In accordance with one or more embodiments, the full-physics NSE (shown above in Equations (1) and (2)) offer the most accurate approach for estimating the hydraulic properties of rock fractures theoretically. The full-physics NSE is created by considering the steady-state of incompressible, Newtonian fluid within the laminar flow regime and with negligible gravity effects. In accordance with one or more embodiments, the flow rate (denoted as Q) may be calculated by integrating the velocity across the outlet, see Equation (3) below.

$\begin{array}{cc}Q={\int }_0^w{\int }_0^a\left({\stackrel{\rightharpoonup }{u}}_{\mathrm{outlet}}▯\stackrel{\rightharpoonup }{n}\right)\mathrm{dwda}& \mathrm{Equation}\left(3\right)\end{array}$

In Equation (3). W is the width of the rock fracture; a is the aperture of rock fracture; is the velocity at the outlet; and is the unit vector normal to the outlet. In accordance with one or more embodiments, fluid flow in the rock fractures may also be described by Darcy's law within the assumption of the laminar flow regime. Combined with Cubic law, fracture hydraulic aperture and corresponding permeability (denoted as k) are determined using Equations (4) and (5) below:

$\begin{array}{cc}{a}_h={\left(\frac{12Q\mu }{w\nabla P}\right)}^{1/3}& \mathrm{Equation}\left(4\right)\end{array}$ $\begin{array}{cc}k=\frac{a_h^2}{12}={\left(12\right)}^{1/3}{\left(\frac{Q\mu }{w\nabla P}\right)}^{2/3}& \mathrm{Equation}\left(5\right)\end{array}$

In Equation (4) and (5) above, a_h is the fracture hydraulic aperture, and ∇P is the pressure gradient across the flow direction. Equation (6), below, shows the stress-dependent permeability model that quantifies the relation between fracture permeability and effective normal stress. In Equation (6), k₀ and a₀ are the fracture permeability (K) and aperture (a) at initial normal stress (σ_n,e,0). ε is the root mean square value of the height distribution.

$\begin{array}{cc}k/k_0={\left[1-\left(\sqrt{2}\epsilon /a_0\right)\ln \left({\sigma }_{n,e}/{\sigma }_{n,e,0}\right)\right]}^3& \mathrm{Equation}\left(6\right)\end{array}$

FIG. 7 shows a stress-dependent permeability (SDK) workflow under coupled flow-normal-shear condition. The SDK workflow includes four major steps. In S700, normal closure under normal stress is calculated using composite topography and a B-S normal closure model (also known as a theoretical normal closure model).

Combined with the effective normal stress, the B-S normal closure model is represented by Equation (7) below, where σ_n,e=effective normal stress as input; σ_n=normal closure displacement as output; η=total number of local maxima per unit area; ψ=mean of tangential stress correction factor; ψ=mean of elastic constant; β^1/2=mean of square root of curvature term; z=height of local maximum in FIG. 4 (b); d₀=distance between reference planes at σ_n,e=0; ϕ(z)=probability density function for the height of local maxima on the composite topography.

$\begin{array}{cc}{\sigma }_{n,e}=\frac{4}{3}\eta \langle \psi \rangle \langle E^{\prime }\rangle \langle {\beta }^{1/2}\rangle {\int }_{d_0-\delta }^{\infty }{\left(z-d_0+{\delta }_n\right)}^{3/2}\varphi \left(z\right)\mathrm{dz}& \mathrm{Equation}\left(7\right)\end{array}$

In S702, normal displacement caused by shearing is calculated using a Joint Roughness Coefficient (JRC) mobilization model and a B-C shear dilation model (also known as a shear dilation model). The B-C shear dilation model characterizes the relationship between normal and shear displacements. Shear-induced normal displacement (closure or dilation) occurs due to the roughness of the fracture surface. In accordance with one or more embodiments, the B-C shear dilation model could be written as Equation (8) below. α_mob in Equation (8) is given by Equation (9) below.

$\begin{array}{cc}{u}_n={\delta }_s\mathrm{▯tan}{\alpha }_{\mathrm{mob}}& \mathrm{Equation}\left(8\right)\end{array}$ $\begin{array}{cc}{\alpha }_{\mathrm{mob}}=\frac{1}{M}▯{\mathrm{JRC}}_{\mathrm{mob}}\mathrm{▯log}\left(\mathrm{JCS}/\sigma \right)& \mathrm{Equation}\left(9\right)\end{array}$

In Equations (8) and (9), un is the normal displacement corresponding to the shear displacement δ_s, α_mob is the mobilization dilation angle, M is the damage coefficient dependent on normal stress (with values of 1 and 2 corresponding to low and high normal stress, respectively), JRC_mob is the mobilized JRC value, and JCS is the fracture wall compression strength.

In S704, the high-fidelity Navier-Stokes simulations are conducted in the updated 3D representation of the modeled fracture (302, 304) using a full physics NSE simulation. In S706, corresponding fracture permeability under specific flow-normal-shear-condition is calculated based on Darcy's law and Cubic's law. In accordance with one or more embodiments, S700-S706 are sequentially iterated with a loop until the final state of stress is reached.

Turning back to FIG. 1, in S106, multi-phase dynamic properties of the modeled fractures (302, 304) are estimated, using the lab scale set of models, based on geometric properties of the modeled fractures and the estimated hydraulic properties. In accordance with one or more embodiments, the FracDPro module (210) uses the geometric properties uploaded/determined from the FracGeo module (206) and the hydraulic properties determined from the FracSPro module (208) to estimate multi-phase dynamic properties of the modeled fractures (302, 304).

In accordance with one or more embodiments, the FracDPro module (210) estimates the multi-phase dynamic properties of the modeled fractures (302, 304) using the fracture geometry, the rock/fluid wettability, the fracture relative permeability, fractional flow, entrapped saturations, capillary pressure, and the Leverett J-function corresponding to different rock-fracture types. The following options are used in the FracDPro module (210) to perform the above estimation: AI-based models (SVM, PCE, ANN, CNN, PINN, PGNN), high-fidelity multi-phase simulations based on Navier-Stokes, lattice Boltzmann method, modified Reynolds equation, and many others.

FIG. 8 shows a multi-phase simulator (800) based on the lattice Boltzmann method in accordance with one or more embodiments. The simulation shown in FIG. 8 represents the dynamic properties of relative permeability, fractional flow, entrapped saturations, capillary pressure, etc.

Turning back to FIG. 1, in S108, characteristics of a flow regime of a fluid flowing through the modeled fractures (302, 304) is determined using the lab scale set of models, the estimated hydraulic properties, and the estimated multi-phase dynamic properties. In accordance with one or more embodiments, the FracMech module (212) determines the characteristics of the flow regime of a fluid flowing through the modeled fractures (302, 304) using the hydraulic properties estimated using the FracSPro module (208) and the multi-phase dynamic properties estimated using the FracDPro module (210).

In accordance with one or more embodiments, the FracMech module (212) enables the functionalities to conduct mechanistic modeling for single and multi-phase flow in fractured systems, such as a naturally fractured reservoir. Simulations can be conducted using unstructured discrete-fracture-modeling (DFM) to study fundamental recovery mechanisms, such as spontaneous/forced imbibition, gravity drainage, gravity re-infiltration, molecular diffusion, among others.

Several numerical methods are used in this module based on DFM, including two-point flux approximation (TPFA), multi-point flux approximation (MPFA), Galerkin finite element (FE), mixed FE, mimetic finite difference (MFD), hybrid approaches (e.g., hybrid of mixed FE and discontinuous Galerkin). On the backend, the FracMech module (212) may be hooked up with a reservoir simulator, such as the MATLAB-based public domain simulator, the MATLAB reservoir simulation toolbox (MRST), or other commercial & in-house simulators, such as Eclipse and GigaPOWERS™.

The FracMech module (212) may also be used to determine and verify the shape factor for dual-porosity-dual-permeability models discussed below in the FracScale module (220). With the presence of uncertainties, the FracUQ module (224) can also be used with the FracMech module (212).

FIG. 9 shows a multi-phase flow simulator in accordance with one or more embodiments. Specifically, FIG. 9 shows simulations of the multi-phase flow in the matrix and the modeled fractures (302, 304). The simulation on the left shows the matrix flow (900) and the simulation on the right shows the fracture flow (902).

In accordance with one or more embodiments, the matrix flow (900) is the flow of fluid through the pores of the reservoir that are not part of the natural fractures. The fracture flow (902) is the flow of the fluid through the naturally occurring fractures in the reservoir that are represented in the NAFT suite (200) as the modeled fractures (302, 304).

In accordance with one or more embodiments, the FracPTA module (214) is used to calibrate the modeled fractures (302, 304) and the properties of the modeled fractures (302, 304) that were determined in the other modules of the NAFT lab (202). The FracPTA module (214) may also show how the modeled fractures (302, 304) contribute to flow dynamics, such as well skin, fracture conductivity, connectivity, and density, given the well testing data.

In accordance with one or more embodiments and in the presence of well-testing data, the FracPTA module (214) can provide useful calibration and insights about the fractures' contribution to flow dynamics, including well skin, fracture conductivity, connectivity, and density. Here the flow dynamics include estimating the well skin, fracture conductivity, and density using the well testing data.

In accordance with one or more embodiments, the pressure transient analkysis performed in the FracPTA module (214) may be used to determine well deliverability, characterize formation damage and other sources of skin effect, identify produced fluids and determine their respective volume ratio, measure reservoir pressure and temperature, obtain representative fluid samples suitable for PVT analysis, evaluate completion efficiency, and evaluate workover or stimulation treatments. In further embodiments, descriptive reservoir tests are conducted to assess reservoir extent and geometry, determine hydraulic communication between wells, characterize reservoir heterogeneities, and evaluate reservoir parameters.

Like the FracMech module (212), the FracPTA module (214) may be driven by a reservoir simulator using single-phase and multi-phase near-wellbore radial and Cartesian models. With the presence of uncertainties, the FracUQ module (224) can also be used with the FracPTA module (214).

Turning back to FIG. 2, the NAFT field (204) is the set of modules targeting the characterization and modeling of fractured reservoirs at the field scale. In accordance with one or more embodiments, the focus in NAFT field (204) is to consolidate the lab-scale mechanistic models developed in NAFT lab (202) to bridge the gap between the lab scale and the field scale.

Fracture geometric properties, such as surface roughness and aperture variations, are considered at the field scale. This consideration results in capturing the anisotropy in the fracture hydraulic properties, which is often ignored in common practices. In accordance with one or more embodiments, the objective of NAFT field (204) is to develop a field-scale workflow to capture the significant governing recovery mechanisms with accuracy and affordable computational resources.

Dual porosity/dual permeability (DPDK) models may be used to model flow at the field scale. Advanced numerical methods based on discrete-fractured-modeling (DFM) with uncertainty quantifications may be used to upscale and guide the DPDK model. In accordance with one or more embodiments, NAFT field (204) consists of several integrated modules that may be integrated within existing reservoir modeling platforms, such as Petrel.

Turning back to FIG. 1, in S110, a discrete fracture network (DFN) of the reservoir is modeled using the field scale set of models, the modeled fractures (302, 304), the estimated hydraulic properties, the estimated multi-phase dynamic properties, and the characteristics of the flow regime using artificial intelligence and stochastic methods to capture an anisotropy and heterogeneity of each individual modeled fracture. In accordance with one or more embodiments, the FracDNet module (216) uses the modeled fractures (302, 304) and their associated hydraulic properties, multi-phase dynamic properties, and characteristics of the flow regime, all of which were determined using the modules in the NAFT lab (202), to model the DFN.

In accordance with one or more embodiments, the FracDNet module (216) provides a DEN approach based on AI & stochastic methods that successfully captures the variation of the geometric properties of the natural fractures in the reservoir. This approach allows capturing the anisotropy and heterogeneity at the individual fracture level. Accordingly, the effective hydraulic properties could be locally upscaled using the FracSPro module (208). FIG. 10 shows a fracture geometric property simulation (1000) upscaled into a fracture permeability map (1002) in accordance with one or more embodiments.

In accordance with one or more embodiments, traditional approaches to DFN do not consider the heterogeneity of fracture aperture. Specifically, traditional approaches assume a constant effective property for the whole individual fracture. This approach is not representative of natural fractures that often exhibit heterogeneity in their petrophysical properties which are, to some extent, similar to the matrix. Thus, the traditional approach that often uses a homogeneous and uniform planar (2D plan or ellipse) to represent a fracture is not adequate. As such, the FracSPro module (208) is based on a stochastic approach to capture the variation of the geometric properties of fractures. The proposed model considers the fracture heterogeneity and represents it with a 2D grid that captures its heterogeneity, porosity distribution, and anisotropy.

In accordance with one or more embodiments, obtaining the detailed information of large-scale fractures in the subsurface, such as roughness, aperture variation, is impossible. The available data is laboratory measurements from fractured cores or wellbore imaging at core-scale. This laboratory data is then populated at the field-scale based on a stochastic approach. Generally, the multi-point statistic (MPS) method is used to estimate local characteristics within the fracture plane. The MPS method can produce realistic and continuous features but most likely fail to condition the “hard data” (measured data). Recent advancements in image-based machine learning models (include CNN, UNet, GAN and its variants) could improve process of conditioning with the “hard data”. Therefore, the present disclosure outlines a hybrid approach. On one hand, the FracSPro module (208) takes advantage of MPS in generating realistic and continuous geological features, on the other hand, the FracSPro module (208) adopts image-based machine learning models to condition the “hard data”.

In accordance with one or more embodiments, the FracDiag module (218) generates multiple realizations of the DEN by quantifying and ranking the static and dynamic properties of the DFN. In accordance with one or more embodiments, the static properties include the modeled fracture's (302, 304) geometric connectivity, the heterogeneity based on the static Lorenz coefficient, dynamic properties (such as dynamic connectivity), and the heterogeneity based on a F-Phi curve.

In accordance with one or more embodiments, the FracDiag module (218) combines mimetic finite difference (MFD) and discrete-fracture model (DFM), for quick diagnostics in the naturally fractured reservoir. The use of MFD significantly improves the general applicability by enabling polygonal grids of any shape with full-tensor permeability. The implementation of streamline tracking significantly enhances the computation efficiency.

FIG. 11 shows a streamlined diagnostic of a DFN, used to assess the fracture connectivity and the volumetric flow velocity, in accordance with one or more embodiments. In accordance with one or more embodiments, FIG. 11 shows the DFN broken up into a plurality of streamlines (1100). The streamlines (1100) represent the velocity and direction of the flid flowing through the modeled fractures (302, 304) of the DFN. The time of flight (ToF) is based on particle tracking within the streamlines (1100) and can be used to rank the different DFN realizations. The dynamic Lorenz coefficient and the F-Phi curve can capture the overall heterogeneity of the model.

Turning back to FIG. 1, in S112, the DFN is upscaled by determining a shape factor using a physics-informed neural network. In accordance with one or more embodiments, the FracScale module (220) upscales the DFN that was created by the FracDNet module (216) and further processed by the FracDiag module (218). Conventional upscaling techniques have major limitations. For example, upscaling of a naturally fractured formation is done using static-based methods, such as Oda's method. Furthermore, transfer functions that are used in the context of dual-porosity/dual-permeability (DPDK) methods require the determination of the shape factor to describe the matrix-fracture interaction.

The shape-factor calculation methods are based on idealized assumptions, such as the sugar-cube configuration and single-phase flow. The shape factor is not well determined for cases with unstructured fracture configuration and variable fracture properties, such as those that exist in naturally fractured carbonate reservoirs. Furthermore, the flow-based upscaling methods that use high-resolution numerical methods are often too complex and time-consuming.

As such, the FracScale module (220) determines the shape factor using an algorithm based on a physics-informed neural network (PINN). The advantage of PINN over the traditional numerical methods is that it does not need to create a simulation mesh (PINN is a mesh-less method). On the other hand, the PINN approach is based on automatic differentiation, which is not convenient to be applied for systems exhibiting discontinuities in the spatial properties, such as permeabilities.

To overcome this challenge, the algorithm has a coupled approach in which the DEN is separated into the matrix of the reservoir and the fractures in the reservoir. FIG. 12 shows the coupled PINN approach used to upscale the DFN in accordance with one or more embodiments.

As can be seen in FIG. 12, the DFN (1200) is separated into the fracture media (1202) and the matrix media (1204). A separate PINN is used to solve Darcy's law (1206) and the continuity equation (1208) within each medium. The coupling between the two media is achieved by imposing the continuity of matrix-fracture flux (1210), referred to as the transfer function. The computed upscaled properties include the matrix permeability, fracture permeability, and the shape factor. In accordance with one or more embodiments, these properties may be exported in formats suitable to be loaded into industrial simulators such as ECLIPSE and GigaPOWERS.

In accordance with one or more embodiments, the FracPilot module (222) is used to capture localized flow mechanisms at the matrix-fracture interface. Homogenized models, such as DPDK, are supported by industrial simulators and convenient to perform simulations at the field scale. However, it is not suitable to use the homogenized models to capture localized flow mechanisms at the matrix-fracture interface (such as capillarity, buoyancy, and diffusion) which could are crucial properties to model when performing shale development, EOR, and CO2 sequestration. The DFM, on the other hand, is more accurate as the fractures are modeled explicitly in the system, but they are computationally prohibitive to be applied at the field scale.

The FracPilot module (222) is designed to leverage the state-of-the-art development in advanced numerical methods to perform high-resolution DFM simulations at the pilot scale, thus enabling accurate assessment for localized mechanisms reflecting the matrix-fracture interaction. The results may then be used to calibrate the DPDK model. Several numerical methods may be used in this module to perform the DFM simulations, such as two-point flux approximation (TPFA), multi-point flux approximation (MPFA), Galerkin finite element (FE), mixed FE, control-volume FE (CVFE), mimetic finite difference (MFD), hybrid approaches (e.g., hybrid of mixed FE and discontinuous Galerkin), etc.

FIG. 13 shows a hybrid numerical model used for DFM modeling in accordance with one or more embodiments. For a fractured reservoir, the DEN is built based on seismic, outcrop, etc. Once DEN is confirmed, grids are generated as shown in FIG. 13. Then, multi-phase flow simulations are performed on these generated finite-element-based grids.

Turning back to FIG. 1, in S114, the DEN is calibrated using machine learning algorithms and the field scale set of models to create a calibrated DFN. In accordance with one or more embodiments, the FracUQ module (224) is used to calibrate the DFN. The FracUQ module (224) combines Long Short-Term Memory (LSTM) and advanced Bayesian Markov Chain Monte Carlo (MCMC) algorithms for assisted reservoir model calibration.

FIG. 14 shows a workflow for calibrating the DFN in accordance with one or more embodiments. In S1400, the uncertainty parameters are determined. In S1402, the prior distributions of the uncertainty parameters is determined. In accordance with one or more embodiments, the uncertainty parameters are assessed based on the data and existing knowledge. Each parameter should be described with its range of variability (upper and lower bounds) and distribution. If there is not enough information, uniform distribution may be assigned. This step is important to provide initial prior distribution for the history-matching process. Revisiting this step may be needed to readjust the range of parameters if the model predictions do not enclose the observed data.

In S1404, the prior distributions of the uncertainty parameters are used to construct a LSTM surrogate module to improve the computational efficiency in which Bayesian optimization is used to automate the tuning of LSTM hyperparameters. In accordance with one or more embodiments, a Latin hypercube design of experiments are used to construct the LSTM surrogate module. The LSTM is subsequently trained and verified.

In S1406, a Bayesian MCMC is run using the LSTM surrogate to obtain the posterior prediction of the uncertainty parameters. In accordance with one or more embodiments, the application of the Bayesian MCMC follows the algorithm of Affine Invariant Ensemble. The likelihood calculations are performed by running the LSTM network and evaluating the mismatch with the data. The initial priors are fed into the LSTM. A good quality matching is achieved when the mean solutions are close to the data.

In accordance with one or more embodiments, a quality assessment is performed to ensure an acceptable surrogate model and good history matching. The process includes the evaluation of the discrepancy between the data and output produced from the posterior distribution. As such, the Root Mean Square Error (RMSE) is determined between the observation and the predictions. A threshold of RMSE=0.003 is used to accept or reject the match. If the accuracy is poor, LSTM network needs to be reevaluated. If the quality is still poor after reevaluation, the workflow recommends revisiting the initial prior distributions in S1402.

In accordance with one or more embodiments, Latin Hypercube Sampling (LHS) is a way of generating random samples of parameter values. LHS may drastically reduce the number of runs necessary to achieve a reasonably accurate result. In other words, the implementation of LHS ensures that the generated datasets are distributed in a space-filling manner instead of clustering distribution. Thus, it could significantly reduce the number of runs required to reach a reasonably accurate result. Low-fidelity model refers to a proxy or surrogate model, which is efficient but may lose some accuracy. For example, the LSTM used herein is the low-fidelity model.

In accordance with one or more embodiments, the MCMC algorithms may include Metropolis-Hastings (MH), Adaptive Metropolis (AMH), Hamiltonian Monte Carlo algorithm (HMC), or Affine invariant ensemble (AIES) algorithms. However, the first three algorithms require extensive tuning to improve the performance of the convergence process. AIES may be used to overcome this issue. The advantage of the AIES algorithm is that the target distribution remains invariant. In accordance with one or more embodiments, the stretch move is used to achieve affine invariance. The new candidate is generated by the following equation:

$\begin{array}{cc}{\theta }_i^{*}={\theta }_i^{\left(t\right)}+Z\left({\theta }_j^{\left(t\right)}-{\theta }_i^{\left(t\right)}\right),\mathrm{where}Z~p\left(z\right)=\{\begin{array}{cc}\frac{1}{\sqrt{z}\left(2\sqrt{a}-\frac{2}{\sqrt{a}}\right)}& \mathrm{if}z\in \left(1/a,a\right),\\ 0& \mathrm{otherwise}.\end{array}& \mathrm{Equation}\left(10\right)\end{array}$

The distribution p(z) may be controlled by only one tuning parameter a>1. This tuning parameter may be set to 2. The following process (Algorithm 1 below) is used to perform the iteration θ^(r)→θ^(t+1), using one stretch move per walker, where L is the number of walkers.

Algorithm 1: Affine invariant ensemble algorithm:
for k = 1, . . . , L
Choose θj ∈ θ[k](t) at random
Generate θ* = θj + Z(θk(t) − θj), all Z choices independent
Accept                    ,                                                                  set                        ⁢                                                                                                                            ⁢                                                  θ                          k                                                      (                                                          t                              +                              1                                                        )                                                                                              =                                              θ                        *                                                              ,                                          with                      ⁢                                            probability                      ⁢                                            min                      ⁢                                              {                                                  1                          ,                                                                                    z                                                                                                N                                                                      i                                    ⁢                                    n                                                                                                  -                                1                                                                                      ⁢                                                                                          P                                ⁡                                (                                                                                                      θ                                    i                                    *                                                                    |                                  Y                                                                )                                                                                            P                                ⁡                                (                                                                                                      θ                                    i                                                                          (                                      t                                      )                                                                                                        |                                  Y                                                                )                                                                                                                                    }
Otherwise reject, set Θk(t+1) = Θk(t)
end

In accordance with one or more embodiments, quality assessment is used to ensure an acceptable surrogate model and good history matching. The process includes the evaluation of the discrepancy between the observation data and output produced from the posterior distribution. In particular, the Root Mean Square Error (RMSE) is calculated between the observation and the predictions. A threshold of RMSE=0.003 is used to accept or reject the match. If the accuracy is poor, we need to re-evaluate the LSTM network. If the quality is still poor after re-evaluation of LSTM proxy or surrogate, the workflow recommends revisiting the initial prior ranges (i.e., S1402).

In S1408, the posterior prediction of the uncertainty parameters are uploaded into a high-fidelity model simulation. Specifically, the high-fidelity models (simulations) are then used to ensure accurate and physically consistent history matching. In accordance with one or more embodiments, the high-fidelity model refers to the full-physics simulation model (herein, the full-physics simulation model is run by CMG). This step includes a loop to repeat the LSTM construction process (i.e., S1404) if the accuracy is below a required threshold.

In S1410, observation mismatch is evaluated. Specifically, the difference between the posterior prediction by the Bayesian inversion and the real data is analyzed. Another inspection may be applied if the history match is not satisfactory, which may require revisiting the distribution of the uncertainty parameters and their ranges (i.e., go back to S1402).

In S1412, the final posterior distribution of the uncertainty parameters is determined. Specifically, after going through S1400-S1410, the posterior distribution of the uncertainty parameters is expected to be with narrow ranges and shows good match with the real data. The one with narrow ranges are referred to as the final posterior distribution.

Turning back to FIG. 1, in S116, an EOR operation is designed and performed on the reservoir using the calibrated discrete fracture network. FIG. 15 shows an example EOR operation being performed on a reservoir (1502). Specifically, a water-alternating-gas (WAG) operation is being performed. A person skilled in the art will appreciate that the operation shown in FIG. 15 is used for example purposed only and any type and design of EOR operation may be used without departing from the scope of the disclosure herein.

Turning to FIG. 15, a fluid pumping system (1512) contains an injection system (1514) and a production system (1516). The injection system (1514) includes an injection well (1506) extending into a formation (1504) under the earth's surface (1522) and penetrating through a hydrocarbon reservoir (1502). The production system (1516) includes a production well (1508) extending into a formation (1504) under the earth's surface (1522) and penetrating through a hydrocarbon reservoir (1502).

In accordance with one or more embodiments, the hydrocarbon reservoir (1502) is naturally fractured. Thus, FIG. 15 shows a plurality of fractures (1501) extending through the hydrocarbon reservoir (1502). The calibrated discrete fracture network modeled in S100 may be used to determine the optimal location of the injection well (1506). Furthermore, the calibrated discrete fracture network modeled in S100 may be used to determine how the injection well (1506) is completed such that the injection of the treatment fluid into the hydrocarbon reservoir (1502) is done in the most optimal way to recover all of the hydrocarbons located in both the matrix of the hydrocarbon reservoir (1502) and the fractures (1501).

In accordance with one or more embodiments, the injection system (1514) is located at a distance across the hydrocarbon reservoir (1502) from the production system (1516). The area of the hydrocarbon reservoir (1502) containing the injection well (1506) and the production well (1508) and a plurality of WAG zones (1550, 1548, 1542, 1538, 1540, 1544) is considered the well environment. The plurality of WAG zones (1550, 1548, 1542, 1538, 1540, 1544) extends between the injection system (1514) and the production system (1516).

When injection well (1506) and production well (1508) are vertical wells, each WAG zone (1550, 1548, 1542, 1538, 1540, 1544) may form a layer, where the layers layer horizontally between the injection well (1506) and (1508). Such a layer is termed herein a “horizontal layer.” When an injection well and a production well are horizontal or lateral wells (not shown), each WAG zone may form a layer, where the layers layer vertically between the injection well and production well. Such a layer is termed herein “vertical layer.”

In an exemplary configuration WAG zones (1550, 1548, 1542, 1538, 1540, 1544) include aqueous zone (1550), first carbon dioxide zone (1548), accumulation zone (1542), second carbon dioxide zone (1538), miscible zone (1540), and oil zone (1544). In accordance with one or more embodiments, aqueous zone (1550) contains a treatment fluid phase that includes the treatment fluid. The treatment fluid phase may be a drive phase that is initiated by a WAG cycle.

In accordance with one or more embodiments, the first carbon dioxide zone (1548) contains a carbon dioxide phase that includes carbon dioxide. The carbon dioxide in zone (1548) may be supercritical at reservoir conditions. Accumulation zone (1542) contains an accumulation phase. The accumulation phase may include fluid with accumulated nanoparticles and trapped carbon dioxide caused by a prior WAG cycle. Second carbon dioxide zone (1538) contains a carbon dioxide phase that includes carbon dioxide. The carbon dioxide in zone (1538) may be supercritical at reservoir conditions.

In accordance with one or more embodiments, the miscible zone (1540) contains a miscible phase. The miscible phase may be caused by interaction of carbon dioxide and previous contact of residual oil with the treatment fluid. The previous contact may specifically be previous contact of residual oil with the nanoparticles in the treatment fluid. Oil zone (1544) contains a mobilized oil phase. The mobilized oil phase may be caused by the WAG cycle. Zone (1544) may also contain portions of the treatment fluid and the carbon dioxide as byproducts. The production system (1516) is configured to extract hydrocarbons from zone (1544).

It will be understood that throughout the present disclosure, where carbon dioxide is described, the carbon dioxide may be mixed with other miscible gases, such as, but not limited to, propane and butane, without departing from the scope of the disclosure.

Referring still to FIG. 15, the production well (1508) may use a production pump (1552) to extract the mobilized hydrocarbons, byproduct treatment fluid, and byproduct carbon dioxide. The production well (1508) may lead to a separation unit (1561) through a production process line (1560) containing a produced amount of hydrocarbons and byproduct amounts of the treatment fluid and carbon dioxide.

The separation unit (1561) may be in fluid communication with one or more hydrocarbon storage tanks (1558), one or more treatment fluid storage tanks (1532), and/or one or more carbon dioxide storage tanks (1520). Alternately or in combination, the treatment fluid and carbon dioxide or portions thereof may be routed directly back to the injection well (1506) as a recycle stream without any storage unit. A combination of direct recycle streams and storage tanks may be used for the treatment fluid and/or the carbon dioxide. There may be a carbon dioxide process line (1528) containing a carbon dioxide pump (1529) from the carbon dioxide storage tank (1520) to the injection well (1506). There may be a treatment fluid process line (1534) containing a treatment fluid pump from the treatment fluid storage tank (1532) to the injection well (1506).

A tanker (1524) may be used to provide carbon dioxide or the treatment fluid to the injection well (1506). Fluid pumping system (1512) may be equipped to allow the tanker (1524) to tie into either the carbon dioxide process line (1528) or the treatment fluid process line (1534). In accordance with one or more embodiments, the system contains a computer (1602) system, described below, that communicates between the production well and the injection well to optimize the extraction of hydrocarbons and maximize a storage amount of carbon dioxide under the earth's surface (1522). In further embodiments, the computer (1602) system, or one similar to it, may have been used to model the calibrated discrete fracture network outlined in S100.

FIG. 16 shows a computer (1602) system in accordance with one or more embodiments. Specifically, FIG. 16 shows a block diagram of a computer (1602) system used to provide computational functionalities associated with described algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure, according to an implementation. The illustrated computer (1602) is intended to encompass any computing device such as a server, desktop computer, laptop/notebook computer, wireless data port, smart phone, personal data assistant (PDA), tablet computing device, one or more processors within these devices, or any other suitable processing device, including both physical or virtual instances (or both) of the computing device.

Additionally, the computer (1602) may include a computer that includes an input device, such as a keypad, keyboard, touch screen, or other device that can accept user information, and an output device that conveys information associated with the operation of the computer (1602), including digital data, visual, or audio information (or a combination of information), or a GUI.

The computer (1602) can serve in a role as a client, network component, a server, a database or other persistency, or any other component (or a combination of roles) of a computer system for performing the subject matter described in the instant disclosure. The illustrated computer (1602) is communicably coupled with a network (1630). In some implementations, one or more components of the computer (1602) may be configured to operate within environments, including cloud-computing-based, local, global, or other environment (or a combination of environments).

At a high level, the computer (1602) is an electronic computing device operable to receive, transmit, process, store, or manage data and information associated with the described subject matter. According to some implementations, the computer (1602) may also include or be communicably coupled with an application server, e-mail server, web server, caching server, streaming data server, business intelligence (BI) server, or other server (or a combination of servers).

The computer (1602) can receive requests over network (1630) from a client application (for example, executing on another computer (1602)) and responding to the received requests by processing the said requests in an appropriate software application. In addition, requests may also be sent to the computer (1602) from internal users (for example, from a command console or by other appropriate access method), external or third-parties, other automated applications, as well as any other appropriate entities, individuals, systems, or computers.

Each of the components of the computer (1602) can communicate using a system bus (1603). In some implementations, any or all of the components of the computer (1602), both hardware or software (or a combination of hardware and software), may interface with each other or the interface (1604) (or a combination of both) over the system bus (1603) using an application programming interface (API) (1612) or a service layer (1613) (or a combination of the API (1612) and service layer (1613). The API (1612) may include specifications for routines, data structures, and object classes. The API (1612) may be either computer-language independent or dependent and refer to a complete interface, a single function, or even a set of APIs. The service layer (1613) provides software services to the computer (1602) or other components (whether or not illustrated) that are communicably coupled to the computer (1602).

The functionality of the computer (1602) may be accessible for all service consumers using this service layer. Software services, such as those provided by the service layer (1613), provide reusable, defined business functionalities through a defined interface. For example, the interface may be software written in JAVA, C++, or other suitable language providing data in extensible markup language (XML) format or other suitable format. While illustrated as an integrated component of the computer (1602), alternative implementations may illustrate the API (1612) or the service layer (1613) as stand-alone components in relation to other components of the computer (1602) or other components (whether or not illustrated) that are communicably coupled to the computer (1602). Moreover, any or all parts of the API (1612) or the service layer (1613) may be implemented as child or sub-modules of another software module, enterprise application, or hardware module without departing from the scope of this disclosure.

The computer (1602) includes an interface (1604). Although illustrated as a single interface (1604) in FIG. 16, two or more interfaces (1604) may be used according to particular needs, desires, or particular implementations of the computer (1602). The interface (1604) is used by the computer (1602) for communicating with other systems in a distributed environment that are connected to the network (1630). Generally, the interface (1604) includes logic encoded in software or hardware (or a combination of software and hardware) and operable to communicate with the network (1630). More specifically, the interface (1604) may include software supporting one or more communication protocols associated with communications such that the network (1630) or interface's hardware is operable to communicate physical signals within and outside of the illustrated computer (1602).

The computer (1602) includes at least one computer processor (1605). Although illustrated as a single computer processor (1605) in FIG. 16, two or more processors may be used according to particular needs, desires, or particular implementations of the computer (1602). Generally, the computer processor (1605) executes instructions and manipulates data to perform the operations of the computer (1602) and any algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure.

The computer (1602) also includes a non-transitory computer (1602) readable medium, or a memory (1606), that holds data for the computer (1602) or other components (or a combination of both) that can be connected to the network (1630). For example, memory (1606) can be a database storing data consistent with this disclosure. Although illustrated as a single memory (1606) in FIG. 16, two or more memories may be used according to particular needs, desires, or particular implementations of the computer (1602) and the described functionality. While memory (1606) is illustrated as an integral component of the computer (1602), in alternative implementations, memory (1606) can be external to the computer (1602).

The application (1607) is an algorithmic software engine providing functionality according to particular needs, desires, or particular implementations of the computer (1602), particularly with respect to functionality described in this disclosure. For example, application (1607) can serve as one or more components, modules, applications, etc. Further, although illustrated as a single application (1607), the application (1607) may be implemented as multiple applications (1607) on the computer (1602). In addition, although illustrated as integral to the computer (1602), in alternative implementations, the application (1607) can be external to the computer (1602).

There may be any number of computers (1602) associated with, or external to, a computer system containing computer (1602), each computer (1602) communicating over network (1630). Further, the term “client,” “user,” and other appropriate terminology may be used interchangeably as appropriate without departing from the scope of this disclosure. Moreover, this disclosure contemplates that many users may use one computer (1602), or that one user may use multiple computers (1602).

Although only a few example embodiments have been described in detail above, those skilled in the art will readily appreciate that many modifications are possible in the example embodiments without materially departing from this invention. Accordingly, all such modifications are intended to be included within the scope of this disclosure as defined in the following claims.

Claims

1. A method comprising: modeling a reservoir having natural fractures using a lab scale set of models and a field scale set of models, wherein modeling the reservoir using the lab scale set of models comprises: scanning a sample of the reservoir into the lab scale set of models to create modeled fractures that represent the natural fractures;estimating hydraulic properties of the modeled fractures using analytical models or artificial intelligence based models;estimating multi-phase dynamic properties of the modeled fractures based on geometric properties of the modeled fractures and the estimated hydraulic properties; anddetermining characteristics of a flow regime of a fluid flowing through the modeled fractures using the estimated hydraulic properties and the estimated multi-phase dynamic properties;wherein modeling the reservoir using the field scale set of models comprises: modeling a discrete fracture network of the reservoir using the modeled fractures, the estimated hydraulic properties, the estimated multi-phase dynamic properties, and the characteristics of the flow regime using artificial intelligence and stochastic methods to capture an anisotropy and heterogeneity of each individual modeled fracture;upscaling the discrete fracture network by determining a shape factor using a physics-informed neural network; andcalibrating the discrete fracture network using machine learning algorithms to create a calibrated discrete fracture network; anddesigning and performing an enhanced oil recovery operation on the reservoir using the calibrated discrete fracture network.

2. The method of claim 1, wherein designing and performing the enhanced oil recovery operation on the reservoir using the calibrated discrete fracture network further comprises pumping a fluid into the reservoir using an injection well.

3. The method of claim 2, wherein designing and performing the enhanced oil recovery operation on the reservoir using the calibrated discrete fracture network further comprises flowing the fluid through the natural fractures of the reservoir towards a production well.

4. The method of claim 3, wherein a location of the injection well is based on the calibrated discrete fracture network.

5. The method of claim 1, wherein modeling the reservoir using the lab scale set of models further comprises calibrating the modeled fractures using a reservoir model using single-phase and multi-phase near-wellbore radial and cartesian models.

6. The method of claim 1, wherein modeling the reservoir using the field scale set of models further comprises quantifying and ranking static and dynamic properties of the discrete fracture network.

7. The method of claim 6, wherein quantifying and ranking the static and dynamic properties of the discrete fracture network further comprises determining a geometric connectivity and the heterogeneity of the modeled fractures based on a static Lorenz coefficient.

8. The method of claim 1, wherein modeling the reservoir using the field scale set of models further comprises modeling localized flow mechanisms in an interface between the modeled fractures and a matrix of the reservoir using numerical methods.

9. The method of claim 1, wherein estimating the estimated hydraulic properties of the modeled fractures using the analytical models or the artificial intelligence based models further comprises estimating a stress-dependent permeability under effective normal stress and coupled flow-normal shear conditions.

10. The method of claim 1, wherein estimating the estimated multi-phase dynamic properties of the modeled fractures based on the geometric properties of the modeled fractures and the estimated hydraulic properties further comprises using a Leverett J-function corresponding to different rock-fracture types.

11. A system comprising: an enhanced oil recovery system having an injection well and a production well drilled into a reservoir having natural fractures;at least one sample of the reservoir; anda computer system configured to model the reservoir using a lab scale set of models and a field scale set of models, wherein the computer system is configured to model the reservoir using the lab scale set of models by: scanning a sample of the reservoir into the lab scale set of models to create modeled fractures that represent the natural fractures;estimating hydraulic properties of the modeled fractures using analytical models or artificial intelligence based models;estimating multi-phase dynamic properties of the modeled fractures based on geometric properties of the modeled fractures and the estimated hydraulic properties; anddetermining characteristics of a flow regime of a fluid flowing through the modeled fractures using the estimated hydraulic properties and the estimated multi-phase dynamic properties;wherein the computer system is configured to model the reservoir using the field scale set of models by: modeling a discrete fracture network of the reservoir using the modeled fractures, the estimated hydraulic properties, the estimated multi-phase dynamic properties, and the characteristics of the flow regime using artificial intelligence and stochastic methods to capture an anisotropy and heterogeneity of each individual modeled fracture;upscaling the discrete fracture network by determining a shape factor using a physics-informed neural network; andcalibrating the discrete fracture network using machine learning algorithms to create a calibrated discrete fracture network.

12. The system of claim 11, wherein the injection well is configured to pump a fluid into the reservoir.

13. The system of claim 12, wherein the fluid is configured to flow through the natural fractures of the reservoir towards the production well.

14. The system of claim 13, wherein a location of the injection well is based on the calibrated discrete fracture network.

15. The system of claim 11, wherein modeling the reservoir using the lab scale set of models further comprises calibrating the modeled fractures using a reservoir model using single-phase and multi-phase near-wellbore radial and cartesian models.

16. The system of claim 11, wherein modeling the reservoir using the field scale set of models further comprises quantifying and ranking static and dynamic properties of the discrete fracture network.

17. The system of claim 16, wherein quantifying and ranking the static and dynamic properties of the discrete fracture network further comprises determining a geometric connectivity and the heterogeneity of the modeled fractures based on a static Lorenz coefficient.

18. The system of claim 11, wherein modeling the reservoir using the field scale set of models further comprises modeling localized flow mechanisms in an interface between the modeled fractures and a matrix of the reservoir using numerical methods.

19. The system of claim 11, wherein estimating the estimated hydraulic properties of the modeled fractures using the analytical models or the artificial intelligence based models further comprises estimating a stress-dependent permeability under effective normal stress and coupled flow-normal shear conditions.

20. The system of claim 11, wherein estimating the estimated multi-phase dynamic properties of the modeled fractures based on the geometric properties of the modeled fractures and the estimated hydraulic properties further comprises using a Leverett J-function corresponding to different rock-fracture types.