NUMERICAL METHODS FOR SOLVING IMPROPER PROBLEMS OF FILTRATION THEORY

This paper is devoted to the development and investigation of methods of mathematical and computer simulation of the process of fluid filtration in a porous medium. The methods of numerical solution of the problems of the filtration theory of build-up of conditions in the catchment and discharge areas boundaries, identification of filtration-capacitive parameters of the effective formation and determination of free (unknown) boundaries and creation of computational algorithms for analysis and forecast of technological indicators of oil and gas fields are considered. Methods and models of continuum mechanics, filtration theories, and methods for solving ill-defined problems, numerical modeling and computer programming were used. Approximate Numerical methods for solving direct and inverse problems of filtration theory, mathematical models for single-phase isothermal filtration of a gas mixture in a horizontal formation at small concentration gradients of components, studying the properties of self-similar solutions, as well as numerical solving the problem of identifying the capacitive parameters of the water-bearing stratum.


INTRODUCTION
Oil and gas deposits are confined to formation water-drive systems. In this case, a natural filtrational flow of water usually takes place. The features of the geological structure of the water-bearing stratum, the change in the area of collecting properties, the presence both catchment area and discharge area has a significant effect on the degree of manifestation of the water-drive regime. These factors determine the nature of the interference of the wells, uncontrolled movements of gas or oil, etc. In recent years, more and more attention has been paid to studying the effect of natural (or artificially created) filtrational flow of water on the configuration of gas-water, water-oil contacts. A special role belongs to the natural filtrational flows and the conservation of gas and oil deposits in hydrodynamic traps. The opposite problem turns out to be significant: establishing, for example, conditions in the catchment and discharge areas, as well as identifying the filtration-capacitive parameters that led to the displacement of the hydrocarbon deposits in the stream of the formation water or the presence of the deposits in the hydrodynamic trap. Suppose a gas deposit is confined to a hydrodynamic trap. It is required to predict its behavior in the process of development and the manifestation of the water-drive regime. This means that it is necessary to set the initial condition in the form of the initial distribution of pressure in the water-bearing stratum, hence the pressure in the catchment and discharge areas boundaries. It turns out that the problem arises here of the inconsistency of the initial condition with the real gas deposit. Thus, if the initial condition is incorrectly specified for gas sampling, for example, with infinitely small rates, a significant defor-mation of the gas-water contact may occur, which contradicts the physics of the process. Elimination of such misunderstandings leads to the need for numerical identification of the parameters of the water-drive basin. This task is very multifaceted. Let us consider only the problem of build-up of pressures in the catchment and discharge areas, as well as identifying the filtration parameters of the water-bearing stratum. Even this, one of the simplest tasks, is characterized by a number of difficulties. First, the location of the catchment and discharge areas is not always known with certainty. Secondly, it is difficult to indicate the outer boundaries of the water-bearing stratum. Thirdly, many features of the geological structure of the basin under consideration are not clear. Therefore, it is necessary to talk about the creation of an equivalent computed model of the water-bearing stratum for a given deposit. A part of the water-bearing stratum is allocated around the deposit, possibly, along which there is a certain amount of initial information. This means that, with some degree of conventionality, certain formation boundaries are characterized as impenetrable and the catchment and discharge areas boundaries are distinguished. There are (basically approximate) maps of equal values of the coefficient of permeability k and the formation thickness h. We also have the values of pressures and filtration parameters in a series of exploratory or piezometric wells, including at the deposit boundary prior to its development. Often, the pressure on the catchment and discharge areas is unknown, as well as the actual distribution of the filtration parameter. Therefore, the following inverse problem is relevant. It is required to find such pressures in the selected catchment area P* n and discharge area P* p boundaries, as well as the actual distribution of the filtration parameter b ф (x, y), which predetermined the observed values of pressures and filtration parameters at different points in the water-bearing stratum and the established configuration of the gas deposit.

LITERATURE REVIEW
This also determines the relevance of the scientific research devoted to solving the problems of identifying of the filtration-capacitive parameters and establishing conditions in the catchment and discharge areas boundaries of the water-bearing stratum. Inverse problems for formation system are in most cases incorrect. They are characterized by instability in the initial data, expressed in the fact that small variations in the initial data can cause significant variations in the solution. The foundations of the theory of the approximate solution of ill-defined problems are laid down in the works of A.N. Tikhonov [1]. Basic principles for solving ill-defined problems are developed by A.N. Tikhonov [1], M.M. Lavrentyev [2], G.I. Marchuk [3], V.G. Romanov [4] et al. First applied to the inverse problems of the filtration theory in the works of W.H. Chen, K. Coats, G.R. Gavalas, J.H. Seinfeld, M.L. Wabermann a variational approach is described. At present, algorithms for solving inverse two-and three-dimensional, single-and multi-phase problems of filtration theory have been created, which on an impartial basis allow refining the collecting properties of the productive stratum with the use of actual development data. For these purposes, the time dependences of flowrates of liquid, gas and formation pressures on wells are needed. Practically there are no similar studies for stationary filtration formations. A wide range of studies devoted to the solution of the problem of identification of formation parameters in two and three-dimensional cases for single-phase and multiphase filtration are presented in the works of Yu.K. Alexeev [5], D.M. Akhmetzyanov [6,7], V.Z. Baishev [8], V.Ya. Bulygin [9], A. Zholdasov [10], B.M. Palatnik [11,12] and other authors . For the first time, methods for refining the collecting properties of the formation were based on conducting of research of individual wells with established and unsteady filtration regimes. In the first case, the indicator lines (diagrams) were removed from the wells, in the second case, the curves, the coefficients of the filtration resistances in the equation of the flow of gas (liquid) to the well in case of violation of the linear filtration law (Darcy's law), the hydraulic conductivity parameter, and the coefficient of the formation-pressure conductivity. This means that on the basis of special well studies, when interpreting the obtained results, the parameters related to the individual wells have been refined. In recent years, the progress made in the formulation and solution of inverse problems, allows to solve multidimensional inverse problems of the theory of development of oil and gas fields. In this case, actual information on changes during the last period of development of well flow rate, formation pressure in them, and coefficients of gas, oil, and water saturation of the formation in separate operational and observing wells is used. In this case, for example, the problem of identifying of the filtration-capacitive parameters (J) of a formation reduces to the task of minimizing of some objective functional of the form: where k, m -the required filtration-capacitance parameters: P(x i , t) -estimated pressure, obtained as a result of solving a direct boundary value problem, describing the filtration, in the formation: P r (t) -real (actual) pressure in the wells; T -time period of observation; x -coordinates of the point at which the pressure measurement is carried out.  [14] to minimize the functional are used second-order methods.

METHODOLOGY
In this study, an algorithm for solving direct and inverse problems of refinement of technological parameters of the effective formation and boundary conditions in the catchment area boundaries of the well is considered. Two approaches to the formulation of inverse problems on the specification of technological parameters of the effective formation and boundary conditions on the catchment area boundaries are also explored. Such tasks are incorrect. For their solution, iterative algorithms are constructed. If from the beginning of field development in the formation water-drive systems to which they are confined presents a natural filtrational flow of water, then it is characterized by pressure drop when water moves from the catchment area to the discharge area. Certainly, the pressures distribution over the area of such a water-drive system is entirely determined by the filtration parameters of the collectors, the tectonic structure of the water-bearing basin, and the presence of oil and gas deposits [15]. Algorithms for solving inverse problems are based on algorithms for solving direct boundary value problems. Therefore, the study deals with statements and algorithms for solving both direct and inverse problems of the filtration theory. Consider the formulation of a direct boundary-value problem.

Figure 1: The isobar map and the layout of exploratory wells in the heterogeneous by its collecting properties water-bearing stratum
Suppose: G -arbitrary area R n with border Г 0 and Г1: more precisely: Г i -is a continuously differentiable manifold of dimensionality n−1, G is on one side of Г i , i = 0,1. Suppose A -elliptic differential operator in G of the second order: where: , Filtration process of water in the heterogeneous by its collecting properties water-bearing stratum with an isolated gas (oil) reservoir is described as follows by a differential equation of elliptic type relating to the reduced pressure P*: To solve the direct problem of interest to us, the integration of equation (1) is carried out under the following boundary conditions: Thus, for a natural filtrational flow in the water-bearing stratum with a separated gas (oil) reservoir, the direct boundary value problem (3) -(5) holds, where -P*= P ± ρ B gl ; P -pressure at the point with the coordinates x and y; ρ B -density of water; g -acceleration of gravity; l -vertical distance from the given point with coordinates and up to datum plane; n -outer to G normal; k -coefficient of formation permeability; h -formation thickness; Г 1 -the boundary of the gas (oil) deposit; Г 2 -impermeable boundaries of the water-bearing stratum; Г 3 -catchment area boundary; Г 4 -discharge area boundary. (Figure 1.) Equation (3) describes the process of stationary filtration. The boundary condition (4) takes into account the impermeability of the outer boundaries of the formation and the gas, oil deposits (impenetrability for the flow of water prior to its development). The condition (5) where m -coefficient of porosity; C -concentration of the component under consideration; λ 1 , λ 2 -respectively longitudinal and transverse scattering parameters among constants (have the dimensions of length); v 1 ,v 2 -the components of the filtration rate respectively along and the axes. If for the equation ∆ϴ and the initial condition C (x,y)=1, t=0 (8) then we obtain a direct boundary value problem for determining the concentration field over the entire area of the water-bearing stratum at different times, including at the time of calculations (to date) Points (wells) -information sources are located randomly across the entire area of the water-bearing stratum (in Figure 2).

Figure 2: Iso-concentration map and the layout of exploratory wells in the water-bearing stratum heterogeneous by its collecting properties
The filtration rates and depend on the pressure distribution P * over the area. Therefore, in order to solve the problem (6) - (8) it is necessary to have a solution of the next problem: There is no inverse dependence P * from C, therefore, the above problems (9) -(11) and (6) -(8) split into two autonomous direct boundary-value problems. Thus, under given collecting properties of the formation k and h, λ 1 , λ 2 under boundary conditions, the solution of the equation makes it possible to determine the distribution of the reduced pressure P * over the entire area of the water-bearing stratum. Further, under the given collecting properties of the water-bearing stratum under the initial condition (8) and the boundary conditions (7), the solution of equation (6) using the function P * from the previous problem allows determining the concentration distributions C over the entire area of the aquifer. Statement of the inverse problem. We assume that the actual values P * (x,y) and C(x,y,T) are known over the entire area of the water-bearing stratum at the end of the time interval [0,T ] at the present. Have P * and C in the wells at the period of geological time T, data on the geometry of the formation and approximate values of its collecting properties. It is necessary to determine (specify) the collecting (filtration-capacitive) properties in all points of the water-bearing stratum, as well as the conditions on the catchment and discharge area boundaries on the basis of available actual data. We will solve the inverse problem as an optimization problem. Construct a functional J equal to the sum J 1 and J 2 . Obviously, this functional depends on the reservoir properties of the reservoir: The inverse problem is as follows: find the values P c * , P d * , b, m, λ 1 , λ 2 which minimize the functional (12). As a result of solving this problem, the collecting properties of some equivalent geological model of the formation are determined, ensuring the best match of the calculated and actual values of the reduced pressures, filtration parameters and concentrations from the point of view of the introduced criterion and the accepted mathematical model of filtration. Method for solving the inverse problem. To minimize the functional, we use the iterative gradient method, for example, to refine the porosity coefficient at different points in the formation, we have the following recurrence relation where S -iteration number. Similarly, gradient procedures for all other required parameters are written. As a zero approximation for the values P c * , P d * , b, m, λ 1 , λ 2 geological maps, data of geophysical and hydrodynamic studies of wells and strata can be used. To implement gradient type procedures it is necessary to be able to calculate the functional derivatives, the derivation of formulas for which is given below. Now we turn to the derivation of formulas for the expressions of functional derivatives.

Derivation of expressions for functional derivatives
Let us proceed to the derivation of formulas for the calculation of functional random variables. The required parameters of the water-bearing stratum are P c * , P d * , b, m, λ 1 , λ 2 . As noted above, the direct problem of finding solutions P * and C and splits into two autonomous tasks: -the problem of finding a function P * as a solution of equation (9), satisfying the boundary conditions (10)- (12) -the problem of finding a function C as a solution of equation (6) satisfying the initial condition (8) and the boundary conditions (7), using the function P * from the previous problem. To solve the inverse problem, the residual functional is chosen in the form (12) In Section 1 of this study, we considered the inverse problem of restoring conditions on the catchment and discharge areas boundaries and identifying the filtration parameters of the water-bearing stratum with the criterion of optimization J 1 . In the same place, formulas were obtained for calculating the functional derivatives .
Therefore, we restrict ourselves to deriving formulas for calculating the functional derivatives from dJ 2 according to required parameters. In obtaining expressions for the functional derivatives, use the technique given in [18][19][20][21][22][23][24]. The definition of the derivatives of the functional J 2 (m, λ 1 , λ 2 ) is based on obtaining the decomposition of the variation dJ 2 in the independent variations dm, dλ 1 , dλ 1 [25][26][27]. We decompose the variation in independent variations.
Coefficients with independent variations and will represent the desired functional derivatives according to refined parameters. We use the property of the Dirac delta function and write the expression for the discrepancy ɛ i where ∫ G C(x, y, T)dG = C ppa i (T) -the calculated concentration corresponding to the locations of the i -th well at the time point; T; ɛ i -the discrepancy between the calculated and actual concentrations in the i -th well at the time point T. Since the calculated concentrations in the wells are the result of the solution of the problem (6) -(8) for the specification of (known) parameters m, λ 1 , λ 2 , then the value ɛ i depends on these parameters. As noted above, the minimum of the following functional of the total quadratic residual, calculated and actual concentrations, is taken as the optimization criterion for finding the parameters of the water-bearing stratum (m, λ 1 , λ 2 ): where γ i -weight function (multiplier).
Let us write down the expression for dJ 2 , using formulas  (15) The validity of the equality obtained follows from the fact that the variation from the actual concentration is zero. We reverse summation and integration. As a result, the expression under integral sign will have a factor We introduce the following boundary-value problem, which is called conjugate, with respect to the function (conjugate function). This function can be interpreted as a result of a perturbation caused by the action of point sources of runoff from the density of selection of which is proportional to ɛ i [29,30]. Thus, it is required to solve the adjoint differential equation under the following boundary and initial conditions: Substituting the expression for the right-hand side of (19) in (21), we obtain: Further transformations (22) require the representation of the derivatives with respect to the spatial and temporal coordinates of the variation δC(x, y, t) through the variations of the sought-for parameters dm, dλ 1 , λ 2 . Therefore, it is necessary to vary equation (4), have When the equations (6) are varied, neglecting terms of higher than the first order of smallness relating to variation.
Variations of the initial and boundary conditions (7), (8) have the form Since the first and second terms on the right-hand side of (25) are of the same form, the detailed calculations are carried out only for the first term. This term is twice Here the indexes "h" and "b" denote the lower and upper limits of the integration region along the corresponding coordinate. By condition (25), the first term in (26) is equal to zero. Consequently Due to condition (24), the first term in (27) is equal to zero. We transform the third term on the right-hand side of (22): Expanding the brackets under the integral, we integrate by parts the second and third terms on the right-hand side of (30): By condition (7), the first term on the right-hand side of (26) is equal to zero. As a result, we obtain from the expression (30) the following T  T  T  2  2  2  1  2  2  2  1  2  2  2  2  2  2  2  1  2  1  2  1  2 T  T  2  2  2  2  1  2  2  1  1  2  2  2  2  2  2  2  2  2  1  2  1  2  1  2  1

RESULTS AND DISCUSSION
The table shows the distributions of P 0 * , P 3 , P 9 , P 24 and P a * by wells, illustrates the change in the values of reduced pressures on the contours of the supply and unloading areas, as well as the values of the target functional , depending on the number of iterations when solving the inverse problem. Average convergence is achieved after 14 iterations, and the solutions at the nodes corresponding to the wells have a small error. Studies of the problem of identifying filtration parameters and restoring conditions on the contours of the supply and discharge areas showed the following main conclusions: 1. When applying the gradient descent method with step splitting, there is a rapid convergence to the exact solution on the average. As the number of iterations increases, the solution converges exactly uniformly. 2. The solution of the problem does not significantly depend on the errors of measurements of the reduced pressures and the filtration parameter on the wells, and if these errors are random, then the convergence on average and uniform convergence across the wells is preserved. It can be concluded that if the model used corresponds to the object under study, then the optimization formulation of the identification and restoration problem is adequate: with a decrease in the residual functional of the calculated and measured reduced pressures and filtration parameters in the wells, the obtained distributions b are approximated to the true ones, as well as the conditions on the contours of the regions supply and unloading.  The above research results should not be absolutized, since they relate to well-defined initial data, to a specific number of wells -sources of information and their placement in the area. Nevertheless, they qualitatively fully reflect the features of the solutions of the studied inverse problems. And most importantly, the performed mathematical experiments quite convincingly, from our point of view, traced the practical suitability of the proposed formulations and algorithms for solving the corresponding inverse problems. Studies of the problem of identifying filtration parameters and restoring conditions on the contours of the supply and discharge areas showed the following main conclusions.

Iteration
1. When applying the gradient descent method with step splitting, there is a rapid convergence to the exact solution on the average. As the number of iterations increases, the solution converges exactly uniformly. 2. The solution of the problem does not significantly depend on the errors of measurements of the reduced pressures and the filtration parameter on the wells, and if these errors are random, then the convergence on average and uniform convergence across the wells is preserved. It can be concluded that if the model used corresponds to the object under study, then the optimization formulation of the identification and restoration problem is adequate: with a decrease in the residual functional of the calculated and measured reduced pressures and filtration parameters in the wells, the obtained distributions b are approximated to the true ones, as well as the conditions on the contours of the regions supply and unloading.
The above research results should not be absolutized, since they relate to well-defined initial data, to a specific number of wells -sources of information and their location on the area. Nevertheless, they qualitatively fully reflect the features of the solutions of the studied inverse problems. And most importantly, the performed mathematical experiments quite convincingly, from our point of view, traced the practical suitability of the proposed formulations and algorithms for solving the corresponding inverse problems. The use of the obtained formulas for functional derivatives in the procedure of the gradient method of minimizing the residual functional makes it possible to refine the collecting properties of the formation. The convergence by functionals and convergence of reduced pressures on the contours of the discharge supply areas to the actual ones depending on the number of iterations when solving the inverse problem. Table 1 illustrates the change in the values of the reduced pressures on the catchment and discharge areas boundaries, as well as the values of the objective functional depending on the number of iterations in solving the inverse problem (also Figure 3). Table 2 shows the convergence by the functionals and and the convergence of the capacitive parameters to the actual one, depending on the number of iterations in solving the inverse problem (also visualized in Figure 4). The convergence by the functionals and the convergence of the capacitive parameters to the actual one, depending on the number of iterations in solving the inverse problem.

CONCLUSION
Let us formulate the conclusions concerning the properties of the numerical solution of the problem of identifying capacitive parameters. The solution of the problem of identifying capacitive parameters does not essentially depend on the choice of the initial approximation. There is rapid convergence to an exact solution. It is possible to draw the following conclusion: when the objective functional of the residual of the calculated and measured concentrations in the wells decreases, the approximation obtained to the true takes place.