Definition of the reservoir permeability field according to pressure measurements on wells with the use of spline function

The problem of reservoir permeability identification based on known well pressures under conditions of single-phase fluid filtration is considered in the article. The permeability field is identified in the spline function class from the solution of the inverse coefficient problem for the filtration equation. The problem of identification is reduced to the problem of minimizing the residual function, having the form of a sum of squares of the difference between the pressure values known from measurements at the wells and obtained with the help of a numerical model. Minimization of the residual function is carried out by the Levenberg-Marquardt method. The solutions of model problems of permeability identification for a two-dimensional reservoir, penetrated by a system of production and injection wells, are presented. The calculated permeability fields are close to the true fields. The example of a problem with errors in pressure measurements shows the stability of the solution.


introduction
In solving tasks of single-phase fluid filtration, field of the reservoir permeability must be known. In practice, permeability is usually determined only at individual points of the reservoir by hydrodynamic methods or by the results of core studies. In this paper, the permeability field is identified in the class of spline functions from the solution of the inverse coefficient tasks for the filtration equation. Methods for solving inverse coefficient tasks are divided into explicit and implicit tasks (Sun, 1994). In explicit methods, the values of the parameters are determined from the solution of the nonlinear system of equations (Golubev et al., 1978;Zinoviev, 1984). In this case, the pressure field must be known. If pressure values are known only at the wells, the parameter values are determined only in the wellbore regions.
In implicit methods, iterative procedures are used to determine the values of parameters of the entire reservoir, in which only pressure values are used in wells (Neuman, Carrera, 1986;Khairullin et al., 2006;Elesin et al., 2009;Khairullin et al., 2017). Unlike the proposed approach, the reservoir parameters are represented in the form of piecewise-constant functions.
The use of spline functions is one of the parametrization methods (Sun, 1994), and in this case for the obtained permeability field, no further processing and correction is required in most cases and it can be easily and unambiguously converted to any grid. To solve inverse tasks, in addition to the data needed to solve a direct hydrodynamic task, additional information is used. In the present work, such information is taken to be the bottomhole pressure, which is known along with well production rates. The flow rate is used to solve a direct task for determining the pressure field. Known wellhead pressures are included in the residual function, in the process of minimization of which the field of reservoir permeability is determined. To minimize the residual function, the classical Levenberg-Marquardt method is used. The stability of the solution to the errors in pressure measurements is investigated.

Formulation of the task
Single-phase stationary filtration in a two-dimensional layer is described by the equation (Aziz, Settari, 1982;Basniev et al., 1986) where s = kh/µ -coefficient of hydroconductivity, kpermeability, h -thickness of the reservoir, µ -viscosity of the liquid, p -pressure, , Q i , (x i , y i ) -flow rate and coordinates of the i-th well, M -number of wells, d(x i , y i ) -delta function. For equation (1), the boundary conditions are given: ( 2) where Г 1 +Г 2 =Г -boundary of the reservoir W, -the normal vector to the reservoir boundary, wnormal component of the filtration rate. Equation (1) with boundary conditions (2) is solved numerically. To approximate the spatial variables, the method of control volumes on a rectangular grid is used. The resulting system of linear algebraic equations is solved by the method of conjugate gradients with preconditioning as an incomplete Kholesky decomposition (Golub, Van Loan, 1999;Hill, 1990;Larabi, De Smedt, 1994).
The determination of the pressure field from solution (1)-(2) is a direct task. The inverse task is to determine permeability values in all control volumes at known pressure points at individual points. To obtain a unique solution of the inverse task, it is necessary that the number of identifiable parameters does not exceed the number of known pressure values. Since in solving practical tasks the number of control volumes covering the design area (layer) is much larger than the number of known pressures, two approaches are commonly used to reduce the number of identifiable parameters. In the first approach, the calculated region is divided into zones, each of which is characterized by a constant value of permeability. The second approach uses different interpolation options.
First, the permeability values at the interpolation nodes are determined, and the remaining values are obtained by interpolation throughout the calculation area. In this paper, unlike the second approach, in the course of solving the inverse task, the field of permeability in the form of a spline function is directly restored (Ashkenazy, 2003;Harder, Desmarais, 1972), the number of determining parameters of which corresponds to the number of wells with known bottomhole pressure.

interpolation by a spline function
Let us assume that the values of a certain quantity a i are known at the points P i (x i ,y i ), i=1, n, of a twodimensional region. The task of interpolation is to construct a spline function φ(x,y) (spline surface) defined on the whole region so that its values at the points P i coincide with the values a i . Points are called interpolation nodes. Interpolation by a spline surface has a simple mechanical meaning. The spline surface is a model of an elastic thin plate bent under the action of external forces applied at points P i . Finding such a spline surface is to solve the variational task of finding the minimum free energy of a thin plate. The splinesurface equation has the form: , ( where . To determine coefficient c i , i=1,n+3, of spline function φ(x,y), it is necessary to solve the system of equations: For n > 3 there is solution to this system, and the solution is unique, if among the points (x i ,y i ), i=1,N, there are at least three points not lying on one line (Ashkenazy, 2003).

Method for solving the inverse task
A frequently used method for solving inverse coefficient tasks is to reduce them to the tasks of minimizing the residual function (Sun, 1994). In this paper, the residual function is constructed from known pressure values and has the form, , where K -the control vector, arguments of which , k i -permeability values at the spline interpolation nodes, p j , p j * -pressure tasks, obtained as a result of solving equation (1), and known from the well measurements, M -the number of known pressure values.
The minimization process is carried out in two stages. At the first stage, the permeability of the entire reservoir is considered constant and is determined in the process of minimizing the residual function by the gradient method (Panteleev, Letova, 2005). At each iteration of the gradient method, the permeability is recalculated using formula: , where g -gradient of the residual function, pitch r is determined by the golden section method. The resulting permeability value is then used as the initial value. At the second stage, the permeability values at the interpolation nodes are determined in the process of minimizing the residual function by the Levenberg-Marquardt method (Aziz, Settari, 1982;Dennis, Schnabel, 1988;Panteleev, Letova, 2005). The new parameter values at each iteration of the Levenberg-Marquardt method are calculated by the formula: , where E -the unit matrix, h = A T A -the approximate matrix of the second derivatives, -the sensitivity matrix, µ n -the Marquardt parameter, n -the iteration number. The initial value of the Marquardt parameter is chosen one order of magnitude larger than the maximum singular number of the matrix h. If the residual function is decreased at the current iteration J(K n ) < J(K n-1 ) the Marquardt parameter is halved, if the descending condition is violated, the Marquardt parameter is doubled until this condition is satisfied. Then, a new iteration is performed. Elements of the sensitivity matrix are computed numerically. The process of minimizing the residual function was stopped by fulfilling one of two criteria: the achievement of a given accuracy with respect to pressure measurements or the slow convergence of the minimization process J n -J n+1 < 0,01J n during 3 iterations.

Model tasks
In the model tasks, the exact solution is always known. This allows us to test the solution methods and evaluate the sufficiency of the initial data to obtain an exact solution. The model tasks of reservoir permeability identification are constructed as follows. First, the formation points (interpolation nodes) are selected in which permeability values are specified. From these values, the spline function (3) is constructed and the permeability field of the entire reservoir is calculated, which is taken as the true permeability field. Then, from the solution of Eq. (1) using the Peaceman formula (Peaceman, 1978), the pressure values at the wells are determined. After that, it is assumed that the permeability values are unknown and it is required to determine them in the process of minimizing the residual function (4) from known well pressures.

The model task 1
A rectangular reservoir with dimensions of 2000 m × 2000 m with a capacity of 10 m, opened by 5 injection wells and 20 producing wells is considered. The radius of the wells is 0.1 m. A pressure of 20 MPa is set at the boundary of the reservoir. The viscosity of the liquid is 10 mPas. The coordinates of the wells, their production rates and the given permeability values ktr are given in Table 1. The coordinates of the interpolation nodes of the spline function coincide with the well coordinates (Fig. 1). To approximate the filtration equation (1) with respect to spatial variables, the layer is covered by a square grid with a step of 40 m (2500 reference volumes).  To study the stability of the solution, random errors εi were introduced into the pressure measurements. The results of solving the model task without error and with errors in pressure measurements are given in Table 1. The solution k 1 was obtained without error in pressure measurements (ε i = 0 MPa), solutions k 2 , k 3, k 4 were obtained in tasks with errors in pressure measurement |ε i |<0,1 MPa, |ε i |<0,01 MPa, |ε i |<0,001 MPa respectively. When solving tasks, both with an error and without an error in measuring pressure, a predetermined accuracy of d = 0,01 МПа is achieved. The calculated permeability values without error in pressure measurements are close to the true ones (Fig. 2).
From the results of the solution of tasks with errors in pressure measurements, it can be seen that with decreasing error, the values of the parameters approach the true values. Note that the maximum relative deviations in model task 1 are observed at the interpolation nodes with the maximum permeability values. This is explained by the fact that at the same values of flow rate, the value of bottomhole pressure at the well is more sensitive to a change in permeability at its small values, which is clearly seen from the Peaceman formula. Therefore, in order to achieve a given pressure accuracy, the small permeability values should be closer to the true values in comparison with the larger permeability values.

The model task 2
In the first model task, the true and calculated permeability fields were determined from the same nodes of interpolation of the spline function. The second model task differs from the first in that the true field of permeability was constructed from the interpolation nodes located at the nodes of the square grid with steps of 1000 m (9 knots). The coordinates of the interpolation nodes in solving the inverse task coincided with the well coordinates. The true permeability values obtained in the result of the solution of the second model task on wells without an error in pressure measurements are given in Table 2.
The corresponding permeability fields are shown in Fig. 3, 4. When solving the task, a predetermined accuracy of d = 0,01 MPa from measurements of pressure at the wells was achieved. The calculated permeability field is close to the true field.  The model tasks of identifying the permeability of a two-dimensional formation, opened by a system of producing and injection wells, were determined based on pressure measurements at wells under conditions of stationary single-phase fluid filtration. The permeability field was approximated by a spline function, constructed from the values at the wells. Unknown permeability values at the wells were determined in the process of minimizing the residual function using the Levenberg-Marquardt method. When solving model tasks without errors in pressure measurements, the calculated permeability fields practically coincide with the given fields. For the task with errors in pressure measurements, the stability of the solution obtained is shown.