Regular Article
A numerical well model considering capillary end effect for the lowpermeability gas reservoir
Department of Thermal Science and Energy Engineering, University of Science and Technology of China, Hefei, Anhui 230026, PR China
^{*} Corresponding author: xhwang@ustc.edu.cn
Received:
18
July
2019
Accepted:
4
November
2019
In the recovery of low permeability reservoir, the capillary pressure has an important effect, which may reduce gas production. Due to the capillary end effect, there are two flow patterns in the vicinity of the production well: both water and gas phases can flow into the production well and only gas phase can enter the production well. Based on the analytical equations of onedimensional radial flow considering the capillary end effect, an alternative numerical well model for low permeability gas reservoir is constructed. Numerical examples show that the proposed model can reflect the dramatic change for saturation and gasphase pressure in the vicinity of the production well both for two flow patterns, and therefore can predict the gas and water productions accurately at different grid scales. In contrast, the original Peaceman well model for multiphase flow, which is directly extended from that for singlephase flow, only can provide good prediction under the condition of that the grid size is enough small. Especially, for the original Peaceman well model, this problem is out of control since it is difficult to estimate a suitable grid size.
© Y.M. Yao et al., published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
B: Capillary pressure constant, MPa
h: Thickness of the reservoir, m
K: Absolute permeability, mD, 10^{−15}m^{2}
K_{rg,eq1} : Equivalent relative permeability
L_{eq} : Equivalent interblock distance, m
m: Pseudopressure, kg/(m^{3} · s)
P_{c} : Capillary pressure, Pa
P_{init} : Initial pressure of gas reservoir, Pa
P^{(e)} : Pressure at equivalent radius
S: Water phase normalized saturation
S_{o} : Outlet saturation in the condition of only water flowing into wellbore
S_{rw} : Water saturation at the entrance of wellbore
T_{eq} : Equivalent transmissibility, m^{4} · s kg^{−1}
λ_{w,eq1} : Equivalent mobility, m · s kg^{−1}
Subscripts and superscripts
(i, j): Physical quantity in wellcontained grid
(adj): Physical quantity in adjacent grid
inlet: Physical quantity at the inlet
outlet: Physical quantity at the outlet
1 Introduction
Numerical reservoir simulation for oil or gas production must account for the presence of wells. It is necessary to use grid blocks whose horizontal dimensions are much larger than the diameter of a well (Williamson and Chappelear, 1981). As a result, the pressure calculated for a well block is greatly different from the bottomhole pressure of the well. In the reservoir simulation, it is necessary to build a numerical well model to calculate the well flow rate associating with the well block pressure and the bottomhole pressure. Peaceman well model, constructed from the analytical solution of OneDimensional (1D) radially symmetric singlephase flow, is widely used in the general reservoir simulator. Initially, Peaceman’s finite difference well model was for the well located at the center of square grid for single phase flow (Peaceman, 1978). Then, it was extended into various applications, including horizontal wells, anisotropic reservoirs, and also multiphase flow, etc. (Babu and Odeh, 1988; Babu et al., 1991; Peaceman, 1983; Shiralkar, 1989; Su, 1995). Chen and Zhang (2009) extended the well model for other numerical methods such as standard finite element, control volume finite element, and mixed finite element methods.
By analogy to singlephase flow, numerical well models are directly extended to multiphase flow without the solid support of mathematical analysis. The singlephase flow in the neighborhood of the well is described by the linear equations, in contrast to that the multiphase flow should be described by the nonlinear equations. Only when the phase saturations remain constant, the direct extension to multiphase flow is accurate. It has been shown that, at the interface of two different type rocks which have the different relative permeability and capillary pressure curves, multiphase flows in porous media face the problem of the discontinuities of fluid properties and transport properties (Ding et al., 2016; Wang et al., 2006). This implies that the physical quantities might vary quickly in the neighborhood of the production wells.
It has been observed that in the core experiment, when twophase fluid flows into the noncapillary force container from the core sample containing capillary force, the wetting phase saturations at the outlet end of the core increase significantly compared with the rest of the area (Richardson et al., 1952). This phenomenon was called as the capillary end effect. Ignoring the capillary end effect, the measurements of relative permeability will produce big errors. Assuming that waterphase saturation is 1 at the outlet end of the waterwet core, a saturation distribution equation near the outlet end was established (Hadley and Handy, 1956). Then the way to correct the errors in the relative permeability measurements, considering the capillary end effect, has been proposed (Huang and Honarpour, 1998; Virnovsky et al., 1995).
In the lowpermeability reservoir, the capillary pressure has an important role on multiphase flows. The capillary pressure can be ignored inside the wellbore while it should be considered for the multiphase flow inside the reservoir. Analogous to the core experiment described above, the capillary end effect may also have a significant impact on fluid flows in the vicinity of production well. Actually, it has already been shown that the capillary end effect in waterwet reservoir creates a water layer near the wellbore with residual saturation of gas or oil. This layer will cause the hydraulic resistance for gas or oil flow (Barenblatt et al., 1989; Bedrikovetsky, 2013). Holditch (1979) found out that water phase will be completely trapped in reservoir if the pressure drop between the formation pressure and the bottomhole pressure cannot overcome the capillary end effect. Even when the pressure drop is large enough to overcome the capillary end effect, severe reduction in gas production may also occur, depending on the level of capillary pressure inside the reservoir. Settari and Aziz (1974) analyzed the incompressible twophase flow in the vicinity of the production well. Ding (1995) presented a scalingup procedure for the radial flow near the well which consists in scaling up transmissibility and numerical PI in the vicinity of the well, then he extended this technique to the nearwell upscaling and formation damage modeling which can give quite satisfactory results for longterm productivity prediction (Ding, 2004, 2009, 2011). Naik et al. (2018) established a semianalytical modeling of steadystate gaswater flow in the rock with constant and piecewiseconstant contact angles accounting for gas compressibility and found out that, water saturation significantly rises at the endeffect zone, creating a hydraulic resistance for gas flow. The capillary end effect may have the significant impact on the production well. As far as we know, this capillary end effect has not been considered in the existing well model for reservoir simulation. In this article, we will construct an alternative numerical well model to describe the gas recovery for the lowpermeability waterwet reservoir.
For the convenience of generalization to nonuniform grids and flexible grids, Ding proposed a more accurate and effective numerical well model for the singlephase flow (Ding and Jeannin, 2004; Ding and Renard, 1994; Ding et al., 1998). The technique proposed by Ding et al. is applied. The article is organized as follows. The numerical well models are introduced in Section 2. In Section 3, the 1D radial steady flow for immiscible twophase fluid toward well is investigated. Section 4 describes the procedure to establish the well model for multiphase flow. In Section 5, the accuracy of the proposed well model is verified by numerical examples and the capillary end effect on the gas and water production is analyzed. Conclusion and discussion are given in Section 6.
2 Brief review of the previous numerical well models
In the neighborhood of the production well, the singlephase flow can be approximated by 1D radially symmetric steady flow:
(1)Assuming the isolated well is at (x_{0}, y_{0}), the pressure distribution is given by:
Consider the single phase flow in two dimensional grids with Δx = Δy (see Fig. 1). For the incompressible singlephase flow, the exchange flow rate between the control volume (i, j) and its neighbors can be obtained from the discretization of Darcy’s equation as follows:
(3)where the variables ρ, μ and K represent the fluid density, fluid viscosity and permeability, respectively. If supposing P^{(i+1,j)} = P^{(i−1,j)} = P^{(i, j+1)} = P^{(i,j1)} = P^{(adj)}, the relation equation (3) becomes:
Fig. 1 Block (i, j) containing a well and its four neighboring blocks. 
To establish the numerical well model, an equivalent radius r_{e} is introduced, where the flowing pressure equals P^{(i,j)}. With the boundary condition and P_{r=Δx} = P^{(adj)}, we can integrate equation (1) to get the pressure distribution. The flow rate is then calculated to be
Combing equations (4) and (5), we can get
Integrating equation (1) with the boundary condition and again, we have
(7)where P_{wf} and r_{w} are respectively the bottomhole pressure and the radius of the production well. The relation equation (7) with the definition of equivalent radius equation (6) is just the wellknown Peaceman well model. The above derivations can be found out in the literature (Peaceman, 1978).
Without the solid mathematical analysis, Peaceman well model equation (7) for singlephase flow is directly extended to multiphase flow and then get
(8)where the subscript α represents the different phase, K_{rα} is the relative permeability and h is the layer thickness. Here, r_{e} is still calculated by the relation equation (6). The Peaceman well models equations (6) and (7) are widely used in reservoir simulators.
As pressure varies rapidly logarithmically near the well, its gradient cannot be approximated by linear relation. To improve the accuracy and the efficiency for the prediction, Ding and Renard (1994) proposed an alternative numerical well model based upon the analytical solution equation (2). The finitevolume method applied to the flow equation for the wellblock yields:
Take the interface Γ1 between wellblock and adjacent block (i + 1, j) as example, the flowterm at the interface is
Ding et al. (1998) defined an equivalent interblock distance L_{eq1} to replace the original interblock distance Δx_{1/2}. The pressure P^{(e)} is located at the equivalent radius r_{e}. The numerical flowterm at the interface can be expressed as
(11)where T_{eq,1} is the equivalent transmissibility to relate P^{(e)} and adjacent grid pressure through equivalent interblock distance. From equations (2), (10) and (11), the equivalent interblock distance and equivalent transmissibilities can be respectively given as:
For the other three directions, equivalent interblock distances and equivalent transmissibilities can be calculated using the same way.
Substituting the flowterms at four directions into equation (9), the pressure P^{(e)} located at equivalent radius r_{e} can be given by:
Substituting the P^{(e)} into the analytical solution equation (2), the equation relating pressure, wellbore pressures and the flow rate of the well can be expressed as:
Equations (13) and (14) constitute the new numerical well model for the incompressible singlephase flow. It should be noted that the equivalent radius r_{e} can be any value in this model.
3 Two flow patterns for the onedimensional radial steady flow for immiscible twophase fluid toward well
In this section, the 1D radial steady flow for immiscible two phase fluid (gas and water) toward well is investigated. A condition of constant temperature is supposed in this paper. And the governing equations can be written as:
Only the gas compressibility is considered here, and it varies with the change of gas phase pressure:
Denote S_{wi} as the irreducible water saturation and S_{gr} as the residual saturation of gas phase, then the water saturation S_{w} can be normalized as:
The relative permeability K_{rα} of the αphase can be expressed as the function of the normalized water saturation. That is:
(18)with K_{rw} being a strictly increasing function, but K_{rg} being a decreasing function. The capillary pressure, defined as the difference between the gas and water pressures, is also a decreasing function of the normalized water saturation S:
Capillary pressure is relatively large inside the reservoir but can be ignored inside the wellbore. There are two different capillary pressure curves on the two sides of the entrance of the production well. Analogous to fluid flows through the interface of the different type rocks, physical quantities may change significantly around the interface. Settari and Aziz (1974) pointed out that both water and gas can flow into the production well when the normalized water saturation reaches the value of 1 on the wall of the wellbore and otherwise water phase will be trapped in the reservoir. That means there are two different flow patterns at the vicinity of the wellbore.
Define P_{inlet} and S_{inlet} as the gasphase pressure and normalized waterphase saturation of the point r = r_{inlet} in the vicinity of the production well. When P_{inlet} ∈ (P_{wf} + P_{c}(S_{inlet}), ∞), both two phase can flow into wellbore, i.e., q_{g} > 0 and q_{w} > 0 . When P_{inlet} ∈ (P_{wf}, P_{wf} + P_{c}(S_{inlet})], only the gas can flow into wellbore, i.e., q_{g} > 0 and q_{w} = 0. Details of these two flow patterns are provided in the following sections.
3.1 Flow pattern 1: both two phases flowing into wellbore
In the case where both water and gas flow into the well, we will have:
The reason that the relation equation (20) holds is as follow. If , then and . Thus, . Since there is no capillary pressure in the well, we have . The water pressure must be continuous at r = r_{w}, so . Combined with , the relation will result in infinite velocity of gas phase, which is impossible. Thus, the relation equation (20) holds. Under the condition , we have . Considering the gas flow rate q_{g} > 0, so . This will lead to the relation (see Fig. 2).
Fig. 2 The typical radial distribution of the saturation and the pressures for the flow pattern 1, where q_{g} > 0 and q_{w} > 0. 
3.2 Flow pattern 2: only gas phase flowing into wellbore
If P_{inlet} ∈ (P_{wf}, P_{wf} + P_{c}(S_{inlet})], then we have and . Thus, the water phase is trapped in the reservoir and only the gas phase can flow into the wellbore. In this article, we only consider the production well case, so we have:
(22)where S_{rw} represents saturation at the entrance of wellbore. The saturation and pressure distributions are shown in Figure 3.
Fig. 3 The typical radial distribution of the saturation and the pressure for the flow pattern 2, where q_{g} > 0 and q_{w} = 0. 
Given the pressure P_{inlet} and saturation S_{inlet}, the pressure and saturation in the vicinity of the production well for the two flow patterns can be determined, which will be given in the next section.
4 Modified well model
In Section 2, numerical well models for the incompressible singlephase flow are reviewed. In the following, we will use Ding’s finitevolume method to derive an alternative numerical well model for multiphase flow based upon the analytical equation of the two flow patterns described in the above section.
The finitevolume method applied to the flow equation for the wellblock yields:
(23)where λ_{α} = K_{rα}/μ_{α} is the mobility of the αphase fluid and m is the pseudopressure defined by:
Take the interface Γ1 between wellblock and adjacent grid (i + 1, j) as example. Considering the axisymmetric properties of the phase pressure distributions, from equation (15), the flowterm between the adjacent grid and wellblock grid can be obtained
To let the equivalent interblock distance be phaseindependent, we choose its value L_{eq1} = [Δyln(Δx_{1/2}/r_{e})]/[2arctan(Δy/Δx)], just the same as that in the singlephase flow given by equation (12). Considering that the different phases have the different pressure distributions, we put the flux corrections into the mobility and relative permeability terms, which is expressed as following:
(26)where λ_{w,eq1} and K_{rg,eq1} are the equivalent mobility and equivalent relative permeability, respectively.
With the help of the boundary condition and , we integrate equation (15) from r = r_{e} to r = Δx_{1/2} and then obtain
From equations (25)–(27), the equivalent mobility (λ_{w,eq1}), equivalent relative permeability (K_{rg,eq1}) and equivalent transmissibilities for the water and gas phases can be obtained:
The other three adjacent grid equivalent mobilities (or relative permeability) and equivalent transmissibilities can be calculated as the same method. Substituting flowterms at four directions into equation (23), the pressure P^{(e)} located at r_{e} can be determined by:
After obtaining the pressure P^{(e)} located at equivalent radius, integrating equation (15) on the interval [r_{w}, r_{e}] again in the wellblock grid, we obtain the calculation formula for twophase flow expressed as:
Equations (29) and (30) constitute the twophase flow well model based on finitevolume method. And the equivalent radius can also be any value.
To get the gas and water production from the formula, the saturation distribution on the interval [r_{w}, r_{e}] should be calculated. In the following, we present the procedure to calculate the saturation distribution for two different flow patterns.
4.1 Both of two phase flowing into wellbore (when )
The phase pressures can be given by the following equation:
(31)with the boundary condition:
4.2 Only gas phase flowing into wellbore (when )
Substituting equation (21) into equation (15), we obtain the following saturation equation:
(33)with the boundary condition:
(34)where the saturation S_{rw} is determined by equation (22) and the phase pressures are given by equation (35):
Equations (31)–(35) can be solved numerically. In order to describe the dramatic changes of saturation and pressure near the wellbore, the logarithmic coordinate grids are adopted. In the following numerical tests, the solution interval [r_{w}, r_{e}] is divided into enough equal parts on the logarithmic coordinate. NewtonRaphson iteration method is applied in the solution of the discretized nonlinear equations for equations (31)–(35) to obtain the saturation and pressure distributions for twophase flow and only gas flow, respectively. Practically, a short subroutine can be added into the numerical simulator. Due to solving the ordinary differential equations, the computation speed is very fast. The additional modification has little influence on the program speed.
5 Numerical results
In this section, numerical tests are performed to validate the accuracy of the proposed model. In the TwoDimensional (2D) examples of Test 1 and Test 2, the layer thickness is assumed to be unit thickness, i.e., h = 1. For comparison, the numerical results calculated from the original Peaceman well model equation (7) are also provided. Based upon the calculations, the capillary end effects on the predictions of the gas recovery for the lowpermeability reservoir are analyzed. In the proposed numerical well model the equivalent radius (r_{e}) can be any value. Without loss of generality, in the following calculations, we choose r_{e} = 0.208Δx, just the same as that in the Peaceman well model.
5.1 Test 1
In this numerical test, the simulation area is a 2D square reservoir with the size 500 × 500 m. There are five production wells in the gas reservoir, whose distribution is shown in Figure 4. The outer boundary is imposed with impermeable boundary conditions. Twophase immiscible flow for water and natural gas is considered here. Water is incompressible with invariable density 10^{3} kg/m^{3} and gas is compressible with its density varying with its pressure. The relevant parameters of the reservoir and fluids are shown in Table 1. The capillary pressure curve takes the form P_{c} = −BlnS, where B = 1 MPa. The capillary pressure curves and relative permeability are shown in Figure 5. In this test, we aim to explore the accuracy of proposed well model for the multiwell problem. Since the grid size is much bigger than wellbore radius, the impermeable boundary has little effect on the well model performance.
Fig. 4 The gas reservoir of five production wells. 
Fig. 5 The capillary pressure curves and relative permeability for test 1 and test 2. 
Relevant parameters of the reservoir and fluids.
To test the accuracy of the original Peaceman and the proposed well models, different grid systems of 50 × 50, 150 × 150 and 250 × 250 are adopted in the computations. Due to that the locations of four production wells in the corner have the geometric symmetry, only the calculations of Well 1 and Well 2 are analyzed. For the calculations of the proposed well model, it is shown that the water and gas productions weakly depend upon the grid size, i.e., the calculated results obtained from the different grid sizes are very close (see Fig. 6). However, if adopting the original Peaceman model, there are the obvious deviations among the results for water and gas productions of the different grid sizes.
Fig. 6 a) Gas production of Well 1; b) water production of Well 1; c) gas production of Well 2; d) water production of Well 2. 
The water saturation increasing phenomenon will cause the hydraulic resistance for gas flow and increase the water production. All the grid sizes in this example are much larger than the size of the water layer near the wellbore. Therefore, without considering the capillary end effects, Peaceman well model overestimates the gas production and underestimates the water production under the coarse grid system. But it is indicated that, under the mesh refinement, the calculations of the Peaceman well model gradually approach to those of the proposed well model.
In Test 2, a further analysis of the capillary end effect on the well productions is presented.
5.2 Test 2
In this numerical test, the simulation area is a 2D square reservoir with the size 500 m × 500 m. Parameters are the same as Table 1. A production well is located at the center of the reservoir. In this example, the twophase flow can be simplified to be 1D radial flow since the wellbore radius is much smaller than the computational domain size. It means that we can get the reference solution by solving the ordinary difference equation of equation (15) for two different flow patterns, respectively. Though actually the simulation area is not radial symmetry, the radially symmetric solution can be seemed as a good reference in the neighborhood of the well. To test the accuracy of the original Peaceman and the proposed well models, the different size grids of 50 × 50, 100 × 100, 200 × 200 and 300 × 300 are adopted in the computations. In the following, first we consider the capillary pressure curve with B = 1 MPa in the calculations for Figures 7 and 8 and then investigated the influence the capillary pressure curve with different values of the constant B in the calculations for Figures 9 and 10.
Fig. 7 a) Gas phase flux profiles; b) water phase flux profiles; c) total flux of gas phase; d) total flux of water phase at different grid size for different well models. 
Fig. 8 a) Saturation distribution and b) pressure distribution near the well calculated by reference solution in 300th day. 
Fig. 9 The comparisons between the calculations from the proposed well model and Peaceman well model for different capillary pressure curves in the size grids of 100 × 100. 
Fig. 10 a) Saturation distributions and b) pressure distributions near the well calculated by reference solution with different capillary pressures in 300th day. 
Figure 7 shows that all the calculations with different grid size are agreed with the reference solution very well when adopting the proposed well model. Nevertheless, adopting Peaceman well model, the calculations of the coarse grid system 50 × 50 deviate from the reference solution even up to 40%. With the mesh subdivision, the result of Peaceman well model is gradually approaching to the reference solution. The calculations indicate that extremely fine grid is needed if using the original Peaceman well model. The capillary end effect has the important influence for the prediction of the lowpermeability reservoir recovery. Due to the capillary end effect, the saturation and the gasphase pressure change rapidly in the vicinity of the production well and the saturation and the gasphase pressure gradients may even diverge at the entrance of the wellbore (see Fig. 8), which is the same as Figure 2. In the original Peaceman well model, these dramatic changes are not considered and therefore may cause a poor prediction for the lowpermeability reservoir recovery.
To investigate the influence of the strength of the capillary pressure upon the predictions of the original Peaceman well model, we perform the calculations with the grid number of 100 × 100 for the different values of the constant B = 0.1 MPa, 1 MPa and 3 MPa. When the capillary pressure is much smaller than the pressure difference between the formation pressure and the well bottomhole pressure, for instance, in the case of B = 0.1 MPa, Peaceman well model can provide a rather accurate result (see Fig. 9). However, along with the increase of the strength of the capillary pressure, zone affected by capillary pressure is also increasing (see Fig. 10) and therefore the gap between the calculations of Peaceman well model and the proposed model becomes larger and larger. When B = 3 MPa, the water phase is blocked inside the reservoir and then cannot flow out of the production well. In this case, the original Peaceman well model yield an incorrect calculation if the grid size is not enough small.
6 Conclusion
In the recovery of low permeability reservoirs, the capillary pressure has an important effect, which may reduce gas production. When pressure difference between the formation pressure and the well bottomhole pressure is not large enough to overcome the capillary end effect, only the gas phase can flow into the wellbore and cause the gas production significantly reduced. When the pressure drop is enough to overcome the capillary end effect, both water and gas phases can flow into the wellbore. In this case, the saturation and gasphase pressure gradients diverge in the vicinity of wellbore. To construct an effective numerical well model for gas reservoir simulation is the main objective of this article. In this article, other nearwellbore processes, such as perforation damage, fines movement, nearwellbore fracturing and wettability changes introduced by stimulation, are not considered. The conclusion can be summarized as follows:

Considering the capillary end effect, there are two flow patterns in the vicinity of the production well: both water and gas phase can flow into the production well and only gas phase can enter the production well. With the help of the solution of the 1D radially symmetric flow for twophase immiscible flow, an alternative well model for low permeability gas reservoir is constructed, where two different flow patterns are considered.

When the capillary pressure is much smaller than the pressure difference between the formation pressure and the well bottomhole pressure, the original Peaceman well model can provide a rather accurate result. The Peaceman well model of multiphase flow is extended directly from that of singlephase flow, where the capillary end effect is not considered. As a result, along with the increase of the strength of the capillary pressure, the error of the predictions of Peaceman well model becomes larger and larger. The saturation and the gasphase pressure changes rapidly in the vicinity of the production well and their gradients may even diverge at the entrance of the production well. As a result, the original Peaceman well model can provide good prediction only when the grid size is enough small. This problem is out of control since it is difficult to estimate the suitable grid size for the original Peaceman well model.

In contrast to the original Peaceman well model, the proposed well model can provide rather accurate prediction for gas and water productions. The prediction weakly depends on the grid size adopting in the computation, since the dramatic change for saturation and gasphase pressure in the vicinity of the wellbore has been already considered in the proposed well model. In the recovery of the lowpermeability gas reservoirs, the capillary end effect is the main reason for waterphase trapping. The increasing water phase saturation reduces gas mobility near wellbore and then affects well productivity and gas recovery. The proposed well model can accurately describe the capillary end effect in the recovery procedure for lowpermeability recovery.
For the singlephase flow, Ding et al. (1998) established the well model based upon the finitevolume method by using the axisymmetric analytical solution. In this article, flow patterns in lowpermeability gas reservoirs are obtained from the near well axisymmetric equation. With the help of the solution equations (31)–(35) of the 1D radially symmetric flow for the twophase flow, the well model for low permeability gas reservoir in cartesian mesh is established, based upon Ding’s method. It is expected that the proposed method can be extended to flexible grids based upon the nearwell axisymmetric solutions. This remains the future research.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (No. 11572315).
References
 Babu D.K., Odeh A.S. (1988) Productivity of a horizontal well appendices A and B, SPE Annual Technical Conference and Exhibition, 2–5 October, Houston, Texas. [Google Scholar]
 Babu D.K., Odeh A.S., AlKhalifa A.J., McCann R.C. (1991) Numerical simulation of horizontal wells, Middle East Oil Show, 1619 November, Bahrain. [Google Scholar]
 Barenblatt G.I., Entov V.M., Ryzhik V.M. (1989) Theory of fluid flows through natural rocks, Kluwer, Dordrecht. [Google Scholar]
 Bedrikovetsky P. (2013) Mathematical theory of oil and gas recovery: with applications to exUSSR oil and gas fields, Vol. 4, Springer Science & Business Media, Dordrecht. [Google Scholar]
 Chen Z., Zhang Y. (2009) Well flow methods for various numerical methods, Int. J. Numer. Anal. Mod. 6, 3, 375–388. [Google Scholar]
 Ding D.Y. (1995) Scalingup in the vicinity of wells in heterogeneous field, 13th SPE Symposium on Reservoir Simulation, 12–15 February, San Antonio, TX. [Google Scholar]
 Ding D.Y. (2004) Nearwell upscaling for reservoir simulations, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 59, 2, 157–165. [CrossRef] [Google Scholar]
 Ding D.Y. (2009) Modeling formation damage for flow simulations at reservoir scale. 8th European Formation Damage Conference, 27–29 May, Scheveningen, The Netherlands. [Google Scholar]
 Ding D.Y. (2011) Coupled simulation of nearwellbore and reservoir models, J. Petrol. Sci. En. 76, 1–2, 21–36. [CrossRef] [Google Scholar]
 Ding D.Y., Jeannin L. (2004) New numerical schemes for nearwell modeling using flexible grid, SPE J. 9, 1, 109–121. [CrossRef] [Google Scholar]
 Ding D.Y., Renard G. (1994) A new representation of wells in numerical reservoir simulation, SPE Reserv. Eng. 9, 2, 140–144. [CrossRef] [Google Scholar]
 Ding D.Y., Renard G., Weill L. (1998) Representation of wells in numerical reservoir simulation, SPE Reserv. Eval. Eng. 1, 1, 18–23. [CrossRef] [Google Scholar]
 Ding Y.X., Shi A.F., Luo H.S., Wang X.H. (2016) Adaptive mesh refinement for nonisothermal multiphase flows in heterogeneous porous media comprising different rock types with tensor permeability, Numer. Heat Transf. A Appl. 69, 1, 31–50. [CrossRef] [Google Scholar]
 Hadley G.F., Handy L.L. (1956) A theoretical and experimental study of the steady state capillary end effect. Fall Meeting of the Petroleum Branch of AIME, 14–17 October, Los Angeles, CA. [Google Scholar]
 Holditch S.A. (1979) Factors affecting water blocking and gas flow from hydraulically fractured gas wells, J. Pet. Technol. 31, 12, 1515–1524. [CrossRef] [Google Scholar]
 Huang D.D., Honarpour M.M. (1998) Capillary end effects in coreflood calculations, J. Petrol. Sci. En. 19, 1–2, 103–117. [CrossRef] [Google Scholar]
 Naik S., You Z., Bedrikovetsky P. (2018) Productivity index enhancement by wettability alteration in twophase compressible flows, J. Nat. Gas Sci. Eng. 50, 101–114. [Google Scholar]
 Peaceman D.W. (1978) Interpretation of wellblock pressures in numerical reservoir simulation (includes associated paper 6988), Soc. Pet. Eng. J. 18, 3, 183–194. [CrossRef] [Google Scholar]
 Peaceman D.W. (1983) Interpretation of wellblock pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability, Soc. Pet. Eng. J. 23, 3, 531–543. [CrossRef] [Google Scholar]
 Richardson J.G., Kerver J.K., Hafford J.A., Osoba J.S. (1952) Laboratory determination of relative permeability, J. Pet. Technol. 4, 8, 187–196. [CrossRef] [Google Scholar]
 Settari A., Aziz K. (1974) A computer model for twophase coning simulation, Soc. Pet. Eng. J. 14, 3, 221–236. [CrossRef] [Google Scholar]
 Shiralkar G.S. (1989) Calculation of flowing well pressures in reservoir simulation using ninepoint differencing, J. Can. Petrol. Technol. 28, 6, 73–82. [CrossRef] [Google Scholar]
 Su H.J. (1995) Modeling of offcenter wells in reservoir simulation, SPE Reserv. Eng. 10, 1, 47–51. [CrossRef] [Google Scholar]
 Virnovsky G.A., Skjaeveland S.M., Surdal J., Ingsoy P. (1995) Steadystate relative permeability measurements corrected for capillary effects. SPE Annual Technical Conference and Exhibition, 22–25 October, Dallas, TX. [Google Scholar]
 Wang X.H., Quintard M., Darche G. (2006) Adaptive mesh refinement for onedimensional threephase flow with phase change in porous media, Numer. Heat Transf. B Fundam. 50, 3, 231–268. [CrossRef] [Google Scholar]
 Williamson A.S., Chappelear J.E. (1981) Representing wells in numerical reservoir simulation: part 1Theory, Soc. Pet. Eng. J. 21, 3, 323–338. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Block (i, j) containing a well and its four neighboring blocks. 

In the text 
Fig. 2 The typical radial distribution of the saturation and the pressures for the flow pattern 1, where q_{g} > 0 and q_{w} > 0. 

In the text 
Fig. 3 The typical radial distribution of the saturation and the pressure for the flow pattern 2, where q_{g} > 0 and q_{w} = 0. 

In the text 
Fig. 4 The gas reservoir of five production wells. 

In the text 
Fig. 5 The capillary pressure curves and relative permeability for test 1 and test 2. 

In the text 
Fig. 6 a) Gas production of Well 1; b) water production of Well 1; c) gas production of Well 2; d) water production of Well 2. 

In the text 
Fig. 7 a) Gas phase flux profiles; b) water phase flux profiles; c) total flux of gas phase; d) total flux of water phase at different grid size for different well models. 

In the text 
Fig. 8 a) Saturation distribution and b) pressure distribution near the well calculated by reference solution in 300th day. 

In the text 
Fig. 9 The comparisons between the calculations from the proposed well model and Peaceman well model for different capillary pressure curves in the size grids of 100 × 100. 

In the text 
Fig. 10 a) Saturation distributions and b) pressure distributions near the well calculated by reference solution with different capillary pressures in 300th day. 

In the text 