CN108779669A

CN108779669A - The continuous fully implicit solution well model with tridiagonal matrix structure for reservoir simulation

Info

Publication number: CN108779669A
Application number: CN201780015177.0A
Authority: CN
Inventors: 阿里·H·多鲁
Original assignee: Saudi Arabian Oil Co
Current assignee: Saudi Arabian Oil Co
Priority date: 2016-03-04
Filing date: 2017-03-02
Publication date: 2018-11-09
Also published as: WO2017151844A1; WO2017151844A8; EP3423676A1; CA3015833A1

Abstract

The subterranean hydrocarbon reservoir with horizontal well or multiple vertical wells is simulated by the continuous solution of reservoir and well equation, to simulate the flowing of the fluid in reservoir and well productivity.Well processing is defined bottom hole pressure well by the continuous solution of reservoir and well equation.This avoids the big matrix caused by the simultaneous solution of reservoir and well equation is solved, it can be computationally very expensive for a large amount of unknown number and need dedicated sparse matrix solver.This continuous solution is related to the reservoir systems solver for the rule supplemented by the minor matrix of the numerical solution of bottom pore pressure force.To the tridiagonal matrix of adjacent with well unit reservoir units at perforated interval；And the vector of the unknown reservoir gesture of adjacent reservoir units executes the solution.

Description

Continuous full-implicit well model with three-diagonal matrix structure for reservoir simulation

Cross reference to related applications

This application is a continuation-in-part application of applicant's pending co-owned U.S. patent application No.14/040,930, filed 2013, 9, 30 and claims priority, and claims priority from U.S. patent application No.13/023,728 filed 2011, 2, 9 (which is now U.S. patent No.9,164,191 issued 2015, 10, 20).

The present application also relates to commonly owned U.S. patent application Ser. No.15/061,572 entitled "Sequential full Implicit Well Model with triple Matrix Structure for Reservoir Simulation" (attorney docket No.004159.005477) and having the same inventor as the present application, filed on even date herewith.

Technical Field

The present application relates to computerized simulation of hydrocarbon reservoirs in the earth, and in particular, to simulation of flow profiles along wells in a reservoir.

Background

Well models have played an important role in numerical reservoir simulation. Well models have been used to calculate oil, water and gas production rates from wells in oil and gas reservoirs. If the well production rates are known, they are used to calculate the flow profile along the perforated section of the well. With the increased ability to measure flow rates along the perforated section of the well, an appropriate numerical well model is necessary to calculate the correct flow profile to match the measurements in the reservoir simulator.

It is generally known that if all reservoirs are vertically connected for any vertical well in the reservoir simulator, a simplified well model (e.g., an explicit or semi-implicit model) may be sufficient. For these models, well productivity is assigned to the perforations in proportion to the zone production index (or total mobility). Thus, the calculation is simple and computationally inexpensive. The structure of the resulting coefficient matrix of reservoir unknowns remains unchanged. In particular, the coefficient matrix maintains a regular sparse structure. Thus, any such sparse matrix solver can be used to solve a linear system of grid block pressures and saturations for each time step of the entire reservoir simulation model.

However, for highly heterogeneous reservoirs with some vertical non-communicating layers, the well model mentioned above does not yield the correct physical understanding. However, they produce incorrect flow profiles and, in some cases, cause simulator convergence problems.

As reservoir models increase in complexity, the number of vertical slices has reached the order of hundreds to represent reservoir heterogeneity. A fully implicit fully coupled well model with simultaneous solutions of reservoir and well equations is already necessary to correctly simulate the flow profile along the well, and is also necessary for the numerical stability of the reservoir simulation. To solve for a fully coupled system, the well equations are typically first eliminated. This produces an unstructured coefficient matrix of reservoir unknowns to be solved. The solution of this type of matrix requires a dedicated solver with a dedicated preconditioner. For wells with many completions and many wells in the simulation model, this approach has become computationally expensive in terms of processor time.

Disclosure of Invention

Briefly, the present invention provides a new and improved computer implemented method for solving well equations as well as reservoir equations in a reservoir simulation model having a formation with vertical fluid flow and a flow barrier without vertical fluid flow. By forming the flow layers having flow between the fluid barriers in the reservoir model only in the vicinity of the wells as a single flow layer, the computer-implemented method constitutes a simplified well model system of horizontal and multiple vertical wells in the reservoir simulation model. For a well of a given production rate, the method then solves the simplified well model system by matrix calculation (using a direct sparse matrix solver) for well and reservoir unknowns of bottom hole pressures for the grid blocks through which the well passes (simplified well model system). The method is repeated for one or more wells in the reservoir simulator. The method then solves the full reservoir simulation model by treating each well as having a determined bottom hole pressure and based on steady state volume balance relationships of formation completion rates (completion rates), formation pressures, and hydraulic conductivity to determine completion rates of one or more wells of the full reservoir model and determines a total rate of the well from the determined completion rates of the one or more wells. The method then records the determined completion rate of the formation and the determined overall rate of the one or more wells. The simulator calculates a pressure distribution of a multiphase flow fluid saturation distribution in the reservoir using the reservoir data assigned to the simulator grid blocks and the calculated perforation rates. Thus, the simulator is able to calculate the pressure and saturation distribution within the reservoir at a given time given the production and injection rate at the well.

The present invention provides a new and improved data processing system for forming a well model for horizontal and multiple vertical wells through reservoir simulation of a subsurface reservoir from a reservoir model having a formation with fluid flow and a flow barrier without fluid flow. The data processing system includes a processor that performs the steps of: the simplified well model system model is solved by matrix calculations to obtain bottom hole pressures for one or more wells and reservoir unknowns around the one or more wells, and the full reservoir model is solved by processing the well(s) to have the determined bottom hole pressures and based on steady state volume balance relationships of formation completion rates, formation pressures, and water conductivity to determine completion rates for each layer of the well(s) of the full reservoir model. The processor also determines a total rate for the well(s) based on the determined completion rates for the layers of the well(s). The data processing system also includes a memory that forms a record of the determined completion rates and the determined one or more aggregate rates for each of the layers.

The present invention also provides a new and improved data storage device having stored in a computer readable medium computer operable instructions for causing a data processor to form a well model of horizontal and multiple vertical wells by reservoir simulation of a subsurface reservoir from a reservoir model, the reservoir simulation having a formation with fluid flow and a flow barrier with no fluid flow, to perform the steps of: the simplified well model system model is solved by matrix calculations to obtain bottom hole pressures and reservoir unknowns for the one or more wells, and the full reservoir model is solved by processing the well(s) to have the determined bottom hole pressures and based on steady state volume balance relationships of formation completion rates, formation pressures, and water conductivity to determine completion rates for the layers of the well(s) of the full reservoir model. The instructions stored in the data storage device further include instructions that cause the data processor to determine a total rate for the well(s) based on the determined completion rates for the zones of the well(s), and form a record of the determined completion rates for the zones and the determined total rate for the well(s).

Drawings

Fig. 1A and 1B are schematic diagrams of a well model in a reservoir simulator of a plurality of subsurface formations formed as a single layer above and below a fluid barrier in a reservoir according to the present invention.

Fig. 2A and 2B are schematic diagrams of a well model in a reservoir simulator of a plurality of subsurface formations formed as a single layer above and below a number of vertically spaced fluid barriers in a reservoir in accordance with the present invention.

FIG. 3 is a schematic diagram of a well model for simulation based on an explicit modeling approach shown in radial geometry.

FIG. 4 is a schematic diagram of a well model for simulation based on a fully implicit fully coupled model approach.

FIG. 5 is a schematic diagram of a well model of a reservoir simulation with a comparison of flow profiles obtained from the models of FIGS. 3 and 4.

FIG. 6 is a schematic diagram of a finite difference grid system for the well model of FIGS. 3 and 4.

Fig. 7A and 7B are schematic diagrams of a well model reservoir showing an unmodified conventional well model and a well model according to the present invention, respectively.

Fig. 8A and 8B are schematic diagrams of flow profiles illustrating a comparison between the models of fig. 7A and 7B, respectively.

FIG. 9 is a schematic illustration of a well formation model with two fracture layers in radial coordinates.

FIG. 10 is a functional block diagram or flow chart of the data processing steps of the method and system for a fully-implicit, fully-coupled well model with simultaneous numerical solutions for reservoir simulation.

FIG. 11 is a functional block diagram or flow diagram of the data processing steps of the method and system for continuous fully-implicit well model for reservoir simulation according to the present invention.

FIG. 12 is a schematic of a linear system of equations with a three diagonal coefficient matrix for an explicit well model for a one-dimensional reservoir simulator with vertical wells.

FIG. 13 is a schematic of a linear system of equations with a three diagonal coefficient matrix for a fully implicit well model for a one dimensional reservoir simulator with vertical wells.

FIG. 14 is a schematic illustration of a linear system of equations with a three diagonal coefficient matrix for a constant bottom hole pressure well model for a one-dimensional reservoir simulator with vertical wells according to the present invention for processing.

FIG. 15 is a functional block diagram or flow chart illustrating the steps of an analysis method of reservoir simulation according to the present invention.

FIG. 16 is a schematic diagram of a computer network of a continuous fully-implicit well model for reservoir simulation according to the present invention.

FIG. 17 is a schematic diagram of a three diagonal coefficient matrix.

FIG. 18 is a schematic of a linear system of equations with a three diagonal coefficient matrix.

FIG. 19 is a schematic diagram of a finite difference grid system of a horizontal well model oriented in the y-axis direction according to the present invention.

Fig. 20A and 20B are schematic diagrams of a horizontal well model in a reservoir simulator of a plurality of subsurface formations in a region of flow impairment in a reservoir before and after formation as a simplified horizontal well model, in accordance with the present invention.

FIG. 20C is a schematic of a linear system of equations of the simplified horizontal well model of FIG. 20B.

Fig. 21 and 22 are functional block diagrams or flowcharts of data processing steps of a simplified horizontal well model method and system according to the present invention.

FIG. 23 is a schematic diagram of a finite difference grid system of a reservoir model with multiple wells according to the present invention.

Fig. 24A and 24B are schematic diagrams illustrating symbols of grid naming in the multiple well reservoir model of fig. 19 and 23.

Fig. 25A, 25B, 25C, and 25D are schematic diagrams illustrating grid numbering of a three-dimensional reservoir model.

FIG. 26 is a functional block diagram or flow chart of data processing steps of the method and system of a modification of the block diagram of FIG. 11 for a plurality of wells in a reservoir having a fully-implicit, fully-coupled well model of a simultaneous numerical solution to the reservoir simulation.

FIG. 27 is a schematic diagram of a simplified linear system of equations for a two-well, three-dimensional reservoir model.

Detailed Description

By way of introduction, the present invention provides a continuous fully-implicit well model for reservoir simulation. Reservoir simulation is the mathematical modeling science of reservoir engineering. Fluid flow within an oil or gas reservoir (porous medium) is described by a set of partial differential equations. These equations describe the pressure (energy) distribution, oil, water and gas velocity distribution, and the volume fractions (saturations) of oil, water and gas at any point in the reservoir at any time during the life of the reservoir from which the oil, gas and water are produced. Fluid flow within the reservoir is described by tracking movement of constituents of the mixture. Expressed as components (e.g. methane, ethane, CO) in unit mass or mole₂Nitrogen gas, H₂S and water).

Since these equations describing fluid flow and the associated thermodynamic and physical laws are complex, they can only be solved on a digital computer to obtain the amount of pressure distribution, velocity distribution, and fluid saturation or component mass or molar distribution within the reservoir at any point in time. This can only be done by solving the equations numerically, rather than analytically. Numerical solutions require that the reservoir be subdivided in regional and vertical directions (x, y, z-three-dimensional space) into computing elements (cells or grid blocks), and time into intervals of days or months. For each element, the unknowns (pressure, velocity, volume fraction, etc.) are determined by solving complex mathematical equations.

In fact, a reservoir simulator model can be thought of as a collection of rectangular prisms (like bricks in a building wall). Variations in the pressure and velocity fields occur due to oil, water and gas production at wells distributed within the reservoir. The simulation is performed over time (t). Typically, the production or injection rate of each well is known during the production history of the reservoir. However, since the well traverses several reservoirs (elements), the contribution of each reservoir element (well perforation) to production is calculated by a different method. The present invention relates to the calculation of how much each well perforation contributes to the total well production. Since these calculations are expensive and very important boundary conditions for the simulator, the proposed method suggests a practical method for correctly calculating the flow profile along the well trajectory. As will be described, it can be shown that: some other methods used will result in incorrect flow profiles, which can cause problems in obtaining the correct numerical solution, and these methods can be computationally expensive.

Name of

Grid size in x, y and z directions

k_x，k_y，k_zPermeability in x, y and z directions

p is pressure

Rho ═ fluid (petroleum) density

g is gravitational constant

z is the vertical depth from the depth of the reference plane

r_oPicleman (Peaceman) radius 0.2 Δ x

r_wRadius of borehole

T_x,T_y,T_zRock hydraulic conductivity in x, y and z directions

An equation describing a generic reservoir simulation model and indicating the well terms of interest relevant to the present invention is set forth below in equation (1):

wherein, Delta_xIs a difference operator in the x-direction of the reservoir, T_x,T_y,T_zIs the rock water conductivity in the x, y and z directions as defined in the following equation), j represents the number of fluid phases, n_pIs the total number of fluid phases, typically 3 (oil, water and natural gas), sigma is the summation term, p_i,jIs the density, λ, of component i in fluid phase j_jIs the fluidity of phase j (equation 6), Φ_jIs the fluid potential (reference surface corrected pressure) of fluid phase j, and similarly, Δ_yIs a difference operator in the y-direction of the reservoir, and Δ_zIs a difference operator in the z-direction of the reservoir, q_i，w，kIs the well term (source or sink), Δ, of component i of grid block (cell) k_tIs a difference operator in the time domain, n_iIs the total number of moles of component i, and n_cIs the total number of components (methane, ethane, propane, CO) in the fluid system₂Etc.).

Equation (1) is a system of partial differential equations describing the nonlinear coupling of fluid flow in the reservoir. In the above equation set, n_iRepresenting the ith component of the fluid. n is_cIs the total number of components of hydrocarbon and water flowing in the reservoir. In this context, by component is meant, for example, methane, ethane, propane, CO₂、H₂S, water and the like. The amount of the component depends on the hydrocarbon water system available to the reservoir of interest. In general, the number of components can vary from 3 to 10. Equation (1) combines the continuity equation and the momentum equation.

In equation (1), q_i，w，kIs the position x of component i_k，y_k，z_kThe well perforation rate, and k is the grid block (cell) number. Again, it is the subject of the present invention to make this calculation based on the specified production rate at the well head.

In addition to the differential equation in equation (1), the pore volume constraint at any point (element) in the reservoir must satisfy:

wherein, V_pIs the cell pore volume, p (x, y, z) is the fluid pressure at point x, y, z, N_jIs the total number of moles in fluid phase j, ρ_jIs the density of fluid phase j.

In equations (1) and (2) there is n_c+1 equations and n_c+1 unknown. For each component i, these equations are solved simultaneously with the thermodynamic phase constraints by equation (3):

wherein f is_iIs the component viscosity, the superscript V represents the vapor phase, L represents the liquid phase, n_iIs the molar total of component i, P is the pressure and T represents the temperature.

In a fluid system, there are typically three fluid phases in the reservoir: petroleum phase, natural gas phase and water phase. Each fluid phase may contain different amounts of the above-described components based on reservoir pressure and temperature. The fluid phase is depicted by the symbol j. The symbol j has a maximum value of 3 (petroleum phase, aqueous phase and natural gas phase). Symbol n_pIs the maximum number of phases (which can sometimes be 1 (oil), 2 (oil and gas or oil and water), or 3 (oil, water and gas)). Number of phases n_pBased on reservoir pressure (p) and temperature (T). Symbol n_iIs the total number of moles of component i in the fluid system. Symbol n_cIs the maximum number of components in the fluid system. Determining the number of phases and each phase n according to equation (3)_i，jFraction of each component in (1) and phase density ρ_jAnd ρ_i，j. In equation (3), V represents the vapor (natural gas) phase and L represents the liquid phase (oil or water).

The total number of moles in fluid phase j is defined by the following formula:

the total component moles are defined by equation (5).

The definitions of the phase fluidity, the relationship between the phases, the fluid potential, and the differential sign in equation (1) are defined in equations (6) to (9).

λ_j＝k_r，j/μ_j(6)

In equation (6), the numerator defines the relative permeability of the phase and the denominator is the viscosity of the phase.

The capillary pressure between the phases relative to the phase pressure is defined by equation (7):

the fluid potential of phase j is defined by:

Φ_j＝P_j-gρ_jz (8)

the discrete differential operators in the x, y and z directions are defined by:

where Δ defines the sign of the discrete difference and U is an arbitrary variable.

Writing equation (1) and the constraints in equations (3) through (9) using a controlled volume finite difference method for each grid block (cell) (some of the grid blocks may include wells) in the reservoir simulatorAnd definition. The resulting equations are solved simultaneously. Finding a given well productivity q for each well after completion_TN of (A) to (B)_i(x, y, x, and t), and P (x, y, z, and t), the component rates in equation (1) are calculated from the given well production rate using the new well model formulation according to the present invention. To solve equation (1) and equation (2), the reservoir boundary in (x, y, z) space, the rock property distribution K (x, y, z), the rock porosity distribution and the fluid properties and saturation related data are entered into the simulation.

In accordance with the present invention and as will be described below, a simplified well model system is formed that yields the same calculated bottom hole pressure determination as the complex, computationally time consuming prior fully coupled well model.

According to the present invention, it has been determined that for a grid block where well trajectories pass through several formations to communicate vertically, the communicating layers can be combined to treat as a single layer, as schematically indicated for the well model in fig. 1 and 2. This is done by: identifying vertical flow barriers in a reservoir of well units, and combining layers above or below various flow barriers of well units. Thus, the full well model system is reduced to a smaller sized well model system with many fewer layers to incorporate into the well model for processing.

As shown in fig. 1, the well model L represents, in simplified schematic form, a complex subsurface reservoir grid block (cell) through which a well passes, which is made up of seven individual earth formations 10, each of the earth formations 10 being in vertical fluid communication with an adjacent layer 10. The model L includes another set of ten formations 12 around the well 12, each of the formations 12 being in vertical fluid communication with an adjacent layer 12. The groups of formations 10 and 12 in the model L that are in fluid communication with other similar adjacent layers are separated as shown in fig. 1 by a fluid-tight barrier layer 14 that forms a barrier to vertical fluid flow.

According to the invention, the well model L is simplified for processing purposes to a simplified or reduced well model R (fig. 1A) by bringing together or combining the layers 10 above the fluid barriers 14 of the well model L into a composite layer 10a in the reduced model R for the purpose of determining potential Φ and completion rate. Similarly, the layers 12 below the fluid barrier 14 of the well model L are combined into a composite layer 12a in the simplified well model R.

Similarly, as indicated in FIG. 2, the reservoir model L-1 is made up of five upper individual formations 20, each of the formations 20 being in vertical fluid communication with an adjacent layer 20. Model L-1 includes another set of seven strata 22, with each of the strata 22 being in vertical fluid communication with an adjacent layer 22. The groups of formations 20 and 22 in the model L-1 in fluid communication with other like adjacent layers are separated as indicated in fig. 2 by a fluid-tight barrier layer 24 that forms a barrier to vertical fluid flow. As indicated in model L-1, another set of nine earth formations 25 in fluid communication with each other below the fluid barrier layer 26 as a vertical fluid flow barrier is spaced from the layer 22. The lowermost set of ten formations 27 in fluid communication with each other is located below the fluid flow barriers 28 in the model L-1.

In accordance with the present invention, the well model L-1 is reduced to a simplified or reduced well model R-1 (FIG. 2A) for disposal purposes by bringing together or combining the layers 20 above the fluid barriers 24 of the well model L-1 into a composite layer 20a in the reduced model R for purposes of determining the well layer potential Φ and completion rate. Similarly, the other layers 22, 25 and 27 below the fluid barriers 24, 26 and 28 of the well model L-1 are combined into composite layers 22a, 25a and 27a in the simplified well model R-1.

A simplified well model system or well model according to the present invention is solved for reservoir unknowns and bottom hole pressure. The wells in the full reservoir simulation model system are then considered to be the specified bottom hole well pressure and the reservoir unknowns are solved implicitly. The diagonal elements and right-hand term vectors of the coefficient matrix of the reservoir model are only the components that are modified in the process according to the invention, and these are only slightly modified. A rule sparse solver technique or method is then used to solve for reservoir unknowns. The perforation rate is calculated by using the reservoir unknowns (pressure and saturation). These rates are then summed to calculate the total well rate. The error between the determined overall well rate and the input well rate according to the present invention will decrease with the non-linear newton iteration of the simulator for each time step.

The flow rates calculated according to the present invention also converge with the flow rates calculated by the fully-coupled simultaneous solution for the entire reservoir simulation model, which includes many wells. Because the method needs to solve the small well system model, the calculation cost is low. It has been found that the method of the present invention converges if the simplified well system is suitably constructed by suitably using upscaling (upscaling) in combining the communicating layers.

It is generally known that if all reservoirs are vertically connected, then a simple well model (e.g., an explicit or semi-implicit model) is usually already sufficient. As shown in FIG. 3, the explicit well model E is comprised of a number Nz of reservoirs 30 in vertical flow communication, each having a permeability k defined as indicated in FIG. 3_x，i(in this context, i denotes the layer number, not the component) and the thickness Δ z_iAnd rate of perforation of layer q_i. Then, in FIG. 3, as indicated in equation (3) in the same figure, the overall productivity q of the explicit model E_TSingle production rate q of Nz layers being an explicit model_iThe sum of (1).

For the explicit model, well productivity is assigned to perforations in proportion to the zone production index (or total mobility). Therefore, the calculation is simple. The resulting coefficient matrix of unknowns remains unchanged, i.e., a regular sparse structure is maintained, as shown in the matrix format in fig. 12. Thus, any sparse matrix solver can be used to solve for a linear system of grid block pressures and saturations for each time step.

Well model

The approach of several vertical well models of reservoir simulators is also presented simply based on a simple fluid system in the form of a flow of slightly compressible single-phase petroleum flow within the reservoir. However, it should be understood that the invention is generally applicable to reservoirs and can be used with any number of wells and fluid phases in a regular reservoir simulation model.

FIG. 6 shows a view of the illustration used in this specificationFinite difference grid G for vertical well models. As shown, the well is located at the center of the central unit in the vertical direction. The model presented below also takes into account the well completed in the vertical Nz direction, and the potentials in adjacent cells: phi_BW,Φ_BE,Φ_BN,Φ_BSIs constant and known from the simulation run (previous time step or iteration value). Here, subscript B refers to "boundary", W refers to the west-side neighboring cell, E refers to the east-side neighboring cell, and N and S refer to north and south, respectively. Also, Φ describes the fluid potential (reference surface corrected pressure).

As shown in fig. 6, the well penetrates several reservoirs in a vertical direction represented by the index I (I ═ 1,2,3 … Nz), where Nz is the total number of vertical layers in the reservoir model. For each layer i, there are four neighboring cells at the same area plane (x, y). These neighboring cells are located east-west (x-direction) and north-south (y-direction).

The potentials Φ of the east, west, north, and south side neighboring cells are known from the simulator time step calculation, and those potentials are set to assume that they do not change with the simulator time step. The vertical potential variation at the well location is then considered, but it is assumed that the adjacent potentials that can vary in the vertical direction are known.

The steady state volume balance equation for cell (i) in FIG. 6 is as follows:

T_wi(Φ_Bw-Φ_i)+T_Ei(Φ_BE-Φ_i)+T_Ni(Φ_BN-Φ_i)+T_Si(Φ_BS-Φ_i)+ (10)

T_Up,i(Φ_i-1-Φ_i)+T_Down,i(Φ_i+1-Φ_i)-q_i＝0.

in equation (10), T represents the water conductivity between cells. Subscripts W, E, N, and S denote the west, east, south, and north directions, and (i) denotes the unit index.

The water conductivity between the cells in three directions is defined by the following equation (11):

in the above equation (11), Δ x_iIs the lattice block size (cell size) in the x direction of the lattice block (cell) number i. Similarly, Δ y_iIs the grid block size (cell size) in the y direction of the grid block (cell) number i, and Δ z_iIs the grid block size (cell size) or the grid layer thickness in the z direction of the grid block (cell) number i. K_z，i-1/2Is the vertical permeability at the intersection of cell i and cell i-1. Similarly, K_z，i+1/2Is the vertical permeability at the intersection of cell i and cell i + 1. As shown in FIG. 6, cell i is at the center and i +1/2 is at the boundary between cell i and cell i + 1. For simplicity, when indicating east-west flows, the subscript j is added. Similarly, (i, j-1) is the neighbor to the north of the central cell (i, j). Thus, the symbol (i, j-1/2) is the interface between the central cell and the north neighboring cell in the y-direction.

The same notation is used for adjacent cells on the south side: (i, j +1/2) denotes a boundary between the central cell and the south-side neighboring cell (i, j +1) in the y direction. For simplicity, in the above equation, in representing the y (or j) direction, the index i is added. Obviously, the hydraulic conductivity in equation (11) is defined in a similar manner.

In fig. 14, the diagonal term is defined by the following equation (17)And the water conductivity between the cells in the three directions is as defined in equation 11 described above.

Item in FIG. 14(right-hand term) is expressed by equation (17b) herein. Extracted by the following formula:

where i is 1,2, … Nz, where Nz is again the total number of cells in the vertical direction, the number of vertical layers. In the extraction of equation (17b), the layer productivity index PI_iDefined by equation (17), and a potential term Φ_BIs a known boundary potential at the boundary of each of the adjacent cells to the west, east, north and south of the central cell. Typically, terms, cells, and gridblocks are used interchangeably.

Conventional well models can be generally classified into three groups: (a) an explicit well model; (b) a well model for bottom hole pressure specification; and a fully implicit well model (Azizz K, Settari A, Petroleum Reservoir Simulination, applied science Publishers Ltd, London 1979). For a better understanding of the invention, each well model is briefly introduced.

Explicit well model

For an explicit well model, the source term q in equation (10)_iAccording to the followingEquation (12) is defined:

wherein q is_iIs the production rate of a cell (grid block) i, where the well passes through and is perforated.

Substituting equation (12) into equation (10) for cell i

T_UpiΦ_i-1+T_C,iΦ_i+T_Down,iΦ_i+1＝b_i(13)

Wherein,

T_c,i＝-(T_Up,i+T_down，i+T_Wi+T_Ei+T_Ni+T_Si) (14a)

and

writing equation (13) for all cells i around a well of a well cell to 1, Nz results in only a linear system of equations with a tri-diagonal coefficient matrix of the type shown in fig. 12, which can be written in matrix vector notation as:

in equation (15), A_RRIs a (Nz x Nz) three diagonal matrix, and phi_RAnd b_RIs the (Nzx 1) vector. Solving equation (15) by computer processing the grid blocks through which the well passes using a tri-diagonal linear system solver (e.g., Thomas algorithm) to obtain the reservoir unknown potential Φ_R。

Three diagonal matrices and systems

A three-diagonal matrix is a matrix with only three diagonals in the middle, with real or complex numbers as entries in the diagonals. These diagonals are referred to as the "lower diagonal", "center diagonal" and "upper diagonal". The remaining elements or entries in the three diagonal matrix are zero. For example, fig. 17 shows a three-diagonal matrix a of 8 × 8 matrix size (or order n ═ 8). In fig. 17, element a; the ith element A, i.e. a₁,a₂,a₃,a_iI is 1,8 denotes the lower diagonal; b₁,b₂,b₃,b_iThe central diagonal line; and c₁,c₂,c₃,c_iThe upper diagonal element.

An example three diagonal system of equations is shown in FIG. 18. The solution to the linear system of equations with the three diagonal matrix as described above and shown in fig. 18 is easily solved by gaussian elimination. An example of a solution to an equation in this manner is contained at https:// en.

Therefore, x shown below and in fig. 18 is performed by solving the matrix relationship of equations (16a) to (16e) as follows_iThe solution of (a):

b′_i＝1；i＝1，2，...，n (16b)

well model with bottom hole pressure regulation

Φ_wIs a uniform potential along The wellbore opening to production using conventional techniques according to The techniques explained in The literature (e.g., in The text of Muskat, "Physical Principles of OilProduction," McGraw-high Book Co. (1949) and "The Flow of Homogeneous fluids through port Media," McGraw-high Book Co. (1937)). For purposes of modeling in this context, the fractional pressure drop along the well is considered negligible. Suppose phi_wAs is known (or specified), the oil production rate (oil rate) from the perforation is calculated by the following equation (17):

wherein, PI_iIs the layer productivity index, [ phi ] w is the specified bottom hole potential (base plane corrected pressure), [ phi ] i is the reservoir grid block pressure of the grid block (cell) i through which the well passes, r_o,iThe well block radius of Pesmann, referred to as grid block i, is defined as 0.2 Δ x, r_wIs the radius of the well.

The variables in equation (17) are explained in the "naming" section above. Substituting equation (17) into equation (10) and collecting the terms for cell i, the following results occur for cell i:

T_UpiΦ_i-1+T_C,iΦ_i+T_Down,iΦ_i+1W＝b_i(18)

order to

T_ci＝-(T_Up,i+T_down,i+T_Wi+T_Ei+T_Ni++T_Si+PI_i) (18a)

b_i＝-(PI_iΦ_W+T_WiΦ_BW+T_EiΦ_BE+T_NiΦ_BN+T_SiΦ_BS) (18b)

When equation (18) is written for all the grid blocks i — 1, Nz, the matrix system shown in fig. 14 results. It can be seen that the matrix for the bottom hole pressure specified well model in fig. 14 is similar to that of fig. 12, and in a comparable manner, equation (18) is similar to equation (13). The bottom hole pressure specified well model can be easily solved by matrix computer processing using three diagonal equation solver methods.

Fully-concealed well model

Specifying a total production rate q for a well according to a fully implicit well model according to equation (19)_T。

Calculating the individual completion rate q from equation (17)_i. For implicit well models, the wellbore potential Φ_WIs assumed to be constant throughout the well but is unknown.

Substituting equation (19) into equation (10) and collecting the terms for cell i, for cell (i), the following expression is obtained:

T_UpiΦ_i-1+T_C,iΦ_i+T_Down,iΦ_i+1W-PI_iΦ_W＝b_i(20)

order to

T_ci＝-(T_Up,i+T_down,i+T_Wi+T_Ei+T_Ni+PI_i) (20A)

b_i＝-(T_WiΦ_BW+T_EiΦ_BE+T_NiΦ_BN+T_SiΦ_BS) (20B)

Writing equation (20) for all cells yields a cell having the form shown in FIG. 13The upper diagonal solid line represents T as defined by equation (11)_Up,iAnd the lower diagonal solid line describes what is referred to as T described in equation (11)_Down,iOf (2) is used. Central item T_C,IDefined by equation (20A), and the right-hand term b_iDefined by equation (20B).

The linear system of the matrix (equation 20) of fig. 13 can be represented using vector matrix notation as follows:

in equation (21), A_RRIs a (Nz x Nz) three diagonal matrix, A_RWIs (Nz x1) vector (of reservoir PI), A_WRIs a (1x Nz) vector (of PI), and A_WWIs a (1x1) scalar. For this example:

write algebraically:

solving phi from equation 22_WObtaining:

into equation (21)

Collecting the terms in equation (25) yields:

coefficient matrix (A) of equation (26)_RR-A_ww ^-1A_wRA_RW) Is a (Nz x Nz) full matrix.

The resulting coefficient matrix can be defined as:

and

equation (26) can be written as:

the matrix of equations (28) can be solved by using a direct solver of the gaussian elimination method or any other suitable conventional solver for a full matrix.

If in the implicit well model the number of zones Nz is large and the well is fully completed in all zones, the solution of equation (28) becomes computationally expensive. To this end and in addition, if many wells are involved, equation (28) is typically solved by an iterative method. This example is described in "A full-Imperial full-Coupled Well Model for Parallel Mega-Cell Reservoir Simulation", SPE Technical Symposium of Saudi Arabia Section,14-16 May 2005.

Flow chart I (fig. 10) indicates the basic computer processing sequence of a fully implicit fully coupled well model simultaneous solution of the matrix type shown in fig. 13. During step 100, the simulation is started by reading reservoir data and production data. Reservoir data includes reservoir geometry information-its size, extent (length) in the x, y and z directions, reservoir properties (e.g., permeability distribution, porosity distribution), thickness of the layer, relative permeability data, capillary pressure data, fluid property data (e.g., fluid density table, formation volume factor table, viscosity table), location of wells within the reservoir, location of well perforations.

The production data includes the measured or specified production rates of oil, water and gas for the well as defined in the previous step. The production data also includes a minimum bottom hole pressure for each well.

In many cases, only oil production rates and water production rates are entered if gas production data is not available. If no water is produced in the field, only oil production rates are read as inputs.

During step 102, the time step is increased by one and the iteration count of the number of non-linear iterations performed during this time step is set to zero. During step 104, a Jacobian matrix of reservoir data is formed. In step 106, the resulting linear system of equation (19) is then solved by an Iterative method using a sparse preprocessor (Youef Saad, Iterative Methods for sparse Linear Systems, Society of Industrial & Applied Mattics (SIAM) Publication, 2003).

During step 108, a convergence step is performed to determine whether the non-linear iteration has converged. The individual residuals (residual) of the equations obtained in step 106 are checked against a user-specified tolerance (tolerance). If these tolerances are met, the non-linear iterative loop is exited, the solution output is written to the file for the current time step, and processing returns to step 102 for the time step to be increased, and processing continues for the increased time step, as shown. If the user-specified tolerance is determined not to be met during step 108, processing according to the non-linear iterative loop returns to step 104 and continues. If the number of non-linear iterations becomes too large, it may be decided to decrease the size of the time step and return to step 102.

However, based on the strength of the preprocessing, this method can also be very expensive in computation time, since there is no precise way to represent the full matrix in equation (19). For difficult equations, the iterative method may not converge with extreme non-uniformities and small layer thicknesses.

Furthermore, for highly heterogeneous reservoirs with some vertically non-communicating layers, the well model described above does not yield the correct physical understanding. Instead, they can produce incorrect flow profiles, and in some cases, they can lead to problems with simulator convergence.

The invention

In the case of some vertical flow barriers, explicit and fully implicit models can produce completely different flow profiles. This is shown in fig. 5. As shown in fig. 5, the reservoir model V is made up of a relatively low permeability upper layer 50 and, in the case of vertical flow, is positioned above an isolated high permeability layer 52 without vertical flow communication with an adjacent layer. A fluid barrier 52 is positioned above a layer 54 of medium permeability and vertical flow communication in the reservoir V. As can be seen in fig. 5, the production rate qT of the well model V is the same for both implicit and explicit well modeling methods.

However, the productivity of the layers 50, 52 and 54 q1, q2 and q3 are significantly different. The production rate q2 of the layer 52 contributes mainly to the total production rate qT shown by the curve 56 of the explicit model. In contrast, for the implicit model shown by curve 58, the production rate q2 for layer 52 is significantly less.

In the fully implicit well model, the productivity of layer 2 is significantly less due to the fact that: the method takes into account internal reservoir heterogeneity (all reservoir properties in the reservoir) as well as heterogeneity around the well block except for the perforation index (or explicit methods use only the properties of layer 2 as the only data to assign a rate score to the perforation). For example, the fully-implicit well method focuses on no fluid being supplied to layer 2 from layers 1 and 3 due to a fluid-tight barrier between layer 2 and layers 1 and 3. Thus, once some fluid is produced from the well model layer 2, the pressure of layer 2 should drop and the layer should not be supplied into the well at a high rate, despite the very high permeability of the layer.

On the other hand, the explicit method allocates the rate of layer 2 based only on the permeability of the layer without considering the layer connection with the upper and lower layers. Based on this, the explicit method will assign a very high rate to this layer and will hold it at the next analog time step. In a later simulation time step this will lead to instability in productivity, i.e. the simulator will reduce the time step size and will take a very long time to complete the simulation. For the same input data, it is not satisfactory for reservoir simulations to have the model provide different results based on the modeling technique selected for use.

According to the present invention, a more comprehensive well model is provided. The well model according to the invention is referred to as a coupled reservoir well model. The associated numerical solution is referred to as a fully implicit fully coupled and simultaneous solution. The fully implicit fully coupled reservoir well model produces the correct fluid profile along the perforated well interval as will be described. As shown in fig. 4, the reservoir model O is composed of i individual layers 1 to Nz of number z (each layer having a permeability k as defined in fig. 4_x,iAnd a thickness Δ z_iAnd potential phi_i) And upper and lower layers 40 of relatively low permeability and, in the case of vertical flow, above and below the isolated high permeability layer 42, respectively, without vertical flow communication with adjacent layers. As can be seen in FIG. 4, the production rate q of a layer i of the well model O_iIndicated by the expressions listed in fig. 4.

The fully-implicit coupled reservoir well model V of fig. 5 is set forth in equation (20) above, and is also set forth herein in a matrix format as further described:

the invention is based on the following facts: the bottom hole pressure of a stratified reservoir with vertical wells is the same as the system according to the invention. The system according to the invention is formed by identifying fluid obstructions and bringing together or merging reservoirs communicating in a vertical direction around the well for treatment purposes. Care should be taken in forming the simplified system according to the present invention. The simplified system must be formed correctly, otherwise errors in the construction of the simplified system can increase the total number of non-linear newton iterations.

The system according to the invention solves for the bottom hole pressure. The solution is performed by treating the well as a prescribed bottom hole pressure well. The process is fully implicit; however, it is not a simultaneous solution. Instead, the solution is continuous. The method of the present invention is convergent in that it is part of the global newton iteration of the simulator. Thus, if the model according to the invention is constructed correctly, any possible error in the rate calculation will be small and will decrease with the simulator newton iterations.

Flow chart F (fig. 11) illustrates the basic computer processing sequence according to the present invention and the computational method that occurs during an exemplary embodiment of a continuous full implicit well model (sequential full implicit well model) utilizing the reservoir simulation of the present invention.

During step 200, the simulation begins by reading reservoir data and production data. The reservoir data and production data read during step 200 are of the types described above. During step 200, the reservoir simulator is also initialized, with the simulation date and time step set to zero. During step 202, the time step is increased by one and the iteration count of the number of non-linear iterations performed during this time step is set to zero.

During step 204, a Jacobian matrix of reservoir data is formed. In step 206, a simplified system is formed according to the invention as defined by the model R described above and in the manner described above for the bottom hole potential Φ_wAnd solving the matrix. In the step ofDuring 208, the modified linear system matrix a is solved in accordance with the present invention in the manner described above.

FIG. 15 illustrates forming a simplified well model system matrix R and solving for the bottom hole potential Φ according to steps 204 and 206 of FIG. 11_wThe method of (1). Vertical flow obstructions in the original well model system are identified, as indicated at step 210. This can be done based on well log data or by specification by the reservoir analyst from data in the original reservoir model.

After step 210, a simplified well model system is then formed by the data processing system D during step 212. Those layers in the well model that are between the flow barriers and have vertical flow are combined together for the purpose of analyzing the model.

Next, during step 214, the resulting simplified well model system is solved by computer processing using the techniques discussed in equations (17), (17a), and (17b) to obtain the bottom hole potential Φ_wAnd reservoir unknowns. Step 216 then solves the full well model system structural matrix of equation (27) using a direct solver or other suitable technique as described above with data processing system D. Completion rate q is then determined based on the results of step 216 using data processing system D_iAnd the total well flow rate q_T。

Referring again to FIG. 11, during step 220, a convergence step is performed to determine whether the non-linear iteration converges. The individual residuals of the equations resulting from step 216 are checked against a user specified tolerance. If these tolerances are met, the non-linear iterative loop is exited, the solution output is written to the file at the current time step, and processing returns to step 202 for the time step to be increased, and processing for the increased time step continues, as shown. If the user-specified tolerance is determined not to be met during step 208, processing according to the non-linear iterative loop returns to step 204 and continues. If the number of non-linear iterations becomes too large, it may be decided to adjust the model.

Horizontal well

Fig. 19 shows horizontal well 300 oriented in the y-direction in three-dimensional reservoir model H. Model H of fig. 19 is similar to fig. 6, except for a finite difference grid of horizontal wells 300. The well 300 is located in the center of a central cell in each of a series of subsurface formation segments 304 extending in the y-direction along a vertical plane or extending in the z-direction. Fig. 24A and 24B show symbols of grid cells in the horizontal well model H of fig. 19. As shown in FIG. 19, well 300 passes N in the horizontal or y-direction_yIndividual cells were completed and potentials in adjacent cells: phi B_Up，ΦB_E，ΦB_Down，ΦB_WIs constant and known from the simulation run (previous time step or iteration value). Here, subscript B refers to "boundary", Up indicates the upper adjacent cell, E indicates the east adjacent cell, down indicates the lower adjacent cell, and W indicates the west adjacent cell, where, again, Φ describes the fluid potential or the pressure of the reference plane correction. In FIG. 19, well 300 extends horizontally in the y-direction along the longitudinal axis through each of a series of grid blocks 302. If the total well rate q is known_TIs input into the reservoir simulator, the reservoir simulator calculates the potential value Φ for each grid block 302 shown in fig. 19 for each time step.

As will be stated, at a rate q forming the sum up to a known sum_TThe perforation rate q of the layer(s)_iIn terms of measurements, the present invention improves the computer performance of the reservoir simulator.

Similar to equation (16), the puncturing rate q for each lattice block_iCan be expressed as:

q_i＝PI_i(Φ_i-Φ_w) (29)

in equation (29), the productivity index is governed by the following equation:

the above expressions of reservoir parameters and their physical relationships are in accordance with Chen and Zhang, "Well Flow Models for variance Numerical Methods". International Journal of Numerical Analysis and modeling, Vol 6, No 3, pages 378-388.

Fig. 21 and 22 collectively represent a flow chart G illustrating the basic computer processing sequence according to the present invention and the computational method that occurs during an exemplary embodiment of a continuous, fully implicit well model utilizing the reservoir simulation of the present invention, such as the horizontal well model shown in fig. 19. During step 400, the simulation begins by initializing the reservoir model H in the data processing system D and reading reservoir data and production data. The reservoir data and production data read during step 400 are of the type described above with respect to the vertical well model of flowchart F and step 200 (fig. 11).

During step 400, the cell potentials Φ for the k cells of the horizontal well model H are formed_kIs estimated. As schematically indicated in fig. 19, k ═ nx ny × nz. During step 402, an estimate of borehole potential is formed according to equation (32) set forth below:

wherein Φ in equation (32) is calculated from the initial potential distribution determined in step 400_i。

During step 404, the puncture rate q for each cell I is determined according to the relationship expressed in equation (29) above_i. In step (b)In step 406, the reservoir simulator is initialized and the simulator iteration count v is set to zero. For the initial iteration, the simulation time step t is also initialized to zero. The iteration count and the time step count are then incremented during subsequent iterations before step 406 is performed, as will be set forth below. The grid block potential for the initial time step is determined during step 406 by solving the three-dimensional potential equation with the reservoir simulator according to equation (1) for the single-phase petroleum stream. During step 408, each perforation Φ is determined for each perforation (i) (i ═ 1,2, Ny)_BW，Φ_B，Up，Φ_BE，Φ_B，DownThe surrounding boundary potentials and the boundary potentials are stored in a memory of the data processing system D.

Step 410 is then performed to form a simplified well model H-1 (FIG. 20B) by bringing together the grid blocks 302 with no flow obstructions between them (as shown at 306 in FIG. 20A). To perform step 410, a horizontal flow obstruction is identified as shown at 303 in the original horizontal well model system H. This can be done based on well log data or as specified by the reservoir analyst from the data in the original reservoir model. The simplified well model system model formed during step 410 would result in those layers 302 in the well model H that are between the flow barrier gridblocks or layers 303 and have horizontal flow between them, which layers 302 are combined together for analytical model purposes.

As can be seen in the example of fig. 20A and 20B, the horizontal well model H of fig. 20A with two sets of four grid blocks 302 on opposite sides of the flow obstacle grid block 303 is transformed into a simplified well model H-1 with two grid blocks 308 on each side of the flow obstacle grid block 303 as a result of step 410.

Thus the simplified well model equation for the simplified well model system H-1 in FIG. 20B becomes as shown in FIG. 20C:

wherein,

then, during step 410 (FIG. 22), the simplified well model is solved, resulting in And phi_wWhereinthe grid block potential of the simplified horizontal well model H-1 is shown.

During step 412,. phi._wThe horizontal well model is converted to a fixed flowing bottom hole potential well. Step 414 involves solving the three diagonal matrix system of FIG. 20C of the simplified horizontal well model H-1 shown in FIG. 20B to determine the potential of each grid 302 shown in FIG. 19.

Referring again to fig. 22, during step 416, a convergence step is performed to determine whether the non-linear iteration has converged. The individual residuals of the equations resulting from step 414 are checked against a user specified tolerance. If these tolerances are met, the non-linear iterative loop is exited, the solution output is written to the file at the current time step and stored in the memory of data processing system D, and processing returns to step 418 for the time step to be increased, and processing for the increased time step continues to step 406, as shown. If the user-specified tolerance is determined not to be met during step 416, processing according to the non-linear iterative loop returns to step 406 and continues. If the number of non-linear iterations becomes too large, it may be decided to adjust the model.

Multiple vertical wells

The 3-dimensional reservoir model M is shown in fig. 23 with multiple vertical wells 500, assuming each has a single phase, slightly compressible petroleum stream. The tissue model M is named according to the symbol schematically illustrated in fig. 25A, 25B, 25C and 25D.

The model M of fig. 23 is similar to fig. 6, except for a finite difference grid of several vertical wells 500. The wells 500 are each located in the center of a central cell in each of a series of subsurface formation segments 504 extending horizontally in the z-direction or extending in the y-direction and z-direction.

As shown in fig. 6, each well 500 is completed in the vertical or z-direction by Nz cells, and the potentials in adjacent cells: phi B_N，ΦB_E，ΦB_S，ΦB_WIs constant and known from the simulation run (previous time step or iteration value). As shown in fig. 6, here, subscript B refers to "boundary", N to the upper adjacent cell, E to the east adjacent cell, down to the lower adjacent cell, and W to the west adjacent cell, where again Φ describes the fluid potential or the reference plane corrected pressure.

Each well 500 extends horizontally in the z-direction along the longitudinal axis through each of a series of grid blocks 502, and each well 500 is segmented by completions 506 at specific formations at various depths. Fig. 25A to 25D schematically show the number symbols designated as well 1 and well 2 in the model M of the well 50 in fig. 23.

If the known total well rate qT for each well 500 is provided as an input parameter into the reservoir simulator, the reservoir simulator calculates the potential value Φ for each grid block 502 shown in FIG. 23 for each time step.

For the general case of a 3-dimensional vertical reservoir model M, n_wThe number of wells is represented and the total well rate per well 500 is represented by q_T(l) Given and expressed, l ═ 1, n_wIs known from production data and is provided as an input parameter into the reservoir simulator.

FIG. 26 represents a flow chart of a process according to the invention in which multiple vertical wells are present, as in the case of model M of FIG. 23. Thus, for 3 dimensionsThe reservoir model M is vertical, and in step 602, the reservoir simulator is initialized and the reservoir and production are read from memory for processing. This is done in a manner comparable to step 200 of fig. 11. The simulation process continues in 602 shown in fig. 26, where an estimate of the wellbore potential for each well, Φ w (l), l 1, 2.. nw, is formed according to equation (32), where the production rate index PI is determined according to the measurements expressed in equation (17). During step 604, based on the estimates of the borehole potential from step 602, a perforation rate q is calculated for each well_i(l)，l＝1，...nw。

Step 606 involves forming a simplified well model of the 3-dimensional multiple well model M for each well l ═ 1, 2.. nw. Reservoir simulation is performed by combining or merging adjacent well cells 502 in a formation, the adjacent well cells 502 having fluid communication with each other and also being located between flow barriers with no flow therebetween. This is done in three-dimensional layers, such as layers adjacent to a single vertical low well as schematically shown in fig. 1A and 1B and layers adjacent to horizontal flow wells as schematically shown in fig. 20A and 20B.

In step 608, a simplified well matrix for each well is formed based on the simplified well model of step 606 in a manner comparable to the simplified well matrix of fig. 20C of the horizontal well model H of fig. 19.

Then, in step 610, the bottom hole potential Φ w (l), l 1, 2.. nw for each well is determined by solving the simplified well matrix obtained in step 608. In step 612, the diagonal of the primary matrix for each well is modified according to the relationships involving equations 18,18A, and 18B aboveAnd the right end item b_iAn item. Fig. 27 is an example schematic diagram of a simplified well matrix of a linear system of equations for a simplified two-well, three-dimensional reservoir model or 3x3x2 reservoir model using the numbering system explained in fig. 24A and 24B and fig. 25A through 25D.

During step 614, unknowns for model MAll of the grid blocks 502 of the number form an overall matrix as in the example of fig. 27. Processing after step 614 of FIG. 26 then continues to a convergence test in the manner of step 416 of FIG. 27, with iterations and time steps increasing as one returns to step 602 for further processing including n_wIteration of the total system for each of the individual wells.

Data processing system

As shown in FIG. 16, a data processing system D in accordance with the present invention includes a computer 240 having a processor 242 and a memory 244 coupled to the processor 242 for storing operating instructions, control information, and database records therein. If desired, the computer 240 may be a portable digital processor, such as a personal computer of the form: a laptop computer, a notebook computer or other suitably programmed or programmable digital data processing apparatus, such as a desktop computer. It should also be appreciated that the computer 240 may be a multi-core processor with nodes (e.g., from Intel corporation or advanced microelectronic device corporation (AMD)) or any conventional type of mainframe computer with appropriate processing capabilities (e.g., international business machines corporation (IBM) available from Armonk, n.y.) or other source.

The computer 240 has a user interface 246 and an output display 248 that are used to display the output data or records of the processing of log data measurements performed in accordance with the present invention to obtain measurements and to form a model of the determined well production of formations in the well(s) of the subsurface formations. The output display 48 includes components, such as a printer and an output display screen, which are capable of providing printed output information or visual output in the form of output records or images: graphics, data tables, graphical images, data graphs, and the like.

The user interface 246 of the computer 240 also includes appropriate user input devices or input/output control units 250 to provide user access to control or access information and database records and to operate the computer 240. Data processing system D also includes a database 252 stored in computer memory, which may be internal memory 244 or external networked or non-networked memory, as indicated at 254 in an associated database server 256.

Data processing system D includes program code 260 stored in non-transitory memory 244 of computer 240. The program code 260 according to the present invention is in the form of computer operable instructions that cause the data processor 242 to form a continuous, fully-implicit well model of a reservoir simulation according to the present invention in the manner already mentioned above.

It should be noted that program code 260 may be in the form of a microcode, program, routine, or symbolic computer operable language that provides a specific set of operations that control the functionality of data processing system D and direct the sequence of its operations. The instructions of program code 260 may be stored in memory 244 of computer 240 or on a computer diskette, magnetic tape, conventional hard disk drive, electronic read-only memory, optical storage device, or other suitable data storage device having a computer usable non-transitory medium stored thereon. Program code 260 may also be embodied on a data storage device (e.g., server 64) as a non-transitory computer readable medium, as shown.

Two illustrative example model equations are presented below: seven homogeneous reservoirs with one fracture flow barrier (fig. 7); and twenty-two homogeneous reservoirs with two fracture flow barriers (fig. 9).

Seven-layer uniform well model

Fig. 7A shows seven reservoirs and attributes of an original model 70 and a simplified model 71. As shown, assume that the reservoir has seven layers. The layer thickness of each layer 72 with vertical flow is 10 ft. There is a layer 73 with a thickness of 1ft indicating no break in the vertical flow. It is further assumed that layer 73 does not communicate with the upper and lower layers 72. Assume that there is a vertical well in the middle, as indicated by the arrow. For the model of fig. 7A, the initial reservoir potential (base-surface corrected pressure) was 3,000 psi. Assuming that each layer 72 has a planar permeability k of 10md_xAnd k_yAnd a vertical permeability k of 1md_z。

Table 1 summarizes the reservoir and grid properties of model 70. Assume that the grid size in the planar direction (square grid) is 840 ft. The oil viscosity is set to 1cp and the oil formation volume factor is assumed to be 1. The total oil production rate of the well was set to 1,000B/D. The zone productivity index PI for each zone completion was calculated by the picosmann method, as described, and is also shown in table 1.

Table 1: equation 1-reservoir Properties

Full-implicit full-coupling simultaneous solution

The coefficient matrices for the solutions to reservoir pressure and bottom hole pressure are formed in a similar manner as described above with respect to equations (18-19) and shown in FIG. 13. It can be appreciated that there are only 8 unknowns (7 potentials or reference plane corrected pressures and one bottom hole potential) and the coefficient matrix is non-sparse. For unknown reservoir (stratum) potential phi_iI ═ 1,7 and other unknowns Φ_WThe linear system of equations can be solved by a direct method (e.g., gaussian elimination).

Results

Table 2 summarizes the calculated interval potentials, wellbore potentials, and interval (completion) flow rates for model 70 of fig. 7A.

TABLE 2 exact solution of original equation

According to the calculation result, the calculated bottom hole potential phi is seen_W＝1257.36psi。

Formation of equations according to the invention

According to the reservoir data in table 1 and as shown in fig. 7A, there is only one layer 73 that is not vertically in communication with the other layers. Thus, as shown in FIG. 7B, in accordance with the present invention, the layers 72 above the fracture layer 73 are combined into a single layer to form a simplified well model. Similarly, the layers 72 below the layer 73 are combined into a single layer. It can now be seen that the simplified model has only three layers. The total number of unknowns is 4 compared to 8 in the full model.

Table 3 summarizes the attributes of the simplified well model 71 formed in accordance with the process of the present invention.

TABLE 3 simplified well model

Layer(s)	Thickness ft	K_x＝K_y,md	K_z,md	PI,b/d/psi
					1	20	10	1	0.24
2	1	10,000	0	12.14
					3	40	10	1	0.40

The linear system of equations (equation 20) of the simplified system also has an unstructured coefficient matrix, but with a number of unknowns below 50%. In a practical reservoir, the well model size reduction according to the invention will be strong in case of hundreds of layers and only some flow barriers, e.g. the simplified well model system model according to the invention may be one hundredth of the size of the whole system. The simplified system is solved for layer potential and bottom hole potential by a direct solver. The results are shown in Table 4.

TABLE 4 results of simplified System

It can be seen that the calculated bottom pore potential: phi_W1257.36psi with full model calculated phi_WAre identical.

The determined well potential is the only information needed for the next step in the development. The well is then processed into a prescribed bottom hole pressure (potential) model. Followed by computer processing according to the process described using the matrix of fig. 14 and equations (16 and 18) to calculate the flow profile (layer velocity) and the total well rate. In fig. 14, the upper diagonal solid line of the matrix represents T defined by equation (11)_Up,iAnd the lower diagonal solid line description of the matrix, also defined by equation (2), is called T_Down,iOf (2) is used. Central item T_C,iIs defined by equation (17a), and the right-end term bi is defined by equation (17 b).

The results are summarized in table 5. Note that the total calculated well rate is exactly the same as the input value of 1,000 b/d.

TABLE 5 results of Using the overall System of the invention

Layer(s)	Potential, psi	Rate/d
			1	2630.61	166.67
2	2630.61	166.67
			3	1257.37	0.0
4	2630.61	166.67
			5	2630.61	166.67
6	2630.61	166.67
			7	2630.61	166.67
Total of		1,000

The results presented in table 5 are the same as table 1 for the fully implicit well model. The difference or error between the calculated and input well rates for the well is zero in this case and no additional iterations are required. This is due to the fact that: the reservoir was homogeneous and no gross errors were made in forming the simplified system. The diagonal elements and right-hand terms of the matrix are the same as in fig. 14, i.e., the lower diagonal solid line represents T defined by equation (11)_up,iAnd the upper diagonal solid line of the matrix depicts what is referred to as T above_Down,iOf (2) is used. Central item T_C,iIs defined by equation (17a), and the right-end term bi is defined by equation (17 b).

The term PI appearing on equation (17) is the perforation productivity index of the square grid, defined by the following formula:

wherein r is_wIs the borehole radius.

Comparison with explicit well models

In several reservoir simulators, either a semi-implicit well model or an explicit well model is used. If the formulation of the well model is semi-implicit but it collapses to explicit under the pressure variable, then the formulation collapses to the explicit well model. The explicit well model is directed to equations (12-14) obtained by a computer process that follows the matrix of FIG. 12. In FIG. 12, the term T appears on the diagonal element_Down,i,T_Up,iIs defined by equation (11), and T_c,i,b_iDefined by equation (14a) and equation (14 b). Fig. 8A shows seven reservoirs and attributes of the implicit well model 80, and fig. 8B shows an explicit well model 81 of the same structure as the models of fig. 7A and 7B. As shown, assume that the reservoir has seven layers. The layers 82 each have a potential Φ of 2630 psi. Layer 83 has a potential Φ of 1257psi, layer 83 indicating a fracture and further assumed not to communicate with the layers 72 above and below. Assume that there is a vertical well in the middle, as indicated by the arrow. Fig. 8A and 8B compare the results of the implicit and explicit models. The calculated perforation (layer) rates are summarized in table 6.

TABLE 6 comparison of perforation (interval) rates for different well models

As shown, the model 81 according to the explicit method is inaccurate; which miscalculates the puncture rate completely. The explicit modeling approach actually distributes all well production from the thin fractured layer 83 as indicated in fig. 8B because that layer has the highest productivity index.

Implicit methods (whether computationally intensive fully implicit models or simplified models according to the present invention) do not make this assignment, but instead determine that layer 83 is not supported by fluid from layers 82 above and below. A fracture layer of the type shown at 83 present in a real reservoir can get the only fluid to support the adjacent cells from its plane. However, since the fracture layer is a very thin layer, the water conductivity in these directions is by its nature small. Therefore, the fracture layer cannot support the fluid at the rate simulated by the explicit model.

In fact, conventional implicit well models show that during transient times, the fracture layer supports most well production, as do explicit methods. However, the formation pressure in layer 83 drops rapidly and assumes a uniform value of wellbore potential (constant bottom hole pressure). After the pressure drops, it reaches steady state and the well production rate is actually produced by contributions from layers 82 above and below the vertical flow barrier 83.

Twenty-two layer uniform reservoir model

Model mesh system 90 includes twenty-two layers as shown in fig. 9. The position of the high permeability fracture layers 6 and 12 is seen as moving down through the layers and is indicated schematically at 91 and 92. Above layer 91 there are 5 layers 93, numbered 1 to 5, each with vertical flow. Also in the model 90 there are five layers 94 with vertical flow between the flow barrier layers 91 and 92 and ten layers 95 with vertical flow below the flow barrier layer 92. Reservoir data for model 90 is shown in table 7.

TABLE 7-reservoir data for 22-layer equation

Permeability of plane adjacent unit 20mD

Total well productivity of 2,500B/D

Wells completed in all zones.

As a result:

full-implicit full-coupling simultaneous solution

Calculated bottom hole potential Φ_W＝1421.247psi

TABLE 8 potential distribution, psi

Layer(s)	Φ_W	Φ_i
			1	1421.25	2901.32
2	1421.25	2900.82
			3	1421.25	2900.38
4	1421.25	2900.38
			5	1421.25	2900.38
6	1421.25	1421.25
			7	1421.25	2864.43
8	1421.25	2864.37
			9	1421.25	2864.10
10	1421.25	2863.88
			11	1421.25	2864.03
12	1421.25	1421.25
			13	1421.25	2841.62
14	1421.25	2841.58
			15	1421.25	2841.48
16	1421.25	2841.33
			17	1421.25	2841.13
18	1421.25	2840.65
			19	1421.25	2840.08
20	1421.25	2839.65
			21	1421.25	2839.37
22	1421.25	2839.24

The model 90 is then subjected to explicit modeling methodology techniques of the type described above and the determined flow data. A comparison of the fully implicitly and explicitly processed flow rate profiles of the reservoir model 90 using the previously described techniques is presented in table 9.

TABLE 9 comparison of flow rates for fully implicit and explicit well methods

Layer(s)	Hidden type	Display type
			1	35.93	14.56
2	89.78	36.39
			3	53.85	21.83
4	179.48	72.78
			5	89.74	36.39
6	0.00	727.80
			7	105.09	43.67
8	52.54	21.83
			9	157.60	65.50
10	210.10	87.34
			11	87.55	36.39
12	0.00	727.80
			13	129.29	54.59
14	129.28	54.59
			15	129.28	54.59
16	129.26	54.59
			17	129.24	54.59
18	158.49	66.96
			19	158.42	66.96
20	158.37	66.96
			21	158.34	66.96
22	158.33	66.96

Simplified model construction

Since there are only two vertical flow barrier layers 91 and 92, the layer 93 above the layer 91 in fig. 9 can be combined into one layer; the layers 94 below the layer 91 are combined into one layer and the layers 95 below the layer 92 are combined into another single layer. The total number of layers according to the invention is therefore 5. The attributes of the simplified model are as follows:

TABLE 10 attributes of simplified well model

Calculated bottom hole potential

Φ_W＝1421.34psi

The simplified model according to the invention continues to determine the bottom hole potential Φ_w. The results for the simplified model with five layers are listed in table 11.

TABLE 11 potential distribution

Layer(s)	Pot wf	Pot
			1	1421.34	2900.53
2	1421.34	1421.35
			3	1421.34	2864.16
4	1421.34	1421.35
			5	1421.34	2740.61

By using the bottom hole potential phi calculated from the simplified model_wAnd a prescribed bottom hole potential phi using a full model_wPotential is calculated and completion interval rate is calculated according to equation (17). The results are indicated below in table 12.

TABLE 12 calculated interval rates

Newly calculated q_t＝2499.84b/d

Error 2,500-2499.8488

＝0.15b/d

The error in the total rate and the calculated bottom hole pressure vanish with the non-linear newton iteration of the simulator. The present invention achieves a simplified model with a productivity of acceptable accuracy, but where the model complexity and computer processing time are significantly reduced, compared to the results obtained by prior art fully implicit fully coupled processing techniques.

The present invention as already described above does not require a dedicated linear solver for the solution of coupled reservoir and well equations. In contrast, the coefficient matrices for the previously used coupled reservoir and well equations do not have a regular sparse structure. Thus, conventional types of coupled reservoir and well equations require specialized solvers that are expensive and can face convergence equations.

It can be seen that the present invention, as described above, does not require any special solvers of the solution of the coupled reservoir and well equations. The same solver for reservoir equations is utilized. Only the modification to the coefficient matrix is in the diagonal entries.

The present invention solves reservoir simulation equations in which a vertical well has many completions (which are often found in a reservoir). In recent simulation studies, wells with more than 100 vertical zones (completions) are common. Fully-coupled fully-implicit well models with simultaneous solutions are very expensive in these cases. The invention can save a great deal of computer time.

The present invention is very useful for wells having hundreds of perforations completed in highly uniform reservoirs. By identifying and advantageously utilizing the physical principles involved, the present invention reduces the large, time-consuming equations of well modeling simulations in a reservoir with a large number of formation (completion) equations to small equations. In the context of the present invention, it has been found that vertically communicating layers can be synthesized as a single layer. The simplified model so formed retains the same bottom hole pressure as the original full model. Once the simplified model is solved for the bottom hole pressure, the wells in the large system are then treated to the specific bottom hole pressure and are easily solved by a conventional linear solver. Thus, the present invention eliminates the need to write and acquire non-structural linear solvers for many wells with hundreds of completions and which can be expensive.

The present invention has been described sufficiently that a person having average knowledge of the subject can reproduce and obtain the results mentioned in the present invention. However, any person skilled in the art of the subject matter of the present invention may carry out modifications not described in the present claims to apply these modifications to a determined structure or to require, during the manufacturing thereof, the subject matter claimed in the appended claims; such structures are intended to be included within the scope of the present invention.

It should be noted and understood that improvements and modifications to the invention described in detail above may be made without departing from the spirit or scope of the invention as set forth in the appended claims.

Claims

1. A computer-implemented method of forming a model of well production rates of constituent fluids from a plurality of wells in a subterranean reservoir determined from measured total well production and determined zone completion rates of perforated intervals in said wells during reservoir simulation of well production using a coupled well reservoir model organized into a reservoir grid subdivided into a plurality of reservoir cells at said perforated intervals during the life of the subterranean reservoir, said perforated intervals in said reservoir being located at a plurality of formations having unknown well potentials and fluid completion rates of constituent fluids at said time steps and said formations comprising a vertical fluid flow layer from which fluids flow vertically and a flow barrier layer from which no fluids flow vertically, said formations further having a permeability, a pressure, and a flow rate, Thickness and stratigraphic potential, the coupled well reservoir model further having a plurality of well units at locations of the wells in formations of the reservoir, the computer-implemented method determining a stratigraphic completion rate of the constituent fluids from the formations of the wells and a well productivity rate of the constituent fluids from the wells, the computer-implemented method comprising the steps of:

(a) forming a fully-computational matrix reservoir model of reservoir data for cells of the model, the fully-computational matrix reservoir model including the reservoir data for the reservoir cells at the perforated interval, the reservoir data including the permeability, thickness, and potential of the formation;

(b) forming a simplified well model system matrix by combining data for vertical fluid flow layers in the reservoir model having vertical fluid flow therebetween and located between flow barriers into a single vertical flow layer in the matrix;

(c) determining a bottom hole pressure of the well;

(d) forming a coupled reservoir well model comprising the full computational matrix reservoir model and the simplified well model system matrix, the wells being treated as bottom hole, pressure designated wells having a determined bottom hole pressure, the coupled reservoir well model being in the form of a matrix of:

wherein A is_RRIs a three diagonal matrix of said reservoir data, A_RWA diagonal matrix that is a productivity index of the formation adjacent to the perforated interval; a. the_WRA diagonal matrix that is a productivity index of the formation from a well to the reservoir; a. the_WWIs an index of the productivity of the wellA matrix;a matrix of unknown reservoir potentials of the cells surrounding the well;is a matrix of unknown well potentials in the wellbore of the well;a matrix of reservoir data constants for the reservoir cells around the well; andis a matrix of the well data constants for the well;

(e) solving the coupled reservoir well model for the fluid flow in the reservoir cells of the formation and the production rates and potentials of the reservoir cells of each of the formation at the time step;

(f) solving the coupled reservoir well model to obtain a productivity index of the well units at the perforated interval of the reservoir at the time step;

(g) determining formation completion rates for the component fluids in the vertical fluid flow and flow barrier layers of the well based on the determined productivity index for the reservoir cells at the perforated interval of the reservoir and the determined productivity index for well cells at the time step;

(h) determining a total well productivity for the well from the determined interval completion rates for the component fluids in the vertical fluid flow zone and the flow barrier zone of the well in the time step; and

(i) forming a record of the determined interval completion rates of the constituent fluids in the vertical fluid flow zone and the flow barrier of the well and the determined total well productivity for the well in the time step.

2. The computer-implemented method of claim 1, wherein the tri-diagonal matrix of the reservoir cells comprises diagonals representing the water conductivity of the cells of the reservoir.

3. The computer-implemented method of claim 1, wherein solving the coupled well reservoir model comprises applying a full matrix solver.

4. The computer-implemented method of claim 3, further comprising the steps of:

reducing residuals from the step of applying the full matrix solver to a specified tolerance; and

if the tolerance is met at the time step and the tolerance is met for each time step during the life of the subsurface reservoir, proceeding to the step of forming a record; and

if the tolerance is not met at the time-step, then for another iteration of processing at the time-step, return to step (a) and repeat steps (b) through (h).

5. The computer-implemented method of claim 3, further comprising the steps of:

increasing the simulator time step if the specified tolerance is met at the time step but not for each time step during the life of the reservoir; and

returning to step (a) and repeating steps (b) through (h) for another iteration of the process with an increased simulator time step.

6. The computer-implemented method of claim 1, wherein solving the coupled well reservoir model comprises: solving for fluid flow based on permeability, thickness and stratigraphic potential of the formation.

7. The computer-implemented method of claim 1, wherein the plurality of wells comprises a plurality of vertical wells.

8. A data processing system for forming a model of the well productivity of constituent fluids from a plurality of wells in a subterranean reservoir determined from measured total well production and the zone completion rates of perforated intervals in said wells determined during reservoir simulation of well production using a coupled well reservoir model organized into a reservoir grid subdivided into a plurality of reservoir cells at said perforated intervals during the life of the subterranean reservoir, said perforated intervals in said reservoir being located at a plurality of formations having unknown well potentials and fluid completion rates of constituent fluids at said time steps and comprising vertical fluid flow layers from which fluids flow vertically and flow barrier layers from which no fluids flow vertically, said formations further having a permeability, a method for producing a model of a subterranean reservoir from a coupled well reservoir model, and a method for producing a model of a subterranean reservoir from a coupled well reservoir model, Thickness and stratigraphic potential, the coupled well reservoir model further having a plurality of well units at locations of the wells in formations of the reservoir, the data processing system determining a stratigraphic completion rate of the constituent fluids of the formations from the wells and a well production rate of the constituent fluids from the wells, and the data processing system comprising:

a processor that performs the steps of:

(b) forming a simplified well model system matrix by combining the data for vertical fluid flow layers in the reservoir model having vertical fluid flow therebetween and located between flow barriers into a single vertical flow layer in the matrix;

(c) determining a bottom hole pressure of the well;

(d) forming a coupled reservoir well model comprising the full computational matrix reservoir model and the simplified well model system matrix, the well being treated as a bottom hole, pressure-designated well having a determined bottom hole pressure, the coupled reservoir well model being in the form of a matrix:

wherein A is_RRIs a three diagonal matrix of said reservoir data, A_RWA diagonal matrix that is a productivity index of the formation adjacent to the perforated interval; a. the_WRA diagonal matrix that is a productivity index of the formation from a well to the reservoir; a. the_WWIs a matrix of productivity indices for the well;a matrix of unknown reservoir potentials of the cells surrounding the well;is a matrix of unknown well potentials in the wellbore of the well;a matrix of reservoir data constants for the reservoir cells around the well; andis a matrix of the well data constants for the well;

(f) solving the coupled reservoir well model for the productivity index of the well units at the perforated interval of the reservoir at the time step;

(g) determining the zone completion rates of the constituent fluids of the vertical fluid flow zone and the flow barrier of the well based on the determined productivity indices of the reservoir cells and well cells at the perforated interval of the reservoir at the time step;

(h) determining a total well productivity for the well from the determined zone completion rates for the vertical fluid flow zone and the constituent fluids of the flow barrier for the well in the time step; and

a memory that performs the steps of: forming a record of the determined interval completion rates of the component fluids of the vertical fluid flow zone and the flow barrier of the well and the determined total well productivity of the well in the time step.

9. The data processing system of claim 8, wherein the tri-diagonal matrix of the reservoir cells comprises diagonals representing water conductivity of the cells of the reservoir.

10. The data processing system of claim 8, wherein the processor in solving the coupled well reservoir model performs the step of applying a full matrix solver.

11. The data processing system of claim 10, further comprising the processor performing the steps of:

if the tolerance is met at the time step and the tolerance is met for each time step during the life of the subsurface reservoir, causing the memory to perform the step of forming a record; and

if the tolerance is not met at the time step, returning to step (a) and repeating steps (b) through (h) for another iteration of the process at the time step.

12. The data processing system of claim 10, further comprising the processor performing the steps of:

returning to step (a) and repeating steps (b) through (h) for another iteration of the process of increasing simulator time step.

13. The data processing system of claim 8, wherein the processor in solving the coupled well reservoir model performs the step of solving for fluid flow based on permeability, thickness, and stratigraphic potential of the formation.

14. The data processing system of claim 8, wherein the plurality of wells comprises a plurality of vertical wells.

15. A data storage device storing computer operable instructions in a non-transitory computer readable medium, the instructions causing a data processor to form a model of well production rates of constituent fluids from a plurality of wells in a subsurface reservoir determined from measured total well production and determined zone completion rates of perforated well sections in the wells at time steps during reservoir simulation of well production using a coupled well reservoir model organized into a reservoir grid subdivided into a plurality of reservoir cells at the perforated well sections, the perforated well sections in the reservoir being located at a plurality of formations having unknown well potentials and fluid completion rates of constituent fluids at the time steps and the formations comprising a fluid flow zone in which fluids flow vertically and a flow barrier zone in which no fluids flow vertically, the formation further having a permeability, a thickness, and a stratigraphic potential, the coupled well reservoir model further having a plurality of well units at locations of the wells in the formations of the reservoir, the stored computer operable instructions causing the data processor to determine a stratigraphic completion rate of the constituent fluids of the formations from the wells and a well productivity rate of the constituent fluids from the wells by performing the steps of:

(c) determining a bottom hole pressure of the well;

(g) determining formation completion rates of the component fluids of the vertical fluid flow zone and the flow barrier of the well based on the determined productivity indices of the reservoir cells and well cells at the perforated interval of the reservoir at the time step;

(h) determining a total well productivity for the well from the determined interval completion rates for the vertical fluid flow zones and the constituent fluids of the flow barriers for the well in the time step; and

(i) forming a record of the determined interval completion rates of the component fluids of the vertical fluid flow zone and the flow barrier of the well and the determined total well productivity of the well in the time step.

16. The data storage device of claim 15, wherein the instructions for causing the processor to perform the step of solving the coupled well reservoir model comprise instructions for applying a full matrix solver.

17. The data storage device of claim 15, wherein the instructions further comprise instructions that cause the processor to:

18. The data storage device of claim 15, wherein the instructions further comprise instructions that cause the processor to:

19. The data storage device of claim 15, wherein the instructions for causing the processor to perform the step of solving the coupled well reservoir model comprise: instructions for resolving fluid flow based on permeability, thickness, and stratigraphic potential of the formation.

20. The data storage device of claim 15, wherein the plurality of wells comprises a plurality of vertical wells.

21. A computer-implemented method of forming a model of well productivity of constituent fluids from horizontal wells in a subterranean reservoir determined from measured total well production and determined completions of layers of perforated intervals in said horizontal wells during reservoir simulation of well production using a coupled well reservoir model organized into a reservoir grid subdivided into a plurality of reservoir cells at said perforated intervals during the life of the subterranean reservoir, said perforated intervals in said reservoir being located at a plurality of formations having unknown well potentials and fluid completions of constituent fluids at said time steps and comprising horizontal fluid flow layers with fluids flowing horizontally therefrom and flow barriers with no fluids flowing horizontally therefrom, said formations further having permeabilities, permeability, and completion rates of wells in a reservoir simulation of well production using a coupled well reservoir model, Thickness and stratigraphic potential, said coupled-well reservoir model further having a plurality of well units at locations of said horizontal wells in a formation of said reservoir, said computer-implemented method determining a stratigraphic completion rate of said constituent fluids from said formation of said horizontal wells and a well productivity rate of said constituent fluids from said horizontal wells, said computer-implemented method comprising the steps of:

(b) forming a simplified well model system matrix by combining the data for horizontal fluid flow layers in the reservoir model having horizontal fluid flow therebetween and located between flow barriers into a single horizontal flow layer in the matrix;

(c) determining the bottom hole pressure of the horizontal well;

(d) forming a coupled reservoir well model comprising the fully computational matrix reservoir model and the simplified well model system matrix, treating the horizontal well as a bottom hole with a determined bottom hole pressure, a pressure-designated well, the coupled reservoir well model in the form of a matrix:

wherein A is_RRIs a three diagonal matrix of said reservoir data, A_RWIs a vector of productivity indices of the formation adjacent to the perforated interval; a. the_WRA vector that is a productivity index of the formation from a well to the reservoir; a. the_WWA linear scalar that is a productivity index of the horizontal well;is a vector of unknown reservoir potentials of the cells surrounding the horizontal well;is a vector of unknown well potential in the wellbore;is a vector of reservoir data constants for the reservoir units surrounding the horizontal well; andis a vector of the well data constants for the horizontal well;

(g) determining formation completion rates of said constituent fluids of said horizontal fluid flow zones and said flow barriers of said horizontal wells based on said determined productivity indices of said reservoir cells and well cells at said perforated interval of said reservoir at said time step;

(h) determining a total well productivity rate for the horizontal well from the determined interval completion rates for the constituent fluids of the horizontal fluid flow interval and the flow barrier interval of the horizontal well at the time step; and

(i) forming a record of the determined interval completion rates of the constituent fluids of the horizontal fluid flow zone and the flow barriers of the horizontal well and the determined total well productivity of the horizontal well at the time step.

22. A modeled data processing system, the model being a model of well production rates of constituent fluids from horizontal wells in a subterranean reservoir determined from measured total well production and determined completions rates of layers of perforated well sections in said horizontal wells at time steps during reservoir simulation of well production using a coupled well reservoir model organized into a reservoir grid subdivided into a plurality of reservoir cells at said perforated well sections, said perforated well sections in said reservoir being located at a plurality of formations having unknown well potentials and fluid completions rates of constituent fluids at said time steps, and said formations comprising horizontal fluid flow layers from which fluids flow horizontally and flow barriers through which no fluids flow horizontally, said formations further having permeabilities, Thickness and stratigraphic potential, said coupled-well reservoir model further having a plurality of well units at locations of said horizontal wells in a formation of said reservoir, said data processing system determining a stratigraphic completion rate of said constituent fluids from said formation of said horizontal wells and a well production rate of said constituent fluids from said horizontal wells, and said data processing system comprising:

a processor that performs the steps of:

(c) determining the bottom hole pressure of the horizontal well;

(h) determining a total well productivity rate for the horizontal well from the determined interval completion rates for the constituent fluids of the horizontal fluid flow zone and the flow barrier zone of the horizontal well at the time step; and

a memory that forms a record of the determined interval completion rates of the constituent fluids of the horizontal fluid flow zone and the flow barriers of the horizontal well and the determined total well productivity of the horizontal well at the time step.

23. A data storage device storing computer operable instructions in a non-transitory computer readable medium that cause a data processor to form a model of well productivity of constituent fluids from horizontal wells in a subsurface reservoir determined from measured total well production at a time step and determined zone completion rates of perforated well sections in said horizontal wells during reservoir simulation of well production using a coupled well reservoir model organized into a reservoir grid subdivided into a plurality of reservoir cells at said perforated well sections during the life of the subsurface reservoir, said perforated well sections in said reservoir being located at a plurality of formations having unknown well potentials and fluid completion rates of constituent fluids at said time step and said formations comprising horizontal fluid flow layers from which fluids flow horizontally and flow barrier layers from which no fluids flow horizontally Said formation further having a permeability, a thickness and a stratigraphic potential, said coupled-well reservoir model further having a plurality of well units at the location of said horizontal well in the formation of said reservoir, said stored computer operable instructions causing said data processor to determine a stratigraphic completion rate of said constituent fluid from said formation of said horizontal well and a well productivity rate of said constituent fluid from said horizontal well by performing the steps of:

(c) determining the bottom hole pressure of the horizontal well;