Yanbin Sun, Zhibin Liu, Yanjie Hu

School of Science, Southwest Petroleum University, Chengdu, 610500, China

ABSTRACT In this paper, we study the state estimation of compressible single phase flow in compressible porous media.The initial pressure distribution is estimated according to discrete adjoint approach based on the collected well pressure data.The first-order Tykhonov regularization method is used to obtain reasonable estimation.By analyzing the optimality condition of estimation problem, the discrete adjoint state equation and discrete adjoint gradient are derived based on the numerical scheme of the continuous equations.A quasi-Newton numerical optimization method related to adjoint gradient is proposed to solve the estimation problem.The estimation results with different regularization coefficients are compared and analyzed by numerical experiments.The deviation between the estimated pressure obtained without regularization and the real pressure is large.Estimation result with smaller deviation and higher smoothness can be obtained through appropriate regularization coefficient.When the observation error is large, the observed values generated by the estimated pressure fit well with the real pressure.

Keywords:Regularization Discrete adjoint approach State estimation Single phase compressible flow Porous media

1.Introduction

Fluid flow in porous media is an important model in many engineering fields.This model uses partial differential equations to describe the variation of state variables such as pressure in porous media.In order to obtain the solution of the equation, the initial condition, boundary condition, parameter value of the equation is necessary, but in many cases, the values of these variables are not easy to achieve.On the other hand, some data can be collected at some point in space, which partly reflects the value of the state variable or parameter on the whole space.We can estimate the distribution of states or parameters across the space utilize these collected data.

There are many ways to use partial observations to estimate state variables or parameters.In the ocean and atmosphere dynamic field, the representer method is a typical method for estimating spatial states or parameters using observed data[1-3].This method performs the optimization in the subspace spanned by collected data and effectively reducing the dimension of the problem.The representer method can also be used to estimate the parameters of porous media based on Darcy flow [4,5] and updating of reservoir parameter [6,7].However, the representer method is typically applicable to the case where the system equation describing the flow is linear.

Alternative method for estimating parameters or state variables using observational data is the ensemble Kalman filter method which is a Monte Carlo random sampling method.This approach uses multiple models to fit observation data and has the advantages of sequence estimation and quantitative analysis of model uncertainty.Ensemble Kalman filtering method is typically used to estimate parameters and states in oceans dynamic models or geological models [8-10].However, this approach requires the generation of multiple system models, which consume a large amount of computing resources.

Another way to estimate the state or parameter of porous media flow model is the adjoint gradient method.In this approach, the estimation problem is expressed as a least squares optimization problem constrained by the system model.The adjoint form of the model is derived according to the variational technique and then the adjoint gradient of the objective function with respect to state or parameter variables is obtained.Adjoint method is a more efficient technique for generating sensitivity of objective function[11].Sun et al.[12,13]research the parameter estimation of groundwater model.The continuous adjoint equation of groundwater model is deduced and the adjoint gradient Algorithm is given.Wu et al.[14]use the discrete adjoint method for generating sensitivity coefficients related to two-phase flow production data.Li et al.[15]developed the adjoint equations for three-phase flow problems and implemented to calculate the sensitivity of production data to permeability fields and well skin factors.Fu et al.[16] use the adjoint method to compute high-resolution sensitivity coefficients for subsurface flow in large-scale heterogeneous geologic formations.Lee et al.[17,18] estimate the absolute permeability in twophase petroleum reservoirs from noisy well pressure data based on the theory of regularization and on spline approximation.Sarma et al.[19] use adjoint method as a means to estimate reservoir parameters and optimize production to realize closed reservoir management.Bukshtynov et al.[20] use an approximate convex programming method based on adjoint gradient to estimate reservoir parameters and optimize production.Oliver et al.[21]summarizes different methods of parameter estimation in porous media.Jansen[22]summarizes the different applications of adjoint gradient method in porous media flow.

However, due to the lack of data in practical application, the direct application of adjoint method may lead to incorrect estimation results.Thus a regularized adjoint gradient method is proposed to estimate the state of nonlinear Darcy flow model in porous media in this paper.The influence of regularization term on estimation results is studied.It is assumed that the system flow model is accurate and the initial pressure variable is needed to estimate.By analyzing the optimality conditions of the estimation problem, the discrete adjoint equation and adjoint gradient are derived.We use the first-order Tykhonov regularization method to regularize the problem to avoid ill-posed of the estimation problem.A quasi-Newton optimization Algorithm is proposed to solve the estimation problem.The feasibility of the algorithm is verified by numerical experiments.Regularization coefficients affect the accuracy and smoothness of estimation results.

This paper proceeds as follows.In section 2, the regularization form of state estimation problem is described.The optimality condition of state estimation problem based on the discrete form of the system equations is analyzed.The adjoint equation and numerical optimization Algorithm for the estimation problem are presented.In section 3, numerical experiments based on homogeneous porous media and heterogeneous porous media are implemented respectively.The influence of regularization term on estimation results is analyzed.Concluding remarks are provided in Section 4.

2.Methods

2.1.System model of compressible flow in porous media

The single phase Darcy flow model of compressible fluids in compressible porous media is often used to describe the flow of groundwater or oil in a rock reservoir.According to the conservation law of matter, the flow process can be represented by the following partial differential equation [23]:

Here ρ is fluid density, φ is porosity of porous media, v is fluid flow velocity,andqis recovery rate which represents the volume of fluid produced per unit of time.The flow velocity of fluid in porous media can be modeled by Darcy’s law.

HereKis the permeability of porous media, μ is the fluid viscosity and is a constant function in this paper,pis the flow pressure.In Equation (1), the density of the fluid and the porosity of the porous media are functions of the flow pressure.The density of compressible fluids related to the flow pressure can be expressed as follows

Here ρcis the fluid density at reference pressurepc,c1is the fluid compressibility.Similarly,the pore volume of porous media related to the flow pressure can be expressed as follows:

Here φcis the porosity at reference pressurepc,c2is the compressibility of porous media.

We assume that there are several wells in the flow field that are used to produce fluids.Each well can be represented by a well model as follows:

whereqwEis the prescribed well production rate of wellwlocated in elementE,pwEis bottom hole flow pressure of wellw,pEis the flow pressure in elementE,WIis the well index which is calculated in advance, ρEis the fluid density computed bypE.The well production rate production rate can be regarded as the recovery rate at the well locationq(E,t) =qwE.

Equation (1) is a continuity equation, which describes the conservation properties of fluids.Equation (2) is the equation of motion,which describes the law of fluid motion.Equations(3)and(4)is a state equation, which describes the nonlinear relationship between the properties of fluid and pressure.Equation(5)is the well equation, which describes the linear relationship between production rate,well pressure and grid pressure.This model is coupled with Equation (1) by variablepEand indicates the linear relationship between well production rate and well bottom hole flow pressure.In many practical applications, well pressure can be measured and measured values include noise.The pressure data collected in the well have the relationship with the actual well pressure aswhereis the observed pressure data at timenin wellw,is the real pressure data of wellwat time stepn.

In order to obtain the solution of nonlinear Equations (1)-(5),the boundary condition and initial pressure conditionp(x,0)=p0need to be specified.The no-flowing boundary condition is specified as v·τ =0,where τ is the outer normal vector of the boundary.

2.2.Regularization form of the estimation problem

When the pressure data are collected at the wells,the system of Equations(1)-(5)becomes over determined,so the distribution of the pressure field needs to be estimated based on the data and the model.This leads to the following optimization problem:

where ωw,nis weighted value, function ϕ is defined as ϕ = (p,pwE),the system Equations (1)-(5) appears as the constraint conditione(ϕ) = 0.Since the system model is assumed to be accurate, we only need to estimate the initial pressure distributionp0.This least square estimation problem is an optimal control problem.The flow pressure function and the bottom hole flow pressure function are state variables and the initial pressure distribution function is the control variable.

However, the estimation problem is ill-posed when the dimension of initial pressure is greater than the number of collected data.Some prior constraints on initial pressure distribution need to be added to the optimization problem.Typically the objective function is modified by adding the penalized term such as Tykhonov regularization.In this paper, we use the first-order Tykhonov regularization term to penalize the objective function,which indicates that a flat initial pressure estimate can be obtained by optimization.

Here λ ＞0 is the regularization parameter which is specified in advance.Symbol ‖·‖2is the 2 norm operator.Appropriate regularization coefficients need to be determined by numerical experiments.

In this paper, we use the discrete adjoint approach to solve the problem(6).In this way,the continuous equation constraint in the problem is replaced by the discrete equation constraint which is the numerical scheme of Equations (1)-(5).In appendix A, the finite volume method is used to discretize this system of equations and the process of solving nonlinear system equations using Newton iteration is described.Based on numerical scheme (A7)-(A10) in appendix A and modified objective function (7), optimization problem can be rewritten in the following form.

whereen(ϕn-1,ϕn)=0 is the discrete form of continuous system Equations (1)-(5), ϕn=(pn,pw,n) is the collection of fluid flow pressure valuespnand bottom hole flow pressure valuespw,nfor each time stepn,Nis the total time step.We note that ϕ0=p0.In order to obtain the solution of the problem,it is necessary to study the optimality condition of the optimal control problem.

2.3.Optimality condition of the estimation problem

Since problem (8) is a constrained optimization problem, the Lagrangian function corresponding to the problem can be written as follows.

We define the adjoint vector ηncorresponding to the residual vectoren(ϕn-1,ϕn)of time stepn.According to optimization theory[24], optimization variables and adjoint variables satisfy the following optimality conditions at the optimal point.

Equation(11)is merely discrete system equations.Equation(10)is the adjoint equation and can be written in detail as follow form.

Here matrixCis a discrete gradient operator.The finite difference scheme can be used to determine the structure of the matrix.The row number of the matrix is the number of internal faces between the cells and the column number of the matrix is the number of cells in the grid.For a 4×4 grid,this grid has four internal faces and four cells.The matrix has the following structure.

The calculation of this partial derivative is adjoint gradient which is denoted bydepends on the solutions of the adjoint Equation (13).Based on the analysis of optimality conditions, the following iterative optimization Algorithm can be used to solve problem (8).

Algorithm.Adjoint gradient algorithm for state estimation

Newton update method can be used to calculate the matrixHk.The line search method is used to choose an appropriate step length α in Step 6.The termination conditions of the Algorithm are mainly determined by two conditions: (1) the difference of the objective function of the adjacent iterations; (2) the maximum number of iterations.This condition is set in advance.

In the process of implementation of the Algorithm, it is necessary not only to solve the nonlinear system Equations (1)-(5), but also to obtain the structure of the Jacobian matrix of the nonlinear equations.Calculating the Jacobian matrix of nonlinear equations manually is boring and error-prone.The more efficient way to generate Jacobian matrix for nonlinear equations is to use automatic differential technique [25-27].This technique requires only the discrete structure of the nonlinear equation to generate the Jacobian matrix at the specified value of the independent variable, which makes discrete adjoint approach more efficient than other methods in solving optimization problem.

3.Results

3.1.Case study 1

To verify the effectiveness of the proposed Algorithm, in this section, we design numerical experiments to study the performance of the algorithm.In the following three numerical experiments, we use the real initial pressure to generate bottom-hole pressure data.These data are disturbed by random noise that satisfies the normal distribution with mean of zero and variance of 1bar and collected at each time step.The data weight ωw,nof each well is set to 1 at each time step.The viscosity of fluid is 5 mPa s and the compressibility is 10-3/bar.The compressibility of porous media is 10-6/bar and the reference pressure is 200 bar.

In the first example, a homogeneous porous media is studied.The field has a 100×100×1 grid with permeability of 200 millidarcy and porosity of 0.3.Each grid covers an area of 1×1×1 m2.The two production wells are located in the northeast corner and the southwest corner, respectively.Fig.1 (a) shows the grid and well placement pattern.The total time to implement the simulation is 200 days and the total time step is 30.The real initial pressure of the model is constant value pr= 200 bar.The flow rate of two production wells is 15m3/day.

Fig.1(b)shows the values of the objective function values with two different λ during the optimization.The pressure value of each element at the initial point of iteration is set to 150 bar.Optimization terminates after the maximum number of iterations 30 is reached.Table 1 summarizes the optimization results with different λ.It can be seen from Fig.1(b)that the convergence rate of the Algorithm is fast and the objective function value can be reduced to a smaller value in only three iterations.

Fig.2 shows the estimated optimal initial pressure distribution after 30 iterations when λ=0 and λ = 1.Fig.3 shows the distribution of deviation values |pr-pe| corresponding to different regularization parameters, wherepris the real pressure andpeis the estimated pressure.

It can be seen from Fig.2 (b) and Fig.3 (b) that the estimated pressure values for most elements are close to the real pressure values except the region near wells.The deviation max|pr-pe| in Table 1 is 25.3870 bar when λ = 0, which indicates that the fluctuation of estimated pressure values is very large without regularization term.The deviation max|pr-pe|from Table 1 is 1.1720 bar when λ = 1, which shows that the optimization result with regularization is reasonable.The quantity‖∇pe‖2denotes the roughness of the estimated pressure.The roughness of the estimated pressure when λ=1 is 2.3032 × 105, which is smaller than the quantity 1.7410×106when λ =0,suggesting that the optimization leads to a more flat pressure distribution.

Table 1 Summary of optimization results in case 1.

Fig.4 shows the variation of the deviation max|pr-pe| for different λ.The deviation decreases first and then increases when λ =0,while the deviation value keeps decreasing when λ =1.This indicates that the regularization term affects the variation of the deviation.Proper regularization coefficient will reduce the deviation gradually in the optimization process and the absence of regularization term will make the deviation uncontrolled.

Fig.5 shows the comparison of observed well pressure, actual well pressure, and estimated well pressure with λ = 1.The actual well pressure is produced by the real initial pressure value and the estimated well pressure is produced by the estimated initial pressure value.The observed well pressure is the data collected in this case, which is produced by real well pressure and random disturbance.It can be seen from Fig.5 that although the observed well pressure values fluctuate greatly, the deviation between the estimated well pressure values and the real well pressure values is very small.This indicates that the estimated initial pressure can produce well pressure consistent with the real well pressure value.

Fig.1.(a) Grid and well placement in case 1; (b) Variation of objective function values with iterations in case 1.

Fig.2.Estimation results of case 1: (a) Estimated initial pressure distribution withλ = 0, bar; (b) Estimated initial pressure distribution withλ = 1, bar.

3.2.Case study 2

In this case study, we use heterogeneous porous media for numerical experiments.The reservoir is also divided into a 100×100×1 grid and each grid covering an area of 1×1×1 m2.The porosity of the reservoir is generated from random variables of normal distribution between 0.1 and 0.25.In this example, five wells were used for production and data collection and located at the four corners and centers of the reservoir respectively.Fig.6(a)shows the porosity and well placement pattern in case 2.

Fig.3.The deviation distribution|pr-pe|of real initial pressure and estimated initial pressure in case 3:(a)Deviation distribution with λ =0,bar;(b)Deviation distribution with λ = 1, bar.

Fig.4.Variation of deviationmax|pr - pe|with iterations in case 1.

Fig.5.Comparison of estimated bottom hole flow pressure, real bottom hole flow pressure and bottom hole flow pressure data in case1.

Fig.6.(a) Porosity and well placement in case 2; (b) The variation of objective function value with iterations in case 2.

The simulation experiment is carried out in the following schedule: the total simulation time is 365 days, which is divided into 50 time steps.The production rate of each production well is 4.5m3/day.The real initial pressure of the reservoir is set to 200 bar in each element.The pressure value of each element at the initial point of iteration is estimated to be 180 bar.Optimization terminates after the maximum number of iterations 30 is reached.

Fig.6 (b) shows the value of the objective function values with different regularization parameters for different iterations.Table 2 summarizes the optimization results corresponding to different λ.Fig.7 show the optimal initial pressure distribution corresponding to different λ after 30 iterations.Fig.8 show the distribution of deviation value |pr-pe| corresponding to different regularization parameters in case 2.

It can be seen from Fig.6 (b) that the descent rate of the objective function corresponding to different λ is different.The larger the regularization coefficient value the slower the decline of the objective function,which indicates that the regularization term may increase the cost of optimization calculation.When λ = 0,although the objective function decreases the fastest,it can be seen from Fig.7 (a) and Table 2 that the deviation max|pr-pe| is 17.9818 bar which is the largest of all values of λ.In contrast, the deviation max|pr-pe|is 0.9302 bar which is the smallest when λ =1.It can be seen from Fig.8 that the deviation value is large when λ = 0, while the deviation value of the other regularization parameters are small.This comparison shows that when the porous media is heterogeneous the ill-posed degree of the estimation problem becomes larger.The regularization term is necessary in order to obtain a reasonable estimation result.

Table 2 Summary of optimization results in case2.

The deviation between the estimated initial pressure and the actual initial pressure is 1.2145 bar when λ=10 which is close to the result obtained when λ = 1.The estimated initial pressure is generally larger than the actual initial pressure when λ = 100,which indicates that the objective function with larger regularization term may require more iterations to obtain reasonable estimation result.With the increase of regularization coefficient,the roughness ‖∇pe‖2of the estimated value in Table 2 decreases gradually.This shows that a large regularization coefficient will make the estimation smoother, but not necessarily reduce the deviation.Among these results, the result is better when λ = 10,which has smaller deviation and higher smoothness.

Fig.9 shows the variation of the deviation max|pr-pe| for different λ in this case.Similar to case 1, the deviation first decreases and then oscillates with iterations when λ = 0.In the iterative process, the deviation only slightly increases locally and the overall trend is downward when λ=1 and λ = 10.When λ =100, the deviation still declines first and then increases.These results show that appropriate regularization coefficients can gradually reduce the estimated deviation,which reducing the cost of optimization calculation.

Fig.7.Estimation results of case 2:(a)Estimated initial pressure distribution withλ =0,bar;(b)Estimated initial pressure distribution withλ =1,bar;(c)Estimated initial pressure distribution withλ = 10, bar; (d) Estimated optimal initial pressure distribution withλ = 100, bar.

Fig.8.The deviation distribution|pr-pe|of real initial pressure and estimated initial pressure in case 2:(a)Deviation distribution with λ =0,bar;(b)Deviation distribution with λ = 1, bar; (c) Deviation distribution with λ = 10, bar; (d) Deviation distribution with λ = 100, bar.

Fig.9.Variation of deviationmax|pr - pe|with iterations in case 2.

Fig.10 shows the comparison of real bottom hole flow pressure,estimated bottom hole flow pressure and collected well pressure data for different production wells with λ = 1.The observed well pressure data are still generated by real pressure data plus random noise with zero mean and variance of 1 bar.It can be seen from Fig.10 that, despite the large value of random disturbances, the deviations between the estimated bottom hole flow pressure values and the real bottom hole flow pressures value except for the first few time steps are very small,which is very similar to Fig.5 in case 1.The lower value of the objective function is also due to the smaller deviation value of the bottom hole flow pressure.This result also shows that the estimation obtained by optimization is reasonable.

3.3.Case study 3

In this section,we use the initial condition of spatial variation to verify the performance of the Algorithm.The reservoir model and well placement are the same as those in case 2 and are shown in Fig.6 (a).The total simulation time is 365 days and is divided into 50 time steps.The production rate of each production well is set to 4.5m3/day.The real initial pressure of the reservoir is shown in Fig.11 (a).

Fig.10.Comparison of estimated bottom hole flow pressure, real bottom hole flow pressure and bottom hole flow pressure data in case 2.

Fig.11(b)shows the value of the objective function values with three differentλ for different iterations.Table 3 summarizes the optimization results corresponding to different λ.Fig.12 show the optimal initial pressure distribution corresponding to different λ.Fig.13 show the deviation distribution |pr-pe| between the real initial pressure and the estimated initial pressure.

In Fig.11 (b), the descent rate of the objective function corresponding to different λ value is different.The descent rate of objective function decreases with the increase of regularization parameters.It can be seen from Table 3 that the deviation value max|pr-pe| corresponding to different parameter values is 4.55,4.98, and 5.31 respectively.There is no significant difference between these deviation values of different parameters.These deviations are larger than the corresponding values in Table 2 of case 2.This comparison shows that when the real pressure distribution is variety in space, the influence of different regularization parameters on the deviation value is small.This may be due to the real pressure of spatial variety makes the high-dimensional optimization more difficult.

Table 3 Summary of optimization results in case 3.

Fig.12.Estimation results of case 3: (a) Estimated initial pressure distribution withλ = 0.01, bar; (b) Estimated initial pressure distribution withλ = 1, bar; (c) Estimated initial pressure distribution withλ = 10, bar.

Table 3 also shows the roughness value ‖∇pe‖2of optimization results corresponding to different parameters.The roughness value of the estimation results decreases with the increase of the regularization parameters.Such a result can also be obtained through Fig.12.This indicates that the regularization parameters have a significant effect on the smoothness of the optimization results.

Fig.13 shows the distribution of deviation value |pr-pe| corresponding to different regularization parameters.Except for a small number of grid points with large deviation values,most of the other grid points have small deviation values.This indicates that the reasonable estimation results can be obtained by optimization.The regularization parameters also affect the smoothness of these deviations.From the comparison of these graphs, it can be seen that the deviation value corresponding to the smaller regularization parameter is rougher.

Fig.14 shows the variation of the deviation max|pr-pe| for different λ in this case.The deviation value always decreases when the regularization parameter is small when λ =0.01.The deviation increases first and then decreases when the regularization parameter λ=1 and λ = 10.These results indicates that appropriate regularization parameter can gradually reduce the estimated deviation, which makes the estimation result obtained by optimization more reasonable.

Fig.15 shows the comparison of real bottom hole flow pressure,estimated bottom hole flow pressure and collected well pressure data for different production wells with λ = 0.01.The observed well pressure data are still generated by real pressure data plus random noise with zero mean and variance of 1barsa.Despite the large value of random disturbances, the deviations between the estimated bottom hole flow pressure values and the real bottom hole pressure values are very small.The estimated well pressure values are slightly smaller than the real well pressure values.

4.Conclusions

Fig.13.The deviation distribution|pr-pe|of real initial pressure and estimated initial pressure case 3:(a)Deviation distribution with λ =0.01,bar;(b)Deviation distribution with λ = 1, bar; (c) Deviation distribution with λ = 10, bar.

Fig.14.Variation of deviation max|pr-pe|with iterations in case 3.

The state estimation problem of compressible flow in porous media is transformed into an optimal control problem.By analyzing the optimality conditions of the optimization problem,the form of discrete adjoint equation and adjoint gradient is obtained.The first-order Tykhonov regularization method is used to reduce the degree of ill-posed estimation problem.Quasi-Newton optimization method based on adjoint gradient is used to solve the problem.

Numerical experiments show that the search direction constructed by the quasi-Newton method can quickly reduce the value of the objective function.Different regularization coefficients will lead to different estimation results.The optimization result without regularization term deviates greatly from the real initial pressure, which is due to the ill-posed estimation problem.The deviation between the real pressure and the estimated pressure is small and the estimated pressure distribution is flat when regularization is adopted.However,large regularization coefficients may require more iteration to get reasonable estimate result.The method of determining the optimal regularization coefficient needs further study.Whether the porous media is homogeneous or heterogeneous, the bottom hole flow pressure produced by the estimated initial pressure is fit well with the true bottom hole flow pressure.

Fig.15.Comparison of estimated bottom hole flow pressure, real bottom hole flow pressure and bottom hole flow pressure data in case 3.

Acknowledgements

This study was financially supported by Science and Technology Project of Sichuan Province under the Grant No.16JC0314.

Nomenclature

φ Porosity,%

ρ Fluid density, kg/m3

v Flow rate, m3/s

qRecovery rate, m3/s

KPermeability, millidarcy

μ Fluid viscosity, mPa.s

pFlow pressure,bar

ρcFluid density at reference pressure, kg/m3

c1Fluid compressibility,%

c2Media compressibility,%

pcReference pressure

φcPorosity at reference pressure,%

qwEWell production rate,m3/s

pwEBot t o m-hole flow pressure,bar

pEPressure at well element, barWIWell index

Pressure data obtained in time step n at well

εnERandom error of pressure data

ωw,nWeight of pressu re data

ϕ Solutions of system equations

JObjective function without regularization term

JmObjective function with regularization term

p0Initial pressure,bar

λ Penalty parameter

LLagrangian function

η Adjoint variable

prReal pressure distribution, bar

peEstimated pressure distribution, bar

nTime step

kIterations

Appendix A

In this appendix, the finite volume method [28] is used to discretize the continuous Equations (1)-(5).Compared with other numerical methods, the finite volume method can preserve the local conservation of matter, thus obtaining a more real physical solution.

The computational area is divided into non-overlapping conformal polygons elements.For each elementEwe denote its area by|E|and its boundary by ∂E.The boundary ∂Eis the union ofmedges ∂E=∂E1∪∂E2∪…∂Emand the length of edge ∂Emis denoted by|∂Em|.The continuous equation is integrated on elementEand Gauss divergence theorem is applied.

Here τdis the outer normal vector of edged.In order to get the integral value of the edge,the velocity Equation(2)is substituted in the pressure Equation (1)and the following approximate estimate is obtained.

whereuE,dis the flux value in edgedof elementE,cE,dis the vector from the center of the elementEto the center of the edged,pEis pressure value of element center and πdis pressure value at edged.

Through the derivation of the above formula, the following formulas can be obtained.

Similarly,assuming that the element adjacent toEon the edgedisF, the following expression can be obtained.

By using Equations(A3)and(A4),the following flux estimates of edgedcan be obtained.

where ρdis the density value on the edged,which is computed by the density weighting of two adjacent elements, Γdis the conductivity coefficient on the edged,which is independent of the flow pressure.The approximate estimate of integral (A1) can be obtained by summing the flux of each edge on elementE.

Next,the backward Euler method is used to obtain the discrete time scheme in space and time.

Discrete Equation(A7)is implicit scheme of system Equation(1).The well index in (A10)is calculated by the Peaceman well model.The well index of a vertical well in a Cartesian grid with dimensions Δx×Δy×Δzcan be calculated as:

wherer0=is the effective well radius andrw=0.2mis the real well radius in this paper.The pressure variable on each element are collected as vectorpnand bottom hole flow pressure on each well is collected as vectorpw,n.We define the discrete state variables of time stepnas ϕn=(pn,pw,n).The discrete form of system Equations (A7)-(A10) is represented bye(ϕn-1,ϕn) = 0,n= 1,2…N.The equations are nonlinear and the following Newton iteration is used to solve the nonlinear equations at time stepn.

Wherekis iteration,ϕn-1is the solution of the system of equations at the previous time step,is the Jacobian matrix of the residual equationen(ϕnk,ϕn-1) with respect to ϕnk.The Jacobian matrix of the system equation can be generated efficiently by using the automatic differential technique.The solution of the system Equations (1)-(5) is obtained by the iterations of each time step.

Appendix B.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.petlm.2020.03.004.