METHOD AND SYSTEM TO ACCELERATE THE NONLINEAR SOLUTION IN NUMERICAL RESERVOIR SIMULATION

Description

BACKGROUND

Numerical reservoir simulation is used in production of hydrocarbons to predict remaining quantities of oil and gas, pressures, and flow-rates in an oil reservoir given certain production scenarios. Such simulations can predict locations and quantities of hydrocarbons in reservoirs over months, years, or even decades, given initial conditions according to chosen production rates and locations.

Numerical reservoir simulation may be performed using a grid of cells that represent the reservoir properties, including rock and fluid properties, over three-dimensional space within the reservoir. To obtain accurate results, fine grid cells may be required, particularly when the reservoir properties vary rapidly in space. Using a large number of fine grid cells may be very computationally expensive when covering any production area of realistic spatial extent. Mass balance equations for each field variable describing the flux of mass in or out of the grid cell and, where applicable a well term, describing injection or extraction, may be solved for each grid cell. Examples of field variables maintained in each grid cell include water, gas, and oil saturations, capillary pressures, relative permeabilities for water, gas, oil/water, oil/gas, and porosities.

The need to speed up reservoir simulations in production scenarios is limitless and never-ending. Thus, there is a clear and pressing need for methods and processes to speed up calculations while maintaining necessary accuracy in calculations of the physical variables.

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.

In general, in one aspect, embodiments disclosed herein relate to a method. The method includes obtaining a reservoir model of a subterranean region of interest, where the reservoir model comprises a plurality of computational cells, and obtaining a production scenario. The method further includes using a reservoir simulator, determining a first set of parameters for a first time from the reservoir model and the production scenario and predicting, using the first set of parameters, a second set of parameters at a second time, wherein the first and second sets of parameters comprise primary variables, wherein the primary variables are used to determine values for a mass, a flux and a well term, wherein the second time is after the first time. The predicting includes, iteratively, or recursively, until a stopping criterion is satisfied, for each of the plurality of computational cells, using the first set of parameters to determine a primary variable update and a residual, determining an updated flux and an updated well term using, at least in part, the primary variable update, predicting the second set of parameters based, at least in part, on the updated flux and the updated well term, and assigning the second set of parameters to the first set of parameters. The method still further includes determining, using an interpretation workstation, a predicted location of unproduced hydrocarbons based, at least in part, on the second set of parameters.

In general, in one aspect, embodiments disclosed herein relate to a system. The system includes a reservoir simulator and an interpretation workstation. The reservoir simulator is configured to receive a reservoir model of a subterranean region of interest, wherein the reservoir model comprises a plurality of computational cells, receive a production scenario, determine a first set of parameters for a first time from the reservoir model and the production scenario, and predict, using the first set of parameters, a second set of parameters at a second time, wherein the first and second sets of parameters comprise primary variables, where the primary variables are used to determine values for a mass, a flux and a well term, and where the second time is after the first time. Predicting includes, iteratively, or recursively, until a stopping criterion is satisfied, for each of the plurality of computational cells using the first set of parameters to determine a primary variable update and a residual, determining an updated flux and an updated well term using, at least in part, the primary variable update, and predicting the second set of parameters based, at least in part, on the updated flux and the updated well term, and assigning the second set of parameters to the first set of parameters. The interpretation workstation is configured to determine a predicted location of unproduced hydrocarbons based, at least in part, on the second set of parameters.

It is intended that the subject matter of any of the embodiments described herein may be combined with other embodiments described separately, except where otherwise contradictory.

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

BRIEF DESCRIPTION OF DRAWINGS

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

FIG. 2A shows a reservoir region in accordance with one or more embodiments.

FIG. 2B shows a discretized reservoir grid model in accordance with one or more embodiments.

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

FIG. 4A shows reservoir parameter curves in accordance with one or more embodiments.

FIG. 4B shows reservoir parameter curves in accordance with one or more embodiments.

FIG. 4C shows reservoir parameter curves in accordance with one or more embodiments.

FIG. 5 shows a reservoir model in accordance with one or more embodiment.

FIG. 6 shows graphs of simulation results in accordance with described simulations.

FIG. 7 shows performance parameter curves in accordance with described simulations.

FIG. 8 is a table of an overview of performance data comparisons for described simulations.

FIG. 9 shows a reservoir model in accordance with one or more embodiments.

FIG. 10 shows graphs of simulation results in accordance with described simulations.

FIG. 11 shows performance parameter curves in accordance with described simulations.

FIG. 12 is a table of an overview of performance data comparisons for described simulations.

FIG. 13 shows a reservoir model in accordance with one or more embodiments.

FIG. 14 shows graphs of simulation results in accordance with described simulations.

FIG. 15 shows performance parameter curves in accordance with described simulations.

FIG. 16 is a table of an overview of performance data comparisons for described simulations.

FIG. 17 displays a computer system in accordance with one or more embodiments.

FIG. 18 shows a flowchart of steps 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.

It is to be understood that the singular forms “a,” “an,” and “the” include plural referents unless the context clearly dictates otherwise. Thus, for example, reference to a “parameter curve” includes reference to one or more of such parameter curves.

Terms such as “approximately,” “substantially,” etc., mean that the recited characteristic, parameter, or value need not be achieved exactly, but that deviations or variations, including for example, tolerances, measurement error, measurement accuracy limitations and other factors known to those of skill in the art, may occur in amounts that do not preclude the effect the characteristic was intended to provide.

It is to be understood that one or more of the steps shown in the method may be omitted, repeated, and/or performed in a different order than the order shown. Accordingly, the scope disclosed herein should not be considered limited to the specific arrangement of steps shown in the method.

Although multiple dependent claims are not introduced, it would be apparent to one of ordinary skill that the subject matter of the dependent claims of one or more embodiments may be combined with other dependent claims.

In the following description of FIGS. 1-18, any component described with regard to a figure, in various embodiments disclosed herein, may be equivalent to one or more like-named components described with regard to any other figure. For brevity, descriptions of these components will not be repeated with regard to each figure. Thus, each and every embodiment of the components of each figure is incorporated by reference and assumed to be optionally present within every other figure having one or more like-named components. Additionally, in accordance with various embodiments disclosed herein, any description of the components of a figure is to be interpreted as an optional embodiment which may be implemented in addition to, in conjunction with, or in place of the embodiments described with regard to a corresponding like-named component in any other figure.

In characterizing subsurface reservoirs, particularly hydrocarbon reservoirs, it is frequently helpful to be able to predict their behavior, such as production rates, in response to operations steps, such as drilling new wells or injecting water around the periphery of the reservoir. These predictions are conventionally made with a reservoir simulator that use numerical descriptions of the subsurface reservoirs (reservoir “models”) and physics-based processes and methods to form the predictions. These reservoir models contain parameters describing the characteristics, particularly the fluid flow characteristics, such as porosity and permeability, of the subsurface reservoir. The accuracy of the predictions, that is the degree to which the prediction matched the actual future behavior of the reservoir, is strongly dependent on the fidelity with which the values of the reservoir model parameters represent the actual reservoir. While conventional methods exist, including well logging and simulating the historical behavior of the reservoir (“history matching”), for determining the values of the reservoir model parameters, the embodiments disclosed herein represent a significant improvement over conventional methods for such a determination.

FIG. 1 shows an example of a typical 3-dimensional reservoir model (100) used to simulate the operation of an oil and/or gas reservoir. FIG. 1 includes a plurality of wells, such as wells (102) and (103). Although the wells shown in FIG. 1 are straight wells or well with minimum curvature and orientated in a vertical or near vertical orientation, typical reservoir models, such as reservoir model (100), may include curved wells, such as deviated, highly deviated, and horizontal wells. Similarly, while FIG. 1 depicts only a few wells, a typical reservoir model may frequently include tens or hundreds of wells. The wells present in a reservoir model may include exploration wells (drilled while searching for the oil and/or gas reservoir depicted in the model), appraisal wells (drilled while determining the extent and quality of the reservoir model) production wells (drilled to extract oil and gas from the reservoir) and injection wells (drilled to inject water or gas, such as nitrogen or air, into the reservoir to maintain reservoir pressure and sweep hydrocarbons towards the production wells). In FIG. 1, the filled cylinders, such as well (102) may depict an injection well, while the open cylinders, such as well (103), may depict a production well. FIG. 1 depicts 2 other injection wells, in addition to well (102), and 5 other production wells, in addition to well (103) . . . . In FIG. 1, the reservoir model (100) is discretized into a plurality of nodes (104) forming a grid mesh (106). In some cases, a grid mesh may be a regular Cartesian grid, in other cases the grid mesh may be a regular hexagonal mesh. In still other cases the grid may be of irregular or non-uniform shape and grid mesh cell size. For example, the grid mesh may be irregular in order to conform to geological features within the reservoir, such as faults, or layer boundaries, or to provide greater resolution and accuracy in the vicinity of wells. In reservoir model (100), the grid cells are of non-uniform size and shape. Thus, although reservoir model (100) depicts a majority of hexagonal grid cells discretizing the model (100), the grid also contains Cartesian cells (cuboids), triangular cells and irregular shapes.

In some embodiments, a reservoir simulator includes functionality for simulating the flow of fluids, including hydrocarbon fluids such as oil and gas, through a hydrocarbon reservoir composed of porous, permeable reservoir rocks in response to natural and anthropogenic pressure gradients. The reservoir simulator may be used to predict changes in fluid flow, including fluid flow into wells penetrating the reservoir, as a result of planned well drilling, and fluid injection and/or extraction. For example, the reservoir simulator may be used to predict changes in hydrocarbon production rate that would result from the injection of water into the reservoir from wells around the reservoir periphery.

The reservoir simulator may use a reservoir model that contains a digital description of the physical properties of the rocks as a function of position within the reservoir and the fluids within the pores of the porous, permeable reservoir rocks at a given time. In some embodiments, the digital description may be in the form of a dense three-dimensional (“3D”) grid with the physical properties of the rocks and fluids defined at each node. In some embodiments, the 3D grid may be a cartesian grid, while in other embodiments the grid may use an alternative coordinate system, such as a hexagonal coordinate system, or may be an irregular grid.

The physical properties of the rocks and fluids within the reservoir may be obtained from a variety of geological and geophysical sources. For example, remote sensing geophysical surveys, such as seismic surveys, gravity surveys, and active and passive source resistivity surveys, may be employed. In addition, data collected from well logs acquired in well penetrating the reservoir may be used to determine physical and petrophysical properties along the segment of the well trajectory traversing the reservoir. For example, porosity, permeability, density, seismic velocity, and resistivity may be measured along these segments of wellbore. In accordance with some embodiments, remote sensing geophysical surveys and physical and petrophysical properties determined from well logs may be combined to estimate physical and petrophysical properties for the entire reservoir simulation model grid.

Reservoir simulators solve sets of mathematical governing equations that represent the physical laws that govern fluid flow in porous, permeable media. For example, the flow of a single-phase slightly compressible oil with a constant viscosity and compressibility obeys the equations captured by Darcy's law, the continuity condition and the equation of state, and may be written as:

∇ 2 p ⁡ ( x , t ) = φ ⁢ μ ⁢ c t k ⁢ ∂ p ⁡ ( x , t ) ∂ t ,

where p represents fluid pressure in the reservoir, x is a vector representing spatial position and t represents time. φ, μ, c_t, and k represent the physical and petrophysical properties of porosity, fluid viscosity, total combined rock and fluid compressibility, and permeability, respectively. ∇² represents the spatial Laplacian operator.

Additional, and more complicated equations are required when more than one fluid, or more than one phase, e.g., liquid and gas, are present in the reservoir. Further, when the physical and petrophysical properties of the rocks and fluids vary as a function of position the governing equations may not be solved analytically and must instead be discretized into a grid of cells or blocks. The governing equations must then be solved by one of a variety of numerical methods, such as, without limitation, explicit or implicit finite-difference methods, explicit or implicit finite element methods, or discrete Galerkin methods.

In some embodiments, reservoir simulation may be performed using the estimated reservoir properties for the unconventional hydrocarbon reservoir. For example, the reservoir simulator may include hardware and/or software with functionality for generating one or more reservoir models regarding the hydrocarbon-bearing formation and/or performing one or more reservoir simulations. For example, the reservoir simulator may store well logs and data regarding core samples for performing simulations. A reservoir simulator may further analyze the well log data, the core sample data, seismic data, and/or other types of data to generate and/or update the one or more reservoir models. In some embodiments, the reservoir simulator may include a computer that is similar to the computer (1700) described below with regard to FIG. 17 and the accompanying description.

Turning back to FIG. 2A, FIG. 2A shows a schematic diagram in accordance with one or more embodiments. As illustrated in FIG. 2A, FIG. 2A shows a geological region (202) that may include one or more reservoir regions (e.g., reservoir region (204)) with various production wells (e.g., production well A (208), production well B (210)). Likewise, a reservoir region may also include one or more injection wells (e.g., injection well C (206)) that include functionality for enhancing production by one or more neighboring production wells. As shown in FIG. 2A, wells may be disposed in the reservoir region (204) above various subsurface layers (e.g., subsurface layer A (212), subsurface layer B (214)), which may include hydrocarbon deposits. Production data and/or injection data may exist for a particular well, where production data may include data that describes production or production operations at a well, such as flow-rates, fluid properties and compositions of samples taken at the wellhead.

In conventional methods to reservoir simulation, the grid of a reservoir model may be refined, to provide enhanced resolution at predetermined spatial locations, such as near wellbores. By way of illustration, FIG. 2B shows a reservoir grid model (216) that corresponds to the geological region (202) from FIG. 2A. More specifically, the reservoir grid model (216) includes grid cells (218) that may refer to an original cell (218) of a reservoir grid model as well as coarse grid blocks (220) that may refer to an amalgamation of original grid cells (218) of the reservoir grid model. For example, a grid cell may be the case of a 1×1 block, where coarse grid blocks may be of sizes 2×2, 4×4, 8×8, etc. Both the grid cells (218) and the coarse grid blocks (220) may correspond to columns for multiple model layers (226) within the reservoir grid model (216).

Prior to performing a reservoir simulation, local grid refinement and coarsening may be used to increase or decrease grid resolution in a certain area of reservoir grid model. For example, various reservoir properties, e.g., permeability, porosity or saturations, may correspond to a discrete value that is associated with a particular grid cell or coarse grid block. However, by using discrete values to represent a portion of a geological region, a discretization error may occur in a reservoir simulation. Thus, finer grids may reduce discretization errors as the numerical approximation of a finer grid is closer to the exact solution, however through a higher computational cost. As shown in FIG. 2B, for example, the reservoir grid model (216) may include various fine-grid models (i.e., fine-grid model A (222), fine-grid model B (224)), that are surrounded by coarse block regions. Likewise, the original reservoir grid model without any coarsening may also be a fine-grid model.

In some embodiments, proxy models or reduced-order models may be generated for performing a reservoir simulation. For example, one way to reduce model dimensionality is to reduce the number of grid blocks and/or grid cells. By averaging reservoir properties into larger blocks while preserving the flow properties of a reservoir model, computational time of a reservoir simulation may be reduced. In general, coarsening may be applied to cells that do not contribute to a total flow within a reservoir region because a slight change on such reservoir properties may not affect the output of a simulation. Accordingly, different levels of coarsening may be used on different regions of the same reservoir model. As such, a coarsening ratio may correspond to a measure of coarsening efficiency, which may be defined as a total number of cells in a coarse reservoir model divided by the original number of cells in the original reservoir model.

Flow properties, such as flux, may be defined as the rate at which a reservoir fluid (e.g., oil or natural-gas) flows between any two grid blocks. Likewise, grid cells or blocks may be upscaled in a method that reduces the computational demand on running simulations using fewer grid cells. However, a grid model may lose accuracy in a reservoir simulation if the underlying properties differ too much from the original fine-grid model.

An example of steps involved in performing one time-step of a reservoir simulation is shown in FIG. 3, in accordance with one or more embodiments. A single time-step uses values of parameters describing each of the grid cells representing the reservoir at a first time to predict the values of parameters at a second, later, time. These steps may include an iterative process in which a conventional Jacobian may be calculated, in step (300), followed, in step (310), by a solution to a linear system of equations. In step (320) the inter-cell fluxes and well terms, may be calculated from the previous values for mass, the updated fluxes, and the updated well flows. In step (330) a new value for mass may be determined from the inter-cell fluxes and well terms. In step (340) the mole fraction and density of each phase are updated using the updated pressures and masses. This is followed by updating the phase state and each phase saturation in step (350). Saturation constraints may be enforced in step (360). In step (370) saturation dependent properties may be determined from the constrained saturations. Candidate nonlinear residuals for the present step of the Newtonian iteration are finally determined in step (380). These candidate non-linear residuals are compared with a convergence criterion, or criteria, the end of the loop over Newtonian iteration in step (390). In some embodiments, the linear solve in step (310) and/or the property calculations in steps (340) and/or (350) may be performed over large arrays of data that represent properties such as, for example, pressure and composition at nodes in the grid mesh (106).

If the convergence criterion or criteria are not satisfied, in step (390) the iterative method may return to step (300) to use the update parameters within the system equations. However, if the convergence criterion or criteria are satisfied, in step (390), then the iterative method may terminate for the current time, in method step (392), and a new time step of the reservoir simulation may be commenced in step (394).

The key feature of the embodiment lies in the approach to updating the primary variables. The embodiment of the inventive method introduces a sequence of steps to accomplish this procedure steps (320) through (370).

In some embodiments, the reservoir simulation steps shown in FIG. 3 is implemented in the following way. A system is solved in step (300) using a fully-implicit method in time and a finite-volume method in space, where each grid cell is a finite-volume. The discrete mass balance equation for fluid component i at each grid cell can be expressed by

V ⁢ ϕ Δ ⁢ t ⁢ Δ [ ρ 0 ⁢ x i ⁢ S o + ρ g ⁢ y i ⁢ S g ] = ∑ j ⁢ T ⁢ ρ o ⁢ x i ⁢ k r ⁢ o μ 0 ⁢ ( Δ ⁢ Φ o ) + ∑ j ⁢ T ⁢ ρ g ⁢ y i ⁢ k r ⁢ g μ g ⁢ ( ΔΦ g ) - q i , ( 1 ⁢ a ) V ⁢ ϕ Δ ⁢ t ⁢ Δ [ ρ w ⁢ S w ] = ∑ j ⁢ T ⁢ ρ w ⁢ k r ⁢ w μ w ⁢ ( Δ ⁢ Φ w ) - q w , ( 1 ⁢ b )

where:

• • k_r=relative permeability
  • S=saturation
  • T=transmissibility
  • V=cell volume
  • Δt=time step size
  • ρ=molar density
  • x=oil mole fraction
  • y=gas mole fraction
  • ϕ=porosity
  • Φ=potential
  • μ=viscosity
  • q_i=well term for component i

Additionally, the saturation constraint, the phase mole fraction constraints, and the capillary pressure relationship are given as follows in equations (2)-(4). The saturation constraint is:

S o + S g + S w = 1 ( 2 )

where S_o, S_g and S_w are the saturations for oil, gas, and water, respectively. The phase mole fraction constraints are the following,

∑ i = 1 n c ⁢ x i = 1 ( 3 ⁢ a ) ∑ i = 1 n c ⁢ y i = 1 ( 3 ⁢ b )

where x_iand y_i are the phase mole fraction for grid cell i of oil and gas, respectively. Capillary pressures are restricted by the following constraints,

P cow = p o - p w = f ⁢ ( S w ) ( 4 ⁢ a ) P cog = p g - p o = f ⁢ ( S g ) ( 4 ⁢ b )

where P_cow and P_cog are the capillary pressures for the oil-water and oil-gas boundaries, respectively, and p_o, p_w and p_g are oil, water and gas pressures in each cell. Each of the capillary pressures of oil-water and oil-gas will decide the water and gas saturations respectively through a function f, that will be empirically determined.

In equations (1a) and (1b), relative permeabilities and capillary pressures are typical input data in the form of two-phase rock tables. A typical set of measured rock curves are shown in FIG. 4A, FIG. 4B, and FIG. 4C, represented by the following parameters:

• • Oil-Water Table: S_w, k_row, k_rw, P_cow
  • Gas-Liquid Table: S_g, k_rog, k_rg, P_cog
    where k_row, k_rog, k_rw and k_rg represent the relative permeability of the rock at the oil-water boundary, oil-gas boundary, single-phase water and single-phase gas, respectively, and:
  • Kr: relative permeability
  • Krow: oil phase relatively permeability in oil-water system
  • Krw: water phase relative permeability in oil-water system
  • Swc: connate water saturation in oil-water system
  • Sgc: connate gas saturation in gas-oil system
  • Sw: water saturation
  • Sg: gas saturation
  • Pc: capillary pressure between oil and water in oil-water system.

Connate water saturation (Swc), also known as irreducible water or bound water, is the water that remains trapped in the pore spaces of a reservoir rock after the initial oil migration and accumulation. Connate water is held in place by capillary forces and cannot be displaced by the flow of oil or gas. The connate water saturation (Swc) represents the fraction of the pore volume occupied by this trapped water.

Irreducible oil saturation (Sorw), also known as residual oil saturation, is the minimum oil saturation that remains trapped in the pore spaces of a reservoir rock after extensive water or gas flooding. Irreducible oil is held in place by capillary forces and cannot be displaced by the flow of water or gas under normal reservoir conditions. The irreducible oil saturation (Sor) represents the fraction of the pore volume occupied by this trapped oil.

Connate gas saturation (Sgc) and Irreducible oil saturation (Sorg) have the same physical meaning as Swc and Sorg, but are applicable to a gas-oil system rather than a oil-water system.

The relative permeability curves, shown in FIGS. 4A and 4B characterizes the ability of a porous medium to conduct one fluid phase in the presence of another phase, when the simultaneous flows of multiple phases pass through the reservoir rock. The phase may be oil phase, water phase or gas phase. The relative permeability is a relative value with reference to the permeability of rock measured for single-phase pure water flow. On microscale, the physical processes govern these curves are based on the interactions between fluids, the pore structure, and the wettability of the rock surface.

A capillary pressure curve, such as shown in FIG. 4C describes the relationship between the capillary pressure and the saturation of the wetting phase in a porous medium. The capillary pressure curve is governed by the physical processes such as interactions between fluids, the pore structure, and the wettability of the rock surface.

In FIG. 4C a negative capillary pressure on the imbibition curve may indicate a state where the wetting phase (usually water) is being spontaneously imbibed or drawn into the porous medium. This phenomenon is driven by the capillary forces that arise from the interaction between the wetting phase, the non-wetting phase (usually oil or gas), and the solid surface of the porous medium. When the capillary pressure is negative, the wetting phase is at a lower pressure than the non-wetting phase. This pressure difference causes the wetting phase to be pulled into the smaller pores of the medium, displacing the non-wetting phase. The more negative the capillary pressure, the stronger the imbibition effect.

Drainage involves the displacement of the wetting phase by the non-wetting phase, while imbibition involves the displacement of the non-wetting phase by the wetting phase. In drainage, the non-wetting phase enters the pore spaces while in imbibition, the wetting phase enters the pore spaces. In a typical capillary pressure curve, the imbibition curve is usually located below the drainage curve.

The permeabilities, saturations and capillary pressures may be experimentally or empirically determined. These functions have discontinuous derivatives at their connate and irreducible saturations Swc, Sorw, Sorg, and Sgc. In simulation practice, each grid cell will have unique connate and irreducible saturation values, which will lead to a different set of rock-fluid tables for each grid cell. The initialization capillary pressure is drainage, which must be consistent with measured water saturation and reservoir original oil in place (“OOIP”) values. As water saturation, Sw, increases during simulation, the capillary pressure, Pcow, will decrease in a water displacing oil scenario, which will be consistent with an imbibition capillarity.

In field-scale simulation, each grid cell can be at one of four different phase states: (1) three-phase oil-gas-water; (2) two-phase oil-water; (3) two-phase gas-water; and (4) single-phase water. The choice of primary variables for the system depends on the phase states. For phase type 1, the primary variables are x₁, x₂, . . . x_nc−1, S_o, S_g and p. For phase type 2, the primary variables are x₁, x₂, . . . x_nc−1, S_o, and p. For phase type 3, the primary variables are y₁, y₂, . . . y_nc−1, S_g and p. The phase-type 4 primary variables are the same as phase-type 1. The reason for this is that single-phase water, phase type 4, becomes singular for solution of the phases, without inserting either a small trace of gas or a small trace of oil within the cell—thereby adding the corresponding mass balance equation for the relevant phase to the governing equations for the grid cell.

In step (300) of FIG. 3, we linearize the coupled nonlinear system of equations via the Newtonian iterative method. In residual form, the mass balance equation for component i for a grid cell has the form

R i n + 1 = 1 Δ ⁢ t ⁢ ( Mass i n + 1 - Mass i n ) + ∑ j ∈ n ⁢ g ⁢ b ⁢ r ⁢ s ⁢ Flux i n + 1 + Q i n + 1 ( 5 )

where n+1 refers to the new time level. For the k^th Newtonian iteration is used

R i n + 1 , k + 1 ∼ R i n + 1 , k + [ ∂ R i ∂ X j ] n + 1 , k * Δ ⁢ X i k = 0 ( 6 )

where

[ ∂ R ∂ X ] n + 1 , k

is the Jacobian matrix J.

In step (310), to find the linear update ΔX_i^k for the primary variables, X_i, for the k^th Newtonian iteration, we solve the linear system as

Δ ⁢ X i k = - J - 1 ⁢ R i n + 1 , k ( 7 )

where ΔX_i^k is the linear solution update for the primary variables. The linear solver method employed here has been disclosed in greater detail in, prior granted, U.S. Pat. Nos. 9,208,268 and 9,069,102 issued to one of the inventors of the current application.

In step (320), the inter-cell fluxes and well layer flow in equation (5) may be updated, using the following formulae,

Flux i k + 1 = Flux i k + [ ∂ Flux i k ∂ X j ] ⁢ { Δ ⁢ X j } , ( 8 ⁢ a ) Q i k + 1 = Q i k + [ ∂ Q i k ∂ X j ] ⁢ { Δ ⁢ X j } . ( 8 ⁢ b )

Using the mass balance Eqn. 5, and absorbing R_i^n+1,k+1 into the new mass updates, in step (330) we calculate the updated mass from the old mass, the updated fluxes and the updated well flows from equations (8a) and (8b) as:

Mass i k + 1 = Mass i k + Δ ⁢ t * ( ∑ Flux i k + 1 + Q i k + 1 ) ( 9 )

If hydrocarbon components do not partition into the aqueous phase, the new water saturation, S_w, can be calculated directly as

S w k + 1 = Mass w k + 1 V ⁢ ∅ ⁢ ρ w k + 1 ( 10 )

where Mass_w^k+1 and ρ_w^k+1 are the updated mass and density of water in the cell, and V and Ø are the volume and water porosity of the cell.

In step (340), for the three-phase grid cell, we use the new pressure to compute the new mole fraction and molar density for the hydrocarbon phases. In step (350), using the new mass, new pressure, and phase properties from step (340), the new phase type of the grid cell is determined from a flash calculation. The hydrocarbon liquid and gas volumes are calculated from the mass and density of each grid cell. The phase saturations are then calculated as

S p ⁢ h k + 1 = Vol p ⁢ h k + 1 V ⁢ ∅ k + 1 ( 11 )

in step (350), where Vol_ph^k+1 is the phase volume and Ø^k+1 is the updated porosity of the grid cell.

During nonlinear iterations, each grid cell can change phase types, and the derivatives of the nonlinear functions are discontinuous at the phase-type boundaries. For a three-phase cell, if pressure increases over time, the saturated oil volume increases until the gas phase disappears, then a further increase in pressure will cause the oil volume to decrease. Similar derivative discontinuities exist for the viscosity functions. It is found that when derivatives are needed across such discontinuity boundaries of physical variables, calculated values can still be used constructively in numerical calculations over a certain distance across the discontinuities. The rules that apply for this process are the subject of experimentation and constitute part of the embodiments of the disclosed invention.

In step (350), if the phase type of a grid cell changes from three-phase (type 1) to two-phase (type 2 or 3) or vice versa, primary variable changes from S_g to x₁ (cell type 1 to type 2) or from S_o to y₁ (cell type 1 to type 3), so the suitable primary variables must be set for the next (k+1) iteration. When switching from type 1 to type 2 or type 1 to type 3, a phase saturation is set to zero and the corresponding phase mole fractions are set to those calculated from step (340), which will be the new primary variables for the (k+1)^th iteration. When switching from two-phase (type 2 or type 3) to three-phase (type 1), the new phase saturations are estimated from equation (11), with the caveat that each accepted saturation change is limited to a maximum fidelity limit of 0.2 of the linear mass updates. A grid cell is considered to be an aquifer cell (type 4) when the new water saturation, S_w^k+1, is greater than 1−ϵ_S.

In step (350), if phase saturation crosses from the mobile region to the immobile region, the saturation is corrected to a suitable point within the immobile region. If it crosses from the immobile region to mobile region, we accept the change to the fidelity limit of the saturation update referred to above.

In the context of two-phase relative permeability curves, the mobile and immobile regions refer to the saturation range in which a particular fluid phase can flow (mobile) or, alternatively, remains trapped in the pore spaces (immobile) of a reservoir rock.

The mobile region represents the saturation range where a fluid phase (oil or water) can flow through the pore spaces of the reservoir rock. In this region, the relative permeability of the fluid phase is greater than zero, indicating that the phase has sufficient connectivity and mobility to flow under the prevailing pressure gradient. For the water phase, the mobile region starts at the connate water saturation (Swc) and extends up to the maximum water saturation (1.0). In this range, water can flow through the pore spaces, and its relative permeability increases as the water saturation increases. For the oil phase, the mobile region starts at the irreducible oil saturation (Sorw) and extends up to the maximum oil saturation (1−Swc). In this range, oil can flow through the pore spaces, and its relative permeability increases as the oil saturation increases.

The immobile region represents the saturation range where a fluid phase (oil or water) remains trapped in the pore spaces and cannot flow under normal reservoir conditions. In this region, the relative permeability of the fluid phase is zero, indicating that the phase lacks sufficient connectivity and mobility to flow. For the water phase, the immobile region is the saturation range below the connate water saturation (Swc). In this range, water is held in place by capillary forces and cannot be displaced by the flow of oil. For the oil phase, the immobile region is the saturation range below the irreducible oil saturation (Sorw). In this range, oil is trapped in the pore spaces by capillary forces and cannot be displaced by the flow of water.

Crossing from a mobile to an immobile region implies a temporal change, due to changes in reservoir conditions, that may happen between two subsequent time steps or iterations.

For two-phase oil-water immiscible system, the new oil saturation, S_o, can be calculated directly, and step (350) is not required:

S o k + 1 = Mass o k + 1 V ⁢ ∅ ⁢ ρ o k + 1 ( 12 )

In step (360), we compute the saturation update error as:

ε S = 1 - S o k + 1 - S w k + 1 - S g k + 1 ( 13 )

If the saturation error is smaller than the specified tolerance, ϵ_S<tol_S, the calculated saturation for the grid cell is accepted as the new saturation, in step (360). In step (370), the saturation dependent properties such as relative permeabilities and capillary pressures are updated. In step (380), the new nonlinear residuals, R_i^n+1,k+1, for each component of each cell may be computed using equation (5). In step (390), the root mean square (RMS) of the nonlinear residual vector for each fluid component and the overall RMS residual vector of all components may be calculated. Any of these RMS values may be compared against input residual preset tolerances, to determine whether acceptable convergence has been reached. The total mass balance error for each fluid component for the time step may also be computed and the error in terms of a fraction of the original mass in the simulation domain may be evaluated. If the estimated total number of timesteps times the fraction is not smaller than an acceptable global accumulated mass balance preset tolerance, convergence is considered not achieved. If convergence has not been achieved on the basis of all the above mentioned criteria, nonlinear iteration is advanced to the (k+1) th Newtonian iteration, and program control goes back to step (300). Otherwise, the current time step has converged, and program control advances to the next time step.

Three test problems are presented below to illustrate the efficiency of the embodiment of the current inventive method—each corresponding to a class of reservoir simulation models: (1) a two-phase oil-water model; (2) a three-phase oil-water-gas black-oil model: and (3) a three-phase wet gas-water-condensate extended black-oil model. In each example, a conventional method from the prior art is compared with an embodiment of the current invention. The prior-art method used is a fully-coupled fully-implicit finite-volume method, but without application of an embodiment of the current inventive method to accelerate nonlinear convergence.

The first model has 50 million cells and consists of a corner-point geometry (CPG) discretized by an irregularly structured grid with dimensions of NX*NY*NZ=214*590*396. A picture of the reservoir model is shown in FIG. 5.

The model of FIG. 5 has several rock-fluid regions as well as multiple equilibrium regions with tilted oil-water contacts. Connate and irreducible saturations are defined on a cell-by-cell basis, allowing for variations of relative permeabilities and capillary pressure. This simulation model has over 5,000 wells.

This reservoir is initialized by honoring the input initial water saturation, S_wi, which is derived from petrophysical interpretation of well logs and/or core samples. This poses difficulties to the specification of the initial vertical equilibrium conditions for the reservoir model. To honor the oil-in-place (OIP) and initial vertical equilibrium, an initial capillary pressure (Pc) value is assigned to each cell. As water saturation increases, the capillary pressure is reduced accordingly to honor imbibition capillarity.

Using the prior art method, this model runs in 11.15 hours with an average time step of 2.95 days. The final mass balance errors are for Oil=4.086E-4, and Water=1.343E-5. The time stepping statistics are as follows: Number of Time Steps=40,914; Number of Newton Iterations=82,133; Number of Solver Iterations=407,201; Number of Time Step Cuts=1,074. In contrast, with the new embodiment of the current inventive method, the model runs in 3.27 hours with an average time-step size of 12.6 days. The final mass balance errors are at Oil=3.153E-5, and Water=1.653E-5. Thus, for this model, we see a speed-up of 3.41 and improvements in mass balance precision. The new run statistics with the embodiment of the current inventive method are Number of Time Steps=9,374; Number of Newton Iterations=21,404; Number of Solver Iterations=184,728; Number of Time-Step Cuts=46. Number of Time Step Cuts means how many time steps were performed which subsequently needed to be re-simulated with a shorter value due to non-convergence of the linear Newtonian iterative steps within the nonlinear time step. Both model runs were conducted on the same HPC clusters using the same amount of hardware resources.

FIG. 6 displays a comparison of predictions of various production parameters generated by the embodiment of the current inventive method and the prior-art method. In each of chart (602), (604), (606), and (608) the dotted lines represent predictions produced by a conventional method and the solid lines represent predictions produced by an embodiment of the current invention. Specifically, chart (602) displays predicted oil production rate on the vertical axis (612), chart (604) displays predicted water rate on the vertical axis (614), chart (606) displays predicted water cut, i.e., the ratio of the water flow rate to the total liquids (oil and water) flow rate, along the vertical axis (616), and chart (608) shows the predicted average pressure of the domain of production on the vertical axis (618). Each of charts (602), (604), (606) and (608) display their predictions as a function of the production time, indicated on their horizontal axes (610).

FIG. 6 indicates that the embodiment of the current inventive method and the conventional method give equivalent predictions of various production parameters. However, FIG. 7 indicates that the numerical performance of the current inventive method is significantly better than that of the conventional method. Specifically, chart (702) shows elapsed simulation time along the vertical axis (712) as a function of production time along the horizontal axis (610). Elapsed time required by the embodiment of the invention, indicated by the solid line, is significantly shorter than that required by the conventional method, represented by the dotted line. For at least this reason, the embodiment of the invention represents a significant improvement over the conventional method.

Similarly, chart (704) displays the number of necessary nonlinear iterations indicated on the vertical axis (714), chart (706) displays the number of necessary linear iterations indicated on the vertical axis (716), and chart (708) indicates the number of time step cuts on the vertical axis (718). In each of charts (704), (706), and (708), the values for the embodiment of the invention are indicated by the solid line, and the values for the conventional method are indicated by the dotted line. In each of charts (704), (706), and (708), the production time (610) is indicated by the horizontal axis. Each of charts (702), (704), (706) and (708) clearly shows an improved performance of the embodiment of the invention over the conventional method.

Further details on the acceleration achieved by an embodiment of the inventive method are listed in the table of FIG. 8. These results demonstrate the efficiency of the new method, which gives the same production results at a better mass balance convergence, and at a fraction of the computational cost of the prior art method. Accordingly, it would be apparent to one of ordinary skill in the art, that the results displayed in FIGS. 6 and 7 demonstrate that the embodiment of the inventive method represents a significant improvement over the prior art.

For another example of embodiment of the method, we have run the 1-billion-cell (639*1821*934) fine grid model and the 9.5-billion-cell (2817*5463*934) fine-grid model with history and prediction period totaling hundreds of years and find the embodiment of the current inventive method to be stable and robust. Prior-art methods have great difficulties running models at these fine scales.

The second model is a full-field three-phase black oil model. A black oil model is a model of a multitude of constituents of different physical properties, in which each phase may be treated as a homogeneous medium. The reservoir has an initial gas cap, oil-water transition zone, and bottom aquifer as illustrated in FIG. 9. Multiple water oil contacts in multiple equilibrium regions are defined, and the capillary pressure curves are prescribed using the Leveret J-function, which is a function characterizing capillary pressure based on permeability and porosity functions. A detailed comparison of the nonlinear convergence behavior for this 50-million-cell model against the prior art method is conducted. By running the model with various time-stepping controls, we demonstrated that the embodiment of the current inventive method has a more consistent convergent behavior and is insensitive to varying these controls. The simulation results are very comparable as shown below in FIG. 10, and mass balance precision is consistently better than the prior-art method. Over the range of control parameters tested, the new method is about 2.61 times faster.

FIG. 10 displays a comparison of predictions of various production parameters generated by the embodiment of the current inventive method and the prior-art method. In each of chart (1002), (1004), (1006), and (1008) the dotted lines represent predictions produced by a conventional method and the solid lines represent predictions produced by an embodiment of the current invention. Specifically, chart (1002) displays predicted oil production rate on the vertical axis (1012), chart (1004) displays predicted water rate on the vertical axis (1014), chart (1006) displays predicted gas production rate along the vertical axis (1016), and chart (1008) shows the predicted average pressure of the domain of production on the vertical axis (1018). Each of charts (1002), (1004), (1006) and (1008) display their predictions as a function of the production time, indicated on their horizontal axes (1010).

FIG. 10 indicates that the embodiment of the current inventive method and the conventional method give equivalent predictions of various production parameters. However, FIG. 11 indicates that the numerical performance of the current inventive method is significantly better than that of the conventional method. Specifically, chart (1102) shows elapsed simulation time along the vertical axis (1112) as a function of production time along the horizontal axis (1010). Elapsed time required by the embodiment of the invention, indicated by the solid line, is significantly shorter than that required by the conventional method, represented by the dotted line. For at least this reason, the embodiment of the invention represents a significant improvement over the conventional method.

Similarly, chart (1104) displays the number of necessary nonlinear iterations indicated on the vertical axis (1114), chart (1106) displays the number of necessary linear iterations indicated on the vertical axis (1116), and chart (1108) indicates the number of time step cuts on the vertical axis (1118). In each of charts (1104), (1106), and (1108), the values for the embodiment of the invention are indicated by the solid line, and the values for the conventional method are indicated by the dotted line. In each of charts (1104), (1106), and (1108), the production time (1010) is indicated by the horizontal axis. Each of charts (1102), (1104), (1106) and (1108) clearly shows an improved performance of the embodiment of the invention over the conventional method.

Evaluation metrics, such as total number of time-steps, total number of nonlinear iterations, total number of linear iterations, the total number of time-step cuts, and the average time-step size are documented in the table of FIG. 12.

The third test case is for a full-field retrograde-condensate gas reservoir as shown in FIG. 13. The initial hydrocarbon in the reservoir is a single-phase wet gas, which becomes two-phase, condensate and wet gas, when the reservoir pressure decreases during production. We have successfully run this 64-million-cell simulation model with a history matching of 21 years in 1.43 hours of wall time, in a 3200-way MPI parallel simulation, using 3200 parallel processors.

FIG. 14 compares the results between the embodiment of the current inventive method and the prior-art method. The new method produces nearly identical results to the prior-art method. However, the prior-art method requires 6.27 hours on the same computer cluster using the same amount of hardware. Thus, for this model, the embodiment of the current inventive method's acceleration is 4.38 times that of the prior-art method.

FIG. 14 displays a comparison of predictions of various production parameters generated by the embodiment of the current inventive method and the prior-art method. In each of chart (1402), (1404), (1406), and (1408) the dotted lines represent predictions produced by a conventional method and the solid lines represent predictions produced by an embodiment of the current invention. Specifically, chart (1402) displays predicted gas production rate on the vertical axis (1412), chart (1404) displays predicted condensate rate on the vertical axis (1414), chart (1406) displays predicted water rate along the vertical axis (1416), and chart (1408) shows the predicted average well pressure of the domain of production on the vertical axis (1418). Each of charts (1402), (1404), (1406) and (1408) display their predictions as a function of the production time, indicated on their horizontal axes (1410).

FIG. 14 indicates that the embodiment of the current inventive method and the conventional method give equivalent predictions of various production parameters. However, FIG. 15 indicates that the numerical performance of the current inventive method is significantly better than that of the conventional method. Specifically, chart (1502) shows elapsed simulation time along the vertical axis (1512) as a function of production time along the horizontal axis (1410). Elapsed time required by the embodiment of the invention, indicated by the solid line, is significantly shorter than that required by the conventional method, represented by the dotted line. For at least this reason, the embodiment of the invention represents a significant improvement over the conventional method.

Similarly, chart (1504) displays the number of necessary nonlinear iterations indicated on the vertical axis (1514), chart (1506) displays the number of necessary linear iterations indicated on the vertical axis (1516), and chart (1508) indicates the number of time step cuts on the vertical axis (1518). In each of charts (1504), (1506), and (1508), the values for the embodiment of the invention are indicated by the solid line, and the values for the conventional method are indicated by the dotted line. In each of charts (1504), (1506), and (1508), the production time (1410) is indicated by the horizontal axis. Each of charts (1502), (1504), (1506) and (1508) clearly shows an improved performance of the embodiment of the invention over the conventional method.

The table of FIG. 16 documents the runtime performance statistics comparison between the two methods. Also, a 1.15-billion-cell fine-grid (2178*2403*220) version of this reservoir model was successfully run using the method of an embodiment of the invention, for which the prior-art method fails to complete.

In some embodiments, the reservoir simulator in FIG. 3 may be implemented in a software executed on a computer system, such as the computer system depicted in FIG. 17. The reservoir simulator may be part of the Application (1716), which may be part of the computer system (1700). Alternatively, the reservoir simulator may be implemented on dedicated computer hardware and software. Such computer hardware may be implemented and/or at a high level include components such as those depicted in the computer system described in FIG. 17.

The computer system in FIG. 17 displays a block diagram (1700) 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 system (1700) is intended to encompass any computing device such as a high-performance computing (HPC) device, 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 system (1700) 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 system (1700), including digital data, visual, or audio information (or a combination of information), or a GUI.

The computer system (1700) 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 system (1700) is communicably coupled with a network (1702). In some implementations, one or more components of the computer system (1700) 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 system (1700) 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 system (1700) 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 system (1700) can receive requests over network (1702) from a client application (for example, executing on another computer system (1700)) 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 system (1700) 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 system (1700) can communicate using a system bus (1704). In some implementations, any or all of the components of the computer system (1700), both hardware or software (or a combination of hardware and software), may interface with each other or the interface (1706) (or a combination of both) over the system bus (1704) using an application programming interface (API) (1708) or a service layer (1710) (or a combination of the API (1708) and service layer (1710). The API (1708) may include specifications for routines, data structures, and object classes. The API (1708) 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 (1710) provides software services to the computer system (1700) or other components (whether or not illustrated) that are communicably coupled to the computer (1700). The functionality of the computer (1700) may be accessible for all service consumers using this service layer. Software services, such as those provided by the service layer (1710), 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 (1700), alternative implementations may illustrate the API (1708) or the service layer (1710) as stand-alone components in relation to other components of the computer (1700) or other components (whether or not illustrated) that are communicably coupled to the computer (1700). Moreover, any or all parts of the API (1708) or the service layer (1710) 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 (1700) includes an interface (1706). Although illustrated as a single interface (1706) in FIG. 17, two or more interfaces (1706) may be used according to particular needs, desires, or particular implementations of the computer (1700). The interface (1706) is used by the computer (1700) for communicating with other systems in a distributed environment that are connected to the network (1702). Generally, the interface (1706) includes logic encoded in software or hardware (or a combination of software and hardware) and operable to communicate with the network (1702). More specifically, the interface (1706) may include software supporting one or more communication protocols associated with communications such that the network (1702) or interface's hardware is operable to communicate physical signals within and outside of the illustrated computer (1700).

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

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

In addition to holding data, the memory may be a non-transitory medium storing computer readable instruction capable of execution by the computer processor (1712) and having the functionality for carrying out manipulation of the data including mathematical computations.

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

There may be any number of computers (1700) associated with, or external to, a computer system containing computer (1700), each computer (1700) communicating over network (1702). 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 (1700), or that one user may use multiple computers (1700).

In some embodiments, the computer (1700) is implemented as part of a cloud computing system. For example, a cloud computing system may include one or more remote servers along with various other cloud components, such as cloud storage units and edge servers. In particular, a cloud computing system may perform one or more computing operations without direct active management by a user device or local computer system. As such, a cloud computing system may have different functions distributed over multiple locations from a central server, which may be performed using one or more Internet connections. More specifically, cloud computing system may operate according to one or more service models, such as infrastructure as a service (IaaS), platform as a service (PaaS), software as a service (Saas), mobile “backend” as a service (MBaaS), serverless computing, artificial intelligence (AI) as a service (AIaaS), and/or function as a service (FaaS).

FIG. 18 shows a flowchart (1800) of steps in accordance with one or more embodiments. In step (1802), a reservoir model of a subterranean region of interest is obtained, wherein the reservoir model includes a plurality of computational cells. In step (1804), we obtain a production scenario, and in step (1806), a reservoir simulator for the remainder of the following method steps is used. These steps (1808) to (1820) are contained within step (1806), see FIG. 18. Firstly, in step (1808), a first set of parameters for a first time is determined from the reservoir model and the production scenario. Secondly, in step (1810), using the first set of parameters, a second set of parameters at a second time is predicted, wherein the first and second sets of parameters include primary variables. The primary variables may include phase mole fractions and saturations of one or more of oil, gas and water, and the primary variables may be used to determine values for a mass, a flux and a well term. Each non-zero porosity grid cell contains either water, oil, and gas, water and gas, water and oil, or only water. The second time is after the first time, and the predicting includes the following steps (1812) to (1820).

In step (1812), iteratively or recursively, until a stopping criterion is satisfied, for each of the plurality of computational cells, perform the following steps from step (1814) to step (1820) inclusive. The stopping criterion may be based, at least in part, upon a root mean square (RMS) of a residual vector of at least one of the primary variables and an input residual preset tolerance. The stopping criterion may additionally or alternatively be based, at least in part, upon a total mass balance error of each primary variable over the plurality of computational cells. Alternatively or additionally, the stopping criterion may be based, at least in part, upon a total accumulated mass balance error of each primary variable over the plurality of computational cells over all time steps. As another additional or alternative choice, the stopping criterion may include a sum of total saturations equalling 1 within a specified tolerance.

In step (1814), use the first set of parameters to determine a primary variable update and a residual. The primary variable update may be determined at least in part by a Jacobian, wherein the Jacobian includes the residual differentiated with respect to the primary variables. In step (1816), determine an updated flux and an updated well term using, at least in part, the primary variable update. In step (1818), predict the second set of parameters based, at least in part, on the updated flux and the updated well term. An updated mass may also be determined, at least in part, by the updated flux and the updated well term. In predicting the second set of parameters, an accepted saturation change may be limited to a maximum of 0.2 times the updated mass, within a preset tolerance. In step (1820), assign the second set of parameters to the first set of parameters.

In step (1822), using an interpretation workstation, a predicted location of unproduced hydrocarbons may be determined based, at least in part, on the second set of parameters. Furthermore, using a wellbore planning system, a planned wellbore trajectory may be planned to intersect the predicted location, and using a drilling system, a wellbore may be drilled, guided by the planned wellbore trajectory.

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.