A comparative study of some well-reservoir coupling models in the numerical simulation of oil reservoirs

In reservoir simulation, the most known well-reservoir coupling technique is based on the Peaceman’s equivalent radius, which is applied on the productivity index computation, that is used to relate flow rate, wellbore pressure, and the mesh block pressure. Original Peaceman’s model considers the assumption of steady-state flow, leading to an artifact on the calculation of the wellbore pressure, called numerical wellbore storage. This artifact is more significant for coarser meshes and initial time instants, and it is also a function of fluid and rock properties. In this context, we have implemented and compared different Peaceman’s technique extensions for productivity index calculation incorporating transient effects to prevent numerical wellbore storage. We have also considered some production scenarios for a vertical well and single-phase oil flow.


I. INTRODUCTION
During the 1960s, the term Numerical Simulation of Reservoirs became common in the oil industry [14], referring to the use of physical-mathematical models and computational tools to predict the performance of hydrocarbon reservoirs under different production scenarios. According to Dumkwu et al. [16], Numerical Reservoir Simulation has become an ally in Reservoir Engineering, and we use it in decision making that involves many financial resources, in estimating underground reserves, and to forecast and improve performance production of reservoirs in the oil and gas industry. This work focuses on the calculation of pressure estimates in wells, using well-reservoir coupling techniques.

Determination of pressure in oil wells
Petroleum is the name given to natural hydrocarbon mixtures that can be found in the solid, liquid, and gas phases, depending on the pressure and temperature conditions. In nature, oil can appear in only one phase, or it can appear as a multiphase mixture in equilibrium [30]. In many cases, the fluid phase in the reservoir depends on the thermodynamics properties of the fluid produced at the surface. Knowledge of the behavior, under pressure and temperature variations, of oil, natural gas, and water, alone or in combination, is of fundamental interest to Petroleum Engineering [14], whether under static or moving conditions, in the rocks of the reservoir or ducts.
Over the years, Petroleum Engineering has recognized the need for accurate information related to physical conditions in the well and inside the reservoir. With the initial progress in oil recovery methods, it became clear that calculations made with information from the surface or the top of the well could often lead to misunderstandings. Sclater et al. [33] described the first pressure value record at the bottom of the well and fluid collection in the well, under pressure, for sampling. These same authors define rock bottom data by referring to pressure, temperature, gas-oil ratio, and the chemical and physical nature of the fluid. Millikan et al. [21] highlighted the need of accurate wellbore pressure measurements when they described the first accurate pressure gauge, and they pointed to the fundamental importance of wellbore pressure knowledge for Petroleum Engineering, in search of more efficient recovery and lifting procedures. Based on this contribution, we were able to measure pressure, which is essential information for reservoir performance calculations [14].
The Well Testing Analysis is the engineering branch which obtains estimates of properties of the well-reservoir system from wellbore data under production/injection conditions, as well as through the application of inverse problems. We are interested in the pressure variation in the well as a function of time. We accomplish this by measuring the pressures at the bottom of the well and the flow at the surface, for example. From the measured pressure response, it is possible to determine reservoir properties useful for planning the completion of the well or the reservoir depletion plan [35]. During a well test, we can create a transient pressure response by a temporary change in the production rate. The well response is usually monitored for a relatively short period, compared to the reservoir life, depending on the test objectives. For estimates involving regions around the wells, we often carried out tests in less than two days. In the case of tests to analyze the reservoir boundaries, it may take several months to register pressure data [20].
It is usual to use analytical solutions to interpret the results of pressure tests in wells. Theis [37] proposed one of the first theoretical solutions for determining the well pressure, depending on the flow rate and the duration of production. Considering the onedimensional flow in the radial direction in a reservoir described in Cylindrical coordinates, producing at a constant flow q sc (defined in standard conditions, negative for production and positive for injection), from the initial time t=0, the analytical solution for the pressure in the reservoir in a radius r is given by [30]: where p i is the initial pressure of the porous medium, µ is the viscosity of the fluid, k is the permeability of the medium, h is the reservoir thickness, E i is the Integral Exponential Function, φ is the porosity of the medium and c t is the total compressibility that includes the contribution of the fluid (c o ) and rock (c φ ) compressibilities (c t = c o + c φ ).
For small values of the argument (X < 0.025), the integral exponential function can be approximated, with an error of less than 1%, by [26], where λ = exp(0.57722), and 0.57722 is the Euler constant. Thus, the pressure in the well (p wf ) is approximated by where r w is the radius of the well. Throughout the years, researchers have developed models for more complex reservoirs. However, there are many situations for which analytical solutions do not exist or are difficult to obtain. Heterogeneous reservoirs with irregular borders and with complex well geometry are examples where numerical simulation is the alternative to determine transient pressure responses.

Numerical reservoir simulation
As already said, we often use reservoir simulation models in the oil and gas industry. The reasons for such acceptance are due to some advances [17] The ultimate goal of a numerical simulation study is to accurately predict flow and pressure in the well and estimate pressure distributions in the reservoir (and saturation, in the multiphase case). When we intend to incorporate well modeling in reservoir simulation, some difficulties appear, which we can separate into three main categories [17]: 1. The cell of the computational mesh that contains the well is usually large compared to the wellbore radius, due to the discrepancy between the spatial scales involved. The average diameter of the production/injection well is about 10 cm, while the dimensions of a reservoir can be kilometers. Therefore, the cell pressure calculated is a poor estimate of the pressure in the well; 2. Coupling is a complex interaction between the reservoir (porous medium) and the well (a duct), particularly in the case of wells that pass through several layers of the porous medium; 3. There are difficulties in describing the wellreservoir coupling in the multiphase flow when the total production flow of the well is specified.
It is also worth noting that there is a variation in the time scale, with phenomena involving strong pressure gradients that may occur in short periods, such as fractions of an hour, and border effects (such as, for example, pressure drop at the reservoir boundaries concerning the initial pressure) that may take months to appear. These times are a function of the characteristics of the well-reservoir system. Other problems can also arise when considering the location of the producer/injector well displaced from the center of the cell, or even when there are several wells in a cell in the computational mesh.
The well-reservoir coupling model appears as a central component in reservoir simulation. Peaceman [26] introduced a well-known and commonly used model, which considers that the flow between the reservoir and the producing well is steady-state. Extensions were also created, for example, for the case of anisotropy [27], the incorporation of transient effects [9], and techniques for horizontal wells [2]. Although this coupling model provides suitable results for reservoirs in a wide range of applications, including their extensions, it introduces a phenomenon called numerical storage (due to the similarity of the results when considering wellbore storage, although it does not correspond to the physical behavior in the well) for the initial instants, and the application of methodologies to circumvent this numerical issue is the main objective of this work.
Therefore, this work aims to implement a model for a transient coupling of the well-reservoir system to minimize the phenomenon of numerical storage originating from the steady-state well-reservoir coupling, using the technique of Peaceman [26]. We consider here a single-phase flow, of a slightly compressible oil, in a reservoir containing a vertical producer well that penetrates the entire oil reservoir. We also use a polynomial function, with low error and wide application range, when calculating the integral exponential E i .

II. FLOW MODELING IN POROUS MEDIA
In this work, we consider the following hypotheses for the flow in the porous medium: 1. the porous medium is homogeneous and anisotropic in terms of its permeability; 2. porosity is a linear function of pressure; 3. the compressibility of the rock is small and constant; 4. the fluid is Newtonian and slightly compressible; 5. the fluid has a constant chemical composition; 6. there are no electrokinetic effects; 7. there are no inertial or turbulent effects; 8. the flow is single-phase and isothermal; 9. there are no chemical reactions.

Governing equations
The continuity equation is a mass balance equation so that the difference in mass entering and leaving a control volume must be equal to the accumulation of mass within the control volume in a given time interval. In oil reservoirs, the control volume is a portion of the porous medium that can contain one, two, or three phases of fluid [17]. It is worth mentioning that we consider the porous medium as a continuum. So, we define the physical properties, at any point, using the concept of Representative Elementary Volume (REV) [41]. For single-phase mass flows in a porous medium, we can write the continuity equation in the form: where ρ is the density of the fluid, v is the apparent fluid velocity (flow rate divided by the cross-sectional area), q m is a source term representing the production or injection of fluid, and V b is the volume of the rock (solid material plus pores). We can also write the continuity equation in terms of the formation volume factor [17], B = ρ sc /ρ, as where the sc subscript indicates the standard conditions. The civil engineer Henry Darcy, in 1856, through experiments of vertical water filtration in columns of homogeneous sand, obtained the first equation used to describe the movement of fluids in porous media. This equation was later called Darcy's law. It is an empirical relationship between the mass flow rate of fluid through a porous medium and the potential gradient, defined in the form [15] where K symbolizes the hydraulic conductivity of the porous medium, which combines fluid and medium properties, and ∇Φ is the potential gradient, defined as [17,30], where γ G = ρg, g is the magnitude of the acceleration due to gravity, and D is the depth (positive in the vertically downward direction). Later, Darcy's law could be deduced mathematically from the balance equations that govern the single-phase flow of a Newtonian fluid [40]. Rewriting Darcy's law in terms of the properties of the fluid, we can express the surface flow velocity as where k is the absolute permeability tensor. From the definition of the compressibility coefficients of the rock and the fluid, it is possible to rewrite the transient term of Eq. (5) in the form dφ dp ∂p ∂t = φ d dp where we considered that φ = φ(p) and B = B(p) and that [17] where the superscript 0 indicates the reference conditions, and c µ is the coefficient of variation of viscosity, and we assume that c o and c φ are small and constants.
We obtain the Hydraulic Diffusivity Equation (HDE) by combining Darcy's law with the mass conservation equation. Now, writing Eq. (8) in terms of its components, In turn, we can write the conservation equation as Next, we multiply Eq. (17) by V b = dxdydz, where A x , A y , and A z are, respectively, the areas of surfaces normal to the x-, y-and z-directions. Then, replacing v x , v y , and v z and using the definition q sc = q m /ρ sc we finally obtain Now, using the result obtained in Eq. (9), (20) where the Γ p coefficient represents the compressibility effects of rock and fluid, given by Incorporating the gravitational effects in the term Γ G , we can rewrite the HDE as follows: where

Initial and boundary conditions
To solve Eq. (22), we need to provide the appropriate initial and boundary conditions. As an initial condition, we take where p ini is the initial pressure distribution before the reservoir is disturbed by fluid production/injection. As external boundary conditions (external borders of the reservoir), we can use prescribed pressure or flow rate. For a sealed rectangular parallelepiped reservoir of edges L x , L y and L z , we have

Source term
In the numerical simulation of reservoirs, wells are considered internal boundaries, and, in the case of the use of Cartesian coordinates, the term source, used to represent wells, is directly related to the use of well-reservoir coupling techniques. We can impose these boundary conditions specifying the pressure in the well (Dirichlet condition) or the flow rate in the well (Neumann condition) [17].
For the representation of wells, it is necessary to obtain an expression that relates the pressure in the porous medium, p, with the pressure in the well, p wf , and the flow rate in the well, q sc . As an example of how we can use expressions of this type, we assume steady-state flow in the region close to the well. Further, we consider the radial flow of an incompressible fluid towards the vertical well of radius r w , in a rock formation with uniform permeability and thickness. So, assuming these conditions, it is possible to write [17]: where k H represents the permeability value in the radial direction. Integrating Eq. (26) between the radius of the well, r w , and an arbitrary radius, r, where r w ≤ r ≤ r e (r e is the external radius), it is possible to obtain an equation for the pressure distribution for the steady-state well-coupling, For r = r e , where we define the pressure as p e , we can rewrite Eq. (27) as so that Eq. (28) provides the well production in terms of external and well pressures.
In general, we express the production rate, under standard conditions, using the pressure of the well and the average pressure of the reservoir [17]. From the development of Van et al. [39], for a reservoir described in Cylindrical coordinates, considering the average volumetric pressure in the oil reservoir, p, between r w and r e , and even though r e r w , we have which differs from Eq. (28) only by the skin factor s, the 1/2 in the denominator, and by replacing p e with p.
The introduction of the skin factor makes Eq. (29) more general, valid for the case where there is an additional pressure drop around the well, for example, due to the formation damage. For the specific case in which the permeability may change, we determine the s factor by [18,30] where k s is the modified permeability, and r s is the radius of the region where the permeability can vary. When k s < k, the value of s is positive, and we have damage to the formation. Otherwise, s is negative, and we have stimulation. On the other hand, when s = 0, there is no influence on the productivity of the well by local permeability changes near the well [30].
Commonly, we rewrite Eq. (29) in a more compact form, introduced by Van et al. [39], where J w is the Productivity Index (PI): or yet, where G w is the geometric factor of the well, defined as The factor F , in Eqs. (32) and (34), assumes the values of 1/2 for steady-state conditions and 3/4 for the pseudo-steady state [17]. The representation of Eq. (32) for the Productivity Index comes from analytical solutions. However, we also apply the concept of the Productivity Index in the context of numerical reservoir simulation, and its calculation is associated, in this case, with the use of well-reservoir coupling techniques.

III. NUMERICAL METHODOLOGY
We describe fluid flow in porous media using a set of partial differential equations, and we can not always obtain analytical solutions due to the non-linear nature of the equations. Therefore, we must use numerical techniques to solve the balance equations that govern the flow. Among the numerical methods that we can apply in its resolution, we choose to use the Finite Difference Method (FDM), which is the most used in the oil industry [17,26].

Discretization
Discretization is the process of converting the PDE, valid in the continuous medium that defines the resolution domain, by a set of algebraic equations defined in a discrete domain, to obtain a numerical solution. The first step in the discretization stage is the choice and construction of the numerical mesh, that is, the partitioning of the resolution domain. In the Finite Difference Method, this implies a computational mesh with a finite number of nodes where we determine the values of the dependent variables. Then, we must approximate the existing derivatives in the governing partial differential equations.
Generally, we use two mesh systems in numerical simulation of reservoirs: the centered block system and the distributed point system [17]. In this work, we use the former.
The mesh system contains, in the x-direction, n x cells that we superimpose on the reservoir. We center each i cell on x i , and we designate its borders by x i−1/2 and x i+1/2 . The cells have dimensions equal to ∆x i , constants or not, satisfying the following relation so that the cells must cover the entire L x length of the reservoir [17]. Similarly, it is also possible to discretize the reservoir in the y-and z-directions. The cells should be small enough to describe the heterogeneous nature of the reservoir and, thus, allowing to represent the flow characteristics adequately. However, it is necessary to carefully determine the total number of cells in the mesh, because the higher the number of them, the higher the number of unknowns in the algebraic system that must be solved, implying an increase in computational effort.
In Fig. 1, a three-dimensional reservoir is discretized into n x blocks (or cells) in the x-direction, n y blocks in the y-direction, and n z blocks in the zdirection, considering that n x =n y =n z =3, for this particular case. In more general situations, the dimensions of the blocks need not be the same, and: where the j and k indexes indicate the blocks in the y-and z-directions, respectively. In this system, we number the border blocks by adding the fraction ±1/2 to one of the i, j, and k indexes, depending on the border that we want to represent. For example, we indicate the boundary between cells i, j, k and i + 1, j, k by i + 1/2, j, k.

Numerical approximation of derivatives
Considering the partitioning illustrated in Fig. 1 and a totally implicit numerical scheme, the first term containing the spatial derivatives of Eq. (22) can be discretized by a central difference approximation of second-order, where with similar equations for y-and z-directions. Still, in Eq. (38), we do not know the pressures at time n + 1, ∆x i,j,k is the spacing of the cell i in the xdirection, and the subscript i ± 1/2 indicates that the variables must be evaluated at the cell boundary in the x-direction. We can also obtain similar forms for the yand z-directions, considering that ∆y i,j,k and ∆z i,j,k are the spacing in these respective directions and that we represent these boundaries by the indexes j ± 1/2 and k ± 1/2. Before proceeding, we introduce transmissibilities as being represented by T x , T y and T z : T n+1 and we calculate the properties at the cell interface using harmonic averages for rock and geometric properties, and arithmetic means for fluid properties. Again, employing central difference approximations, and we can proceed similarly to the y-and z-directions. Now, we apply conservative expansions to the terms of accumulation to preserve the mass balance. It is worth mentioning that the use of non-conservative schemes in the finite difference equations does not necessarily produce results that are not correct [17]. Therefore, the final discretized form is given by: where V b,i,j,k = ∆x i ∆y j ∆z k . We approximate the time derivative by a backward Euler scheme Finally, considering the use of an iterative resolution aiming to linearize the system of equations we obtain where the level v indicates the known values while the level v + 1 the values to be determined. Next, we move on to the discretization stage of the well-reservoir coupling, Through the term (J w ) n+1,v i,j,k the well-reservoir coupling occurs, following the methodologies based on the work of Peaceman [26], which is still the most used to date in numerical reservoir simulation studies.
We note that we use the term q sc to represent a source term for the computational mesh cell. In general, the well passes through a set of cells in the computational mesh. Therefore, the expression to the total flow of the producing well is given by: if we assume a vertical well that contains the cells W i up to W f (Fig. 2) and we do not consider the frictional and convective effects inside the well. We adopted as a reference level k = W for calculating the pressure in the well. Considering only the gravitational effects inside the well, we can write: where and we find the pressure in the well, in each layer, in the computational cell in terms of wellbore pressure in the reference level layer: and ρ n+1,v i,j,k+1/2 is the arithmetic mean between ρ n+1,v i,j,k and ρ n+1,v i,j,W .

Numerical solution of the algebraic equation system
The set composed of Eqs. (47) and (50) forms a linear system of equations whose dependent variables are the pressures in the reservoir and the production well.
We choose the Conjugate Gradient Method [19] to solve the resulting algebraic equation system. Initially, it was a direct method [19]. However, it became known for its properties as an iterative method, especially after the development of sophisticated preconditioning techniques [32].
We can increase the rate of convergence in the iterative process by using preconditioning techniques to reduce the conditioning number of the matrix by changing the original system of equations and, thus, their eigenvalues. Thus, the method used here is the Preconditioned Conjugate Gradient (PCG), using a diagonal preconditioning technique [7,17,32]. If we consider an iterative procedure, we attain the numerical convergence when where χ represents the pressure in the reservoir, and well and tol is a numerical tolerance. We use a twostep iteration procedure [17]: 1. in a set of internal iterations, we apply the PCG method to obtain the pressures, and 2. in a set of external iterations, we update the coefficients, such as the transmissibilities.

IV. WELL-RESERVOIR COUPLING
The modeling of wells in the numerical simulation of reservoirs presents some difficulties as, for example, the discrepancies in spatial scales when comparing the wellbore radius with the dimensions of the oil reservoir, the geometry of complex wells, the flow transition from the porous medium to a duct, among other factors. In this section, we investigate the treatment of the term source, review some wellreservoir coupling models, focusing on the Peaceman [26] model and its extensions, and discuss the methodology that incorporates the transient effects, both in the productivity index and in the equivalent radius.

Some well-reservoir coupling models
One of the difficulties in modeling wells in fieldscale simulations is the fact that the region where the highest pressure gradients typically occur is close to the well, which is much smaller than the usual dimensions of a cell in a computational mesh [13]. Except for simulations that somehow make use of Cylindrical coordinates, source/sink terms are, in general, used for the implementation of internal boundary conditions involving wells [3,17]. Well-reservoir coupling techniques aim to relate the pressure in the reservoir cells to the pressure in the well, by using, for example, equations that contain source/sink terms and the productivity index [35].
Among the first works to related the pressures in the cell and the well, we found that of Van et al. [39]. They stated that the pressure of the cell containing the well (for a wellbore centered in it) must be equal to the average volumetric pressure of the cell. Thus, for the hypotheses of steady-state flow, the average pressure p (or the pressure in the cell), located in the reservoir between the radii r w and r b (r b r w ) is given by: where r b is called the block radius and determined, for a vertical well, as We observe, then, that a methodology to deal with the problem of the representation of wells in Cartesian coordinates consists of initially considering the onedimensional and radial mass flow around the well [3]. The idea is to reconcile the analytical solutions, which we can obtain for the one-dimensional flow in Cylindrical coordinates, with the discretized equations, written in Cartesian coordinates.

Peaceman's equivalent radius
We can obtain the first models for well-reservoir coupling by considering a two-dimensional flow governed by the Darcy equation written in Cartesian coordinates and a one-dimensional flow in a cylindrical reservoir. Peaceman [26], to obtain its well-reservoir coupling model, used the discrete equation for the pressure p(x, y) in the two-dimensional flow together with the analytical solution for the pressure p(r) in the one-dimensional flow, both for a porous medium of the same length ∆z.
From the hypotheses of steady-state flow, homogeneous and isotropic porous medium (k x = k y = k H ), assuming that ∆x = ∆y, and symmetry, we obtain from the discrete form of HDE [17] On the other hand, we know that the analytical solution of the one-dimensional flow in the radial direction, in a reservoir described in Cylindrical coordinates, for a steady-state flow, is given by [12] Peaceman [26] considered the pressure p(r) coinciding with the pressure p i+1,j . So, for r = ∆x, we obtain Combining Eqs. (56) and (58), that is, considering the discrete form in Cartesian coordinates and the analytical expression for the Cylindrical geometry around the well, we can arrive in an equation for the flow in the well, written using the productivity index, where r eq is Peaceman's equivalent radius, r eq ≈ 0.198∆x [17].
Since then, researchers have created several extensions assuming a steady-state flow, from the Peaceman model for the equivalent radius [26], to improve the numerical approximation considering other hypotheses. For example, the concept of the equivalent radius was investigated for the case of k x = k y and ∆x = ∆y, thus allowing us to address more comprehensive situations in numerical simulations [27]. Using coordinate transformations, Peaceman [27] obtained the following expression for the equivalent radius which we must use in conjunction with the equation It is worth mentioning that this derivation for the equivalent radius, for an anisotropic medium, considers that ∆x = ∆y but that the mesh is uniform and that the well is far from the reservoir borders.
Other extensions have been developed, such as the one used by Al-Mohannadi et al. [2] J w = 2π k y k z ∆x Bµ Another well-coupling model, also for a horizontal well, is given by [4,5] where r eq = 0.14 (k y k x ) for a well parallel to the x-axis, where α H = (∆y/∆x) k z /k y . As examples of other methodologies we can quote Peaceman [28], Peaceman [29], and Babu et al. [6]. Researchers have even studied inclined wells [22].
So that we can apply these well-reservoir coupling models, some restrictions must be considered, such as, for example, the existence of a minimum distance between the wells and the reservoir borders, and between the wells themselves. Despite its limitations [23], we extensively use the technique introduced by Peaceman [26] in the numerical simulation of reservoirs.

Extensions for the transient productivity index
As already said, the results obtained with the technique suggested by Peaceman [26] are subject to a numerical artifact, which is not related to the phenomenon of physical wellbore storage. Fundamentally, this happens due to the use of the hypothesis of steady-state flow near the well. In reality, we should have considered the fluid flow as being transient. Therefore, other models have been developed, including the well-reservoir coupling for the transient regime. For flow in the transient and pseudo-steady state regime, in the vicinity of the well, Peaceman [26] suggested that we can obtain the equivalent radius from the following equation: is the dimensionless time and ∆L (∆L = ∆x = ∆y) is the spatial increment of the computational mesh, and is the dimensionless pressure, and p i,j,k is the pressure in the cell that contains the well. Blanc et al. [9] also concluded that when t D > 1, the equivalent radius calculated using Eq. (66) is close enough to its value determined by Peaceman's model for steadystate flow. Over time, we have been using analytical solutions capturing the transient regime, in the porous medium, to interpret the results of pressure tests in wells. Theis [37], for example, offered one of the first solutions to the problem of pressure drop in a well.
Blanc et al. [9] presented a correction to the original Peaceman model and introduced the Transient Well Index: where E 1 (u) = −E i (−u) and r eq = 0.198 ∆L (steadystate flow hypothesis). This model offers better solutions than those obtained with the conventional Peaceman technique for the initial times, but it can present deviations when the transient regime appears in the porous medium [9]. Therefore, instead of using the equivalent radius as r eq = 0.198 ∆x, it is possible to calculate its transient version. The pressure variation, for radial flow in the transient regime and a vertical well, is given by [2]: From Eq. (70), we can use the Newton-Raphson method to determine the r eq as a function of time, where the subscripts n and n + 1 represent the iterations in the Newton-Raphson method and while its derivative is given by In Eq. (72), ∆p = p ini − p i,j,k , where p ini is the initial pressure of the reservoir, and p i,j,k is the pressure in the cell containing the well. From these equations, we can obtain the expression for the determination of the equivalent transient radius and In this work, p i,j,k is evaluated at time t n+1 . When the boundary effects start to affect the wellbore pressure behavior, r eq (t) becomes too large to represent the flow dynamics. In this case, we calculate r eq (t) from Eqs. (70), (71) and (72) while t ≤ t * , where t * is the characteristic time associated with the boundary effects [2]. For t > t * , we must use r eq = r eq (t * ), and we keep the transient equivalent radius, calculated at t * , frozen until the end of the simulations.
Here, a criterion to shift the equivalent transient radius calculation is considered, based on the time required to reach the pseudo-steady regime. Therefore, we use the concept of dimensionless time considering A (the area normal to the z-axis) as a reference [30] so that and we assume that the producing well is in the center of the reservoir and that it has a square drainage area. For t DA < 0.09, the reservoir still behaves as being infinite (transient regime). On the other hand, for t DA = 0.1, we reach the pseudo-steady flow regime [30].

V. NUMERICAL RESULTS
In the study carried out in this work, the wellreservoir coupling models, including transient effects, were incorporated into our simulator, developed in C programming language.
We calculate the well pressure values considering a constant production flow rate under standard conditions. We use two types of graphs in the analysis of the results: 1. specialized plot: well pressure curve as a function of time; 2. diagnostic plot: pressure drop curves in the well, ∆p wf , and the Bourdet derivative [11] ∆p wf = d∆p wf d ln ∆t = ∆t d dt (∆p wf ) .
as a function of time.
Our approach is in the context of reservoir simulation and the transient analysis of pressure tests. We know that from the pressure values, it is possible to determine some reservoir properties, useful for planning the completion of the well or the reservoir depletion plan. In addition to the well-reservoir coupling model of Peaceman [26], considering a steady-state flow, three extensions were implemented, assuming a transient flow regime. In all simulations, we stipulate that ∆x = ∆y and k x = k y = k z = k.

Model 1
It is the conventional model of Peaceman [26], where we assume a steady-state flow and whose productivity index and the equivalent radius considered are, respectively, given by Bµ ln req rw and r eq = 0.198∆x.

Model 2
Peaceman [26] also proposed this model. However, it considers that the flow is transient in the cells containing the well. The productivity index is the same as for Model 1, although the equivalent radius incorporates the transient effects: Bµ ln req rw where and we use here the criterion provided, based on t DA , to shift the calculation of the equivalent radius (we do not update r eq when border effects occur).

Model 3
Blanc et al. [9] originally proposed this model, and it employs a transient productivity index. Nevertheless, the equivalent radius proposed by Peaceman [26], for the steady-state flow (Model 1), is applied where r eq = 0.198∆x.

Model 4
In this well-reservoir coupling model, we consider the incorporation of transient effects in the productivity index and the equivalent radius [2,9]. The productivity index is the same as that defined in Model 3, and we calculate the equivalent transient radius using the Newton-Raphson method. Thus, we have: knowing that and The criterion based on Eq. (77) is also applied to interrupt the update of the transient equivalent radius. Many analytical approaches to evaluate E 1 (u) have appeared in the literature. Tseng and Lee [38] provide an excellent survey on the subject. Note that there is no single analytical approximation that is valid for the entire range for u > 0. Tseng and Lee [38] claim that different methods are suitable for approximating E 1 (u) for a wide range of u values. Nevertheless, in practical applications, we often do not have only a single calculation method that covers the entire range of interest [8].
In this work, as we are interested in the shortest times and flows that may have low permeabilities and high viscosities, using power series to approximate the exponential integrals E i is not suitable for us. Therefore, we chose the approximation of the exponential integral E 1 (x) suggested by Segletes [34], based on a polynomial fit [31]. It is not our knowledge that others have done this previously.

Analytical solution and mesh refinement
Now, we compare the results obtained with the different models of well-reservoir coupling. Table 1 contains the parameters used in the construction of the standard simulation case. Unless when explicitly mentioned, the data in this table are those used in all simulations performed.

Parameter
Value 10,000 ft L y 10,000 ft L z 80 ft L wf 80 ft n x 81n y 81n z 5p ini 8,000 psi p 0 8,000 psi Q sc -400 STB/day In this table, L wf is the length of the producing well, t max represents the maximum production time, δ ∆t is the factor used to vary the time increment (∆t n+1 = δ ∆t ∆t n ), and ∆t ini and ∆t max are the initial and maximum values of the time increment, respectively.
Here, the results show the pressure curves in the well, the pressure drop, and the Bourdet derivative [11,10] as a function of the elapsed production time.
The first test performed was done to compare results from well-reservoir coupling models with an analytical solution for the particular case of constant viscosity and formation volume factor. Figures 3 and 4 show comparisons of the results obtained with Models 1 and 4, using the standard set of parameters for the simulations, with the results obtained with the analytical solution given by Ozkan [25].
It should be made clear that the analytical solution, as already said, makes use of simplifying assumptions, such as B and µ constants. Besides, we must apply a numerical procedure to obtain the analytical solution, which depends on the evaluation of functions of Bessel [24] and numerical inversion of Laplace transform [25], performed using the Stehfest algorithm [36]. Therefore, the results for the analytical solution are not available in the entire range of numerical results. Even so, it is possible to observe that the results are in good agreement with the numerical results obtained using Model 4. On the other hand, we note a deviation when there is a comparison with Model 1 (Fig. 3). When the border effects are present, the three solutions have very similar behaviors (about 250 days of production). Figure 4 corroborates the discussions about Fig. 3, and we see that the Bourdet derivative highlights the difference in the behavior of the numerical solution and the effect of the numerical artifact due to the use of a constant equivalent radius (r eq ).
Also noteworthy is the qualitative similarity of the result obtained with the original technique of Peaceman [26] (Model 1) with the behavior of the pressure when there is physical storage in the well. Also, there is a discontinuity in the derivative curves corresponding to the transient and the pseudo-permanent regimes. It is due to the strategy employed to change the calculation of the equivalent radius. From these results, it is possible to conclude that the results of Model 4 present a behavior compatible with that expected for the real problem.
Aiming to compare the four well-reservoir coupling models that we implemented in this work, in Fig. 5, we show the set of results for the pressure in the well. We can note the existence of a plateau in the curve obtained with Model 1 for the initial production times. For the standard case, this plateau (numerical storage) lasts approximately one day of production, interfering in the pressure response, reaching a difference of about 33 psi for the same time when compared to the Model 4 (in principle, the most correct). Regarding the models that incorporate the transient effects in the equivalent radius, Models 2 [26] and 4, we observe that we obtained better results with these models (superimposed in Fig. 5) in comparison with those based on the steady-state flow assumption. Besides that, in the initial instants of production, the numerical artifact no longer appears, and we correctly capture the expected flow regimes. However, according to Blanc et al. [9], the results of Model 2 are not always as favorable as those obtained here. For example, it can happen if we vary the reservoir and fluid properties.
Given the results obtained with Model 3 (Fig. 5), it is possible to observe the numerical artifact occurring. However, with less intensity than that of Model 1, and we do not perceive it when the production begins. It appears approximately 15 minutes after the beginning of the production. Nevertheless, after about one day, the results show behavior consistent with those of Models 2 and 4. We should remark that the results of Model 3 are in line with that reported by Blanc et al. [9]. On the other hand, in Model 4, when we incorporate transient effects in both productivity index and equivalent radius, the numerical artifact does not exist anymore, and the results are physically correct for the entire production duration. For the four models, we capture the border effect when we reach a time of approximately 250 days of production. Among the four models, we can say that Models 2 and 4 were the ones that presented the best results considering the transient flow regime. Finally, we should stress that the behavior of Model 4 results is also in line with those reported by Blanc et al. [9]. Furthermore, we also carried out a mesh refinement study using Models 1 and 4 that take into account the transient effects of flow. Table 2 shows the number of cells that we used in the generation of the different meshes that we employed in the mesh refinement study. Here, n x , n y , and n z are the number of cells used to discretize the oil reservoir in the x-, yand z-directions, respectively.   Mesh n x n y n z  1  11  11  5  2  21  21  5  3  41  41  5  4  81  81  5  5  161 161 5 For Model 1, it was possible to observe a numerical convergence as we refine the mesh. Except for the initial moments, where the solutions are quite different. We can also realize that the mesh refinement influences the magnitude of the numerical artifact, since the more refined the mesh, the smaller the storage phenomenon. For example, for Mesh 1 (n x =n y =11), we noted that the numerical artifact lasted approximately 30 days, while for Mesh 5 (n x =n y =161) the duration was about 0.1 days.
On the other hand, we did not observe the same behavior when we used Model 4 (Fig. 7), which incor-porates the transient effects. In this case, the pressure variation presents a behavior compatible with that of real fluid flow in an oil reservoir. Although the results are practically overlapping, for all meshes and time, including the border effects, numerical convergence did not occur in a similar way to that of Model 1. However, the difference between the values is small, and we cannot see it in the graph. It is worth pointing out that the same issue was reported in the literature [1,2], indicating a possible loss of accuracy when we utilize very refined meshes. However, this is still an open problem. Another important conclusion that we can draw is the fact that we can use less refined meshes, because the difference between the pressure values, obtained with more refined meshes, is almost imperceptible. Therefore, this fact implies less computational effort.
We also analyzed the influence, on well pressure variation, of the growth rate of the time step in the well-reservoir coupling for Model 4. We obtained the results considering three different values of the growth rate of the time step, δ ∆t =1.15, 1.20, and 1.25. We show the results in Figs. 8 and 9, for the pressure in the well, the pressure drop, and the derivative of the

www.ijaers.com
pressure drop. For all the values, we noticed that the results are practically the same. Furthermore, they coincide throughout all simulation time and present the behavior expected for this type of flow. That is, for the single-phase flow of oil in a reservoir containing a vertical producer well.

Low permeability and high viscosity
We also performed other simulations varying the permeability and viscosity. However, we employed lower values for permeability and higher for viscosity. In this study, we only considered Models 1 and 4 and specialized plots. We did a more detailed analysis to compare the two well-reservoir coupling models, for these permeability and viscosity values, to understand how the choice of the model may impact the results when the numerical artifact is more significant. We considered a total production time of 80 days and a prescribed flow (q sc ) of -100 STB/day. This approach aims to deal with situations closer to those of unconventional reservoirs. This type of study has become a trend more recently, and, in general, the effect of numerical storage tends to be more pronounced.
For the model of Peaceman [26] (Model 1), it is possible to observe that the permeability value directly influences the magnitude of the numerical artifact in the initial moments. The higher the permeability, the smaller the size of the artifact (Fig. 10). It is physically consistent because when the permeability is higher, the lower the flow resistance in the porous medium. Therefore, it reaches the transient regime in the porous medium more quickly, as well as for border effects. It is worth noting that since permeability has a value considered low for this model of a well-reservoir system, the numerical artifact has a longer duration. Thus, if we intend to study the phenomena that occur in the initial moments (or even for longer times) with this model, the results will not be accurate.
As an illustration, in Fig. 10, for the lowest permeability (k=0.25×10 −3 Darcy), the numerical artifact impacted almost all results. The pressure difference, comparing the two models, was 376 psi at t=2 days of production and 54.6 psi at t=70 days. For k=1.00×10 −3 Darcy, the magnitude of the numerical artifact was smaller, but not negligible since, after 0.6 and 19 days of production, there was a pressure difference of about 94 and 10 psi, respectively.
It was also possible to observe that for these low values of permeability, besides the duration of the numerical artifact being longer, in some cases, it was not possible to capture the border effects. It is relevant in the context of the analysis of pressure tests, as it indicates the real need for long-duration tests when it is necessary to determine boundary effects. We know that it is due to the great difficulty of the fluid to flow through the porous medium and, therefore, it will take longer to perceive the reservoir boundary effects (Fig. 10). Therefore, we can conclude that the transient well-reservoir coupling model presented the best results since we were able to prevent the appearance of the numerical artifact. Besides that, the use of the expression of Segletes [34] to calculate E 1 allows us to obtain results in a range in which other formulas failed. Contrary to what happens for the permeability variation, for the conventional model [26] (Model 1), the numerical artifact has a longer duration and magnitude as the viscosity increases (Fig. 11). Therefore, the time of occurrence and the magnitude of the artifact vary depending on the viscosity. Outside the region associated with the numerical artifact, the behavior of the results of the conventional model qualitatively tends to be physically correct. Also, we know that border effects occur more quickly for lower viscosity values. Regardless, we did not detect the boundary effects for any of the viscosity values proposed in our test.

VI. CONCLUSION
The objective of this work was to implement a model for the transient well-reservoir coupling in a reservoir simulator to correctly describe the pressure behavior in the well in the initial production times. We also investigated the effects of the mesh refinement and the size of the time step on the results.
As well known, the traditional well-reservoir coupling model of Peaceman [26] is still widely used nowadays. However, we saw that it has a flaw that leads to numerical storage at the beginning of the simulation. Nevertheless, it is simple, and we can correctly determine the wellbore pressure in a wide range of applications except for the initial instants of production. We must also emphasize that we were able to detect the border effects with this technique, without the need for a specific criterion to change the calculation of the equivalent radius.
In addition to physical properties, the numerical artifact is also affected by mesh refinement. The refinement can be used, in some cases, to mitigate the nonphysical storage, but it leads to an increase in computational cost. On the other hand, despite the qualitatively and quantitatively correct results that we obtained with Model 4, its numerical convergence must be better understood, as already pointed out by other authors. However, this is not a big problem, because we have not to use refined computational meshes when considering Model 4, as we could see in this work. Nevertheless, this issue deserves further studies.
We must also stress that the use of the expression proposed by Segletes [34] to calculate the exponential integral, allowed us to consider flows with low permeabilities, high viscosities, and short production time. Sometimes, this is not possible when we apply other expressions available in the literature.