PARALLEL PRODUCTION SIMULATION OF A MATURE HYDROCARBON FIELD

Abstract

Embodiments for parallel production simulation of an oil field are disclosed. A plurality of production parameter data sets corresponding respectively to a plurality of scenarios are received and aligned in a memory of a vectorized computation hardware. In SIMT-based architecture, each well of each scenario is associated with a thread of a plurality of threads. In a SIMD-based architecture, each well of each scenario is associated with a respective lane of a multi-lane data pipeline.

Description

BACKGROUND

1. Technical Field

The technical domain of the present disclosure is the exploitation of hydrocarbon fields. More particularly the present disclosure concerns a method of building a reliable simulator able to forecast quantities produced vs. production parameters, in the particular case of mature fields.

2. Background Art

Mature hydrocarbon fields represent a special challenge both in terms of investment and allocation of human resources, because the net present value of any new investment diminishes with the degree of maturity. Therefore, less and less time and effort can be invested in reservoir studies to support field exploitation. Still, there remain opportunities to improve the production over a so-called “baseline” or “business as usual” behavior of an entire mature field, even with little investment. Past strategic choices in the way to operate the hydrocarbon field have created some heterogeneity in pressure and saturation. These can be drastically revisited and production parameters reshuffled accordingly. With respect to a mature hydrocarbon field, many production avenues have been explored in the past, and a learning process can be applied: reshuffled parameters can be implemented with low risk.

Two prior art approaches are currently known to model the behavior of a hydrocarbon field and to forecast an expected quantity produced in response to a given set of applied production parameters.

A first approach, called “meshed model” or “finite element modeling” parts a reservoir into more than 100,000s of elements (cells, flowlines . . . ), each cell carrying several parameters (permeability, porosity, initial saturation . . . ), and applies physical laws over each of said cells in order to model the behavior of fluids in the hydrocarbon field. In that case the so-called Vapnik-Chervonenkis, VC-dimension h of the space of solutions S, from which the simulator is selected, is very large. Therefore, the available number m of measured data in history data remains comparatively small, even for mature fields, and the ratio

$\frac{h}{m}$

appears to be very large compared to 1. As a result of the Vapnik learning theory, which is further mentioned later, the forecast expected risk R is not properly bounded (due to the Φ term), and such a simulator cannot be considered to be reliable, even if it presents a very good match with history data. In practice, it is widely recognized that for such meshed models, a good history match does not guarantee a good forecast: there are billions of ways to match the past, leaving large uncertainty on which one provides a good forecast.

A second approach, by contrast, uses over-simplified models, such as, for instance, decline curves or material balance. However, this is too simplified to properly take into account the relevant physics and geology of the reservoir, in particular complex interaction and phenomena. In such case, the forecast expected risk R is not minimized, because no good match can be reached (the empirical Risk R_emp term remains large).

SUMMARY

In summary, the present disclosure represents a proper compromise between too complex and too simplistic modeling approaches. It is valid only for mature fields, which deliver enough past information, in the form of history data HD, to allow building a space S of candidates for becoming a field production simulator, large enough to take into account all key phenomena at stake in the field, without becoming too complex and hence requiring too many history data to be calibrated.

The Vapnik statistical learning theory defines under which conditions such a simulator can be devised. Such a simulator may be designed in such a way it complies with the conditions for a reliable forecasting capacity.

The object of the present disclosure is a production simulator, for simulating a mature hydrocarbon field, providing quantity produced per phase, per well, per layer (or group of layers) and per time as a function of production parameters, wherein said production simulator matches sufficiently well history data of said mature hydrocarbon field and verifies a Vapnik condition. This said Vapnik condition ensures that the quantities calculated to be produced per phase, per well, per layer (or group of layers) and per time to be accurate. This will allow the users of this simulator to play different production scenarios, according to different production parameters, each of them delivering reliable quantities, so that such scenarios can be properly compared to each other and an appropriate one can be selected according to specific criteria. As all production scenarios will deliver reliable quantities, the selected one will also deliver a reliable forecast of the production, and therefore it will become a low risk and preferred path for producing the field.

According to another feature of the present disclosure, a match with history data is obtained when:

$\frac{{Q_{\varphi \mathrm{ktb}}-Q_{\varphi \mathrm{ktbHD}}}_{\left[T_1-\mathrm{Xy},T_1\right]}}{{Q_{\varphi \mathrm{ktbHD}}}_{\left[T_1-\mathrm{Xy},T_1\right]}}\le {\varepsilon }_1,$

for more than p % of existing Q_φktbHD values and the cumulated oil produced over the same [T₁−Xy, T₁] time interval is accurate up to ϵ₂, where

Q_φktb are quantities produced per phase φ∈{oil, water, gas}, per well k, per layer b (or group of layers) and per time t, determined by the production simulator,

Q_φktbHD are the same quantities produced per phase, per well, per layer (or group of layers) and per time, found in history data (HD),

[T₁−Xy, T₁] is the time interval comprising the most recent X years before time T₁,

T₁ being the last date for which history data is available,

∥Z∥_[T1,T2] denotes a norm of Z over time interval [T₁,T₂],

ϵ₁ being a positive number small with respect to 1, and

ϵ₂ being a positive number small with respect to 1, and p being a positive number close to 100%.

Such feature allows to practically defining the situation when a simulator can be considered as delivering a satisfactory history match, which is a necessary condition, to be a candidate for showing satisfactory forecasting capacities.

According to another feature of the present disclosure, X=5, ϵ₁=0.2, and ϵ₂=0.15 and p=90%.

Such feature allows framing the conditions under which the history match can be considered as accurate; it gives an order of magnitude of the accuracy, which can be expected for satisfactory forecasted quantities.

According to another feature of the present disclosure, the Vapnik condition is expressed as

$\frac{h}{m}\le 0,$

where:

h is the Vapnik-Chervonenkis dimension of the space S of solutions, from which the simulator is selected, and

m is the number of independent measures available in history data. Such feature allows framing the conditions, when they can be calculated, for achieving a satisfactory forecasting capacity.

According to another feature of the present disclosure, the Vapnik condition is verified by realizing a satisfactory blind test over N years, as defined by the following steps:

• • determining a production simulator matching history data over a time interval preceding T−Ny,
  • forecasting quantities produced over a time interval [T−Ny, T] given production parameters over the same time interval [T−Ny, T],
  • the blind test being satisfactory when:

$\frac{{Q_{\varphi \mathrm{ktb}}-Q_{\varphi \mathrm{ktbHD}}}_{\left[T-\mathrm{Ny},T\right]}}{{Q_{\varphi \mathrm{ktbHD}}}_{\left[T-\mathrm{Ny},T\right]}}\le {\varepsilon }_1,$

for more than p % of existing Q_φktbHD values and the cumulated oil produced over the same [T−Ny, T] time interval is accurate up to ϵ₂, where

Q_φktb are quantities produced per phase, per well, per layer (or group of layers) and per time, determined by the production simulator,

Q_φktbHD are the same quantities produced per phase, per well, per layer (or group of layers) and per time, found in history data,

[T−Ny, T] is the time interval comprising the most recent N years before time T,

T being the last date for which history data is made available,

∥X∥_[T1,T2] denotes a norm of X over time interval [T₁,T₂],

ϵ₁ being a positive number small with respect to 1, and

ϵ₂ being a positive number small with respect to 1, and p being a positive number close to 100%.

Such feature allows defining a discriminating test for deciding whether a simulator exhibits a satisfactory forecasting capacity, without looking at the way it has been set up, but rather at the actual results it can achieve. Such feature allows carrying out the test without the need to wait for years and measure the results, but rather to use data already available to the third party conducting the test, and not available to the party setting up the simulator (the latter party is “blind”).

According to another feature of the present disclosure, ϵ₁=0.1, and ϵ₂=0.1, p=90% and N=3.

Such feature allows framing practical conditions for the blind test, which can be accepted for properly discriminating a satisfactory simulator from a non-satisfactory simulator.

According to another feature of the present disclosure, the Vapnik condition is a forecast stability property verified when if

$\frac{{\mathrm{PP}-{\mathrm{PP}}^{\prime }}_{\left[T_0,T\right]}}{{\mathrm{PP}}_{\left[T_0,T\right]}}\le \varepsilon \mathrm{and}$ $\frac{{Q_{\varphi \mathrm{ktbHD}}-Q_{\varphi \mathrm{ktbHD}}^{\prime }}_{\left[T_0,T\right]}}{{Q_{\varphi \mathrm{ktbHD}}}_{\left[T_0,T\right]}}\le \varepsilon \mathrm{then}\frac{{Q_{\varphi \mathrm{ktb}}-Q_{\varphi \mathrm{ktb}}^{\prime }}_{\left[T,T+\mathrm{Ny}\right]}}{{Q_{\varphi \mathrm{ktb}}}_{\left[T,T+\mathrm{Ny}\right]}}\le n\varepsilon$

where

PP are production parameters,

PP′ are slightly varying different production parameters,

Q_φktbHD are quantities produced per phase, per well, per layer (or group of layers) and per time, found in history data,

Q_φktbHD′ are slightly varying different quantities produced per phase, per well, per layer (or group of layers) and per time, found in a slightly varying different history data,

Q_φktb are quantities produced per phase, per well, per layer (or group of layers) and per time, determined by the production simulator,

Q_φktb′ are quantities produced per phase, per well, per layer (or group of layers) and per time, determined by another close enough production simulator determined out of slightly varying different history data,

[T₀, T] is the time interval from time T₀ to time T, being respectively the initial and the last date for which history data is available,

[T, T+Ny] is the time interval comprising the Ny years after time T,

∥Z∥_[T1,T2] denotes a norm of Z over time interval [T₁,T₂],

ϵ being a positive number small with respect to 1, and

n being a small integer less than 5.

Such feature allows defining whether a simulator exhibits properties, which are a pre-requisite to reliably calculate forecasted quantities. This test does not require comparing with another existing simulator, nor does it require knowing about the way this simulator has been set up, as only its resulting properties are tested.

According to another feature of the present disclosure, ϵ=0.05, n=2 and N=3.

This feature allows practically framing the forecast stability, over a period of three years.

According to another feature of the present disclosure, the production simulator is built following the steps of:

• • defining initial detailed reservoir partition, rock properties, laws of reservoir physics and laws of well physics,
  • upscaling said reservoir partition, rock properties, laws of reservoir physics and laws of well physics until said Vapnik condition is verified, and
  • calibrating said production simulator by choosing among the production simulator candidate solutions, the candidate solution, defined by α and minimizing an expected so called “empirical risk” R_emp(α).

Upscaling means reducing complexity. Upscaling is carried out in such a way that the space of solutions candidates for becoming a simulator is such as

$\Phi \left(\frac{h}{m},\delta \right)$

is as small as possible, where:

δ is a positive number close to zero, 1-δ defining a probability and

Φ is a positive function defined by:

$\Phi \left(\frac{h}{m},\delta \right)=2\sqrt{\frac{2}{m}\left(h\log \left(\frac{2\mathrm{em}}{h}\right)+\log \left(\frac{2}{\delta }\right)\right)},$

where

h is the Vapnik-Chervonenkis dimension of the space of solutions, and

m is the number of independent measures available in history data,

e is equal to exp(1).

Calibrating the production simulator among the candidates of the space of solutions means making the “empirical risk” R_emp(α) as small as possible, while choosing the appropriate parameters α that completely defines a solution within the space of solutions. This empirical risk R_emp(α) measures a (positive) distance between actual past data and the corresponding data calculated by the simulator.

This upscaling process is a way to reach a trade-off between

$\Phi \left(\frac{h}{m},\delta \right)$

and R_emp(α) values, which need to be together as small as possible, as the objective is to minimize their sum

$R_{\mathrm{emp}}\left(\alpha \right)+\Phi \left(\frac{h}{m},\delta \right)$

with a given probability 1-δ.

In minimizing the above-mentioned sum, this minimizes the expected risk of the forecast, R(α), according to the Vapnik inequality:

$R\left(\alpha \right)\le R_{\mathrm{emp}}\left(\alpha \right)+\Phi \left(\frac{h}{m},\delta \right)$

Such feature allows minimizing the expected risk of the forecast R(α), while starting from large values of

$\Phi \left(\frac{h}{m},\delta \right)$

(as starting models are complex) and small empirical risks R_emp (α) (as such complex models can properly picture past data) and gradually decrease the complexity, and therefore

$\Phi \left(\frac{h}{m},\delta \right)$

while keeping a good match of past data, therefore keeping the empirical risk R_emp(α) small.

According to another feature of the present disclosure, said reservoir partition is upscaled following the steps of:

• • partitioning a reservoir G into elementary parts G_ab, such that

$G=\bigcup _{a=1}^A\bigcup _{b=1}^BG_{\mathrm{ab}}$

with G_ab∩G_a′b′=Ø for (a,b)≠(a′,b′), where a∈{1 . . . A} describes a x-y area, and b∈{1 . . . B} describes one or several layers,

• • grouping adjacent elementary parts exhibiting homogeneous rock properties, into sub-geologies G_z where z∈{1, . . . , Z}.

Such feature gives a practical way of upscaling the field geology in identifying a reasonably small number of sub-geologies G_z.

According to another feature of the present disclosure, rocks properties are upscaled by following a step of averaging out rock properties over each sub-geology, according to formula:

${\mathrm{RP}}_z=\frac{1}{V_z}\underset{G_z}{\int \int \int }\mathrm{RP}\left(x,y,z\right)\mathrm{dxdydz},$

where V_z is the volume of sub-geology G_z.

This feature explains a way to practically define the properties to be used in a given sub-geology G_z.

According to another feature of the present disclosure, laws of reservoir physics are upscaled in such a way they apply with functioning parameters of the sub-geology and wherein space and times scales associated with the sub-geology are determined in such a way that the associated space of solutions is consistent with the complexity of history data at the well level.

Such feature gives the rules to look for upscaled reservoir physics laws, in focusing on their behavior measured at well level.

According to another feature of the present disclosure, the production simulator is built following the steps of:

• • defining initial coarse reservoir partition, rock properties, laws of reservoir physics and laws of well physics,
  • downscaling said reservoir partition, rock properties, laws of reservoir physics and laws of well physics while keeping said Vapnik condition verified, until the production simulator matches history data and
  • calibrating said production simulator by choosing among the production simulator candidate solutions, the candidate solution minimizing a forecast expected risk.

Downscaling means increasing complexity, starting from a simplistic description of the field, and adding relevant reservoir and/or well phenomena, which will properly picture the behavior of the entire field, well by well. Downscaling is carried out in such a way that the space of solutions candidates for becoming a simulator remains such as

$\Phi \left(\frac{h}{m},\delta \right)$

is as small as possible, where:

δ is a positive number close to zero, 1-δ defining a probability and

Φ is a positive function defined by:

$\Phi \left(\frac{h}{m},\delta \right)=2\sqrt{\frac{2}{m}\left(h\log \left(\frac{2\mathrm{em}}{h}\right)+\log \left(\frac{2}{\delta }\right)\right)},$

where

h is the Vapnik-Chervonenkis dimension of the space of solutions, and

m is the number of independent measures available in history data,

e is equal to exp(1).

Calibrating the production simulator among the candidates of the space of solutions means making the “empirical risk” R_emp(α) as small as possible, while choosing the appropriate parameters “α” that completely defines a solution within the space of solutions. This empirical risk measures a distance between actual past data and the corresponding data calculated by the simulator.

This downscaling process is a way to reach a trade-off between

$\Phi \left(\frac{h}{m},\delta \right)$

and R_emp(α) values, which need to be together as small as possible, as the objective is to minimize their sum

$R_{\mathrm{emp}}\left(\alpha \right)+\Phi \left(\frac{h}{m},\delta \right)$

with a given probability 1-δ.

In minimizing the above-mentioned sum, this minimizes the expected risk of the forecast, R(α), according to the Vapnik inequality:

$R\left(\alpha \right)\le R_{\mathrm{emp}}\left(\alpha \right)+\Phi \left(\frac{h}{m},\delta \right)$

Such feature allows minimizing the expected risk of the forecast R(α), while starting from small values of

$\Phi \left(\frac{h}{m},\delta \right)$

(as starting models are simplistic) and large empirical risks R_emp(α) (as such simplistic model do not properly picture past data) and gradually increasing the complexity, and therefore increasing

$\Phi \left(\frac{h}{m},\delta \right)$

in order to better match past data, therefore decreasing the empirical risk R_emp(α).

According to another feature of the present disclosure, said reservoir partition is downscaled following the steps of:

• • starting from the whole reservoir,
  • parting said reservoir into sub-geologies where a substantial change of reservoir properties is present at around the boundary between said sub-geologies.

Such feature describes the process of going from global field description to a refined understanding, where local behaviors are identified at a finer scale G_z.

According to another feature of the present disclosure, wherein rocks properties are downscaled by defining new separate rock field properties over each sub-geology.

Such feature describes the principle of allocating different physical properties at the G_z level.

According to another feature of the present disclosure, laws of reservoir physics are downscaled in such a way they apply with functioning parameters of the sub-geology and wherein space and times scales associated with the sub-geology are determined in such a way that the associated space of solutions is consistent with the complexity of history data at the well level.

Such feature gives the rules to look for downscaled reservoir physics laws, in focusing on their behavior measured at well level.

According to another feature of the present disclosure, the laws of reservoir physics are derived from the (Navier-)Stokes equations of conservation of momentum and conservation of mass for a fluid or gas evolving in a rock modeled as a porous medium, only characterized by its averaged porosity and relative/permeability per time per phase per sub-geology.

This feature explains how general laws of reservoir physics can be simplified at the sub-geology level, while averaging out some key parameters over the entire sub-geology G_z.

According to another feature of the present disclosure, the laws of reservoir physics comprise formulas:

$u_{\varphi k}=-\frac{k}{{\mu }_b}k_{r\varphi }\left(\nabla P_{\varphi }-{\rho }_{\varphi }g\right),$

where:

u_φkz is the velocity of phase φ at well k, at sub-geology z,

μ_φc is the viscosity of phase φ, at sub-geology z,

P_φz is the pressure of phase φ at sub-geology z,

ρ_φc is the density of phase φ, at sub-geology z,

g is the gravity vector,

k is a permeability coefficient, and

k_rφz is a relative permeability coefficient of phase φ, at sub-geology z. k_rφz is a function of time, through its dependency on phase saturation.

This feature explains what law is considered in practice for fluid transport and how it is parameterized.

According to another feature of the present disclosure, the laws of reservoir physics further comprise heat transfer laws between a fluid or a gas and a rock given by the following formulas:

Φ∂_t(ρ_wS_w+ρ_gS_g)+∇·(ρ_wu_w+ρ_gu_g)=0

Φ∂_t(ρ_oS_o)+∇·(ρ_ou_o)=0,

S _o +S _w +S _g=1,

∂_t(E_s+E_f)+∇·(U_f−λ(T)∇T)=0, where

T=T(x,y,z) is the temperature and can vary over the reservoir,

• • φ is the rock porosity,
  • S_φ is the saturation of phase φ,
  • ρ_φ is the density of phase φ,
  • u_φ is the velocity of phase φ,
  • E is the volumetric internal energy (indices s and f corresponding respectively to solid and fluid phases),

U _f =U _g +U _o +U _w is the enthalpy flux, where U_φ=ρ_φh_φu_φ,

h_φ is the specific enthalpy of phase φ,

h_g=h_w+Λ, Λ being the specific heat of water vaporization,

g, w, o indexes being respectively for gas, water and oil,

λ(T) is the reservoir conductivity coefficient.

This feature describes which heat transfer laws need to be considered and how they can be parameterized.

According to another feature of the present disclosure, the laws of reservoir physics further comprise phase laws given by:

μ_φz=μ_φz(P,T) the viscosity of phase φ at sub-geology z, function of the local pressure P and the local temperature T;

ρ_φz=ρ_φz(P,T) the density of phase φ at sub-geology z, function of the local pressure P and the local temperature T;

k_rφz=k_rφz(P,T) the relative permeability coefficient of phase φ at sub-geology z, function of the local pressure P and the local temperature T.

This feature explains how laws of fluid propagation and laws of heat transfer are coupled, through dependency of pressure and temperature.

According to another feature of the present disclosure, laws of well physics comprise formulas:

Q _φktz =T _kz(PP_ktz,u_φkz) where:

Q_φktz is the quantity produced of phase φ at well k in sub-geology z at time t,

T_kz is a transfer function of well k, in sub-geology z

PP_ktz are the production parameters applied to well k in sub-geology z at time t,

u_φkz is the velocity of phase φ at well k in sub-geology z.

This feature explains how wells behaviors can be modeled, for instance through transfer functions. Such functions have the same “fineness” than the modeling of the reservoir: if a sub-geology G_z includes several layers b, these will be averaged out together at the well level.

The present disclosure also concerns a method for optimally exploiting a mature hydrocarbon field, comprising the steps of:

• • building a production simulator according to any one of the preceding embodiments,
  • iterating several runs of said production simulator in order to find the optimal production parameters optimizing a gain value derived from said quantity produced,
  • applying said optimal production parameters so obtained to exploit the hydrocarbon field.

This feature describes a generic optimization process, aimed at maximizing a gain function on the field. It fully uses an aforementioned feature of the production simulator, namely its relative low-complexity, which allows very low computing time for calculating a given production scenario, defined by its corresponding production parameters. This low computing time renders possible playing a vast number of scenarios, typically more than millions of them in 24 hours.

According to another feature of the present disclosure, said optimized gain value is a net present value or reserves of said hydrocarbon field.

This feature illustrates a practical use of the present disclosure, which turns out to generate large financial gains, against non-optimized scenarios.

According to another feature of the present disclosure, said net present value can be determined using the formula:

$\mathrm{NPV}=\sum _i\left(\sum _k\sum _zP_{\mathrm{kz}}*\left(1-R_{\mathrm{ik}}\right)\right)*\frac{S_i}{{\left(1\text{'}0d\right)}^i}-\sum _i\left(\sum _kI_{\mathrm{ik}}+\sum _k{\mathrm{OC}}_{\mathrm{ik}}+\sum _k\sum _z\left({\mathrm{TO}}_i*P_{\mathrm{kz}}+{\mathrm{TL}}_i*L_{\mathrm{kz}}\right)\right)*\frac{1}{{\left(1+d\right)}^i}$

where:

P_kz is oil production (in barrels) for well k and sub-geology z,

R_ik is tax and royalties for well k and year i,

S_i is the oil sale price (per barrel) for year i,

d is the percentage discount rate,

I_ik is investment made on well k during year i,

OC_ik is operating costs for well k during year i,

L_kz is the liquid production (in barrels) for well k and sub-geology z,

TO_i is treatment cost (per barrel of oil), for year i.

TL_i is treatment cost (per barrel of liquid), for year i.

This feature gives a practical way to define a gain function, which can be completely calculated using the production of wells, per phase per time and per sub-geology, which is itself calculated by the production simulator for any given set of production parameters.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated herein and form a part of the specification, illustrate the present disclosure and, together with the description, further serve to explain the principles of the disclosure and to enable a person skilled in the pertinent art to make and use the disclosure.

FIG. 1 shows an illustrative hydrocarbon field.

FIG. 2 is a block diagram of a simulator according to an embodiment.

FIG. 3 is a block diagram detailing a production simulator according to an embodiment.

FIG. 4 is a schematic view of an example oil field.

FIG. 5 illustrates example scenarios of an oil field.

FIG. 6 illustrates an example memory representation of multiple scenarios according to an embodiment.

FIG. 7 illustrates a multiple scenarios simulation run using a Single Instruction, Multiple Threads (SIMT) architecture according to an embodiment.

FIG. 8 illustrates a multiple scenarios simulation run using a Single Instruction, Multiple Data (SIMD) architecture according to an embodiment.

FIG. 9 illustrates a process for simulating the production of an oil field according to an embodiment.

FIG. 10 illustrates an example computer system which can be used to implement embodiments.

DETAILED DESCRIPTION OF EMBODIMENTS

According to FIG. 1, a hydrocarbon field is a subterranean (or submarine) reservoir G comprising rock and gas and fluids of which oil. Said reservoir G is drilled from the surface 10 by several wells 11, in a known fashion. Said wells 11 may be injectors in charge of injecting a fluid, generally water, gas or steam, into the reservoir in order to change its properties, for instance to get its pressure or its temperature modified. Said wells 11 may instead be producers in charge of extracting gas, oil and associated fluids (and maybe other components) from the reservoir. Appropriate fluids or gas can also be injected into these producers, in order to enhance their production effectiveness or to get local reservoir properties to be modified.

The aim of the present disclosure is to calculate future hydrocarbon production per well, per phase and per sub-geology over time.

Another aim of the present disclosure is to determine and apply the optimal production parameters PP, that is, the production parameters PP that optimize a gain value e.g. the production, the NPV, reserves, of a mature hydrocarbon field.

The present disclosure concerns a method, for optimally exploiting a mature hydrocarbon field, based on an optimizer 1. Said optimizer 1 uses a production simulator 2 module that provides quantities produced Q_φzkt per phase φ, per well k, 11, per sub-geology z, and per time t (said time is generally expressed in month) as a function of production parameters PP. The phase φ can be oil, gas, water or other relevant fluids like chemical additives (oil may be considered as a single phase or split into many different hydrocarbon phases). There can be several oil or gas phases, depending of hydrocarbons chains. Gas can be hydrocarbon, vapor or any other element.

Production parameters PP comprise all the driving parameters impacting the production. These production parameters include: injection rates that define the behavior of injectors, production rates, e.g. in terms of pumping rates, that define the behavior of producers but also well work-over, like perforation, completion, acidification, fracturation . . . , or some other impacting choices such as a possible conversion of a producer into an injector or else a drilling of an new additional well (be it a producer and/or an injector). PP can be sub-geology-dependent (e.g. for wells produced with multiple strings).

Such an optimizer 1 allows the test of new production strategies that would imply modifying production parameters against the baseline (“business as usual” production strategy).

Since said production simulator 2 is simple enough to quickly compute the result, the optimizer 1 may iterate many runs of said production simulator 2, thus trying a lot of strategies or scenarios, in order to find the optimal production parameters PP that optimize a gain value 5 derived from said quantity produced Q_φzkt.

As figured by block diagram of FIG. 2, the essential brick of optimizer 1 is production simulator 2. For a given set of production parameters PP, it provides a forecast of quantities produced Q_φzkt.

From said quantities produced Q_φzkt a gain value 5, to be optimized, can be computed by a module 4.

In optimizer 1 a module 3 is in charge of proposing said scenarios in the form of sets of production parameters PP. Several heuristics/algorithms may be used in said module 3, ranging from random Monte-Carlo methods to more efficient context aware methods to propose relevant candidate sets of production parameters PP.

As is figured by a feedback arrow, said module 3 may take advantage of the feedback value 5, to determine new sets, and to drive toward new scenarios, according to known optimization methods.

By so iterating, optimizer 1 may determine an optimal set of production parameters that optimize the gain value 5.

The method also comprises applying said optimal production parameters PP so obtained to drive the exploitation of the hydrocarbon field. As simulated, the expected gain or revenue can be optimal and at least better than the baseline.

To be efficient, according to an essential feature of the present disclosure, without reproducing the drawback of prior art solutions, said production simulator 2 must first match history data HD, that is provides a result that rather well compare to quantities produced Q_φzkt for the past, when applying the production parameters PP historically applied, thus differentiating from too simplistic models. Also, said production simulator 2 must verify a Vapnik condition.

The production simulator 2 is intended to produce estimated data for the future, starting from instant T, as accurate as possible. One necessary condition to pretend to such an accuracy is first that said production simulator matches history data HD, that is its reproduces known data in the past, that is so-called history data HD. Matching history data thus means the production simulator 2 is accurate enough to reproduce data produced such as quantities produced Q_φzkt per phase φ, per well k, per sub-geology z and per time t, and total cumulated oil produced, when known, that is before T. The quality of said match is appreciated by relatively comparing known quantities produced Q_φzktHD recorded in history data HD to quantities produced Q_φzkt as computed by mean of production simulator 2, to which past production parameters are applied. Cumulated oil produced is also compared between the sum of quantities produced Q_φzktHD recorded in history data HD, and the sum of quantities produced Q_φzkt as computed by mean of production simulator 2. Such a comparison can be expressed by

$\frac{{Q_{\varphi \mathrm{ckt}}-Q_{\varphi \mathrm{zktHD}}}_{\left[T_1-\mathrm{Xy},T_1\right]}}{{Q_{\varphi \mathrm{cktHD}}}_{\left[T_1-\mathrm{Xy},T_1\right]}}\le {\varepsilon }_1$

for p % of the quantities Q_φzktHD and

$\underset{k,z,t\in \left[T_1-\mathrm{Zy},T_1\right]}{\sum Q_{\varphi \mathrm{ckt}}}-\underset{k,z,t\in \left[T_1-\mathrm{Xy},T_1\right]}{\sum Q_{\varphi \mathrm{cktHD}}}\le {\varepsilon }_1$

where

Q_φzkt are quantities produced per phase φ, per well k, per sub-geology z and per time t, determined by the production simulator 2,

Q_φzktHD are the same quantities produced per phase φ, per well k, per sub-geology z and per time t, found in history data (HD),

[T₁−Xy, T₁] is the time interval comprising the most recent X years before time T₁,

T₁ being the last date for which history data HD is available,

∥Z∥_[T1,T2] denotes a norm of Z over time interval [T₁,T₂],

ϵ₁ being a positive number small with respect to 1, and

ϵ₂ being a positive number small with respect to 1,

p being a number close to 100%.

The length of the interval considered for the matching condition may be adapted to the available length of history data. However if possible, length of five years is considered to be satisfactory. Such a length is consistent with a comparable horizon of expected accurate forecast for the production simulator 2, of five years.

The particular values of X=5, p=90%, ϵ₁=0.2 and ϵ₂=0.15 have been found satisfactory to detect a good match.

The norm used herein for comparison may be any one among: least square or weighted least square.

In addition to the matching condition, an accurate production simulator 2 must also verify a Vapnik condition. Said Vapnik condition may be expressed in different ways.

One of the results of the learning theory as developed by Vapnik shows that, from a given history data HD set of dimension m, that is, comprising m independent measures, one can derive, by learning, a model within a space of VC-dimension h, with a reliability forecasting expected risk R(α) associated to parameters α, bounded from above by the sum of two additive terms R_emp(α)+Φ. The first term R_emp(α) is an empirical risk associated with parameters α, indicative of the quality of the matching to history data HD provided by said parameters α. The second term Φ is characteristic of the model and can be expressed by

$\Phi \left(\frac{h}{m},\delta \right)=2\sqrt{\frac{2}{m}\left(h\log \left(\frac{2\mathrm{em}}{h}\right)+\log \left(\frac{2}{\delta }\right)\right)},$

wherein δ is a positive number close to zero, 1-δ defining a probability of said forecast expected risk R(α) is indeed bounded by R_emp(α)+Φ, where h is the aforementioned VC-dimension of the space S of solutions, also named Vapnik-Chervonenkis or VC dimension, and in is the aforementioned number of independent measures available in history data HD.

According to said result, prior art oversimplified models all suffer from a high R_emp(α) due to the over simplification of the model. Instead, prior art complex meshed models may provide a small R_emp(α), but suffer from a high Φ value, due to an high

$\frac{h}{m}$

ratio, since the VC-dimension h of the space of solutions is too large with respect to the number m of independent measures available in history data HD.

A trade-off has then to be reached. m value is constant and given by history data HD available for a given hydrocarbon field. The size h of the space of solutions has then to be adapted so that the

$\frac{h}{m}$

ratio remains small. With an average objective of a five years forecast with a reliability of around +/−5%, a value

$\frac{h}{m}\le 0,$

has been found to be compatible with both the objective reliability and the possibility of scaling, as will be detailed further, to obtain a model.

Another way to express and check said Vapnik condition is by realizing a blind test. A blind test over N years may be realized by parting history data into two time intervals. A first “past” interval runs from an initial time T₀ for which history data are available, to a time T−Ny preceding the final time T by N years. A second “blind” interval runs from T−Ny to time T corresponding to the last time for which history data is available. T is generally the actual time. The data of “past” interval are considered as known, and are used to build a matching production simulator 2, by learning on said data. The data of “blind” interval are considered unknown or at least are concealed during the building of production simulator. They are then compared to data forecast by production simulator 2 over said “blind” interval given production parameters over the “blind” interval. A blind test is considered satisfactory when the forecast data reproduce with enough accuracy the history data over said concealed interval.

A blind test then comprises the following steps:

• • determining a production simulator 2 matching history data HD over a time interval preceding T−Ny,
  • forecasting quantities produced over a time interval [T−Ny, T],
  • the blind test being satisfactory when:

$\frac{{Q_{\varphi \mathrm{ktb}}-Q_{\varphi \mathrm{ktbHD}}}_{\left[T-\mathrm{Ny},T\right]}}{{Q_{\varphi \mathrm{ktbHD}}}_{\left[T-\mathrm{Ny},T\right]}}\le {\varepsilon }_1,$

for more than p % of existing Q_φktbHD values and the cumulated oil produced over the same [T−Ny, T] time interval is accurate up to ϵ₂, where

Q_φktb are quantities produced per phase, per well, per layer (or group of layers) and per time, determined by the production simulator,

Q_φktbHD are the same quantities produced per phase, per well, per layer (or group of layers) and per time, found in history data,

[T−Ny, T] is the time interval comprising the most recent N years before time T,

T being the last date for which history data is made available,

∥Z∥_[T1,T2] denotes a norm of Z over time interval [T₁,T₂],

ϵ₁ being a positive number small with respect to 1, and

ϵ₂ being a positive number small with respect to 1, and p being a positive number close to 100%.

The particular values of p=90%, ϵ₁=0.1 and ϵ₂=0.1 have been found satisfactory to assert a good blind test. A period of N=3 years is also satisfactory.

The norm used herein for comparison may be any one, for instance: least square or weighted least square.

Another way to check the Vapnik condition is expressed by measuring a forecast stability property. Such property test the stability of the production simulator 2 toward its output data in response of a slight variation of the input data. Input data are history data HD comprising production parameters PP and quantities produced Q_φktbHD in the past, that is over a [T₀, T] time interval. Based on these nominal data PP and Q_φktbHD, a nominal production simulator 2 is built. These data are then slightly varied to obtain corresponding data PP′ and Q_φktbHD′. Based on these slightly varied input data another production simulator is built that is expected to be close enough to nominal production simulator 2.

Output data, that is quantities produced Q_φktb and Q_φktb′ are forecast respectively by nominal production simulator 2 and by other production simulator, over a future time interval [T, T+My], and are then compared.

The forecast stability property is verified when

$\frac{{\mathrm{PP}-{\mathrm{PP}}^{\prime }}_{\left[T_0,T\right]}}{{\mathrm{PP}}_{\left[T_0,T\right]}}\le \varepsilon \mathrm{and}\frac{{Q_{\varphi \mathrm{ktbHD}}-Q_{\varphi \mathrm{ktbHD}}^{\prime }}_{\left[T_0,T\right]}}{{Q_{\varphi \mathrm{ktbHD}}}_{\left[T_0,T\right]}}\le \varepsilon$

implies that

$\frac{{Q_{\varphi \mathrm{ktb}}-Q_{\varphi \mathrm{ktb}}^{\prime }}_{\left[T,T+\mathrm{My}\right]}}{{Q_{\varphi \mathrm{ktb}}}_{\left[T,T+\mathrm{My}\right]}}\le n\varepsilon$

Where

PP are production parameters,

PP′ are slightly varying production parameters,

Q_φktbHD are quantities produced per phase φ, per well k, per layer (or group of layers) b and per time t, found in history data HD,

Q_φktbHD′ are slightly varying quantities produced per phase φ, per well k, per layer (or group of layers) b and per time t, found in a slightly varying history data,

Q_φktb are quantities produced per phase φ, per well k, per layer (or group of layers) b and per time t, determined by the production simulator 2,

Q_φktb′ are quantities produced per phase φ, per well k, per layer (or group of layers) b and per time t, determined by another close enough production simulator determined out of slightly varying history data,

[T₀, T] is the time interval from time T₀ to time T, being respectively the initial and the last date for which history data (HD) is available,

[T, T+My] is the time interval comprising the M years after time T,

∥Z∥_[T1,T2] denotes a norm of Z over time interval [T₁,T₂],

ϵ being a positive number small with respect to 1, and

n being a small integer lesser than 5.

The important value here is the ratio n between the input variation allowed ϵ, and the output variation nϵ obtained. The values of ϵ=0.05, n=2 have been found satisfactory. A horizon of M=3 years is also satisfactory when considering the accuracy objectives.

The norm used herein for comparison may be any one, for instance: least square or weighted least square.

The three described ways to check the Vapnik condition can be use either separately or concurrently.

As already mentioned, the heart of the present disclosure is the production simulator 2. Said production simulator 2 may use behavioral physics laws such as reservoir physics laws 6 and well physics laws 7 in order to be sufficiently precise to accurately match the history data HD. However the VC-dimension h of the space of solutions S must remain small enough so as to verify said Vapnik condition.

Two approaches are proposed here to reach such result. Both approaches comprise using physics laws and applying them to a reservoir model.

In a first approach, also called upscaling, both said reservoir model and physics laws start with detailed formulations that are gradually simplified/scaled by approximating when and where some homogeneity can be observed, until a stop condition is reached. Said condition is reached by upscaling a reservoir partition, rock properties RP, laws of reservoir physics 6 and laws of well physics 7 until said Vapnik condition is verified.

According to a second approach, also called downscaling, both said reservoir model and physics laws start with coarse formulations that are gradually complexified/scaled by detailing when and where some substantial variation can be observed, until a stop condition is reached. Said condition is reached by downscaling a reservoir partition, rock properties RP, laws of reservoir physics 6 and laws of well physics 7 will keeping said Vapnik condition verified, until the production simulator 2 matches history data HD.

Both approaches may be used alone separately or may be used alternately. For instance, one may upscale the geology from a fine model by merging parts, and then locally downscale back one sub-geology that cannot be matched, thus revealing a substantial change of properties inside.

At this step, whatever the approach used, the structure of the production simulator 2 is defined. However, said production simulator 2 depends on several parameters α, that must be further trimmed in order for the production simulator 2 to be operative. Depending on parameters α, several production simulators 2 may be candidate solutions. These several candidate solutions are checked against a forecast expected risk R(α) for each possible set of parameters α until finding the optimal set of parameters α minimizing said expected risk R(α).

According to FIG. 1, a reservoir G may be modeled by considering its geology, that is, the closure of the rock volume containing oil, and parting it using a reservoir partition.

When using an upscaling approach, a first fine partition may be defined by crossing x-y areas indexed by letter a, each e.g. comprising at least one well 11, and z areas indexed by letter b. In the example of FIG. 1, four such x-y areas (a=1 . . . 4) are used and five (b=1 . . . 5) 1 areas are used to initially part the reservoir G. The reservoir G is thus initially parted into elementary parts G_ab such that

$G=\bigcup _{a=1}^A\bigcup _{b=1}^BG_{\mathrm{ab}}$

with G_ab∩G_a′b′=Ø for (a,b)≠(a′,b′) (two different parts are disjoined), where a{1 . . . A} describes a x-y area, and b{1 . . . B} describes a 1 layer.

Starting from this initial reservoir partition, said partition is upscaled by considering elementary parts G_ab, and gradually merging adjacent ones together when they exhibit homogeneous rock properties RP. The result is a new coarser partition, comprising sub-geologies G_z where z∈{1 . . . Z}. In FIG. 1, three such sub-geologies are shown. G₁ gathers G₁₁, G₁₂, G₁₃, G₁₄ and G₁₅, exhibiting a homogenous behavior. G₂ gathers G₂₁, G₂₂, G₂₃, G₂₄, G₂₅, G₃₁, G₃₂, G₃₃, G₃₄ and G₃₅. G₃ gathers G₃₁, G₃₂, G₃₃, G₃₄ and G₃₅. In real cases a sub-geology typically includes from 3 to 50 wells. Said upscaling of the geology partition is a first way to reduce the dimension h of the space of solutions S. The size of the space of solutions then depends on the number of geologies which is of an order of few units instead of over 100 000s cells as in a meshed prior art model.

When using a downscaling approach, the process starts from a coarse partition, e.g. composed of a single part corresponding to the whole reservoir G. Said coarse partition s then parte into sub-geologies G_z where z∈{1 . . . Z} where and every time a substantial change of properties is present around the boundary between said sub-geologies. So doing, a homogeneous behavior may be expected in every sub-geology, while keeping the number of such sub-geologies as reduced as possible.

Rock properties RP comprises mainly the rock porosity, permeability k and relative permeability k_rφz, used to determine the dynamic of fluid through rock, but also other property such as net pay, heat capacity or conductivity, as relevant. Since the sub-geology becomes the new elementary volume unit, all of these rock properties RP, are considered to be fairly constant over a given sub-geology G_z.

When using a upscaling approach rock properties RP are homogenized over each sub-geology G_z, according to averaging formula:

${\mathrm{RP}}_z=\frac{1}{V_z}\underset{G_z}{\int \int \int }\mathrm{RP}\left(x,y,z\right)\mathrm{dxdydz},$

where V_z is the volume of sub-geology G_z.

When instead using a downscaling approach, rocks properties (RP) are downscaled by defining new separate rock properties RP_z over each new sub-geology (G_z) obtained by parting along a discontinuity.

With respect to FIG. 3, is now described in more details the content of production simulator 2. Said production simulator 2 module may be parted into two main modules 6, 7. A first module implements laws of reservoir physics 6. Given a set of production parameters PP, said module provides dynamic characteristics of the fluids across the reservoir, at least at the entrance of each well 11 and at each sub-geology b. Said dynamic characteristics are e.g. expressed as a velocity field.

Said laws of reservoir physics 6 are either upscaled or downscaled in such a way they apply with functioning parameters of the sub-geology G_z. Also, space and times scales associated with the sub-geology G_z are determined in such a way that the associated space of solutions is consistent with the complexity of history data HD at the well level. For instance, if insufficient data richness is available on layer-related production (m too small), well production will be summed up, partially (some layers are grouped) or totally (all layers are grouped):

$Q_{\varphi \mathrm{kt}}=\sum _{b\in \left\{1\dots B\right\}}Q_{\varphi \mathrm{ktb}}$

In such instance, the production simulator 2 only describes production per well and not production per well and per layer (or group of layers).

The upscaling/downscaling process also applies to the laws. A variation in number and complexity of laws of physics and of well may be obtained by using more simple or instead more detailed behavior models and formulas. For instance Darcy's law is a simplification of the more general Navier-Stokes's law. The order of the formulas can also be changed either increased or decreased according to the Vapnik condition intended. Some parameters may either be constant or refined by a varying function of another parameter or variable. Depending on the way the field has been exploited so far, heat transfer can be considered or not, thus employing or disregarding heat transfer laws.

In a similar fashion, well transfer functions may sometimes be linearized or approximated around their functioning point.

Such concordance is to be considered on partial time intervals or on partial spatial areas of the hydrocarbon field in accordance with the richness (dimension m) of detail of data available in history data HD.

The laws of reservoir physics 6 are, in a known manner, derived from the (Navier)-Stokes equations of conservation of momentum and conservation of mass for a fluid, said fluid successively corresponding to phases, evolving in a rock modeled as a porous medium only characterized by its averaged porosity Φ, permeability k and relative permeability k_rφz.

Usual laws of reservoir physics 6 are derived by using general flow model such as the Navier-Stokes equations coupled with multiphase transport considered at the microscopic scale. Fluids are evolving in a domain obtained as the complement of rock, which may contain a large variety of obstacle topologies and geometries. The upscaling/downscaling process consists in determining the effective dynamics of fluids in the porous medium in which the rock structure is only described by the average porosity, permeability and relative permeability at a given time t.

Fluid propagation in the reservoir G is derived from the general Navier-Stokes equations, which describe conservation of momentum for a fluid:

$\rho \left(\frac{\partial u}{\partial t}+\left(u·\nabla \right)u\right)=\nabla ·T+f$

where u is the fluid velocity, ρ is the fluid density, T the stress tensor and f body forces acting on the fluid (for oil in a reservoir, f is the gravity: f=ρ·g).

These are complemented by the conservation of fluid mass:

$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho u\right)=0$

And for each phase φ∈{oil, water, gas}:

$\frac{\partial \left(\rho S_{\varphi }\right)}{\partial t}+\nabla ·\left(\rho S_{\varphi }u\right)=\nabla ·J_{\varphi }$

Where S_φ denotes phase φ saturation and J_φ at diffusion flux. Of course, one has:

$\sum _{\varphi }S_{\varphi }=1\mathrm{and}\sum _{\varphi }J_{\varphi }=0$

Since reservoir flows are usually characterized by small Reynolds numbers, i.e. creeping flows, the Navier-Stokes equations can generally be simplified into Stokes equations, assuming stationary, incompressible flow:

μ∇²u_i+ρg_i−∂_ip=0 and ∇·u=0 where u is the viscosity, u_i is the velocity in the i direction, g_i is the gravity component in direction i and p is the pressure.

Using a volume averaging procedure, the Stokes equations may be homogenized so that the effective viscous resisting force is proportional to the velocity, and opposite in direction. One may therefore write in the case of isotropic porous media

$-\frac{\mathrm{\mu \Phi }}{k}u_i+\rho g_i-{\partial }_ip=0,$

where k denotes the permeability of the porous medium and Φ is the porosity. This gives the velocity in terms of pressure gradient

$u_i=-\frac{k}{\mathrm{\Phi \mu }}\left({\partial }_ip-\rho g_i\right),$

which gives Darcy's law. In the case of anisotropic porous media, one ends up with

$u_i=-\frac{1}{\mathrm{\Phi \mu }}\sum _jK_{\mathrm{ij}}\left({\partial }_jp-\rho g_j\right),$

where K denotes the symmetric permeability tensor and the sum is over the three directions j. The resulting equation of mass conservation of fluid is given by

$\frac{\partial \left(\mathrm{\rho \Phi }\right)}{\partial t}+\nabla ·\left(\rho u\right)=0.$

Actually, the fluid is generally made of gas (g), oil (o) and water (w), whose composition might be complex, depending for instance on their salinity and hydrocarbon chains. Mass conservation of phase φ∈{o,w,g} at the upscaled level may be written as:

$\frac{\partial \left({\rho }_{\varphi }S_{\varphi }\Phi \right)}{\partial t}+\nabla ·\left({\rho }_{\varphi }u_{\varphi }\right)=0$

where u_φ, ρ_φ, S_φ respectively denote the velocity, density, saturation and diffusion flux of phase φ.

In any case, one has to determine closure laws for phase velocity u_φ in terms of average flow characteristics like velocity u and pressure p and phase φ saturation. Consistency relations have to be satisfied such as

$\sum _{\varphi }S_{\varphi }=1.$

Closure of phase φ velocity involves the so-called relative permeability coefficient k_rφ accounting for differential mobility of phases

$u_{\varphi k}=-\frac{k}{{\mathrm{\Phi \mu }}_{\varphi }}k_{r\varphi }\left(\nabla P_{\varphi }-{\rho }_{\varphi }g\right)$

where μ_φ denotes the viscosity of phase φ.

Finally module 6 implements the laws of reservoir physics by following formulas:

$u_{\varphi \mathrm{kz}}=-\frac{k}{{\mathrm{\Phi \mu }}_{\varphi c}}k_{r\varphi c}\left(\nabla P_{\varphi z}-{\rho }_{\varphi z}g\right),$

where:

u_φkz is the velocity of phase φ at well k, in sub-geology z,

μ_φc is the viscosity of phase φ, in sub-geology z,

ρ_φc is the density of phase φ, in sub-geology z,

P_φz is the pressure of phase φ, in sub-geology z,

g is the gravity vector,

k is a permeability coefficient, and

k_rφz is a relative permeability coefficient of phase φ, in sub-geology z.

The laws of reservoir physics 6 may also comprise, if corresponding data are available in history data HD, heat transfer laws between a fluid or a gas and a rock given by the following formulas:

Φ∂_t(ρ_wS_w+ρ_gS_g)+∇·(ρ_wu_w+ρ_gu_g)=0

Φ∂_t(ρ_oS_o)+∇·(ρ_ou_o)=0,

S _o +S _w +S _g=1,

∂_t(E_s+E_f)+∇·(U_f−λ(T)∇T)=0, where

T=T(x,y,z) is the temperature and can vary over the reservoir,

• • Φ is the rock porosity,
  • S_φ is the saturation of phase φ,
  • ρ_φ is the density of phase φ,
  • u_φ is the velocity of phase φ,
  • E is the volumetric internal energy (indexes s and f corresponding respectively to solid and fluid phases),

U _f =U _g +U _o +U _w is the enthalpy flux, where U_φ=ρ_φh_φu_φ,

h_φ is the specific enthalpy of phase φ,

h_g=h_w+Λ, Λ being the specific heat of water vaporization,

g, w, o indexes being respectively for gas, water and oil,

λ(T) is the reservoir conductivity coefficient.

Additionally, some variables in the production simulator 2 may be constant or may be detailed by considering them as variables e.g. of local pressure P and the local temperature T. So the laws of reservoir physics 6 may further comprise phase laws given by:

μ_φz=μ_φz(P,T) the viscosity of phase φ at sub-geology z, function of the local pressure P and the local temperature T;

ρ_φz=ρ_φz(P,T) the density of phase φ at sub-geology z, function of the local pressure P and the local temperature T;

k_rφz=k_rφz(P,T) the relative permeability coefficient of phase φ at sub-geology z, function of the local pressure P and the local temperature T.

From said velocities u_φkz, known at least at a perforation in layers of the sub-geology z of each well k, second module 7, implementing well models in the form of laws of well physics 7, determines the quantities produced Q_φktz of phase φ in sub-geology z, at well k at time t, in function of the production parameters PP_ktz applied to each well k, in sub-geology z, at time t, and of the velocity u_φhb of phase φ at layer b at well k, using formulas: Q_φktz=T_kz(PP_ktz,u_φkz) where T_kz is a transfer function defined for each well k and each sub-geology z. Such transfer functions are functions depending, inter alia, on local reservoir pressure, geology, skin effect, well/pump design and completion. They can vary over time, as a result of ageing or damage. They are usually tabulated functions, which are derived from past production data.

The space of solutions candidate for becoming a production simulator 2 of the hydrocarbon field, is the combination of T (well transfer physics) and u (propagation physics of the reservoir). Each such candidate is associated to and defined by a set of parameters α. These parameters α are those defining the candidate simulator that will link input production parameters PP_ktz to the quantities produced Q_φktz. The space of candidate simulators in which they are defined is the result of the upscaling/downscaling process of the reservoir, wells and laws. The simulator, defined by such parameters α will fully characterize the field behavior, well by well, sub-geology by sub-geology, as a response to given PP_ktz. These are valid for a given space of allowed values for input production parameters. Said space of allowed values is typically not greater than the envelope of all past production parameters, with an additional variation of +/−30%.

Based on the history data HD, each production simulator candidate solution is compared against past data. The candidate solution whose set of parameters minimizes the forecast expected risk R(α) is determined by iteration over said sets of parameters α.

In the case of mature fields, the amount m of past production or history data HD is large enough to allow defining a space of solutions S embedding the appropriate physical richness able to capture all key phenomena that (have and) will matter for the future production of the same wells. Such situation allows a trade-off for which the VC-dimension h of space of solutions is suitable with respect to the complexity of available data, therefore minimizing the

$\frac{h}{m}$

ratio and hence the Φ and the R(α) values. The forecast expected risk R(α) is minimized and the forecast provided by the production simulator 2 may be considered reliable.

For each candidate solution, the forecast expected risk R(α) is bounded from above using the Vapnik formula:

$R\left(\alpha \right)\le R_{\mathrm{emp}}\left(\alpha \right)+\Phi \left(\frac{h}{m},\delta \right)$

R(α) is the forecast expected risk associated with parameters α,

R_emp(α) is an empirical risk associated with parameters α, determined by a matching process with history data HD,

δ is a positive number close to zero, 1-δ defining a probability that the inequality holds, and

Φ is a function defined by:

$\Phi \left(\frac{h}{m},\delta \right)=2\sqrt{\frac{2}{m}\left(h\log \left(\frac{2\mathrm{em}}{h}\right)+\log \left(\frac{2}{\delta }\right)\right)},$

where

h is the Vapnik-Chervonenkis dimension of the space of solutions, and

m is the number of independent measures available in history data HD,

e is equal to exp(1).

The best retained production simulator 2 is determined and defined by a set of optimal parameters α_opt that minimizes the forecast expected risk R(α).

Once a production simulator 2 is thus determined and optimized so as to respect the matching condition, the Vapnik condition and to minimize the forecast expected risk R(α), it can be used inside an optimizer 1. The optimizer 1 generates scenarios of production parameters PP and applies the production simulator 2 over these scenarios. Many such runs are then iterated over different sets of production parameters PP. At each iteration, a gain value 5 is calculated derived from the quantity produced Q_φktz forecast by production simulator 2. The resulting gain value 5 may be used to select the next scenario. So doing an optimal set of production parameters PP can be obtained that optimize said gain value 5.

Said optimized gain value 5 can be a net present value NPV or reserves RES of the hydrocarbon field.

Any gain values 5 can be determined by module 4 from quantities produced Q_φktz output by optimizer 1, by taking into account the necessary economic parameters or indexes, as is well known by the one skilled in the art.

Said net present value NPV may e.g. be determined using the formula:

$\mathrm{NPV}=\sum _i\left(\sum _k\sum _zP_{\mathrm{kz}}*\left(1-R_{\mathrm{ik}}\right)\right)*\frac{S_i}{{\left(1+d\right)}^i}-\sum _i\left(\sum _kI_{\mathrm{ik}}+\sum _k{\mathrm{OC}}_{\mathrm{ik}}+\sum _k\sum _z\left({\mathrm{TO}}_i*P_{\mathrm{kz}}+{\mathrm{TL}}_i*L_{\mathrm{kz}}\right)\right)*\frac{1}{{\left(1+d\right)}^i}$

where:

P_kz is oil production (in barrels) for well k and sub-geology z,

R_ik is tax and royalties for well k and year i,

S_i is the oil sale price (per barrel) for year i,

d is the percentage discount rate,

I_ik is investment made on well k during year i,

OC_ik is operating costs for well k during year i,

L_kz is the liquid production (in barrels) for well k and sub-geology z,

TO_i is treatment cost (per barrel of oil), for year i.

TL_i is treatment cost (per barrel of liquid), for year i

An alternate choice is reserves RES value of the hydrocarbon field, defined as the cumulative oil produced over a given time period. Others choices are possible.

Parallel Production Simulator

In the above, embodiments were described for building a reliable production simulator for a mature oil or natural gas field. The production simulator allows to model the behavior of the oil or natural gas field and to forecast produced quantities for a given scenario of production parameters. However, in practice, substantial variations in operation can be envisaged such that billions of billions of scenarios may exist. Traditional simulators, such as meshed-based or behavior-based simulators, iterate timewise on the model and as such require an excessive amount of computation time to run multiple scenarios associated with multiple hypotheses.

In the following, embodiments for implementing the production simulator described above to allow for running multiple scenarios in parallel are provided. The embodiments, as further described below, make use of vectorized computation hardware to provide a parallel production simulator over multiple scenarios. By leveraging the computation power of the vectorized computation hardware, multiple scenarios can be computed in the time that it takes to compute a single scenario in a state of the art simulator. Specifically, the production of the oil field is predicted, well by well, in parallel, for each scenario of the multiple scenarios. In an embodiment, the production parameters of the multiple scenarios are aligned in the memory of the vectorized computation hardware, which allows leveraging the hardware power to compute multiple scenarios in a single pass.

For the purpose of illustration, embodiments will be described with reference to an example oil field 400 illustrated in FIG. 4. However, as would be understood by a person of skilled in the art based on the teachings herein, embodiments are not limited to example oil field 400. As shown in FIG. 4, example oil field 400 includes a plurality of wells 404a-g (W1-W7). According to embodiments, wells 404a-g may include horizontal and vertical wells. In addition, alternatively or additionally, wells 404a-g may include, without limitation, field elements like perforations, multiple strings, or any controlled device that can affect production of the oil field.

According to embodiments, oil field 400 can form a single zone or can be divided into multiple zones. In the embodiment illustrated in FIG. 4, oil field 400 is divided into a plurality of zones 402a-c (Z1-Z3). Each well 404 of the plurality of wells 404a-g belongs to a respective zone 402 of the plurality of zones 402a-c. Specifically, wells 404a-c belong to zone 402a, wells 404d-e belong to zone 402b, and wells 404f-g belong to zone 402c. A zone may defined as an arbitrary group of wells, a physical region of the oil field, an administrative partition of the oil field, a vertical separation of layers, a part of the oil field separated by a fault from other parts, a surface installation, a sub-geology or a local region exhibiting a specific behavior. A zone may be independent or connected to other zones of the oil field. For example, in example oil field 400, zones 402a-c may or may not be connected to each other by having fluid transfer between them. Further, a zone may be divided into a number of subzones, which themselves may also contain subzones.

As noted above, embodiments enable scaling the production simulator described above to run multiple scenarios at the same time. A scenario according to embodiments is defined by a data set that includes production parameters of existing wells and, if any, the locations and production parameters of one or more new wells envisioned by the scenario. The scenario may also include production parameters of existing injectors and, if any, the locations and production parameters of new injectors. According to embodiments, the production parameters may include any measured or calculated data, such as, Tubing Head Pressure (THP), Bottom Hole Flowing Pressure (BHFP), well head pressure, gas lift rate, pump frequency, work-over, change of completion, etc., and for new wells, the drilling time or completion, for example.

Different scenarios can be tested using the production simulator according to embodiments. The scenarios are executed on the same model of the oil field, as generated by the production simulator, but using different production parameter data sets. For example and without limitation, scenarios for changing the current production parameters of the oil field, converting an existing producer to an injector, adding one or more new wells (producers, injectors, etc.), perforating into a new layer, and/or isolating an existing layer can be tested.

For the purpose of illustration, FIG. 5 shows example scenarios S1, S2, and S3 that may be tested for example oil field 400 discussed above with reference to FIG. 4. Scenario S1 describes operating oil field 400 with normal production parameters, for example, corresponding to current production parameters. Scenario S1 may be used as a baseline for evaluating other scenarios. Scenario S2 describes making a conversion 502 at well W7 of the oil field. For example, well W7 may be converted from a producer well to an injector well. Scenario S3 describes the addition of a new well 504 (W8) in the zone Z3 of the oil field.

For each scenario, as further described below, the production simulator predicts the production of the oil field, well by well, and in the case of a multi-layered geology, by layer or group of layers. In an embodiment, for each scenario, the production simulator computes the quantities Q_φktz=f(PP) per phase φ (oil, water, and gas), per well ω and per time step t as a function of production parameters PP corresponding to the scenario. In another embodiment, as discussed above, the computed quantities may be of the form Q_φktz, i.e., the quantities produced of phase φ at well k in sub-geology z at time step t.

In an embodiment, to fully leverage the computation power of the vectorized computation hardware, the simulation data is organized in memory to facilitate the simulation of multiple scenarios in parallel. This is illustrated, for example, in FIG. 6 which shows a memory representation of the scenarios S1, S2, and S3 described above. As shown, the simulation data is organized by data type (e.g., production parameter, equation result, output) with each data type being represented by a single one-dimensional array at each time step t of the simulation. Each well of each scenario is assigned the same index in each of the one-dimensional arrays corresponding to the different data types.

At the beginning of the simulation (e.g., time step t=1), production parameter data sets 608a, 608b, and 608c corresponding respectively to the plurality of scenarios being tested are received. The production parameter data sets are aligned into corresponding one-dimensional arrays 602 formed in memory. For example, as shown in FIG. 6, the production parameter data set of each scenario S1, S2, and S3 may include any measured or calculated data, such as, a THP value and a Choke value for each well in the scenario. The production parameter data sets are aligned into a THP one-dimensional array 602a and a Choke one-dimensional array 602b. The THP/Choke one-dimensional array 602a/b includes the THP/Choke parameter value of each well of each scenario. A well that appears in multiple scenarios, such as well 404a for example, will be associated with multiple entries in the one-dimensional array, each entry corresponding to one of the multiple scenarios. The alignment of production parameter data sets into one-dimensional arrays as described above is performed repeatedly for each time step t of the simulation.

In an embodiment, the one-dimensional array containing the production parameters for each scenario at each time step t can be reduced to only include changes across time steps of said production parameter during the duration of the simulation. In another embodiment, the one-dimensional array containing the production parameters for each scenario is only created for time steps of the simulation that induce a change for the corresponding production parameter. This will reduce the memory taken by the simulation.

In an embodiment, the production parameter data sets are grouped per scenario when aligned in the one-dimensional arrays such that the production parameter data set values of a scenario occupy consecutive cells of the arrays. The starting cell index of each scenario is stored during the simulation.

In another embodiment, when a scenario includes new wells or zones, the same new wells or zones may be added to the other scenarios in order to match the scenario with the highest number of wells. The results of the new wells or zones are ignored for the other scenarios. This approach ensures that the production simulator works the same way for all wells without adding special case processing and branching (if-then-else for instance).

The actual locations of wells in memory may vary. For example, depending on the underlying hardware, wells of the same zone and/or scenario may be packed together for better memory locality or dispatched in different memory locations to enhance memory bandwidth and layout. Alternatively, wells of the same zone and/or scenario may be partitioned without memory locality, i.e., without them being stored physically near each other.

Intermediate/equation results (e.g., BHFP, phase velocities, etc.) and outputs (e.g., produced Oil, Water, and Gas volumes and NPV) of the simulation are organized in a similar fashion as the production parameters in one-dimensional arrays 604 and 606 respectively. Arrays 602, 604, and 606 have all the same size.

In an embodiment, arrays 604 and 606 are generated in time step t′=t+1 based on the arrays 602 associated with time step t. In an embodiment, intermediate/equation results and outputs may be stored for each time step of the simulation.

In an embodiment, each well of each scenario is associated with the same index in each of the one-dimensional arrays 602, 604, and 606 (i.e., production parameters arrays, equation results arrays, and outputs arrays). For example, well 404a of the Scenario S1 is associated with the first cell in all of the one-dimensional arrays 602, 604, and 606. Accordingly, a given well of a given scenario can be readily assigned to a data pipeline (e.g., processing thread or lane) by associating its corresponding index with the data pipeline. In an embodiment, once the data is organized in the memory of the vectorized computation hardware, the production simulator is executed in parallel on all scenarios in memory assuming that the computing power is available.

In an embodiment, illustrated in FIG. 7, the vectorized computation hardware supports Single Instruction, Multiple Threads (SIMT) processing. For example, the vectorized computation hardware may be a General-Purpose Graphics Processing Unit (GPGPU), i.e., a graphics processing unit that performs non-specialized functions typically performed on a Central Processing Unit (CPU). In this embodiment, each well of each scenario is assigned a thread of a plurality of available threads based on the cell index of each well.

Using this approach, the production simulator can be scaled up to the maximum number of threads that can be supported by the underlying hardware. Typically, the number of available threads is higher than or equal to the aggregate number of wells in the tested scenarios. If the aggregate number of wells is greater than the number of available threads, the vectorized computation hardware can be leveraged sequentially, testing as many scenarios as can be fitted in a single run and testing additional scenarios in subsequent runs.

As shown in FIG. 7, in a SIMT architecture, a simulation run begins in step 702, which includes, in the first step of the simulation, loading the data into the memory of the vectorized computation hardware. In an embodiment, step 702 includes receiving a plurality of production parameter data sets corresponding respectively to a plurality of scenarios. As discussed above, the plurality of scenarios each includes one or more zones and a set of wells of the oil field, and the plurality of production parameter data sets each includes one or more production parameter values for each well of the corresponding scenario. The one or more production parameter values may include, without limitation, for example, a Tubing Head Pressure (THP), a Choke value, and any other measured or calculated Production Parameter relevant to the model.

In an embodiment, step 702 may further include assigning each well of each scenario a cell index; and assigning each zone of each scenario a zone cell index. For example, with reference to FIG. 7, the well W₁ of the first scenario is assigned the index 0 as a cell index, and the well W₂ of the first scenario is assigned the index 1 as a cell index. The zone comprising the wells W₁ and W₂ in the first scenario is assigned the index 0 as a zone index. Accordingly, the zone cell index assigned to the zone is identical to the cell index assigned to the well W₁ but different than the cell index assigned to the well W₂. According to embodiments, the cell index assigned to a well may be based on or independent from a zone and/or a scenario comprising the well.

Subsequently, step 702 may further include aligning the plurality of production parameter data sets in at least a first one-dimensional array formed in the memory such that a production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, where the respective cell has an index equal to the cell index assigned to the well. For example, with reference to FIG. 6, the first one-dimensional array may be THP array 602a. Accordingly, in this example, the THP production parameter value for the well W₁ of the first scenario is stored in the cell of THP array 602a having the index 0, which is equal to the cell index assigned to the well W₁ of the first scenario.

Next, thread assignment is performed. Specifically, in an embodiment, the cell index, for each well of each scenario, is associated with a first thread of a plurality of thread; and the zone cell index, for each zone of each scenario, is associated with a second thread of the plurality of threads. The first thread may be identical to or different than the second thread. In the example of FIG. 7, the cell index assigned to the well W₁ of the first scenario is associated with a thread 722 of the plurality of threads. The zone cell index assigned to the zone comprising the well W₁ of the first scenario is also associated with thread 722 in this example, though it may be associated with a different thread of the plurality of threads.

Once the data has been loaded and thread assignment has been performed, an iteration of the simulation run is executed. In an embodiment, an iteration of the simulation run includes the steps 704, 706, 708, 710, 712, 714, 716, and 718 further described below. As would be understood by a person of skill in the art based on the teachings herein, these steps may be performed in a different order than shown in FIG. 7. For example, the order of steps 712, 714, and 716 may be different in another embodiment.

In each iteration, steps 704, 706, 708, 710, 712, 714, 716, and 718 are performed in parallel for each thread of the plurality of threads, i.e., for each well of each scenario. For the purpose of clarity, steps 704, 706, 708, 710, 712, 714, 716, and 718 are described below with respect to a given well of a given scenario.

In an embodiment, step 704 includes computing, using the first thread associated with the cell index assigned to the well, an intermediate result for the well based on the one or more production parameter values for the well. In an embodiment, the intermediate result includes a Bottom Hole Flowing Pressure (BHFP) at the well, computed according to the equation BHFP(w,t)=F_(BHFP,w)(THP(w,t),Choke(w,t),Q(φ,w,t−1)), where F_(BHFP,w) is a function taking as input the production parameters THP and Choke and the quantity produced during the previous time step t−1 at the well w, to forecast the local BHFP of well w at time step t. In an embodiment, the intermediate result for the well is stored in a cell of a one-dimensional array having an index that is equal to the cell index assigned to the well.

Step 706 includes synchronizing the plurality of threads and dispatching the intermediate result, computed in step 704, to a cell of a second one-dimensional array, where the cell of the second one-dimensional array has an index equal to the zone cell index assigned to the zone comprising the well. For example, as shown in FIG. 7, the intermediate result associated with well W₁ (as well as the intermediate results associated with the wells W₂ and W₃ belonging to the same zone Z₁ in the first scenario) is dispatched to a cell of a second one-dimensional array 708 having an index (0) equal to the zone cell index assigned to the zone comprising the well W₁ in the first scenario.

Step 708 includes computing, using the second thread associated with the zone comprising the well, a zone-based intermediate result for the well based on the intermediate result of the well. In an embodiment, the zone-based intermediate result comprises one or more velocities of one of more phases for each well in the zone. In an embodiment, the velocity of a phase φ in zone Z at time step t is given by the equation:

$u_{\varphi t}=-\frac{k}{{\mathrm{\phi \mu }}_b}k_{r\varphi t}\left(\nabla P_{\varphi }-{\rho }_{\varphi }g\right)$

where μ_b is the viscosity of phase φ, ϕ the rock porosity, P_φ is the pressure of phase φ, ρ_φ is the density of phase φ, g is the gravity vector, k is a permeability coefficient, and k_rφt is a relative permeability coefficient of phase φ, at time step t through its dependency on phase saturation. In an embodiment, step 708 further includes the computation of related intermediate zone-based parameters. In an embodiment, the computation in step 708 is done in parallel for all zones of each scenario. As noted above, the second thread may be identical to or different from the first thread associated with the cell index assigned to the well. In another embodiment, the second thread can be independent of any threads associated with any wells of the corresponding zone.

Step 710 includes synchronizing the plurality of threads and dispatching the zone-based intermediate result to a cell of a third one-dimensional array having an index equal to the cell index assigned to the well. For example, for the wells W₁, W₂, and W₃ of the first scenario, the zone-based intermediate result computed on the second thread associated with the zone comprising the wells (which, in this example, is identical to the first thread associated with the well W₁) is dispatched to respective cells (of the third one-dimensional array) having indices that correspond respectively to the cell indices assigned to the wells W1, W₂ and W₃ of the first scenario.

Subsequently, steps 712, 714, and 716 include computing, using the first thread, quantity results for the well based on the zone-based intermediate result computed in step 708 and/or the intermediate result computed in step 704. The quantity results computed in steps 712, 714, and 716 may be computed in parallel in an embodiment. In an embodiment, the quantity results comprise one or more volumes of one or more phases, produced at the well during the simulation time step. In an example embodiment, the quantities at the current time step t for each zone Z are given by the equation:

Q _φwt =T _w(BHFP_wt,u_φt)

where Q_φwt is the quantity produced of phase φ (Oil, Water, Gas) at well w in zone Z at time step t, T_w is a transfer function of well w in zone Z, and u_φwt is the velocity of phase φ at well w in zone Z at time step t.

In an embodiment, the simulation run optionally includes in step 718 calculating a Net Present Value (NPV) for each zone of each scenario based on the quantity results computed in steps 712, 714, and 716. In an embodiment, the NPV is computed as a running sum up to the current time step t according to the equation:

$\mathrm{NPV}=\sum _t\left(\sum _w\sum _zP_{\mathrm{wzt}}*\left(1-R_{\mathrm{tw}}\right)\right)*\frac{S_t}{{\left(1+d\right)}^t}-\sum _t\left(\sum _wI_{\mathrm{tw}}+\sum _w{\mathrm{OC}}_{\mathrm{tw}}+\sum _w\sum _z\left({\mathrm{TO}}_{\mathrm{tw}}*P_{\mathrm{wzt}}+{\mathrm{TL}}_{\mathrm{tw}}*L_{\mathrm{wzt}}\right)\right)*\frac{1}{{\left(1+d\right)}^t}$

where P_wzt is the oil quantity (in barrel) produced at time step t for well w∈Z, R_tw is the tax and royalty at time step t for well w∈Z in %, S_t is the oil sale price (per barrel) at time step t, d is the percentage discount rate, I_tw is the investment made at time step t for well w∈Z, OC_tw is the operating cost at time step t for well w∈Z, TO_tw is the treatment cost (per barrel of oil) at time step t for well w∈Z, TL_nt is the treatment cost (per barrel of liquid) at time step t for well w∈Z, and L_wzt is the liquid quantity (in barrel) produced at time step t for well w∈Z.

At the end of each iteration, a check is made to determine if the simulation run final time step t_final has been reached. If not, the current time step value is incremented by one and the above described process is repeated beginning at step 704. Otherwise, the simulation run is terminated in step 720. The produced quantities (Oil volume, Water volume, Gas volume, NPV) for each well of each scenario over the life of the simulation are output to the user.

In another embodiment, illustrated in FIG. 8, the vectorized computation hardware supports Single Instruction, Multiple Data (SIMD) processing. For example, the vectorized computation hardware may be a multi-core processor capable of performing the same instruction on multiple data points simultaneously. In this embodiment, a group of wells corresponding to one or more scenarios are associated with a first data pipeline having a plurality of lanes, and each well of the group of wells is associated with a respective lane of the plurality of lanes of the first data pipeline. Further, a group of zones corresponding to one or more scenarios are associated with a second data pipeline having a plurality of lanes, and each zone of the group of zones is associated with a respective lane of the plurality of lanes of the second data pipeline.

In an embodiment, each lane of the plurality of lanes of the first or second data pipeline runs on a processor core (or logical core if it uses hyper threading) and makes use of SIMD instructions to parallelize data computation. In other words, vectorization takes place at the instruction level of the CPU.

As shown in FIG. 8, an SIMD simulation run process is similar to the SIMT process described above with reference to FIG. 7. However, where in SIMT a well is assigned to a thread based on their cell index, in SIMD, a group of wells corresponding to one or more scenarios are assigned to a data pipeline having a plurality of lanes. Within the data pipeline, multiple processing lanes are provided by virtue of the SIMD processing. Specifically, each individual well of the group of wells is assigned a SIMD lane within the data pipeline based on their cell index. For example, in FIG. 8, a data pipeline 822 is assigned to wells W₁, W₂, and W₃ of the first scenario, and each well is associated with a SIMD lane within data pipeline 822.

In an embodiment, the SIMD simulation run may include steps 802, 804, 806, 808, 810, 812, 814, 816, 818, and 820.

Step 802 is similar to step 702 of the SIMT process described above. In an embodiment, as in step 702, step 802 includes receiving a plurality of production parameter data sets corresponding respectively to a plurality of scenarios. Step 802 further includes assigning each well of each scenario a cell index; and assigning each zone of each scenario a zone cell index. The cell index assigned to a well may be identical to or different than a zone cell index assigned to a zone comprising the well. According to embodiments, the cell index assigned to a well may be based on or independent from a zone and/or a scenario comprising the well. Additionally, step 802 includes aligning the plurality of production parameter data sets in at least a first one-dimensional array formed in the memory such that a production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, where the respective cell has an index equal to the cell index assigned to the well.

In an embodiment, step 802 further includes associating groups of wells and groups of zones with data pipelines. Specifically, in an embodiment, step 802 includes associating contiguous cell indices, corresponding to a group of wells corresponding to one or more scenarios, with a first data pipeline having a plurality of lanes, where each cell index of the contiguous cell indices is associated with a respective lane of the plurality of lanes of the first data pipeline. For example, with reference to FIG. 8, the contiguous cell indices assigned to wells W₁, W₂, and W₃ of the first scenario are associated with a first data pipeline 822 having a plurality of lanes. Each cell index of the contiguous cell indices is associated with a respective lane of first data pipeline 822.

Step 802 also includes associating contiguous zone cell indices, corresponding to a group of zones corresponding to one or more scenarios, with a second data pipeline having a plurality of lanes, where each zone cell index of the contiguous zone cell indices is associated with a respective lane of the plurality of lanes of the second data pipeline.

In an embodiment, the first or second data pipeline is assigned to a thread of a plurality of threads.

Steps 804 through 816 are computed in parallel for each well of the group of wells. For the purpose of clarity, the steps are describe below with respect to a given well of a given group of wells. As would be understood by a person of skill in the art based on the teachings herein, the simulation may involve one or more groups of wells and one or more groups of zones.

Step 804 is similar to step 704 described above except that the computation of the intermediate result is performed using the lane of the first data pipeline associated with the well, instead of a thread associated with the well.

Step 806 is similar to step 706 of the SIMT process, in that it includes dispatching the intermediate result to a cell of a second one-dimensional array having an index equal to the zone cell index assigned to the zone comprising the well. In an embodiment, the synchronization in step 806 may be optional if the size of a SIMD register is higher than or equal to the number of wells in the simulated scenarios.

Step 808 is similar to step 708 described above, except that the computation of the zone-based intermediate result is performed using the respective lane of the second data pipeline corresponding to the zone comprising the well, instead of a thread associated with the zone. As noted above, the second data pipeline may be identical to or different than the first data pipeline.

Step 810 is identical to step 710 of the SIMT process, except that the synchronization may be optional in the SIMD process if the size of a SIMD register is higher than or equal to the number of wells in the simulated scenarios.

Steps 812, 814, 816, and 818 are similar to steps 712, 714, 716, and 718 respectively, except that the computations are performed on the respective lane of the first data pipeline associated with the well, instead of a thread associated with the well. Finally, step 820 is identical to step 720 described above with respect to the SIMT process.

FIG. 9 illustrates an example process 900 for simulating the production of an oil field according to an embodiment. Example process 900 is provided for the purpose of illustration only and is not limiting of embodiments. Process 900 may be performed by a vectorized computation hardware such as an SIMT or SIMD processor.

As shown in FIG. 9, process 900 begins in step 902, which includes receiving a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, where each of the plurality of scenarios includes one or more zones and a set of wells of the oil field and the plurality of production parameter data sets each includes a production parameter value for each well of the corresponding scenario. The production parameter value may include, without limitation, a THP, a Choke value, any other relevant measured or calculated production parameters, and/or a previous time step quantity result associated with the well.

In step 904, process 900 includes assigning each well of each scenario a cell index.

In step 906, process 900 includes assigning each zone of each scenario a zone cell index.

Subsequently, in step 908, process 900 includes aligning the plurality of production parameter data sets in a first one-dimensional array formed in the memory such that the production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, the respective cell having an index equal to the cell index assigned to the well. In an embodiment, the first one-dimensional array is formed in a memory of the vectorized computation hardware.

Next, in step 910, process 900 includes associating the cell index, for each well of each scenario with a first dedicated computation unit of a plurality of computation units. In embodiments, the computation unit may be a thread of a plurality of threads in the case of a SIMT architecture or a lane of a multi-lane data pipeline in the case of a SIMD architecture.

Process 900 terminates in step 912, which includes associating the zone cell index, for each zone of each scenario, with a second dedicated computation unit of the plurality of computation units.

In an embodiment, process 900 may further include in parallel, for each well of each scenario, computing an intermediate result for the well based on the production parameter value for the well; computing a zone-based intermediate result for the well based on previous results computed for each well of the zone comprising the well; and computing a quantity result for the well based on the zone-based intermediate result. The zone-based intermediate result may comprise one or more velocities of one or more phases at the well. The quantity result may comprise one or more volumes of one or more phases, produced at the well during a simulation time step.

Example Computer System

Various embodiments can be implemented, for example, using one or more computer systems, such as computer system 1000 shown in FIG. 10.

Computer system 1000 includes one or more processors (also called central processing units, or CPUs), such as a processor 1004. Processor 1004 is connected to a communication infrastructure or bus 1006.

One or more processors 1004 may each be a graphics processing unit (GPU). In an embodiment, a GPU is a processor that is a specialized electronic circuit designed to rapidly process mathematically intensive applications on electronic devices. The GPU may have a highly parallel structure that is efficient for parallel processing of large blocks of data, such as mathematically intensive data common to computer graphics applications, images and videos. In an embodiment, each of the one or more processors 1004 may implement an SIMT or an SIMD architecture.

Computer system 1000 also includes user input/output device(s) 1016, such as monitors, keyboards, pointing devices, etc., which communicate with communication infrastructure 1006 through user input/output interface(s) 1002.

Computer system 1000 also includes a main or primary memory 1008, such as random access memory (RAM). Main memory 1008 may include one or more levels of cache. Main memory 1008 has stored therein control logic (i.e., computer software) and/or data that when executed by the one or more processors 1008 perform the above described embodiments, including the processes described in FIGS. 7-9.

Computer system 1000 may also include one or more secondary storage devices or memory 1010. Secondary memory 1010 may include, for example, a hard disk drive 1012 and/or a removable storage device or drive 1014. Removable storage drive 1014 may be a floppy disk drive, a magnetic tape drive, a compact disk drive, an optical storage device, a tape backup device, and/or any other storage device/drive.

Removable storage drive 1014 may interact with a removable storage unit 1018. Removable storage unit 1018 includes a computer usable or readable storage device having stored thereon computer software (control logic) and/or data. Removable storage unit 1018 may be a floppy disk, magnetic tape, compact disk, DVD, optical storage disk, and/any other computer data storage device. Removable storage drive 1014 reads from and/or writes to removable storage unit 1018 in a well-known manner.

According to an exemplary embodiment, secondary memory 1010 may include other means, instrumentalities or other approaches for allowing computer programs and/or other instructions and/or data to be accessed by computer system 1000. Such means, instrumentalities or other approaches may include, for example, a removable storage unit 1022 and an interface 1020. Examples of the removable storage unit 1022 and the interface 1020 may include a program cartridge and cartridge interface (such as that found in video game devices), a removable memory chip (such as an EPROM or PROM) and associated socket, a memory stick and USB port, a memory card and associated memory card slot, and/or any other removable storage unit and associated interface.

Computer system 1000 may further include a communications or network interface 1024. Communications interface 1024 enables computer system 1000 to communicate and interact with any combination of remote devices, remote networks, remote entities, etc. (individually and collectively referenced by reference number 1028). For example, communications interface 1024 may allow computer system 1000 to communicate with remote devices 1028 over a communications path 1026, which may be wired and/or wireless, and which may include any combination of LANs, WANs, the Internet, etc. Control logic and/or data may be transmitted to and from computer system 1000 via communication path 1026.

In an embodiment, a tangible apparatus or article of manufacture comprising a tangible computer useable or readable medium having control logic (software) stored thereon is also referred to herein as a computer program product or program storage device. This includes, but is not limited to, computer system 1000, main memory 1008, secondary memory 1010, and removable storage units 1018 and 1022, as well as tangible articles of manufacture embodying any combination of the foregoing. Such control logic, when executed by one or more data processing devices (such as computer system 1000), causes such data processing devices to operate as described herein.

Based on the teachings contained in this disclosure, it will be apparent to persons skilled in the relevant art(s) how to make and use the invention using data processing devices, computer systems and/or computer architectures other than that shown in FIG. 10.

CONCLUSION

The foregoing description of the specific embodiments will so fully reveal the general nature of the disclosure that others can, by applying knowledge within the skill of the art, readily modify and/or adapt for various applications such specific embodiments, without undue experimentation, without departing from the general concept of the present disclosure. Therefore, such adaptations and modifications are intended to be within the meaning and range of equivalents of the disclosed embodiments, based on the teaching and guidance presented herein. It is to be understood that the phraseology or terminology herein is for the purpose of description and not of limitation, such that the terminology or phraseology of the present specification is to be interpreted by the skilled artisan in light of the teachings and guidance.

The breadth and scope of embodiments of the present disclosure should not be limited by any of the above-described exemplary embodiments, but should be defined only in accordance with the following claims and their equivalents.

For the purpose of illustration and not limitation, various aspects of the present disclosure are set out in the following variants.

In a first variant, the present disclosure provides a parallel production simulator for simulating the production of an oil field, comprising: a processor; and a memory coupled to the processor. The processor is configured to: receive a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, the plurality of scenarios each including one or more zones and a set of wells of the oil field, and the plurality of production parameter data sets each including a production parameter value for each well of the corresponding scenario; assign each well of each scenario a cell index; assign each zone of each scenario a zone cell index; align the plurality of production parameter data sets in a first one-dimensional array formed in the memory such that the production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, the respective cell having an index equal to the cell index assigned to the well; associate the cell index, for each well of each scenario, with a first thread of a plurality of threads; associate the zone cell index, for each zone of each scenario, with a second thread of the plurality of threads; and in parallel, for each well of each scenario, using the first thread associated with the corresponding cell index, an intermediate result for the well based on the production parameter value for the well; dispatch the intermediate result to a cell of a second one-dimensional array, the cell of the second one-dimensional array having an index equal to the zone cell index assigned to the zone comprising the well; compute, using the second thread associated with the zone comprising the well, a zone-based intermediate result for the well based on the intermediate result; dispatch the zone-based intermediate result to a cell of a third one-dimensional array, the cell of the third one-dimensional array having an index equal to the cell index assigned to the well; and compute, using the first thread, a quantity result for the well based on the zone-based intermediate result.

In a second variant, which is a variant of the first variant, the production parameter value includes a Tubing Head Pressure (THP) or a Choke value.

In a third variant, which is a variant of any of the preceding variants, the intermediate result comprises a Bottom Hole Flowing Pressure (BHFP) at the well.

In a fourth variant, which is a variant of any of the preceding variants, the zone-based intermediate result comprises one or more velocities of one or more phases at the well.

In a fifth variant, which is a variant of any of the preceding variants, the quantity result comprises one or more volumes of one or more phases, produced at the well during a simulation time step.

In a sixth variant, which is a variant of any of the preceding variants, the zone cell index for a zone comprising the well is identical or different than the cell index assigned to the well.

In a seventh variant, which is a variant of any of the preceding variants, the first thread is identical to or different than the second thread.

In an eighth variant, which is a variant of any of the preceding variants, the processor is further configured to synchronize the plurality of threads before dispatching the intermediate result to the cell of the second one-dimensional array.

In a ninth variant, which is a variant of any of the preceding variants, the processor is further configured to synchronize the plurality of threads before dispatching the zone-based intermediate result to the cell of the third one-dimensional array.

In a tenth variant, which is a variant of any of the preceding variants, the cell index assigned to the well is based on or independent from a zone and/or a scenario of the well.

In an eleventh variant, which is a variant of any of the preceding variants, the processor comprises a Single Instruction, Multiple Thread (SIMT) processor.

In a twelfth variant, which is a variant of any of the preceding variants, the processor is further configured, in parallel for each well of each scenario, to compute, using the first thread, a NPV result for the well based on the quantity result of the well.

In a thirteenth variant, the present disclosure provides a parallel production simulator for simulating the production of an oil field, comprising a processor; and a memory coupled to the processor. The processor is configured to: receive a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, the plurality of scenarios each including one or more zones and a set of wells of the oil field, and the plurality of production parameter data sets each including a production parameter value for each well of the corresponding scenario; assign each well of each scenario a cell index; assign each zone of each scenario a zone cell index; align the plurality of production parameter data sets in a first one-dimensional array formed in the memory such that the production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, the respective cell having an index equal to the cell index assigned to the well; associate contiguous cell indices, corresponding to a group of wells corresponding to one or more scenarios, with a first data pipeline having a plurality of lanes, wherein each cell index of the contiguous cell indices is associated with a respective lane of the plurality of lanes of the first data pipeline; associate contiguous zone cell indices, corresponding to a group of zones corresponding to one or more scenarios, with a second data pipeline having a plurality of lanes, wherein each zone cell index of the contiguous zone cell indices is associated with a respective lane of the plurality of lanes of the second data pipeline; and in parallel, for each well of the group of wells, compute, using respective lane of the first data pipeline, an intermediate result for the well based on the production parameter value for the well; dispatch the intermediate result to a cell of a second one-dimensional array, the cell of the second one-dimensional array having an index equal to the zone cell index assigned to the zone comprising the well; compute, using the respective lane of the second data pipeline corresponding to the zone comprising the well, a zone-based intermediate result for the well based on the intermediate result; dispatch the zone-based intermediate result to a cell of a third one-dimensional array, the cell of the third one-dimensional array having an index equal to the cell index assigned to the well; and compute, using the respective lane of the first data pipeline, a quantity result for the well based on the zone-based intermediate result.

In a fourteenth variant, which is a variant of the thirteenth variant, the production parameter value includes a Tubing Head Pressure (THP) or a Choke value.

In a fifteenth variant, which is a variant of the thirteenth or fourteenth variant, the intermediate result comprises a Bottom Hole Flowing Pressure (BHFP) at the well.

In a sixteenth variant, which is a variant of any of the thirteenth to fifteenth variants, the zone-based intermediate result comprises one or more velocities of one of more phases at the well.

In a seventeenth variant, which is a variant of any of the thirteenth to sixteenth variants, the quantity result comprises one or more volumes of one or more phases, produced at the well during a simulation time step.

In an eighteenth variant, which is a variant of any of the thirteenth to seventeenth variants, the zone cell index for a zone comprising the well is identical or different than the cell index assigned to the well.

In a nineteenth variant, which is a variant of any of the thirteenth to eighteenth variants, the first data pipeline is identical to or different than the second data pipeline.

In a twentieth variant, which is a variant of any of the thirteenth to nineteenth variants, the processor comprises a Single Instruction, Multiple Data (SIMD) processor.

In a twenty-first variant, which is a variant of any of the thirteenth to twentieth variants, the processor is further configured, in parallel for each well of the group of wells, to compute, using the respective lane of the first data pipeline, a NPV result for the well based on the quantity result of the well.

In a twenty-second variant, which is a variant of any of the thirteenth to twenty-first variants, the first or second data pipeline is assigned to a thread of a plurality of threads.

In a twenty-third variant, the present disclosure provides a method for simulating the production of an oil field, comprising: receiving a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, the plurality of scenarios each including one or more zones and a set of wells of the oil field, and the plurality of production parameter data sets each including a production parameter value for each well of the corresponding scenario; assigning each well of each scenario a cell index; assigning each zone of each scenario a zone cell index; aligning the plurality of production parameter data sets in a first one-dimensional array such that the production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, the respective cell having an index equal to the cell index assigned to the well; associating the cell index, for each well of each scenario, with a first dedicated computation unit of a plurality of computation units, wherein a computation unit of the plurality of computation units is a thread of a plurality of threads or a lane of a data pipeline; associating the zone cell index, for each zone of each scenario, with a second dedicated computation unit of the plurality of computation units; and in parallel, for each well of each scenario, computing an intermediate result for the well based on the production parameter value for the well; computing a zone-based intermediate result for the well based on previous results computed for each well of the zone; and computing a quantity result for the well based on the zone-based intermediate result.

The twenty-third variant can be combined with features of any of the second to twelfth and fourteenth to twenty-second variants described above as would be understood by a person skilled in the art based on the teachings herein.

Claims

1. A parallel production simulator for simulating the production of an oil field, comprising: a processor; anda memory coupled to the processor,wherein the processor is configured to:receive a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, the plurality of scenarios each including one or more zones and comprising a set of wells of the oil field, and the plurality of production parameter data sets each including a production parameter value for each well of the corresponding scenario;assign each well of each scenario a cell index;assign each zone of each scenario a zone cell index;align the plurality of production parameter data sets in a first one-dimensional array formed in the memory such that the production parameter value for each well of each scenario is stored in a respective cell of the first one-dimensional array, the respective cell having an index equal to the cell index assigned to the well;associate the cell index, for each well of each scenario, with a first thread of a plurality of threads;associate the zone cell index, for each zone of each scenario, with a second thread of the plurality of threads; andin parallel, for each well of each scenario, compute, using the first thread associated with the corresponding cell index, an intermediate result for the well based on the production parameter value for the well;dispatch the intermediate result to a cell of a second one-dimensional array, the cell of the second one-dimensional array having an index equal to the zone cell index assigned to the zone comprising the well;compute, using the second thread associated with the zone comprising the well, a zone-based intermediate result for the well based on the intermediate result;dispatch the zone-based intermediate result to a cell of a third one-dimensional array, the cell of the third one-dimensional array having an index equal to the cell index assigned to the well; andcompute, using the first thread, a quantity result for the well based on the zone-based intermediate result.

2. The parallel production simulator of claim 1, wherein the production parameter value includes a Tubing Head Pressure (THP) or a Choke value.

3. The parallel production simulator of claim 1, wherein the intermediate result comprises a Bottom Hole Flowing Pressure (BHFP) at the well.

4. The parallel production simulator of claim 1, wherein the zone-based intermediate result comprises one or more velocities of one of more phases at the well.

5. The parallel production simulator of claim 1, wherein the quantity result comprises one or more volumes of one or more phases, produced at the well during a simulation time step.

6. The parallel production simulator of claim 1, wherein the zone cell index for a zone comprising the well is identical or different than the cell index assigned to the well.

7. The parallel production simulator of claim 1, wherein the first thread is identical to or different than the second thread.

8. The parallel production simulator of claim 1, wherein the processor is further configured to synchronize the plurality of threads before dispatching the intermediate result to the cell of the second one-dimensional array.

9. The parallel production simulator of claim 1, wherein the processor is further configured to synchronize the plurality of threads before dispatching the zone-based intermediate result to the cell of the third one-dimensional array.

10. The parallel production simulator of claim 1, wherein the cell index assigned to the well is based on or independent from a zone and/or a scenario of the well.

11. The parallel production simulator of claim 1, wherein the processor comprises a Single Instruction, Multiple Thread (SIMT) processor.

12. The parallel production simulator of claim 1, wherein the processor is further configured, in parallel for each well of each scenario, to compute, using the first thread, a Net Present Value (NPV) result for the well based on the quantity result of the well.

13. A parallel production simulator for simulating the production of an oil field, comprising: a processor; anda memory coupled to the processor,wherein the processor is configured to:receive a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, the plurality of scenarios each including one or more zones and comprising a set of wells of the oil field, and the plurality of production parameter data sets each including a production parameter value for each well of the corresponding scenario;assign each well of each scenario a cell index;assign each zone of each scenario a zone cell index;align the plurality of production parameter data sets in a first one-dimensional array formed in the memory such that the production parameter value for each well of each scenario is stored in a respective cell of the one-dimensional array, the respective cell having an index equal to the cell index assigned to the well;associate contiguous cell indices, corresponding to a group of wells corresponding to one or more scenarios, with a first data pipeline having a plurality of lanes, wherein each cell index of the contiguous cell indices is associated with a respective lane of the plurality of lanes of the first data pipeline;associate contiguous zone cell indices, corresponding to a group of zones corresponding to one or more scenarios, with a second data pipeline having a plurality of lanes, wherein each zone cell index of the contiguous zone cell indices is associated with a respective lane of the plurality of lanes of the second data pipeline; andin parallel, for each well of the group of wells, compute, using the respective lane of the first data pipeline, an intermediate result for the well based on the production parameter value for the well;dispatch the intermediate result to a cell of a second one-dimensional array, the cell of the second one-dimensional array having an index equal to the zone cell index assigned to the zone comprising the well;compute, using the respective lane of the second data pipeline corresponding to the zone comprising the well, a zone-based intermediate result for the well based on the intermediate result;dispatch the zone-based intermediate result to a cell of a third one-dimensional array, the cell of the third one-dimensional array having an index equal to the cell index assigned to the well; andcompute, using the respective lane of the first data pipeline, a quantity result for the well based on the zone-based intermediate result.

14. The parallel production simulator of claim 13, wherein the production parameter value includes a Tubing Head Pressure (THP) or a Choke value.

15. The parallel production simulator of claim 13, wherein the intermediate result comprises a Bottom Hole Flowing Pressure (BHFP) at the well.

16. The parallel production simulator of claim 13, wherein the zone-based intermediate result comprises one or more velocities of one of more phases at the well.

17. The parallel production simulator of claim 13, wherein the quantity result comprises one or more volumes of one or more phases, produced at the well during a simulation time step.

18. The parallel production simulator of claim 13, wherein the zone cell index for a zone comprising the well is identical or different than the cell index assigned to the well.

19. The parallel production simulator of claim 13, wherein the first data pipeline is identical to or different than the second data pipeline.

20. The parallel production simulator of claim 13, wherein the processor comprises a Single Instruction, Multiple Data (SIMD) processor.

21. The parallel production simulator of claim 13, wherein the processor is further configured, in parallel for each well of the group of wells, to compute, using the respective lane of the first data pipeline, a Net Present Value (NPV) result for the well based on the quantity result of the well.

22. The parallel production simulator of claim 13, wherein the first or second data pipeline is assigned to a thread of a plurality of threads.

23. A method for simulating the production of an oil field, comprising: receiving a plurality of production parameter data sets corresponding respectively to a plurality of scenarios, the plurality of scenarios each including one or more zones and comprising a set of wells of the oil field, and the plurality of production parameter data sets each including a production parameter value for each well of the corresponding scenario;assigning each well of each scenario a cell index;assigning each zone of each scenario a zone cell index;aligning the plurality of production parameter data sets in a first one-dimensional array formed in the memory such that the production parameter value for each well of each scenario is stored in a respective cell of the one-dimensional array, the respective cell having an index equal to the cell index assigned to the well;associating the cell index, for each well of each scenario, with a first dedicated computation unit of a plurality of computation units, wherein a computation unit of the plurality of computation units is a thread of a plurality of threads or a lane of a data pipeline;associating the zone cell index, for each zone of each scenario, with a second dedicated computation unit of the plurality of computation units; andin parallel, for each well of each scenario, computing an intermediate result for the well based on the production parameter value for the well;computing a zone-based intermediate result for the well based on previous results computed for each well of the zone; andcomputing a quantity result for the well based on the zone-based intermediate result.