Next Article in Journal
Review of Cybersecurity Analysis in Smart Distribution Systems and Future Directions for Using Unsupervised Learning Methods for Cyber Detection
Next Article in Special Issue
A Coupled Poro-Elastic Fluid Flow Simulator for Naturally Fractured Reservoirs
Previous Article in Journal
Self-Healing Concrete: Concepts, Energy Saving and Sustainability
Previous Article in Special Issue
A Fully Coupled Hydro-Mechanical Approach for Multi-Fracture Propagation Simulations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Testing the INSIM-FT Proxy Simulation Method

1
Petroleum Engineering, Skolkovo Institute of Science and Technologies, 121205 Moscow, Russia
2
Novosibirsk R&D Center LLC, 630090 Novosibirsk, Russia
3
Department of Physics, Novosibirsk State University, 630090 Novosibirsk, Russia
*
Author to whom correspondence should be addressed.
Energies 2023, 16(4), 1648; https://doi.org/10.3390/en16041648
Submission received: 26 December 2022 / Revised: 19 January 2023 / Accepted: 30 January 2023 / Published: 7 February 2023
(This article belongs to the Special Issue Advances in Reservoir Simulation)

Abstract

:
This paper describes testing of the INSIM-FT proxy simulation method (interwell-numerical-simulation model improved with front-tracking method) to assess the dependencies between production and injection wells, as well as to assess the forecast of oil/liquid production by wells depending on their operation parameters. The paper proposes the approach of taking into account the influence of various production enhancement operations. The method was tested on a synthetic hydrodynamic model and on a sector of a real field. The results show a good match between historical data and simulation results and indicate significant computational efficiency compared to classical reservoir simulators.

1. Introduction

Waterflooding of oil fields is used to displace oil with water from the formations and maintain the formation pressure at a given level by injecting water. Water injection through injection wells is the main and most common way to maintain reservoir pressure. In Russia over 90% of oil fields are developed using this technology [1].
Mature oil and gas fields are characterized by a rapid increase in the water cut of the produced liquid and a decrease in the oil flow rate. This stage is associated with complications of the production process. The sweep efficiency and the oil recovery factor can be improved by application of efficient waterflood systems, primarily, focal ones as well as by changing the filtration flows [2]. Moreover, with the increase in the share of hard-to-recover reserves under development, the task of increasing oil production rates becomes relevant.
With the development of the oil industry, geological and hydrodynamic simulation of hydrocarbon recovery has become one of the main tools for both predicting future performance of the reservoir in either the long-term or short-term period and selecting a reasonable development strategy for oil and gas fields. Moreover, it is an essential tool for the decision-making process in field development and petroleum production optimization where different operating conditions and actions are tested to increase productivity from existing wells or well patterns [3].
Full-scale hydrodynamic models are usually built on poor, averaged, and roughened geological models. Therefore, the insufficient accuracy of calculations performed on such models is associated with a high degree of uncertainty of initial information, especially of reservoir properties in the interwell space [4]. Moreover, to assess the efficiency of development methods, it is necessary to be able to carry out operational multivariate calculations. The use of three-dimensional hydrodynamic models in such cases is inefficient due to the following reasons [5]:
  • Long duration of the simulation;
  • Need to use large computational resources;
  • Complexity of modifying the model when making adjustments or changing various parameters.
The use of machine-learning algorithms or simplified models (proxy models) based on material balance methods and various analytical dependencies is practical, making it possible to account for the most important factors affecting the calculated properties. Such models are less time-consuming, adapt faster to actual data, and allow for quick calculations while maintaining the required accuracy [6]. The advantage of proxy models is historical data are used as the input, namely, production and injection data.
It is obvious that in the development of hydrocarbon fields, interconnection and mutual influence of all wells must be taken into account. Under reservoir pressure maintenance conditions, understanding the interaction between production and injection wells, assessing the degree of injection impact on oil production, and plotting the dependencies of recovery on water injection plays a key role in selecting the most optimal waterflood strategy. The use of proxy modeling makes it possible to solve the problems of accounting for the relationship and mutual influence of all wells.
Thus, a promising area of simulation of hydrocarbon recovery is the improvement of approaches to the creation and implementation of proxy models, as well as the development of methods for automated optimization of the oil field development system.
One of the important stages in proxy simulation is history matching, which solves the inverse task of hydrodynamics, i.e., the main reservoir properties (porosity, permeability, net pay thickness, heterogeneity, etc.) are adjusted to fit the actual field data with certain accuracy [7]. These algorithms are used in predicting well performance and in redistributing the volume of injected fluid between the injection wells with minimal costs and high payback.
Currently, one of the most common proxy modeling approaches are based on reduced-order models, data-driven models, physics-based models, and flow-network models. These models provide for solving problems of waterflooding optimization, identifying the efficiency of injection and production, assessing the interaction of wells, forecast fluid, and oil production in a specific period of time [8].
Reduced-order models are used to replace large reservoir models by reducing high-fidelity simulator models or by creating new reduced-order, data-driven models [9]. There are different model reduction techniques in the literature [10,11,12,13,14]. The main disadvantage of such models is that geologic data is required.
As opposed to reduced-order models, data-driven models require no prior knowledge of petrophysical properties or other specific geological information. Recurrent neural networks are common data-driven models used in the petroleum industry [15,16]. Because these are data-driven models, the accuracy of these models solely depends on the quality of the data used for training and can be affected by data noise. Moreover, these models are difficult to train and have pitfalls, such as overtraining, extrapolation, or lack of validation.
One example of the physics-based approach is the capacitance-resistance model (CRM), which uses production and injection rate data and bottomhole pressure to match the model against a particular reservoir [17,18,19]. CRM is analogous to electrical circuits where the compressibility and transmissibility, respectively, are analogous to capacitance and resistance. It is based on a model of nonlinear signal processing in which injection rates are treated as input signals and production rates are output signals [20]. The disadvantages of the CRM method are the productivity of each production well is considered constant throughout the entire history of production, which is not quite correct physically in the conditions of multiphase filtration. Moreover, the saturation distribution in this model is not calculated, and the oil production rate is adjusted indirectly through the displacement characteristics [17].
Flow-network models represent reservoir flow by a coupled-network model in which each pair of wells is connected with a one-dimensional (1D) finite-difference reservoir-simulation model. Each reservoir is defined by two parameters: absolute permeability and pore volume. The advantages of this approach are the construction of this model does not require a high-fidelity geological data and that it follows the physical laws of the realistic multiphase system. In this approach, relative permeabilities are assumed to be known. This model uses the true relative permeability curves, which is physically correct but requires a priori knowledge of the relative permeability curves [21].
This paper focuses on the interwell-numerical-simulation model (INSIM-FT), which is somewhat similar to the CRM and flow-network model. In the INSIM assumptions, a reservoir is viewed as a series of units connecting well pairs, but instead of discretizing those connections as in the flow-network model, INSIM only defines a pair of parameters for each connection. That is a significant reduction in the number of parameters compared to a set of 1D finite-difference reservoir models. The model uses the correct front-tracking procedure to calculate water saturation [22]. A more detailed description and methods used in the work are described in the next chapters.
This article presents software implementation of the method [4] and modification of the method to account for the influence of geological and technical measures, and on this basis, conducts a series of numerical experiments to assess the accuracy, convergence, and sensitivity of the method to parameterization of the historical performance of the hydrocarbon field.

2. Model Description

2.1. INSIM-FT: Theoretical Background

The interwell-numerical-simulation model (INSIM) is used as a tool to calculate the approximation of the well production rate under waterflooding [23]. In INSIM, the reservoir is characterized as a rough model consisting of several interwell control units (flow tubes), where each unit has two specific properties: conductivity, Tij, and control pore volume, Vpij (Figure 1).
By solving the mass balance equation and front motion equations for each of the units, it is possible to obtain the velocities and saturation of interwell fluids for further prediction of the production rate. INSIM is used to adapt the model based on the available data to estimate parameters and to determine interwell correlation and geological characteristics. INSIM has the following advantages:
  • The model parameters estimated based on historical data provide a relative characteristic of the reservoir properties between the wells. The model can handle changes in flow direction caused by well flow rate changes, including well shut-in or conversion of production wells to injection wells;
  • INSIM is able to calculate the oil and water flow rate and the adapted water cut data;
  • It can be used to optimize waterflooding but with considerably less computational effort.
In the described model, only a two-phase flow rate (liquid oil and water) is considered, and the conductivity, Tij, is set as the average total permeability between the i-th and j-th wells. INSIM-FT solves the material balance equation for the j-th well (without accounting for capillary pressure and gravity):
j = 1 n w T i , j t · P j t P i t q i t = c i t · V p , i · d P i d t
where nw is the total number of wells; qi(t) is the flow rate of the i-th well at time, t, and is positive for injection and negative for production; ci(t) is full compressibility, Vp,i. If well j is not interlinked with well i, then Ti,j = 0 in the equation.
The physical meaning of the above-mentioned equation is that the change in the control pore volume due to compressibility is equal to the difference between the injection into the volume, ij, and the production from the volume, i. The equation is a combined pressure equation. Approximation of the equation by the implicit finite difference scheme used in reservoir simulation gives the following expression:
P i n P i n 1 = Δ t n c t , i n V p , i n j = 1 n w T i , j n P i n + j = 1 n w T i , j n P j n + q i n
for i = 1, 2, …, n w ; where Δ t n = t n t n 1 and t 0 = 0 . Throughput (capacity) ( T i , j n ) well drainage pore volumes ( V p , i n -th) and compressibility ( c t , i n -th) may change with time.
In INSIM, similar to the pressure equation in reservoir simulation, non-linear terms are estimated at the previous time level, that is, terms that depend on the pressure and water saturation estimated at the moment in time, t n 1 , instead of t n . Accordingly, in Equation (2) it is necessary to use:
T i , j n 1 = α k i , j A i , j λ t , i , j n 1 L i j = T i , j 0 λ t , i , j n 1 λ t , i , j 0
V p , i n 1 = V p , i 0 1 + c r p i n 1 p i 0
c t , i n 1 = S o , i n 1 c o + S w , i n 1 c w + c r
instead of T i , j n ; V p , i n ; c t , i n , respectively. Here, S o , i and S w , i are the corresponding volumes of oil and water saturation of the well i ; c o , c w , and c r , respectively, represent the compressibility of oil, water, and rock; λ t , i , j is the total mobility, which is calculated by upstream weighting. If p i n > p i n 1 , then λ t , i , j is replaced by the overall fluidity of the well i , λ t , i . It is assumed that the viscosity of oil and water is constant. Weighing with an upward flow means p i n < p i n 1 :
λ t , i , j n 1 = λ t , i n 1 = k r o S w , i n 1 μ o + k r w S w , i n 1 μ w
Otherwise:
λ t , i , j n 1 = λ t , i n 1 = k r o S w , i n 1 μ o + k r w S w , i n 1 μ w
where μ o and μ w are the viscosities of oil and water, respectively. As V p , i , j n 1 changes, we use the following equation:
V p , i , j n 1 = V p , i , j 0 1 + c r 0.5 p i n 1 + p j n 1 p i 0
We designate:
E i n = Δ t n c t , i n 1 V p , i n 1
G i n = E i n j = 1 n w T i , j n 1
M i n = Δ t n q i n 1 c t , i n 1 V p , i n 1
Then, at each step t, the following system of equations is solved:
G 1 n + 1 E 1 n T 1 , 2 n E 1 n T 1 , n w n E 2 n T 2 , n w n G 2 n + 1 E 1 n T 1 , n w n E n w n T n w , 1 n E 1 n T 1 , n w n G n w n + 1 p 1 n p 2 n p n w n = p 1 n 1 p 2 n 1 p n w n 1 + M 1 n M 2 n M n w n
where M i n is positive for the injection well and negative for the production well.
The system of Equation (12) can be solved for p i n 1 on the assumption the necessary saturations at t n 1 are known. After finding p i n , it is necessary to solve the ‘fractional flow’ equations to obtain the water saturation at the moment in time, t n .
In Formula (12), the number of pressure equations is equal to the number of downhole assemblies, while in traditional reservoir modeling, the number of pressure equations is equal to the number of grid blocks. Therefore, the INSIM method requires considerably less computational resources than 3D direct simulation of dynamic processes in the reservoir. After the pressures are calculated, the saturation values are estimated. Further, it is easy to calculate the flow rate between the pairs of wells, which can be used for further diagnostics of downhole dynamics. Since the number of downhole assemblies is limited, the methodology may only give an approximate distribution of the formation fluid. However, the results of the calculations to date show the saturation or pressure distributions calculated by the INSIM model compared with historical data are sufficient to obtain a successful production forecast.
The INSIM model was improved to become the INSIM-FT model—an interwell simulation model with fluid front tracking. The key difference is the calculation of the water cut distribution between the wells. Unlike INSIM, INSIM-FT also includes parameters that define relative permeability power curves.
The distribution of water saturation between two wells is described by the following form of the one-dimensional Buckley–Leverett equation [24]:
S w ( x , t ) t + q t , i , j ( t ) ϕ i , j A i , j f w ( x , t ) x = 0 ,   0 x L i , j , t n 1 t t n
f S w = k r w μ w k r w μ w + k r o μ o
The curves of relative phase permeabilities are set by analytical relations. In this paper, we used the Corey correlation:
k r w S w n = α · S w n n w
k r o S w n = S w n n o
S w n = S w S i w 1 S i w S o r
To solve Equation (13), the front-tracking method is used, which is an adaptation of the method [25]. This method is stable and has a relatively low variance.
The general idea of the solution is as follows:
  • Each pair of wells is represented as a quasi-one-dimensional model. This model is divided into cells in which the water cut value is set. The water cut function is approximated as a sequence of constant values (Figure 2).
  • In each interval, the task of calculating the velocity of the shock wave cluster is solved.
Determination of the right and left boundary values of the water cut, as well as the accuracy of the water cut determination is S w l , S w r , δ S w
Calculation of the Rankine–Hugoniot condition is σ t r i a l = f w S w l f w S w r S w l S w r
If f w S w l > σ t r i a l > f w S w r , then the result is asingle shock wave.
If f w S w l < f w S w r and f w S w l f w S w r > 0 , then the result is a depression wave. It is simulated as a sequence of shock waves connecting values S w l , S w r .
In other cases, a sequence of shock waves to approximate a depression wave connecting values S w l , S w * that follows the shock wave connecting the values S w * , S w r is used.
The value calculation, S w * , is based on the equation f w S w * = f w S w * f w S w r S w * S w r .
For the calculation of the first intersections, in the x-t diagram, the calculated velocities are straight lines with different slopes. As seen on Figure 3, the first intersections appear at the boundaries of the primary areas.
  • The point of intersection becomes a new boundary, and it is necessary to calculate new velocities of the shock waves in new intervals. Then, the process continues until the time value, Δt, is reached.
  • Go to the calculation in the next time interval.

2.2. Model Adaptation

As the model is adapted to history, interwell conductivity and interwell volumes are selected, as well as parameters that determine the power functions of relative phase permeabilities used in the calculation of saturation and other parameters required for solving the previously mentioned equations. The complete list of parameters used in this paper is presented in Appendix A. The INSIM-FT model itself consists of at least 14 adaptation parameters.
Adaptation parameters are distribution functions having the following parameters: minimum value, maximum value, mean value, distribution type, and value of dispersion parameter, σ (Figure 4). During the adaptation period, the parameters variation intervals may change. For correct calibration, it is necessary to determine the intervals of these parameters’ variations with high accuracy.
In adaptation, we use the ES-MDA method [26] (smoothing assembly with multiple data assimilation). The task of adapting historical data in matrix form looks like:
O m = 1 2 g m d o b s T C D 1 g m d o b s
where C D is the covariance matrix that adjusts the sensitivity of the model. The general view of the matrix is as follows:
C D = c o v d c a l c 1 d o b s 1 d c a l c 1 d o b s 1 0 0 c o v d c a l c n d o b s n d c a l c n d o b s n
Volumetric flow rates and water cut values are used as calculated and observed parameters.
In the task of historical data adaptation, it is necessary to minimize the discrepancy between the observed and calculated data. The calculated data are calculated as a direct hydrodynamic task based on the model parameters (capacity, drainage pores volumes, compressibility, etc.). At the same time, it is necessary to periodically check and adjust the model parameters for compliance with historical information. The model parameters are restored from the inverse optimization problem, reducing the error with the observed data, according to the following formula:
m j u = m j p + C M D C D D + α i C D 1 d u c d j p
i = 1 N i 1 α i = 1 ,   i = 1 N
The time interval of using the historical and calculated data is called the data assimilation step N.

2.3. Forecast of Oil and Liquid Flow Rate

The flow rate forecast is based on the Dupuis formula:
Q r = 2 π k h P k P c μ ln R k R c
where k is the permeability coefficient; h —reservoir capacity; P k and P c —pressure on the feed loop and in the well, respectively; R k and R c —the radii of the feed loop and the well, respectively; μ —the fluid viscosity; Q r —well flow rate. Assume all parameters are constants except for the pressure difference. Then, the Dupuis formula can be rewritten as follows:
η = Q Δ P
where Δ P is pressure differential, Q —the flow rate, η —productivity index.
Then, the task of forecasting the flow rate is reduced to finding the productivity index based on the known values of the flow rate and pressure difference. To solve the problem, it is proposed to use the linear regression method.
Knowing the pressure, it is possible to predict the flow rate. Figure 5 shows the forecast results: the blue curve is the initial flow rate; yellow is the initial pressure; the red section of the pressure curve is the time interval, the points of which were used to plot the regression; the green curve is the forecasted flow rate for the entire interval.

3. Testing the Method on Synthetic Test Data

The main purpose of this section is to verify the correctness of the model adaptation on synthetic data. In this paper, a hydrodynamic model in tNavigator format was used. It represents a sector with geomechanical and filtration-capacitive properties. To assess the quality of adaptation on synthetic data, four production and three injection wells were added to the model, and calculations were performed with a period of 25 years.
Figure 6 shows the distribution of water saturation in the field, as well as the layout of production and injection wells.
Below are the results of the developed INSIM-FT prototype (Figure 7 and Figure 8).
The fluid flow rate and water cut were adapted using synthetic data for the time period from 2020 to 2045. Based on the results of the adaptation, it can be concluded the fluid flow rate is reproduced well, and the water cut is reproduced satisfactorily.
To improve the results of water cut adaptation at well-4, sensitivity analysis was carried out to the number of assimilation steps (Figure 9).
It was determined that the number of assimilation steps significantly affects the final result. Increasing the number of steps to 16 reduced the average error from 17.51% to 11.65%.
It should also be noted that using the approach described in the article [27], the adaptation rate directly depends on the number of assimilation steps. The increase in speed is due to the decrease in the number of data assimilation steps since it is multiplied by the sum of the reciprocal values of the MDA coefficients. Thus, it is necessary to select such a number of assimilation steps that will maintain the quality of adaptation at the required level but with minimal adaptation time. These algorithms require a more detailed study.

4. Methods of Accounting for Geological and Technical Measures

4.1. Hydraulic Fracturing

For the INSIM-FT model, we proposed the following solution: according to the source [27], during the injection of water into the injection well, as well as during the formation and propagation of cracks, hydraulic auto-fracturing takes place due to several factors:
  • Equality of pressure in the injection well and horizontal stress, which can cause the formation of a developed pattern of cracks causing a high rate of pressure drop of the injection fluid;
  • Change in effective horizontal stress due to temperature change at the bottomhole;
  • Change in effective horizontal stresses due to changes in pore pressure;
  • Stress contrast in the rock between different geological layers (clays and sandstones);
  • Difference between vertical and horizontal stresses (the stress anisotropy coefficient is higher than 1.15).
It is assumed that during the hydraulic fracturing procedure or the occurrence of auto-fracturing, there is a sharp increase in the absolute permeability parameter k, i, j of the flow tube. However, over time, the parameter should return to the “shelf” (due to natural contamination of the bottom hole zone). That is, it is necessary to multiply it by a time-dependent monotonically decreasing function. The possible form of the additional function will be exponential, and the following parameters may be present in the exponent index:
  • Borehole radius across the bit;
  • Effective well radius;
  • Reservoir permeability;
  • Fluid viscosity;
  • Geometric characteristics of the producing reservoir;
The following type of additional function is proposed as the first approximation:
g s , t = a e s t b + c
where a , b , c are adaptable parameters that depend on the type of geological and technical measures performed, s is skin factor, which is expressed by dependence s = ln R r , r is the wellbore radius across the bit, and R is the effective radius of the wellbore.
An illustration of this dependence is shown in Figure 10.
Thus, it is necessary to create a correspondence table “type of geological and technical measures”—adaptable coefficients a , b , c .
In the case of hydraulic fracturing and hydraulic auto-fracturing, it is possible to consider the formula for the skin factor, which takes into account the dimensionless conductivity of the fracture, C F D , the value depending on the difference in permeability of the proppant and the formation:
C F D = k f w k x f
where k f is proppant permeability (mD), k is formation permeability (mD), w is fracture width (m), x f is the crack half-length (m).
The skin factor is calculated according to the following algorithm.
  • Option 1:
C F D = k f w k x f
r e f = r w + x f 2 1 + C F D 1.7 1.01
S = ln r e f r w
  • Option 2:
C F D = k f w k x f
u = ln ( C F D )
f = 1.65 0.328 u + 0.116 u 2 1 + 0.18 u + 0.064 u 2 + 0.005 u 3
S = f + ln r w x f

4.2. Hydraulic Auto-Fracturing

Assume that auto-fracturing occurs when the injection pressure exceeds the sum of the formation pressure and horizontal stress of the formation:
P h f r > P f o r m + σ h o r
The horizontal rock stress can be expressed through the vertical stress as follows:
σ h o r = 1 3 ÷ 1 2 σ v e r t
In turn, vertical intensity is related to geostatic pressure:
P g e o s t = P f o r m + σ v e r t
Therefore:
σ v e r t = P g e o s t P f o r m
Geostatic pressure is calculated according to the following relationship:
P g e o s t = ρ p r o c k g H
The correction introduced during hydraulic auto-fracturing has the same structure as in the case of hydraulic fracturing, but the point in time when this multiplicative additive is introduced is determined by the condition recorded above.
The next step was to regulate the multiplicative additive: setting the exponent degree, amplitude, and shift. As the first approximation, it is proposed to determine these parameters manually for each type of geological and technical measures. In the INSIM-FT code, these values were used as averages, after which they are edited and adapted.

5. Testing the Model on Real Data While Accounting for the Geological and Technical Measures

Initially, calculations in the INSIM-FT model are carried out without accounting for geological and technical measures. As mentioned in the previous paragraph, it is necessary to introduce a multiplicative exponential additive, which, in fact, shows the pollution of the bottomhole zone and decrease in water saturation. Figure 11 shows an illustration of the real field sector. The wells for which the adjustment for geological and technical measures have been introduced are shown in red.
Let us compare the influence of geological and technical measures on the calculation of water cut. Let us consider two cases: neutral and positive effect of the amendment. The results are presented below (Figure 12 and Figure 13).
Despite a correction introduced on well 1 (hydraulic fracturing was performed on it), the algorithm could not significantly reduce the error. This is probably due to the mutual influence of various types of geological and technical measures. It may be necessary to introduce additional corrections. Adaptation of water cut accounting for the geological and technical measures at well 2 yielded a positive effect. It is clear to the naked eye (Figure 13) that the calculated curve significantly approached the actual one and repeated its shape, accounting for the features.
Calculations have demonstrated the efficiency of this approach. Further, a software realization of the algorithm is required, which will take into account all types of geological and technical measures according to the input parameters. It is also necessary to create a library of multiplicative correction parameters for various types of geological and technical measures or to select coefficients in automatic mode.

6. Discussion

Our results showed the introduction of the geological and technical measures into the approach described in the paper allows in several cases to significantly increase the design parameters approximation accuracy. We associate the partial low sensitivity of the solution to corrections for geological and technical measures with insufficient knowledge and elaboration of the parameterization of certain types of geological and technical measures. It is necessary to carry out a complex analysis of the mutual influence of various types of geological and technical measures and analysis of certain geological and technical measure zones of influence.
The choice of assembly and assimilation parameters, such as: expectation, dispersion, number of iterations of the ES-MDA algorithm, and assembly size, plays a key role for fine settings the model and predictive calculations. The search for balance between computational costs and procedure for approximating historical indicators makes it possible to find sufficiently accurate and computationally efficient approaches to assessing the operating parameters of the hydrocarbon field.
The results of the experiments along with the high computational efficiency of the method showed good potential for parallel implementations since the operations of multithreaded computations on the CPU were used in the paper; it is reasonable to further implement the method for GpGPU platforms using low-level languages, for example, C ++. It is also necessary to carefully consider the choice of a framework for the implementation of calculations on graphic cards; the largest modern vendors of graphic equipment created software solutions only for their “hardware” components without the support of competitors’ hardware and software.
The authors of the paper are convinced that the presented approach has rich opportunities for “tuning”, for example, building proxy connections of wells based on more advanced algorithms, Voronoi diagrams, PEBI grids, etc., using optimized solvers or more efficiently solving differential equations and SLAEs computational approaches.

7. Conclusions

The method of the field hydrodynamic model based on the data-driven (setting to historical data) approach and the proxy model of fluid migration dynamic processes in hydrocarbon reservoir was implemented. Based on efficient averaging of reservoir hydrodynamic parameters for a finite set of production and injection wells, the method shows significant computational efficiency compared to classical reservoir simulators (Eclipse, tNavigator). The proposed additional adaptation of the model to the geological and technical measures parameters has shown its efficiency in predicting the operating conditions of the hydrocarbon field.
The use of multithreaded calculation technologies made it possible to significantly reduce the simulation time for a particular case; in the future, it is planned to transfer the algorithm to the GpGPU platform, which, according to preliminary experiments, will increase the performance by at least an order of magnitude.
For the correct settings of water saturation and reservoir pressures, it is necessary to fine-tune the model parameters: selection of the optimal interval of parameter variation, selection of the parameter distribution function, selection of the initial geological and technical measures adaptation coefficients while considering the interference of different types of geological and technical measures.

Author Contributions

M.O.—investigation, writing—original draft preparation; E.L.—investigation, writing—original draft preparation; A.C.—writing—original draft preparation, writing—review and editing, conceptualization; S.F.— writing—review and editing; R.K.—formal analysis, writing—original draft preparation; E.U.—methodology, resources; V.U.—resources, supervision; D.T.—conceptualization, writing—review and editing; N.K.—writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Funding

The research was carried out in accordance with the FSUS-2022-0020 project.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1. List of adaptation parameters.
Table A1. List of adaptation parameters.
ParameterSymbolUnits
Initial reservoir pressurePPa
Oil viscosityμ0Pa·s
Water viscosityμwPa·s
Rock compressibilitycrPa−1
Water compressibilitycwPa−1
Oil compressibilitycoPa−1
Residual water saturationSiwfraction
Residual oil saturationSrofraction
Parameter in Corey correlationnw
Parameter in Corey correlationno
Parameter in Corey correlationα
Porosity of flow-tubeφi,jfraction
Permeability of flow-tubeki,jmD
Cross-section area of flow-tubeAi,jm2

References

  1. Yaskin, S.A.; Mukhametshin, V.V.; Andreev, V.E.; Dubinskiy, G.S. Forecasting the Parameters of Non-Stationary Water Flooding of Oil Deposits. IOP Conf. Ser. Earth Environ. Sci. 2018, 194, 062037. [Google Scholar] [CrossRef]
  2. Mukhametshin, V.S.; Tyncherov, K.T.; Rakhimov, N. Geological and Technological Substantiation of Waterflooding Systems in Deposits with Hard-to-Recover Reserves. IOP Conf. Ser. Mater. Sci. Eng. 2021, 1064, 012068. [Google Scholar] [CrossRef]
  3. Aziz, K.; Settari, A. Petroleum Reservoir Simulation; Applied Science Publishers Ltd.: London, UK, 1979. [Google Scholar] [CrossRef]
  4. Khalimov, E.M. Detailed Geological Models and Three-Dimensional Simulation. Neft. Geologiâ. Teor. I Pract. 2012, 7, 17–27. [Google Scholar]
  5. Chudin, Y.S.; Chumakov, G.N.; Fedorov, I.A.; Pyatakova, O.A.; Gilev, D.V.; Petrukhin, V.A.; Tarasov, G.V.; Arkhipov, Y.A. Using Proxy Models of Gas Deposits for Production Optimization. Gazov. Promishlennost 2020, 4, 30–36. [Google Scholar]
  6. Bahrami, P.; Sahari, M.F.; James, L.A. A Review of Proxy Modeling Highlighting Applications for Reservoir Engineering. Energies 2022, 15, 5247. [Google Scholar] [CrossRef]
  7. Nehoroshkova, A.A.; Danko, M.Y.; Zavyalov, A.S.; Elisheva, A.O. Critical Analysis of the INSIM-FT Proxy Modeling Method (Interwell Numerical Simulation Front Tracking Models) at Synthetic Models and Real Field. Neft. Gas. Novacii 2019, 12, 49–55. [Google Scholar]
  8. Udy, J.; Hansen, B.; Maddux, S.; Petersen, D.; Heilner, S.; Stevens, K.; Lignell, D.; Hedengren, J. Review of Field Development Optimization of Waterflooding, EOR, and Well Placement Focusing on History Matching and Optimization Algorithms. Processes 2017, 5, 34. [Google Scholar] [CrossRef]
  9. Jafroodi, N.; Zhang, D. New method for reservoir characterization and optimization using CRM–EnOpt approach. J. Pet. Sci. Eng. 2011, 77, 155–171. [Google Scholar] [CrossRef]
  10. Jansen, J.D.; Durlofsky, L.J. Use of Reduced-Order Models in Well Control Optimization. Optim. Eng. 2017, 18, 105–132. [Google Scholar] [CrossRef]
  11. Cardoso, M.A.; Durlofsky, L.J.; Sarma, P. Development and Application of Reduced-Order Modeling Procedures for Subsurface Flow Simulation. Int. J. Numer. Meth. Eng. 2009, 77, 1322–1350. [Google Scholar] [CrossRef]
  12. Gildin, E.; Klie, H.; Rodriguez, A.; Wheeler, M.F.; Bishop, R.H. Projection-Based Approximation Methods for the Optimal Control of Smart Oil Fields. In Proceedings of the ECMOR X—10th European Conference on the Mathematics of Oil Recovery, Amsterdam, The Netherlands, 4–7 September 2006; European Association of Geoscientists & Engineers: Amsterdam, The Netherlands, 2006. [Google Scholar] [CrossRef]
  13. Gildin, E.; Klie, H.; Rodriguez, A.; Wheeler, M.F.; Bishop, R.H. Development of Low-Order Controllers for High-Order Reservoir Models and Smart Wells. In All Days; SPE: San Antonio, TX, USA, 2006; p. SPE-102214-MS. [Google Scholar] [CrossRef]
  14. Liu, F.; Mendel, J.M. Forecasting Injector-Producer Relationships From Production and Injection Rates Using an Extended Kalman Filter. In All Days; SPE: Anaheim, CA, USA, 2007; p. SPE-110520-MS. [Google Scholar] [CrossRef]
  15. Jansen, F.E.; Kelkar, M.G. Non-Stationary Estimation of Reservoir Properties Using Production Data. In All Days; SPE: San Antonio, TX, USA, 1997; p. SPE-38729-MS. [Google Scholar] [CrossRef]
  16. Bailer-Jones, C.A.L.; MacKay, D.J.C.; Withers, P.J. A Recurrent Neural Network for Modelling Dynamical Systems. Netw. Comput. Neural Syst. 1998, 9, 531–547. [Google Scholar] [CrossRef]
  17. Sayarpour, M. Development and Application of Capacitance-Resistive Models to Water/CO2 Floods. Ph.D. Thesis, The University of Texas at Austin, Austin, TX, USA, August 2008. [Google Scholar] [CrossRef]
  18. Yousef, A.A.; Gentil, P.; Jensen, J.L.; Lake, L.W. A Capacitance Model to Infer Interwell Connectivity from Production and Injection Rate Fluctuations. In All Days; SPE: Dallas, TX, USA, 2005; p. SPE-95322-MS. [Google Scholar] [CrossRef]
  19. Cao, F.; Luo, H.; Lake, L.W. Oil-Rate Forecast by Inferring Fractional-Flow Models From Field Data With Koval Method Combined With the Capacitance/Resistance Model. SPE Reserv. Eval. Eng. 2015, 18, 534–553. [Google Scholar] [CrossRef]
  20. Sayarpour, M.; Zuluaga, E.; Kabir, C.S.; Lake, L.W. The Use of Capacitance–Resistance Models for Rapid Estimation of Waterflood Performance and Optimization. J. Pet. Sci. Eng. 2009, 69, 227–238. [Google Scholar] [CrossRef]
  21. Lerlertpakdee, P.; Jafarpour, B.; Gildin, E. Efficient Production Optimization With Flow-Network Models. SPE J. 2014, 19, 1083–1095. [Google Scholar] [CrossRef]
  22. Guo, Z.; Reynolds, A.C.; Zhao, H. Waterflooding Optimization with the INSIM-FT Data-Driven Model. Comput. Geosci. 2018, 22, 745–761. [Google Scholar] [CrossRef]
  23. Guo, Z.; Reynolds, A.C.; Zhao, H. A Physics-Based Data-Driven Model for History Matching, Prediction, and Characterization of Waterflooding Performance. SPE J. 2017, 23, 367–395. [Google Scholar] [CrossRef]
  24. Buckley, S.E.; Leverett, M.C. Mechanism of Fluid Displacement in Sands. Trans. AIME 1942, 146, 107–116. [Google Scholar] [CrossRef]
  25. Godunov, S.K.; Bohachevsky, I. Finite Difference Method for Numerical Computation of Discontinuous Solutions of the Equations of Fluid Dynamics. Mat. Sb. 1959, 47, 271–306. [Google Scholar]
  26. Emerick, A.A.; Reynolds, A.C. History-Matching Production and Seismic Data in a Real Field Case Using the Ensemble Smoother With Multiple Data Assimilation. In All Days; SPE: The Woodlands, TX, USA, 2013; p. SPE-163675-MS. [Google Scholar] [CrossRef]
  27. Baykov, V.A.; Burakov, I.M.; Latypov, I.D.; Yakovlev, A.; Asmandiyarov, R.N. Waterflood Induced Hydraulic Fracturing Control under Reservoir Pressure Maintenance Conditions on RN-Yuganskneftegas Oilfields. Oil Ind. 2012, 11, 30–33. [Google Scholar]
Figure 1. Schematic illustration of the pore channels network model in the INSIM model.
Figure 1. Schematic illustration of the pore channels network model in the INSIM model.
Energies 16 01648 g001
Figure 2. Partitioning of the quasi-one-dimensional model.
Figure 2. Partitioning of the quasi-one-dimensional model.
Energies 16 01648 g002
Figure 3. Intersections of shock wave clusters.
Figure 3. Intersections of shock wave clusters.
Energies 16 01648 g003
Figure 4. Schematic distribution of the parameter in the model.
Figure 4. Schematic distribution of the parameter in the model.
Energies 16 01648 g004
Figure 5. Forecast results.
Figure 5. Forecast results.
Energies 16 01648 g005
Figure 6. Hydrodynamic model of the field (red vertical lines—production wells, blue—injection wells).
Figure 6. Hydrodynamic model of the field (red vertical lines—production wells, blue—injection wells).
Energies 16 01648 g006
Figure 7. Displaying the field map in the INSIM-FT prototype. The studied production wells 1 and 4 are marked with a red rectangle.
Figure 7. Displaying the field map in the INSIM-FT prototype. The studied production wells 1 and 4 are marked with a red rectangle.
Energies 16 01648 g007
Figure 8. Approbation of the model on synthetic data.
Figure 8. Approbation of the model on synthetic data.
Energies 16 01648 g008
Figure 9. Analysis of the forecast sensitivity to the number of assimilation steps.
Figure 9. Analysis of the forecast sensitivity to the number of assimilation steps.
Energies 16 01648 g009
Figure 10. Effect of geological and technical measures on the hydraulic conductivity of the flow tube.
Figure 10. Effect of geological and technical measures on the hydraulic conductivity of the flow tube.
Energies 16 01648 g010
Figure 11. Illustration of the field sector where the geological and technical measures effect was tested. Production wells 1 and 2 are in red rectangles.
Figure 11. Illustration of the field sector where the geological and technical measures effect was tested. Production wells 1 and 2 are in red rectangles.
Energies 16 01648 g011
Figure 12. Comparison of water saturation adaptation for well 1 (on the left—without accounting for geological and technical measures, on the right—accounting for geological and technical measures).
Figure 12. Comparison of water saturation adaptation for well 1 (on the left—without accounting for geological and technical measures, on the right—accounting for geological and technical measures).
Energies 16 01648 g012
Figure 13. Comparison of water saturation adaptation for well 2 (on the left—without accounting for geological and technical measures, on the right—accounting for geological and technical measures).
Figure 13. Comparison of water saturation adaptation for well 2 (on the left—without accounting for geological and technical measures, on the right—accounting for geological and technical measures).
Energies 16 01648 g013
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Ovsepian, M.; Lys, E.; Cheremisin, A.; Frolov, S.; Kurmangaliev, R.; Usov, E.; Ulyanov, V.; Tailakov, D.; Kayurov, N. Testing the INSIM-FT Proxy Simulation Method. Energies 2023, 16, 1648. https://doi.org/10.3390/en16041648

AMA Style

Ovsepian M, Lys E, Cheremisin A, Frolov S, Kurmangaliev R, Usov E, Ulyanov V, Tailakov D, Kayurov N. Testing the INSIM-FT Proxy Simulation Method. Energies. 2023; 16(4):1648. https://doi.org/10.3390/en16041648

Chicago/Turabian Style

Ovsepian, Mkhitar, Egor Lys, Alexander Cheremisin, Stanislav Frolov, Rustam Kurmangaliev, Eduard Usov, Vladimir Ulyanov, Dmitry Tailakov, and Nikita Kayurov. 2023. "Testing the INSIM-FT Proxy Simulation Method" Energies 16, no. 4: 1648. https://doi.org/10.3390/en16041648

APA Style

Ovsepian, M., Lys, E., Cheremisin, A., Frolov, S., Kurmangaliev, R., Usov, E., Ulyanov, V., Tailakov, D., & Kayurov, N. (2023). Testing the INSIM-FT Proxy Simulation Method. Energies, 16(4), 1648. https://doi.org/10.3390/en16041648

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop