Introduction

The idea of optimization of oil and gas production using mathematical modeling originated in the 1960s. However, its implementation was very much restricted by such factors as poor capacity of RAM and the computing performance characteristics. In the 1990s, when this barrier was finally overcome and it became possible to calculate natural gas reservoirs the mathematical models of which comprised a large number of calculation points and within a timeframe that was acceptable to the user, the market of the software products for reservoir modelling saw the appearance of competing companies.

The leading positions in this respect have for years been held by the Schlumberger’s software products which are used both as stand-alone software packages and in combination with other applications, as can be seen, for example, in the study of Gerogiorgis et al. (2006). This study discusses strategy for interfacing reservoir simulation (ECLIPSE®) with equation-oriented process optimization (gPROMS®) and presents a relevant application.

Alongside with the Schlumberger company, such companies as Smedvig Technologies, Roxar Software Solutions, Western Atlas, Landmark Graphics, Paradigm Geophysical, CogniSeis, CGG Petrosystems, PGS Tigress, Seismic Microtechnology, GeoMatic, Quick look, Tigress, Western Atlas have also been involved in the development of software products for reservoir simulation.

However, most of these software products calculate only existing modes of reservoir development or their change according to a predetermined option. The search for the best possible mode for reservoir development is usually carried out on the basis of a certain number of options, while it is absolutely not obvious that there will be an optimal one among them. This results in a decrease of oil and gas production than it would be if a really optimal option was chosen. It is determined by the redistribution of individual well flow rates, as proposed in this study or by optimizing the initial placement of wells as proposed by Forouzanfar and Reynolds 2013, Forouzanfar et al. (2012).

When developing gas reservoirs, three stages (periods) are distinguished. The first one is a period of incrementally increasing production (a period of drilling up the reservoir, construction of related facilities, and transition of the reservoir to constant production level). Duration of the stage is up to 7–10 years, and gas extraction rate reaches 20–25% of initial reserves. The second period is characterized by constant production level (drilling of the reservoir continues to maintain a constant production capacity and a booster compressor station is being built or its capacity is increased). During this period, up to half of the initial gas reserves is extracted and a gas recovery factor reaches 60–70%. The period of constant production depends on a gas extraction rate achieved: the higher it is, the shorter is the period. The third period is characterized by declining production, a decrease in the number of producing wells, their production rates, the appearance of water in wells and a significant decrease in reservoir pressure. When a reservoir enters the period of declining production, there comes an issue of an optimal choice and the most effective use of the technologies for additional development of the reservoir that makes it possible to conduct cost-effective production of low-pressure gas (Ping et al. 2014; Curtis 2003). The problem of daily production optimization in the exploration and production of oil and gas was also studied by Samuelson (2008). The main objective of his research is to find the optimal solutions that utilize the production system efficiently and maximize revenue. Mathematical models are used to find optimal operations in such processes.

Currently, there are many different ways to increase the reservoir productivity which is at the third stage of development. Among the main methods, the following ones can be singled out (Krishnamoorthy et al. 2016; Feng et al. 2012; Xiaogang et al. 2010; Ren et al. 2007; Zhu et al. 2009; Huang et al. 2012; Liao et al. 2003; Biagi et al. 2014):

  • drilling of horizontal wells;

  • hydraulic fracturing (the method is widely used in most reservoirs with declining production in conjunction with other measures to increase reservoir productivity);

  • gas production with simultaneous water pumping out of the reservoir;

  • gas production with simultaneous injection of water, nitrogen or carbon dioxide;

  • pressure-dependent separate piping to exclude the impact of wells with high pressure on the regular performance of low-pressure wells.

  • optimization of reservoir development, for example, due to acid dissolution of deposits in wells in some places of the reservoir and temporary blockage of wells in others; through the infill drilling; through the use of drilling technology that reduces hydrostatic pressure in the wellbore, etc.

A rather complete review of optimization methods for increasing natural gas production such as planned production, wells placement, design of a gas production system is presented in the study of Zheng et al. (2010).

All of the above-mentioned methods require additional and often very significant financial costs. Redistribution of gas extraction rates in different wells to evenly reduce reservoir pressure and to shorten the final stage of reservoir development can be introduced at minimum expenditures. This approach proposed in this study is of practical interest in the development of low-permeability gas reservoirs.

To enhance the gas extraction, when the reservoir is at its final stage of development, it is necessary to redistribute gas flow rates in the wells during the stage of continuous production (while preserving the planned level of annual extraction rate) in such a way as to enable all producing wells have the same recharge area taking into account permeability, porosity, and thickness. In this case, optimal flow rates in wells will be ensured during the entire final period of the reservoir development, the duration of the final stage will be reduced, and it will be possible to bring the gas recovery factor up to 95–97%.

In some cases, when natural conditions do not impose restrictions on well flow rates, gas extraction rates from wells are set based on the cost-effectiveness analysis or consumer needs. In one way or another, technological regimes contain some restrictions that must be taken into account when operating gas wells.

The approach proposed in our study can be used as a base for the proactive control of gas recovery. Contrary to the reactive control that responds quickly to unexpected events like water breakthrough, proactive control, that is based on mathematical simulation of reservoir, allows to find the optimum well flow rates and thus to increase the total gas recovery.

Mathematical and physical foundations of the model

Formulation of the filtration problem

In our previous study (Kalugin et al. 2015), we have formulated and solved the problem of optimizing the distribution of well flow rates to maximize the recovery of gas condensate C5+ taking into account phase transitions for both the cycling and depletion modes. Based on the obtained solution, an algorithm was developed that allows to optimally manage the process of reservoir developing taking into account the constraints imposed on flow rates of certain wells, as well as their shutoffs for remedial maintenance.

Initial equations for the isothermal filtration of a multicomponent mixture in gas-bearing strata while operating the system of producing and injection wells were recorded using an assumption of local thermodynamic phase balance, the validity of the generalized Darcy law neglecting the influence of capillary, diffusion forces, and gravity:

$$h\frac{{\partial \left( {m{z_k}{F_k}} \right)}}{{\partial t}} - div\left( {{k_0}h{z_k}{\beta _k}grad{\kern 1pt} P} \right)=\sum\limits_{{n=1}}^{N} {\rho _{{g,k}}^{{std}} \cdot {Q_{k,n}}\left( t \right)\delta \left( {x - {\eta _n},\;y - {\xi _n}} \right)} , \quad \left( {k=\overline {{1,{N_C}}} } \right),$$
(1)

where \({F_k}=\frac{{\left( {1 - S} \right){\rho _g}{k_k}+S{\rho _w}}}{{1+W\left( {{k_k} - 1} \right)}},\;{\beta _k}=\left( {\frac{{{\rho _g}\,{f_g}{k_k}}}{{{\mu _g}}}+\frac{{{\rho _w}\,{f_w}}}{{{\mu _w}}}} \right)\frac{1}{{1+W\left( {{k_k} - 1} \right)}}\), \({\beta _k}=\left( {\frac{{{\rho _g}\,{f_g}{k_k}}}{{{\mu _g}}}+\frac{{{\rho _w}\,{f_w}}}{{{\mu _w}}}} \right)\frac{1}{{1+W\left( {{k_k} - 1} \right)}}.\)

The unknown functions of coordinates \(x\), \(y\), and time \(t\) in the above equations are the pressure P and molar fractions of the hydrocarbon components in the mixture \({z_k}\). Each equation from the system (1) is a continuity equation for one of the components of the multicomponent mixture in its differential form.

The well-known Peng–Robinson equation was applied to close the system of equations for multicomponent filtration with regard to phase transitions. In the case of enhanced gas extraction, we used the set of Eq. (1) for the mixture components containing \({C_i}\left( {i=\overline {{1,4}} } \right)\), i.e., for those that remain in the gaseous state during extraction.

Since no phase transitions occur for these components at the pressure change, we have used the Clapeyron equation instead of Peng–Robinson one while calculating the state of the actual gas, corrected for the deviation of actual gas behavior from that of perfect gas via a compressibility factor Z:

$$P={\rho _g}\frac{{RT}}{Z}.$$
(2)

The formulation of the filtration problem given above is a special case of the filtration problem for a multicomponent mixture with phase transitions and is equivalent to the «black oil» model. As a result of some transformations (Kalugin et al. 2015), the original set of Eq.(1) can be reduced to a form convenient for solving it through numerical methods:

$$mh\frac{{\partial {\rho _g}}}{{\partial t}} - \frac{\partial }{{\partial x}}\left( {{k_0}h\frac{{{\rho _g}}}{{{\mu _g}}}\frac{{\partial P}}{{\partial x}}} \right) - \frac{\partial }{{\partial y}}\left( {{k_0}h\frac{{{\rho _g}}}{{{\mu _g}}}\frac{{\partial P}}{{\partial y}}} \right) = \rho _{g}^{{std}}\sum\limits_{{n=1}}^{{{N_w}}} {{Q_n}\left( t \right)\delta \left( {x - {\eta _n},\;y - {\xi _n}} \right)} ,$$
(3)
$$m{F_k}\frac{{\partial {z_k}}}{{\partial t}}={k_0}{\beta _k}\frac{{\partial P}}{{\partial x}}\frac{{\partial {z_k}}}{{\partial x}}+{k_0}{\beta _k}\frac{{\partial P}}{{\partial y}}\frac{{\partial {z_k}}}{{\partial y}},\;\left( {k=\overline {{1,{N_C}}} } \right).$$
(4)

Initial conditions determine the value of unknown functions prior to the external actions on the stratum. In the case of 2D filtration, these conditions are written in the form:

$$P\left( {x,y,0} \right)={P^0}\left( {x,y} \right),\quad {z_k}\left( {x,y,0} \right)=z_{k}^{0}\left( {x,y} \right),\quad \left( {x,y} \right) \in G,\quad \left( {k=\overline {{1,{N_C}}} } \right),$$
(5)

where \({P^0}\left( {x,y} \right)\) and \(z_{k}^{0}\left( {x,y} \right)\) are given values of the unknown functions at the initial moment \(t=0\).

When modeling the enhanced development of the reservoir at its boundary\(\Gamma \left( {x,y} \right)\), the conditions of impenetrability are set

$${\left. {\frac{{\partial P\left( {x,y,t} \right)}}{{\partial n}}} \right|_{\left( {x,y} \right) \in \Gamma }}=0,\quad {\left. {\frac{{\partial {z_k}\left( {x,y,t} \right)}}{{\partial n}}} \right|_{\left( {x,y} \right) \in \Gamma }}=0,\quad t>0.$$
(6)

Thus, under adopted boundary conditions, only the operation of the wells causes all changes in unknown functions within the area. The set of Eqs. (2)–(4) given above, as well as initial and boundary conditions (5), (6) allows finding the values of simulated pressure function \(P\left( {x,y,t} \right)\), density \({\rho _g}\left( {x,y,t} \right)\), and molar fractions of the mixture components \({z_k}\left( {x,y,t} \right)\), \(\left( {k=\overline {{1,{N_C}}} } \right)\).

Algorithm for numerical solution of the filtration problem

The solution of the set of equations described above for practical filtration problems in inhomogeneous domains is obtained by a well-known, widely used and reliable approximate numerical method of finite differences. An absolutely stable implicit conservative difference scheme was used to approximate the initial differential equations. This scheme leads to the necessity at each step of time to solve a set of algebraic equations in which the number of unknowns corresponds to the number of internal nodes of the reticular domain (for example, for the canonical rectangular region, the number of equations \({N_e}\) is equal to \({N_e}={N_x}{N_y}\)).

The direct methods used for solving sets of algebraic equations requires for each time step the following time-consuming operations: calculation of the equation coefficients to populate the matrix having an order \(N_{e}^{2}\) and the vector of absolute terms consisting of \({N_e}\) components; calculation of the main determinant of the system and then the matrix inversion, when doing the calculation of each unknown.

Therefore, the solution of sets of algebraic equations was obtained with the help of a new iterative procedure developed by the authors that ensures high convergence and the required accuracy of the computational process and is convenient for the development of algorithms and programs for a PC. The technique of the numerical solution is described in detail by Kalugin et al. (2007).

To implement the described algorithm, a computer program in Delphi has been developed for any number of layers, wells, and mixture components. When solving specific problems, most of the restrictions are related to the hardware capabilities of PCs at the disposal of specialists, completeness, and quality of the initial information and the amount of time allotted for the research. It is obvious that with the further development of hardware, it will be efficient to solve such problems on multi-processor PCs with parallel processing of computations.

Formulation of the optimization problem

Let a reservoir have \(N\) producing wells. Let \(Q\) be the planned production level from the whole reservoir, and \({q_i}\), \(\left( {i=\overline {{1,N}} } \right)\) be extraction rates in individual wells. It is required to optimize extraction rates for each well in such a way that pressure drop in the recharge area of each well would be the same. Considering that the recharge area depends on the reservoir thickness, its porosity and permeability, this condition can be formulated as maximizing the weight-average reservoir pressure for the group of producing wells (Kalugin et al. 2015).

Gas production capabilities of a well are determined by the maximum possible flow rate \({q_{i\hbox{max} }}\) which corresponds to the difference between formation pressure in an area of the \(i\)th production well and minimum pressure at the inlet to a gas treatment plant. When the optimization algorithm is running, \({q_{i\hbox{max} }}\) is used as an upper restriction on the allowable production rate of the well.

Thus, the problem of distribution of gas extraction to reach its maximum production for the defined period of time \({T_{{\text{opt}}}}\) is formulated as an optimization problem in the following way: it is required to find such a distribution of extraction rate within producing well stock that ensures maximum value of the objective function (corresponding to the weighted-average formation pressure in the group of producing wells) with appropriate restrictions at each instants of time:

$$\left\{ \begin{gathered} \Phi (\vec {x}) \to \mathop {\hbox{max} }\limits_{{\vec {x} \in X}} \hfill \\ \sum\limits_{{i=1}}^{N} {{q_i}} =Q \hfill \\ {q_{i\hbox{min} }} \leq {q_i} \leq {q_{i\hbox{max} }},\quad \left( {i=\overline {{1,N}} } \right) \hfill \\ \end{gathered} \right.,$$
(7)

where \(\Phi (\vec {x})=\sum\limits_{{i=1}}^{N} {{P_i}{q_i}}\)—is an objective function;

$$\vec {x}=\left( {{q_1},{q_2},...,{q_N}} \right),$$

\({q_i}\)—well flow rate of an \(i\)-well,

\(X=\left\{ {\left. {\vec {x}} \right|{q_i} \geqslant 0,i=1,2,...,N} \right\} \in {R^N}\)—allowable set,

\(N\)—number of wells working at the time of development;

\({P_i}\)—reservoir pressure in the area of the \(i\)th production well;

\({q_{i\hbox{min} }}\)—minimum allowable production rate, determined by pressure level at the inlet to a gas treatment plant. At \({q_i}<{q_{i\hbox{min} }}\) \(i\)th well is shut down.

It should be pointed out that the optimization to the maximum \(\sum\nolimits_{{i=1}}^{N} {{P_i}{q_i}}\) results in a higher load of the well with the highest formation pressure.

The above-mentioned problem is a conditional optimization problem, i.e., additional conditions are imposed on the required solution (restrictions on wells). This problem is a set of \(N\) equations relative to the desired flow rates \({q_1},{q_2},...,{q_N}\) of production wells for the solution of which and based on the steepest descent method with an alterable step \({\lambda ^r}\), an iterative scheme for a certain grid area was developed. This approach enables to achieve higher accuracy results compared with a constant step, as well as faster to achieve convergence. Here, we assume that the \(\Phi (\vec {x})\)function is everywhere differentiable in the \(N\)-dimensional Euclidean space and the direction of descent in gradient methods of optimization coincides with the direction of an antigradient of the minimized function \(\Phi (\vec {x})\). An iterative formula of gradient optimization methods has the following form:

$${\vec {x}^{i+1}}={\vec {x}^i}+{\lambda ^i}{\vec {S}^i},$$

where \({\lambda ^i}\)—a step length at the \(i\)th iteration in the direction \({\vec {S}^i}\),

\({\vec {S}^i}= - \frac{{\nabla {\Phi ^i}}}{{\left\| {\nabla {\Phi ^i}} \right\|}}\)—unit vector of the direction of an antigradient function \(\Phi (\vec {x})\) at the point \({\vec {x}^i}\),

\(\left\| * \right\|\)—some vector norm, for example, Euclidean.

It is known that the gradient of the function \(\Phi (\vec {x})\) at the point \({\vec {x}^i}\) is the value of the partial derivatives vector of this function at the point \({\vec {x}^i}\):

$$\nabla {\Phi ^i}=\nabla \Phi \left( {{{\vec {x}}^i}} \right)={\left. {\left( \begin{gathered} \frac{{\partial \Phi }}{{\partial {x_1}}} \hfill \\ \frac{{\partial \Phi }}{{\partial {x_2}}} \hfill \\ ...... \hfill \\ \frac{{\partial \Phi }}{{\partial {x_N}}} \hfill \\ \end{gathered} \right)} \right|_{\vec {x}={{\vec {x}}^i}}}.$$

The scheme of the steepest descent method with an alterable step is as follows:

  1. 1.

    Choose an initial (zero) approximation \({\vec {x}^0}=\left( {q_{1}^{0},q_{2}^{0},...q_{N}^{0}} \right)\) that has no impact on the result but affects the number of iterations. The values should be chosen as the zero approximation, which were obtained from the solution for the previous point in time. Then assign an initial step size \({\lambda ^0}\) and the step fragmentation factor \(\nu\). Assume that the iteration counter is \(i=0\). Assign the step for flow rates \(\Delta q\). The indicated values ​​for \(\varepsilon\) and \(\Delta q\) were determined as a result of a numerical experiment and provide a sufficiently high accuracy of results with a relatively short calculation time.

  2. 2.

    Calculate the gradient \(\nabla {\Phi ^i}=\nabla \Phi \left( {{{\vec {x}}^i}} \right)={\left. {\left( \begin{gathered} \frac{{\partial \Phi }}{{\partial {q_1}}} \hfill \\ \frac{{\partial \Phi }}{{\partial {q_2}}} \hfill \\ ...... \hfill \\ \frac{{\partial \Phi }}{{\partial {q_N}}} \hfill \\ \end{gathered} \right)} \right|_{\vec {x}={{\vec {x}}^i}}}={\left. {\left( \begin{gathered} \frac{{{\Phi ^0} - \Phi _{1}^{0}}}{{\Delta q}} \hfill \\ \frac{{{\Phi ^0} - \Phi _{2}^{0}}}{{\Delta q}} \hfill \\ ............. \hfill \\ \frac{{{\Phi ^0} - \Phi _{N}^{0}}}{{\Delta q}} \hfill \\ \end{gathered} \right)} \right|_{\vec {x}={{\vec {x}}^i}}}\),

where \(\Phi _{k}^{0}\) is total volume of gas produced when the gas extraction is increased by the \(k\)th well at \(\Delta q\) (at the same time, it is necessary to observe the condition that the overall production is constant \({Q_{sum}}=\sum\nolimits_{{k=1}}^{N} {{q_k}}\)), \({\Phi ^0}={\Phi ^0}\left( {{q_1},{q_2},...,{q_N}} \right),\)

$$\Phi _{k}^{0}=\Phi _{k}^{0}\left( {{{q^{\prime}}_1},{{q^{\prime}}_2},...,{{q^{\prime}}_N}} \right),$$
$${q^{\prime}_i}=\left\{ \begin{gathered} {q_i} - \frac{{\Delta q \cdot {q_i}}}{{{Q_{sum}} - {q_k}}}\quad i \ne k \hfill \\ {q_k}+\Delta q\quad \quad \quad \quad i=k \hfill \\ \end{gathered} \right..$$

Calculate the vector norm \(\left\| {\nabla {\Phi ^i}} \right\|=\sqrt {\sum\limits_{{i=1}}^{n} {{{\left( {{\Phi ^0} - \Phi _{i}^{0}} \right)}^2}} }\) and the unit vector \({\vec {S}^i}\)of the antigradient direction of function \(\Phi \left( {\vec {x}} \right)\)at the point \(\vec {x}\):

$${\vec {S}^i}= - \frac{{\nabla {\Phi ^i}}}{{\left\| {\nabla {\Phi ^i}} \right\|}}.$$
  1. 3.

    By the formula \({\vec {x}^{i+1}}={\vec {x}^i}+{\lambda ^i}{\vec {S}^i}\), calculate the vector components \({\vec {x}^{i+1}}\).

  2. 4.

    Recalculate the vector components \({\vec {x}^{i+1}}\) taking into account constraints on the wells.

  3. 5.

    Calculate the value \(\Phi \left( {{{\vec {x}}^{i+1}}} \right)\) that is a value of the function \(\Phi \left( {\vec {x}} \right)\) at the point \({\vec {x}^{i+1}}\).

  4. 6.

    If the condition \(\Phi \left( {{{\vec {x}}^i}} \right) - \Phi \left( {{{\vec {x}}^{i+1}}} \right) \geqslant \nu \cdot {\lambda ^i} \cdot \left\| {\nabla {\Phi ^i}} \right\|\) is met then proceed to the next step. Otherwise, go to item 8.

  5. 7.

    Assume that \({\lambda ^i}=\nu {\lambda ^i}\) and proceed to item 3.

  6. 8.

    Check for the completion condition of calculation: \(\left| {\Phi \left( {{{\vec {x}}^{i+1}}} \right) - \Phi \left( {{{\vec {x}}^i}} \right)} \right|<\varepsilon\). If it is fulfilled, assume that the calculation \({\vec {x}^{i+1}}={\vec {x}^i}\) is completed. Otherwise, assume that \(i=i+1\) and proceed to item 3.

The solution for the optimization problem of distribution of well flow rates can be found on the basis of the above-mentioned two-dimensional hydrodynamic model with the corresponding boundary conditions at the current moment of a time period \(\Delta {T_{{\text{opt}}}}\). The received values of flow rates are given as the pre-defined mode of well stock operation for hydrodynamic calculations. The desired extremum of the sum of target functions was determined using the steepest descent method. The hydrodynamic modeling results in the values of formation pressure \(P\left( {x,y,t} \right)\)as well as the concentrations of mixture components \({z_k}\left( {x,y,t} \right)\) \(\left( {k=\overline {{1,{N_C}}} } \right)\) at each point of the stratum (the computational mesh), correspond to the optimal mode of development.

It should be noted that for natural gas reservoirs with a large number of wells and small time steps (\(\Delta t<<1\) month), the number of computational operations in solving the optimization problem is out of PC capabilities. In this connection, a special computational procedure was developed to significantly reduce the calculation time (see Fig. 1).

Fig. 1
figure 1

Procedure of conducting numerical experiments

Full size image

According to the proposed procedure, the solving for the optimum of the gas production problem for the entire period of reservoir development was presented as a sequence of optimization tasks for a period of time \(\Delta {T_{{\text{opt}}}}\) (optimization step). An optimization step is characterized by a constant level of gas production of the overal production during the whole period \(\Delta {T_{{\text{opt}}}}\). The results of numerical experiments indicate that the optimization step can be assigned as \(\Delta {T_{{\text{opt}}}}=1\) month without considerable losses in computational accuracy.

Results and discussion

Numerical experiments were carried out based on the 2D model of the natural gas reservoir Kotelevske in Ukraine. The scheme of the reservoir and the location of all wells are shown in Fig. 2.

Fig. 2
figure 2

Scheme of the Kotelevske natural gas reservoir

Full size image

Petrophysical properties of the reservoirf, such as porosity, permeability, and thickness, are inhomogeneous and are assumed to be predetermined functions of \(\left( {x,y} \right)\). Their distribution over the whole territory of the reservoir was received from the retrospective (historical) data gathered during the gas extraction (in cycling mode) from the gas reservoir Kotelevske and then processed with the specially developed computer application designed for solving the problem of multicomponent filtration.

For this purpose, a geological model of the reservoir was created based on the available data, in which the data on the effective thickness, porosity and permeability were averaged over the stratum thickness and further transferred to a finite-difference mesh of the reservoir. Then, the filtration problem of a multicomponent mixture was numerically solved for this reservoir with allowance for phase transitions (1), (5), (6) (Kalugin et al. 2015) and with given full-scale data on extraction of «wet» gas from producing wells and injection of «dry» gas into intakes wells.

Having solved the problem, we obtained the gas-condensate factor (GCF) within the entire reservoir and in each producing well. The calculated values of GCF at each production well were compared with full-scale GCF data for this well. In primary calculations, the error for different wells ranged from 5 to 30% which was caused, in our opinion, by an error in setting the geological model, mainly relative to the effective thickness and permeability.

Having analyzed the variations of GCF in these wells and taking into account the time of «dry» gas breakthrough into production wells, we calibrated the model. During the calibration, the effective thickness and/or permeability in the well area was increased/decreased based on physical considerations. Such changes were made in the places, where the greatest difference was observed. Then, the calculation of the filtration of a multicomponent mixture was carried out again. Thus, after several iterations, the rms deviation for GCF values in each producing well did not exceed 3–5%. After that it was concluded that the geological model of the reservoir was calibrated, and therefore, the optimization of the reservoir development could be carried out.

A similar calibration can be done for dry gas reservoirs by altering the wellhead pressure at each well during the retrospective period of gas development.

The source data describing the gas reservoir are presented in Table 1.

Table 1 Source data for the Kotelevske gas reservoir
Full size table

The following gas composition is adopted in the calculation: 83.63% methane, 8.10% ethane, 4.46% propane, 1.11% of butane, 0.96% nitrogen, and 1.74% of carbon dioxide.

The gas reservoir development was carried out with 18 wells during the research period.

When calculating the number of nodes, it was chosen that \({N_x}=141\), \({N_y}=67\) and the steps along the axes were equal \({\text{d}}x=58.6\) m and \({\text{d}}y=59.0\) m on the coordinate axis \(Ox\) and \(Oy\), respectively. The calculations were performed without any restrictions on the operation of the wells. In some moments of time, up to 15 wells were operating. The adopted calculation accuracy was \(\varepsilon ={10^{ - 4}}\).

To assess the impact of uniformity of the reservoir pressure drop on the gas extraction, two different variants of the reservoir development were calculated: optimal and «actual».

When calculating the optimum development regime for the reservoir, the flow rates variation affected all 18 production wells (No. 15, No. 24, No. 26, No. 76, No. 79, No. 81, No. 82, No. 85, No. 87, No. 89, No. 96, No. 98, No. 100, No. 157, No. 163, No. 164, No. 166, and No. 167), involved in the development of the reservoir. The overall calculation period embraced 170 months. The forecast calculation was carried out to select the optimal mode in which the overall gas production would be maximum, with the optimization step \(\Delta {T_{{\text{opt}}}}=1\) month. That means that the calculation of flow rate redistribution for the whole operating well stock was carried out for each month.

When calculating the «actual» mode of reservoir development, the flow rates of production wells were adopted according to the actual data of the Kotelevske reservoir. When the gas reservoir was coming into the third stage of its development (the period of significant reduction in reservoir pressure), it was impossible to keep planned flow rates of gas production for all of the operational wells, due to the fact that at \(P<{P_{\hbox{min} }}=2{\text{MPa}}\), and the pressure at the inlet to the complex gas preparation plant exceeds the bottomhole pressure. That is why some of the wells were deactivated until their pressure restored to the level that allows going on with the gas extraction in the planned mode. As can be seen from the curves in Fig. 3, it was so, e.g., with the well №15 at the time points of 123, 130, and 132 months, with the well №26 during the 87th and 94th months of development, etc. As a result of the unscheduled outages, the real gas production was significantly reduced compared to the planned one (see Figs. 4, 5).

Fig. 3
figure 3

Pressure dependence on time for some wells in the Kotelevske reservoir under optimal and «actual» regimes of development

Full size image
Fig. 4
figure 4

Total gas production from the Kotelevske gas reservoir with optimal and «actual» modes of development

Full size image
Fig. 5
figure 5

Monthly gas production from the Kotelevske gas reservoir with optimal and «actual» modes of development

Full size image

The results of the calculations are presented in Figs. 3, 4, 5, and testify that it is possible to choose the optimum alternative of gas reservoir development at its final stage with regulating the gas production from individual wells.

Fig. 6
figure 6

Pressure profiles at different points of time in well range 89-164-163-157-82-167-166-76-79-85-87-96 in the Kotelevske gas reservoir under optimal and the «actual» development modes

Full size image

The obtained optimum alternative of reservoir development allowed to increase significantly the gas recovery factor during 170 months and as a result to increase the gas production on 6.18% or 333 million m3 (5418 million m3 and 5751 million m3 in the «actual» and optimal development options, respectively).

The second stage of the reservoir development (permanent gas production) was extended for 61 months in the optimum alternative, and the beginning of the final development stage was delayed from the 64th months to the 123rd (when the pressure in the 157th well got less than the pressure at the inlet to the complex gas preparation plant).

Figures 4 and 5 present the curves of the total and monthly gas production from the Kotelevske gas reservoir for different development alternative. In fact, the monthly gas production for every time point is higher in the optimum alternative than in the «actual» one. It explains the permanent additional increase in total gas production in the optimum alternative.

Figure 6 illustrates the impact of gas extraction regulation on pressure field when the flow rates of production wells are specified in the course of solving the optimization problem. As an example, this figure presents pressure profiles of several sectors of the reservoir, situated between the wells No. 89, No. 164, No. 163, No. 15, No. 82, No. 100, No. 98, No. 83, No. 87, and No. 96 (30, 60, 90, 130, and 160 months after the reservoir development started).

Fig. 7
figure 7

Minimum pressure in the Kotelevske gas reservoir for optimal and «actual» alternatives

Full size image

It is obvious that the pressure profile for the optimum alternative is smoother than for the «actual» one. In this case, the cones of depression of production wells are smaller and, as a result, they work for a longer time at the final stage of reservoir development. The range of pressure differential in the whole reservoir is less and the pressure over the most area of the reservoir is higher for the optimum alternative than for «actual» one (see Fig. 7). The optimization of reservoir development thus leads to the pressure equalizing inside the area and, as a result, to the stabilization of recharge area of each well.

As an example, Fig. 8 presents the plots of monthly gas production for 6 (No. 15, No. 26, No. 81, No. 82, No. 85, and No. 167) of the 18 operating wells that are involved in development process. The results of the calculations reveal the fact that most operating wells work in nonoptimal mode during the «actual» gas development. In most wells, «actual» extraction exceeds significantly the optimal one. As a consequence, the pressure differential between untouched reservoir region and production wells rises which results in an increase of interformational flows and then leads to the temporary deactivation of some wells and their withdrawal from the development process. Just this fact suggests broad potentials for application to increasing gas production via optimization by well flow rates control.

Fig. 8
figure 8

Monthly gas production for development wells No. 15, No. 26, No. 81, No. 82, No. 85, and No. 167 for «actual» and optimal alternatives of reservoir development

Full size image

Consequently, the obtained results of numerical modelling offers great opportunities for using the proposed model and optimization programs to manage gas development from natural reservoirs for the purpose of increase in gas production.

Conclusions

An optimization method has been proposed to increase gas production from natural gas reservoirs. The method is based on the mathematical model for 2D filtration of a two-phase multicomponent hydrocarbon mixture and has for an object the redistribution of flow rates in production wells. An own computer program was composed to do necessary calculations for this method. The basis for the computer program is the own algorithm that allows optimum managing the gas extraction process and takes into account the restrictions on the production rates in some operating wells, as well as their shutoffs for remedial maintenance.

Efficiency of the proposed optimization method has been illustrated by data from the natural gas reservoir Kotelevske (Ukraine). It is demonstrated that the pressure differential inside the gas reservoir decreases owing to the well flow rates control. This resulted to the 6% increase in gas production over the short time interval.