Nonlinear flow model for well production in an underground formation

Fluid flow in underground formations is a nonlinear process. In this article we modelled the nonlinear transient flow behaviour of well production in an underground formation. Based on Darcy’s law and material balance equations, we used quadratic pressure gradients to deduce diffusion equations and discuss the origins of nonlinear flow issues. By introducing an effective-well-radius approach that considers skin factor, we established a nonlinear flow model for both gas and liquid (oil or water). The liquid flow model was solved using a semi-analytical method, while the ga flow model was solved using numerical simulations because the diffusion equation of gas flow is a stealth function of pressure. For liquid flow, a series of standard log-log type curves of pressure transients were plotted and nonlinear transient flow characteristics were analyzed. Qualitative and quantitative analyses were used to compare the solutions of the linear and nonlinear models. The effect of nonlinearity upon pressure transients should not be ignored. For gas flow, pressure transients were simulated and compared with oil flow under the same formation and well conditions, re ulting in the c clusion that, under the same volume rate production, oil wells demand larger pressure drops than gas wells. Comparisons between theoretical data and field data show that nonlinear models will describe fluid flow in underground formations realistically and accurately.


Introduction
Among the many nonlinear geophysical processes, transient fluid flow through porous media is of particular interest (Cao et al., 2004;Finjord, 1990;Finjord and Aadnoy, 1989;Giachetti and Maroscia, 2008;Liang et al., 2001).This process, which is of practical importance, is governed by the diffu-sivity equation, an equation describing the nonlinearities resulting from the dependence of fluid and medium properties on pressure.When porosity, permeability and fluid density depend exponentially on pressure, the diffusivity equation reduces to a diffusion equation containing a squared gradient term.Many published articles have described analytical solutions for this equation through variable modifications (Chakrabarty et al., 1993a, b;Jelmert and Vik, 1996;Odeh and Babu, 1998;Wang and Dusseault, 1991), which are special cases of the Hopf-Cole transformation (Marshall, 2009).Applications in dual-porosity (Bai et al., 1994;Bai and Roegiers, 1994) and fractal (Tong and Wang, 2005) media have also been described.Pressure transients are graphic plots of theoretical solutions for diffusion equations under specific initial and boundary conditions for the well-test model that represents a reservoir-well system (Chaudhry, 2004).Pressure transients were plotted according to these analytical solutions and were used in the interpretation of well-test transients (Braeuning et al., 1998).This research showed that the quadratic pressure gradient term influenced the pressure transient solutions.If the nonlinear term is ignored, significant errors in predicted pressures during certain live oil and low permeability reservoir operations, such as hydraulic fracturing, large-drawdown flows, slug testing, drill-stem testing and large-pressure pulse testing will occur.Even though the importance of nonlinearity has been heavily emphasised, no application of modern well-test interpretation (Onur et al., 2003) on field data has been found; in fact, a standard set of log-log type curves for modern well-test interpretation has not been developed, except for type curves for nonlinear spherical flow (Nie and Ding, 2010).In addition, there has not been any research or discussion on the effects of quadratic pressure gradients on gas flow.
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Therefore, the main tasks of this article is to (1) deduce the nonlinear governing equations with the quadratic pressure gradient term for both liquid and gas flow in porous media and discuss the origins of the nonlinear flow issue; (2) establish and solve a nonlinear flow model for well production; (3) plot a series of standard type curves for wellbore pressure-transient analysis for liquid flow in an underground formation; (4) analyze the nonlinear flow characteristics of liquid flow using type curves; (5) make qualitative and quantitative comparisons of pressure and its derivative transients between nonlinear and linear liquid flow models; (6) simulate wellbore pressure transients for nonlinear gas flow using numerical simulations and compare these pressure transients with nonlinear oil flow under the same formation and well conditions; and (7) match theoretical data from nonlinear and linear models against field data to assess the differences between actual and theoretical applications.The results of this research suggest the use of a nonlinear flow model in actual studies.
Compared with previous publications, such as Cao et al. (2004), Chakrabarty et al. (1993b), Marshall (2009), this research highlights (1) newly considered use of real parameters for wells with skin factor in nonlinear flow models; (2) application of modern standard type curves to intuitively observe nonlinear transient flow behaviour; (3) recognition of flow regimes from type curves, including recognition of different external boundary responses; (4) thorough analysis of parameter sensitivity to type curves; (5) use of quantitative methods of "DV" and "RDV" to describe solution differences between nonlinear and linear models; (6) establishment of numerical modelling of nonlinear gas flow and simulation and comparison of nonlinear gas and oil flow pressure transients; and (7) consideration of real world applications through comparisons of theoretical data and field data.

Liquid flow
For vertical well production in a homogeneous formation, flow in a vertical plane is not significant and a radial cylindrical coordinate system without a z coordinate can be employed to describe the diffusion equation: where p is pressure, MPa; r is radial cylindrical coordinate, cm; t is time, s; ϕ is rock porosity, fraction; k is radial permeability, µm 2 ; µ is viscosity, mPa s; C t is total compressibility of rock and liquid, MPa −1 ; C ρ is liquid compressibility, MPa −1 .The governing partial differential equation is nonlinear because of the quadratic pressure gradient term (i.e. the second power of the pressure gradient in Eq. 1).
The appendix deduces the diffusivity equation containing the quadratic gradient term when porosity and fluid density depend exponentially on pressure (Marshall, 2009;Nie and Ding, 2010): where ρ is density, g cm −3 ; C f is rock compressibility, MPa −1 ; ρ 0 , ϕ 0 , p 0 are reference values, which are usually used in standard conditions.
The function e x using Maclaurin series expansion can be written by If we use Maclaurin series expansion for Eqs. ( 2) and (3) and neglect second and higher order items, Eqs. ( 5) and ( 6) can replace Eqs. ( 2) and (3), respectively Appearance of the quadratic pressure gradient term is due to a lack of simplification in the state equations (Eqs. 2 and 3) when deducing the diffusion equation.If we use Eqs. 5 and (6), we obtain the conventional linear flow equation without the quadratic pressure gradient term.Therefore, the linear flow equation is an approximation and simplification of the nonlinear flow equation which includes the quadratic pressure gradient term.In fact, flow of liquid in porous media is a complex nonlinear process and the nonlinear flow law is equal to the flow law of liquid in porous media.

Gas flow
Gas flow in porous media is different from that of liquid due to a different state equation for fluid (Nie et al., 2012a): where V is gas volume, cm 3 ; Z is a compressibility factor, fraction; n is number of gas moles, mol; R is a universal gas constant, J (mol K) −1 ; T is gas temperature, K. Gas volume is a function of mass and density: where m is gas mass, g; ρ is density, g cm where M is average molecular weight of gas, g mol −1 .Substitute Eqs. ( 8) and ( 9) into Eq.( 7): We consider isothermal equations and Darcy's flow and substitute Eq. ( 10) into Eq.(A3): where x, y, and z are Cartesian coordinates; k x is permeability in the x direction, µm 2 ; k y is permeability in the y direction, µm 2 ; k z is permeability in the z direction, µm 2 .We assume an elastic and slightly compressible rock and isotropic and constant permeability in both the horizontal and vertical plane: where k h is horizontal permeability that is equal to k x and k y is isotropic formation in the horizontal plane, µm 2 .The partial differential terms in the left of Eq. ( 12) are expressed by The partial differential term in the right of Eq. ( 12) is expressed by Gas isothermal compressibility is a function of pressure and compressibility factor: Substitute Eqs.(13-15) into Eq.( 12): where C g is gas compressibility, MPa −1 ; C t is total compressibility of rock and gas, MPa −1 .Equation ( 16) is a diffusion equation containing the quadratic derivative term of pressure square in Cartesian coordinates for gas flow in a homogeneous formation.Equation (16) shows that flow of gas in porous media is a nonlinear process.Compared with the governing differential equation of liquid (Eq.1), Eq. ( 16) shows a more complex nonlinear properties.It is difficult to solve this diffusion equation because the term of ∂ ln (µZ) /∂p 2 is a stealth function of pressure and not a constant.Usually the pseudo-pressure (or "potential") function (Ertekin and Sung, 1989;King and Ertekin, 1988;Nie et al., 2012a) can be used to describe the governing equation of gas flow: where ψ is gas pseudo-pressure, MPa 2 (mPa s −1 ) −1 ; p sc is pressure at standard conditions, MPa.Derivatives of pseudo-pressure to coordinates can be expressed by Derivative of pseudo-pressure to time can be expressed by Substitute Eq. ( 20) into the right of Eq. ( 12): Substitute Eq. ( 19) into the left of Eq. ( 12): www.nonlin-processes-geophys.net/20/311/2013/ Nonlin.Processes Geophys., 20, 311-327, 2013 Then, the diffusion equation expressed by gas pseudopressure can be obtained by Equation ( 23) is a diffusion equation without a quadratic derivative term.According to Eqs. ( 15) and ( 17), total compressibility of rock and gas is a stealth function of pressure because gas compressibility is a stealth function of pressure.
In addition, gas viscosity is also a stealth function of pressure.Therefore, Eq. ( 23) is still a nonlinear equation.This nonlinear equation is very difficult to solve using analytical methods; however, a numerical approach will be implemented.
3 Description of the nonlinear flow model

Physical model
Physical model assumptions: 1. Single vertical well production at a constant rate in a homogeneous and isotropic formation saturated by a single-phase fluid (gas, oil or water) and the external boundary of formation may be infinite, closed or of constant pressure.
2. Slightly compressible rock and liquid (oil or water) with a constant compressibility are considered, while the compressibility of gas changes with depletion of pressure.
3. Isothermal equations and Darcy's flow, ignoring the impact of gravity and capillary forces.
4. Wellbore storage effect is considered when the well is opened, while fluid stored in the wellbore starts to flow and when fluid in formation does not.
5. Skin effect is considered near the wellbore where the formation could be damaged by drilling and completion operations (there could be an additional pressure drop during production, with the "skin" being a reflection of additional pressure drop).
6.At time t = 0, pressure is uniformly distributed in formation, equal to initial pressure (p i ).

Mathematical model
For the convenience of well-test analysis, the mathematical model was established using a set of applied engineering units.

Liquid flow in a formation (1) Establishment of mathematical model
The governing differential equation in a radial cylindrical system where C ρ is liquid compressibility, MPa −1 ; C t is total compressibility of rock and liquid, MPa −1 ; k h is radial formation permeability, µm 2 ; p is formation pressure, MPa; r is radial radius from the centre of wellbore, m; t is well production time, h.Initial condition where p i is initial formation pressure, MPa.
Well production condition based on effective radius kh µ r ∂p ∂r r=r wa = 1.842 where B is oil volume factor, dimensionless; C s is wellbore storage coefficient, m 3 MPa −1 ; p w is wellbore pressure, MPa; q is well rate at wellhead, m 3 d −1 ; r wa is effective wellbore radius, m.The effective wellbore radius r wa is defined as (Agarwal et al., 1970;Chaudhry, 2004) where r w is real wellbore radius, m; S is skin factor, dimensionless.
External boundary conditions: where r e is external boundary radius, m.
The following dimensionless definitions are introduced to solve the mathematical model: Dimensionless wellbore storage coefficient: Dimensionless time: t D = 3.6kt/ ϕµC t r 2 w .
Dimensionless coefficient of nonlinear term: The dimensionless mathematical model is as follows: Governing differential equation in a radial cylindrical system equals Initial condition: Well production condition: where p wD is dimensionless wellbore pressure.External boundary conditions: (2) Linearization of dimensionless mathematical model Equation ( 32) is a nonlinear partial differential equation.
The following variable modifications are introduced to solve the dimensionless mathematical model (Nie and Ding, 2010;Odeh and Babu, 1998): Then Substitute Eqs. ( 38)-( 42) into Eqs.( 31)-( 37), the model can be converted to ξ r D =r eD = 0 (constant pressure) , (47) (3) Solution to dimensionless mathematical model Introduce the Laplace transform based on T D : The dimensionless mathematical model in Laplace space is as follows: The general solution of the Eq. ( 50) is Substitute Eq. ( 55) into Eqs.( 50) and ( 51): Substitute Eq. ( 55) into Eqs.( 52)-( 54): where A and B are undetermined coefficients; I 0 ( ) is a modified Bessel function of the first kind, zero order; I 1 ( ) is a modified Bessel function of the first kind, first order; K 0 ( ) is a modified Bessel function of the second kind, zero order; K 1 ( ) is a modified Bessel function of the second kind, first order.In Eqs. ( 57)-( 61), there are three unknown numbers (A, B, ξ w ) and three equations, solutions to the model in Laplace space can be easily obtained by using linear algebra (Nie et al., 2011a, b), such as a Gauss-Jordan reduction.
In real space, ξ w and the derivative (dξ w /dT D ) can be obtained using a Stehfest numerical inversion (Stehfest, 1970) to convert ξ w back to ξ w , and then dimensionless wellbore pressure (p wD ) and the derivative (dp wD /dT D ) can be obtained by substituting ξ w into Eq.( 38).The standard log-log type curves of well-test analysis (Nie et al., 2012b, c) of p wD and (p wD • t D /C D ) vs. t D /C D can then be obtained.

Gas flow in a finite formation (1) Establishment of mathematical model
Governing differential equation in a radial cylindrical system where C t is total compressibility of rock and gas, MPa −1 ; k is radial formation permeability, µm 2 ; ψ is gas pseudopressure, MPa 2 (mPa s) −1 ; r is radial radius from the centre of wellbore, m; t is well production time, h.Initial condition where the subscript "i" means "initial".
Well production condition based on effective radius where C s is a wellbore storage coefficient, m 3 MPa −1 ; ψ w is wellbore pseudo-pressure, MPa 2 (mPa s) −1 ; q g is well rate at bottom hole, m 3 d −1 ; r wa is effective wellbore radius, m.At the external boundary for a finite formation: where r e is the external boundary radius, m. (

2) Solution of the nonlinear mathematical model
Because well production pressure is largely depleted near the wellbore, we use logarithmic-uniform radial meshes in space to discretize the equation and obtain relatively dense grids near the wellbore.A new space variable is taken by The mathematical model can be converted to The central difference method is employed to solve the model: where n is the previous time level; (n + 1) is the current time level; 0 is the initial time; j is the grid node; J is the number of grid node at external boundary; t is the time step size; R is the grid step size.Note that Eqs. ( 68) and ( 69) are nonlinear equations because µ, Z, C t are stealth functions of pressure and pseudopressure, therefore, in order to solve the model, we used the previous time level (µC t ) n j in Eq. ( 73), and (1) Type curves of the linear flow model The standard log-log type curves of the linear flow model (see Fig. 1) are well known.An entire flow process can be discerned from the type curves: i. Regime I, pure wellbore storage regime, slope of pressure and the pressure derivative both equal one.
ii. Regime II, wellbore storage and skin effect regime, shape of the derivative curve is a "hump".
iii.Regime III, radial flow regime, slope of the pressure derivative curve equals zero, and all the pressure derivative curves converge to the "0.5 line", indicating the logarithmic value of the pressure derivative is 0.5.iii.Regime III, radial flow regime, different type curves are associated with the two flow models.Pressure derivative curves of the nonlinear flow model do not conform to the "0.5 line" law.The curves are located along the "0.5 line" and are inclined instead of horizontal (see Figs. 2-5).As time elapses, the pressure derivative curves gradually deviate from the "0.5 line".
where DV is the differential value between linear and nonlinear models and RDV is the relative differential value between linear and nonlinear models.
Tables 1 and 2 show the quantitative differences of nonlinear influence on type curves for "β = 0.01" and "β = 0.1", respectively.Dimensionless pressure values and their derivative values in Tables 1 and 2 were calculated by setting C D e 2S as 10 3 , the corresponding type curves are shown in Fig. 3.The tables show that dimensionless pressure and its derivative differ between linear and nonlinear models.For "β = 0.01" and "t D /C D = 10 2 ", DV and RDV of pressure are 0.1711 and 2.81 %, respectively; when "t D /C D = 10 5 ", DV and RDV of pressure are 0.4311 and 4.48 %, respectively; when "t D /C D = 10 2 ", DV and RDV of pressure derivative are 0.0348 and 6.17 %, respectively; when "t D /C D = 10 5 ", DV and RDV of pressure derivative are 0.0439 and 8.78 %, respectively.For "β = 0.1" and "t D /C D = 10 2 ", DV and RDV of pressure are 1.3030 and 21.38 %, respectively; when "t D /C D = 10 5 ", DV and RDV of pressure are 2.8667 and 29.82 %, respectively; when "t D /C D = 10 2 ", DV and RDV of pressure derivative are 0.2212 and 39.24 %, respectively; when "t D /C D = 10 5 ", DV and RDV of pressure derivative are 0.2457 and 49.14 %, respectively.DV and RDV increase over time (Tables 1 and 2).RDV of the pressure derivative is greater than that of pressure at a fixed time, such as when "t D /C D = 10 3 " for "β = 0.01" (Table 1).The RDV of the pressure derivative is 6.91 %, which is greater than that of pressure, 3.43 %.It can be observed from the tables that DV and RDV increase in parallel with β, such as when "t D /C D = 10 4 ", RDV of pressure for "β = 0.01" is 3.96 % and RDV of pressure for "β = 0.1" is 27.39 %.
According to the equation, β = (1.842× 10 −3 qBµC ρ )/(kh), and probable values of β (Table 3), it is demonstrated that β is proportional to liquid viscosity, µ, and inversely proportional to formation permeability, k, and formation thickness, h.Formations with low permeability, heavy oil or thin thickness have a larger β, causing the influence of nonlinearity to be more intense for these formations.
In general, flow of fluid in porous media is a nonlinear process and the nonlinear quadratic pressure gradient term should be retained in diffusion equations.

Comparision to Chakrabarty's model
Homogenous reservoir models are among recent pressuretransient models containing the quadratic pressure gradient term.However, these models do not study pressure-transient curves of well-test analysis, such as Marshall (2009), Giachetti and Maroscia (2008), Liang et al. (2001); while others do not model homogenous reservoirs, such as Tong and Wang (2005), Bai et al. (1994).The solutions of these models cannot be compared to the model discussed here.Chakrabarty et al. (1993b)    Explanations: k is formation permeability; µ is liquid viscosity; β is dimensionless coefficient of the nonlinear term; β was calculated according to its definition under a set of fixed parameters: set liquid rate q = 25 m 3 d −1 ; liquid volume factor B = 1.004; formation thickness h = 10 m and liquid compressibility C ρ = 0.001 MPa −1 .
The model of this article considered skin factor and regulated skin factor (S) and wellbore storage coefficient (C D ) to a parameter group (C D e 2S ).S would need to be set as zero to make an effective comparison with Chakrabarty's model.In addition, the log-log pressure-transient curves were plotted based on (t D /C D ), therefore the abscissa would need to be converted to t D .Please note that the type curves based on (t D /C D ) will cross over the origin of coordinates (10 −2 , 10 −2 ) and the type curves based on t D do not.Figures 8 and 9 of Chakrabarty et al. (1993b) are log-log pressure curves and log-log pressure derivative curves, respectively.In order to have a convenient comparison, the same parameter values and the same range of coordinate scales as those shown Figs. 8 and 9 were used.Equation (6a) of Chakrabarty et al. (1993b) defined the dimensionless quadratic gradient coeffi-cient as where α is the dimensionless quadratic gradient coefficient; q is rate at wellbore, cm 3 s −1 ; µ is viscosity, mPa s; c is fluid compressibility, atm −1 ; k is permeability, D; h is formation thickness, cm.
Comparing the definition of α in Chakrabarty et al. (1993b) with the definition of β in this article, they used the rate at wellbore and we used the rate at wellhead, and β is the opposite number of α.Figures 6 and 7 show the comparisons of pressure curve and pressure derivative curve, respectively.The solutions of our model are completely the same as those of the model in Chakrabarty et al. (1993b).
In conclusion, our model can be reduced to the model of   (1993b).β is the opposite number of α.For convenience, the same parameter values are set.The solutions of this article are the same as those of Chakrabarty et al. (1993b).Chakrabarty et al. (1993b) by setting skin factor (S) as zero.
However, our research is different from that of Chakrabarty et al. (1993b).Significant improvements are (1) our model considers well skin factor (S) using the effective well radius method (Agarwal et al., 1970;Chaudhry, 2004), while Chakrabarty's model did not, therefore our model is more realistic; (2) we used modern standard type curves (Corbett et al., 2012) based on (t D /C D ) to regulate pressure and its derivative curves for different external boundaries making it convenient to use these curves to observe nonlinear transient flow behaviour.The analysis curves of Chakrabarty et al. (1993b) were based on t D and did not use modern standard type curves; (3) different flow regimes were recognised from the researched type curves and nonlinear flow characteristics in every flow regime were analyzed.Chakrabarty et al. (1993b) did not follow this procedure; (4) we thoroughly analyzed parameter sensitivities to type curves, including the nonlinear coefficient (β) and the parameter group (C D e 2S ); Chakrabarty et al. (1993b) did not; (5) we used "DV" and "RDV", clearly showing the difference in quantitative data between nonlinear and linear solutions; Chakrabarty et al. (1993b) did not; and (6) we researched real world applications by matching the theoretical nonlinear model against field data, but Chakrabarty et al. (1993b) did not.The above analysis methods used to research nonlinear flow issues of liquid in an underground formation are different from previous publications, such as Cao et al. (2004), Chakrabarty et al. (1993b), Marshall (2009).

Data preparation
PVT parameters of gas can be calculated from gas composition (Hagoort, 1988).Using an example gas well of sandstone formation, basic data of formation and gas composition are shown in Table 4 (formation at mid-depth is 1380 m; formation thickness is 9.4 m; formation porosity is 0.08; formation permeability is 0.014 µm 2 ; compressibility of formation rock is 0.00038 MPa −1 ; the well radius is 0.062 m; formation temperature is 338.5 K; formation pressure is 20.5 MPa; and mole composition of methane, ethane, propane and nitrogen are 99.711%, 0.092 %, 0.03 % and 0.167 %, respectively).Formation parameter data are from well logging interpretations and mole composition data are from laboratory tests.Gas property parameters are a function of pressure and temperature and were calculated using mole composition data (such as gas compressibility, viscosity, compressibility factor and pseudo-pressure) in the commercial software, Saphir.Figure 8 shows the relationship curves between gas property parameters and pressure.The gas compressibility factor decreases and then increases with an increase in pressure (Fig. 8a).Both gas viscosity and pseudopressure increase with an increase in pressure (Fig. 8b and c).
Gas compressibility decreases with an increase in pressure (Fig. 8d).Compressibility is large at low pressures, especially in critical isotherm conditions (Marshall, 2009).Critical pressures calculated using Saphir are 4.62 MPa with gas compressibility of the example well = 0.024 MPa −1 .Gas flow in formation is a nonlinear process as shown by changes in property parameters with pressure.For convenience, the following formulations were regressed using polynomial or   .Relationship curves of gas viscosity (µ), gas compressibility (C g ), compressibility factor (Z) and pseudo-pressure (ψ) with pressure (p) for the example well.Parameters in the (3 ∼ 25) MPa pressure range were calculated using gas compositions found in the commercial software, Saphir.In order to conveniently calculate these parameters in numerical simulations, relationship equations were regressed with pressure using polynomial or power functions in Microsoft Excel.
power functions in Microsoft Excel: In numerical simulations, each value of a property parameter in the current time level was approximated using a value associated with a previous time level.Each value of a property parameter in the previous time level can be calculated using Eqs.( 81)-( 83).Pressure transients can be simulated simulta-neously according to the relationship formulation of pseudopressure with pressure (Eq.84).
In order to compare gas and oil flow, numerical simulations must be done under the same formation and well production conditions.Wellhead gas rate of the example well is 4.3 × 10 4 m 3 d −1 , and the bottom hole gas rate is equal to the product of the gas volume factor B g (0.006) with a wellhead rate of 259 m 3 d −1 .If one was to simulate oil flow in such a sandstone formation with 9.4 m thickness and 0.08 porosity, the oil rate at the bottom hole would not reach 259 m 3 d −1 .Therefore, we simulated pressure transients by setting a relatively small bottom hole rate (20 m 3 d −1 ) for both gas and oil production.With similar formation and well production conditions, property parameters of oil must be prepared before Nonlin.Processes Geophys., 20, 311-327, 2013 www.nonlin-processes-geophys.net/20/311/2013/ numerical simulations.We set oil viscosity and compressibility as 2 mPa s and 0.0034 MPa −1 , respectively and ignored changes in oil viscosity and compressibility with depletion of pressure.In addition, the wellbore storage coefficient and skin factor were set as 0.001 m 3 MPa −1 and 0.5, respectively, with a 3000 m closed external boundary for the formation.Substituting the above data into Eqs.( 57)-( 61) for oil wells and Eqs. ( 73)-( 77) for gas wells, wellbore pressure transients were calculated under the above conditions.

Numerical simulation results
Figure 9 shows the numerical simulation results for both gas and oil well production.The horizontal coordinate represents production time (t) and the vertical coordinate represents production pressure drop ( p), or production pseudopressure drop ( ψ).Simulated production time is 60 h. Figure 9 shows (1) the location of pseudo-pressure and its derivative curves (curve 1 ) are at their highest because the value of pseudo-pressure is larger than that of pressure (see Table 5); and (2) the location of pressure and its derivative curves for oil well production (curve 2 ) is higher than that of gas well production because oil wells need a larger pressure drop for the same rate of production.The model assumes well production is at a constant volume rate.The very high dilatability of gas with the depletion of pressure makes the volume supply of gas to well rate easier than that of oil; therefore, the pressure drop of gas well production is lower than that of oil well production.Figure 9 also shows the difference between external boundary response time for oil and gas well production.A larger pressure drop for oil well production results in a faster propagation of pressure waves, therefore propagating time for a pressure wave to reach the closed external boundary for oil well production (t eo ) is shorter than that of gas well production (t eg ).According to the simulated results, external boundary response time of oil and gas flow is 0.6 h and 4.3 h, respectively.As shown in Fig. 9, after a pressure wave reaches the closed external boundary, transient flow will reach a pseudo-steady state flow, in which the derivative curves converge to a straight line (Chaudhry, 2004;Nie et al., 2011a).Table 5 lists partial data associated with the numerical simulation results and shows quantitative differences in numerical values among wellbore pseudo-pressure, wellbore pseudo-pressure drop, wellbore pressure and wellbore drop for both gas and oil well production under the same conditions.Differences in quantitative analyses can be described from the table : (1) pseudo-pressure values are near 1500 times that of pressure values.Therefore, comparisons between pseudo-pressure and pressure are not significant and comparisons of real wellbore pressure between gas and oil well production need to be calculated; (2) at the beginning of production, pressure drop in oil well production is about 19 times that of gas well production because a rate of 20 m 3 d −1 is small for gas well production and high for oil well pro- duction in a formation with 10 m thickness; (3) at a time of "t = 0.6 h" when the pressure wave of oil well production reaches the external boundary, the pressure drop in the oil well is 7.21 MPa, while the pressure drop in the gas well is 0.346 MPa; (4) at the time of "t = 4.3 h" when the pressure wave of gas well production reaches the external boundary, the pressure drop of the gas well is 0.352 MPa, while the pressure drop in the oil well is 9.57 MPa; and (5) when the wellbore pressure of oil well production decreases to atmospheric pressure (0.1 MPa), which is an absolute open flow status, the production time is 21.5 h which is the limit time of oil well at 20 m 3 d −1 rate production.If production time exceeds limit time, the formation could not support an oil rate of 20 m 3 d −1 and the rate must decline.However, at the limit time of oil well production the pressure drop of gas well production is only 0.358 MPa, which means the formation can easily support a gas rate of 20 m 3 d −1 for an extended period of time.
In conclusion, numerical simulations show that there are obvious differences in wellbore pressure transients for oil and gas flow under the same formation and well conditions owing to the difference of fluid properties.Compared with gas wells, oil wells demand a relatively larger wellbore pressure drop and a relatively faster depletion.Numerical solutions comparisons for oil and gas flows.For convenience, the same formation and well ction conditions were used in the numerical simulations.Pseudo-pressure (see curve ①) and pressure ents (see curves ③) for gas well production, and pressure transients for oil well production (see curves ②) simulated.Pressure depletion of the oil well is greater than that of a gas well under the same conditions.il well production, propagating time of the pressure wave to the closed external boundary (t eo ) is ximately 0.6h.For gas well production, propagating time to the boundary (t eg ) is about 4.3h.Build-up pressure curve of well-testing at wellbore of the example well.The relationship between up pressure and shutting-down time is shown.Wellbore pressure when shut down and at testing ation were 13.79MPa and 15.23MPa, respectively.Shutting-down time of well-testing was 326h.

Field application
This study used a pressure buildup testing well in a sandstone reservoir with an edge water drive.The curve of wellbore shutting-down pressure, p ws , with shutting-down time, t, is shown in Fig. 10.Formation and well parameters are shown in Table 6.The researched nonlinear flow model from this article was used to make a well-test interpretation.Matching curves from a well-testing interpretation using the nonlinear model with a constant pressure boundary are shown in Fig. 11.The conventional linear flow model with a constant pressure boundary was used in a well-test interpretation.Matching curves are shown in Fig. 11.Matching effects between the two interpretations are desirable.Differences between the theoretical curves and real testing data are noticeable (Fig. 11).Using the matching derivative curves, matching effects for the nonlinear model are more desirable than for the linear model, especially in periods II and IV.Matching curve differences between the two models were inconspicuous; therefore, the dimensionless matching curves of well-test interpretation were used to show the differences between the matching effects associated with the two models (Fig. 12).
The main interpretation parameters are shown in Table 7.The formation permeability (k) of well-test interpretation using the nonlinear model is 0.03925 µm 2 and the dimensionless coefficient of nonlinear term (β) of well-test interpretation using the nonlinear model is 0.0287.We substituted the formation permeability value (0.03925 µm 2 ) into the equation β = (1.842×10−3 qBµC ρ )/(kh) and calculated β.The calculated value of β using interpretation permeability = 0.0285, which is slightly smaller than the interpretation β value (0.0287).The differential value and the relative differential value are 0.0003 and 0.70 %, respectively.Therefore, interpretation results using a nonlinear flow model are credible and there is no need to re-match field testing data by changing theoretical parameters.
Table 7 shows the obvious differences between well-test interpretation parameter values for the two models.The differential value of the wellbore storage coefficient is zero because wellbore pressure of the example well is not affected by nonlinearity of oil flow in the wellbore storage regime (Table 7).The values of parameters S, k, and r e using nonlinear model interpretation are smaller than those using nonlinear model interpretations.The relative differential value of S, k and r e are 29.63 %, 18.53 %, and 10.38 %, respectively.
Interpretation using the conventional linear flow model would enlarge the parameter value.The nonlinear flow model is recommended because it is a useful tool for evaluation of formation properties and prediction of engineering conditions.

Conclusions
Nonlinear diffusion equations of liquid and gas in porous media were deduced, nonlinear flow models were established and solved and nonlinear flow behaviour was simulated and analyzed.The findings are as follows:  We assume the permeability in both the horizontal and vertical planes is isotropic and constant: Dimensionless pressure:p D = kh (p i − p) / 1.842 × 10 −3 qBµ .Skin factor: S = kh p s / 1.842 × 10 −3 qBµ , p s is additional pressure drop near wellbore.Dimensionless radius based on effective well radius:r D = r/ r w e −S .Dimensionless radius of external boundary: r eD = r e / r w e −S .
Type curves of pressure transientsType curves reflect properties of underground formations.Type curves graphically show the process and characteristics of fluid flow in reservoirs.
iv. Regime IV, external boundary response regime.For a constant pressure boundary, where the pressure derivative curve decreases, transient flow ultimately becomes steady state.For closed boundary, as the pressure derivative curve increases, transient flow ultimately becomes pseudo-steady state, where the type curves converge to a straight line with unit slope.

Figure 1
Figure 1 shows type curve characteristics as controlled by different values of parameter group C D e 2S .A larger C D e 2S leads to a higher location of the dimensionless pressure curve.(2) Type curves of nonlinear flow model Figures 2-5 show the typical nonlinear flow characteristics and flow processes of vertical well production in a homogenous formation.Figure 2 contains the type curves of a nonlinear flow model with an infinite boundary.The type curves are obviously controlled by the dimensionless coefficient of nonlinear term β.Pressure transients were simulated by setting β = 0, 0.01, 0.05 and 0.1.A larger β leads to smaller dimensionless pressure curves and associated derivative curves.The "β = 0" curves are associated with the linear flow model.Three main flow regimes can be easily discerned from the type curves:i.Regime I, pure wellbore storage regime, there are no differences in type curves between the two flow models because liquid in formation does not start to flow and the influence of the nonlinear quadratic pressure gradient term is only produced for flow in formation.Wellbore pressure transients are not affected by the nonlinearity of oil flow in this regime.ii.Regime II, wellbore storage and skin effect regime, there are obvious differences in type curves between the two flow models.The nonlinearity of liquid flow positively influences the pressure transients.A larger β means a stronger nonlinear effect of liquid flow on the type curves.

Fig. 1 .
Fig. 1.Type curves for linear flow models used for comparison with those from nonlinear flow models.Type curve characteristics controlled by parameter group C D e 2S and external boundary conditions are shown.Four main flow regimes were recognised: pure wellbore storage regime (I), wellbore storage and skin effect regime (II), radial flow regime (III), and external boundary response regime (IV).Derivative curves converged to the "0.5" line in the radial flow regime for infinite formation, decreased for constant pressure boundary and increased for closed boundary.Pressure curves and their derivative curves ultimately converged to a straight line with a slope of one.

Fig. 2 .
Fig. 2. Type curves for the nonlinear well-test model with an infinite boundary used in the analysis of nonlinear flow processes.The type curve characteristics are controlled by different values of the dimensionless coefficient of the nonlinear term β.Pressure transients were simulated by setting β = 0, 0.01, 0.05, and 0.1.Pressure curves and their derivative curves were positively affected by β.Derivative curves did not converge to the "0.5" line in the radial flow regime.The "β = 0" curves were linear flow model curves.

Fig. 3 .
Fig.3 Type curves for the nonlinear well-test model with a constant pressure boundary used for analysis o 687

Figure 3
Figure 3 contains the nonlinear flow model type curves with constant pressure boundaries.Pressure transients were simulated by setting β = 0, 0.02, 0.1, and 0.2.Nonlinear flow characteristics in regimes I-III are similar to the nonlinear flow model type curves with an infinite boundary.The type curves in regime IV show that the pressure derivative curves decrease and ultimately converge at a point.Figure 4 contains the nonlinear flow model type curves with closed boundaries.Pressure transients were simulated by setting β = 0, 0.02, 0.1, and 0.2.Nonlinear flow characteristics in Regimes I-III are similar to infinite boundary type curves.The type curves in Regime IV show that the pressure curves and their derivative curves increase and ultimately converge to a straight line whose slope is smaller than one, instead of one, which is different from the linear flow model.Figure 5 shows the type curve characteristics that are controlled by values associated with the parameter group C D e 2S .Pressure transients were simulated by setting C D e 2S = 10 5 , 10 3 and 10.A larger C D e 2S leads to larger dimensionless pressure curves, which is similar to the linear flow model.Dimensionless pressure curve characteristics are completely different from those of the linear flow model.Derivative curves cross at a point (see point A on the caption of Fig. 5) before the radial flow regime.A larger C D e 2S leads to larger derivative curves in Regime II and smaller derivative curves in Regime III.4.1.2Quantitative analysis of nonlinear influencePressure curves and the derivative curves of nonlinear flow models deviate gradually from those of linear flow models with time."DV" and "RDV" show the quantitative

Fig. 4 .
Fig. 4. Type curves of nonlinear well-test model with a closed boundary used for the analysis of nonlinear flow processes.The figure shows the type curve characteristics controlled by different values of β.Pressure transients were simulated by setting β = 0, 0.02, 0.1 and 0.2.Pressure curves and the associated derivative curves increase and ultimately converged to a straight line whose slope is smaller than one.

Fig. 6 Fig. 5 .
Fig.5 Type curves of the nonlinear well-test model with an infinite boundary used for analysis of nonlinear 699 curves of the nonlinear well-test model with an infinite boundary used for analysis of nonlinear ow processes.The type curve characteristics are controlled by different values of the parameter group C D e 2S .ressure transients were simulated by setting C D e 2S = 10 5 , 10 3 and 10.Derivative curves crossed at a point (see int A on the caption) before the radial flow regime.ig.6 Pressure curve comparisons with Chakrabarty et al. (1993b) for homogenous infinite formation.Full nes are pressure curves from this article and dotted lines are pressure curves associated with Chakrabarty et al. 993b).β is the opposite number of α.For convenience, the same parameter values are set.The solutions of this ticle are the same as those of Chakrabarty et al. (1993b).

Fig. 6 .
Fig.6.Pressure curve comparisons withChakrabarty et al. (1993b) for homogenous infinite formation.Full lines are pressure curves from this article and dotted lines are pressure curves associated withChakrabarty et al. (1993b).β is the opposite number of α.For convenience, the same parameter values are set.The solutions of this article are the same as those ofChakrabarty et al. (1993b).

Fig. 10 .
Fig. 10.Build-up pressure curve of well testing at wellbore of the example well.The relationship between build-up pressure and shutting-down time is shown.Wellbore pressure when shut down and at testing termination were 13.79 MPa and 15.23 MPa, respectively.Shutting-down time of well testing was 326 h.

738Fig. 11 Fig. 11 .
Fig.11 Matching curves with units of well-test interpretation for the example well.Log-log curve 739

Fig. 12 .
Fig.11 Matching curves with units of well-test interpretation for the example well.Log-log curve 739

Table 1 .
studied a nonlinear pressuretransient model which does not consider skin factor and plotted a group of log-log pressure-transient curves based on t D .Theoretical offset of type curves caused by the nonlinear term (β = 0.01).: β is the dimensionless coefficient of the nonlinear term; t D is dimensionless time; C D is dimensionless wellbore storage coefficient; p wD is dimensionless pressure; p wD is dimensionless pressure derivative; DV is differential value; RDV is relative differential value.The quantitative analysis of nonlinear influence was calculated by setting C D e 2S= 10 3 and β = 0.01.Corresponding type curves are shown in Fig.2.

Table 2 .
Theoretical offset of type curves caused by the nonlinear term (β = 0.1).
Explanations: β is the dimensionless coefficient of the nonlinear term; t D is dimensionless time; C D is dimensionless wellbore storage coefficient; p wD is dimensionless pressure; p wD is dimensionless pressure derivative; DV is differential value; RDV is relative differential value.The quantitative analysis of nonlinear influence was calculated by setting C D e 2S= 10 3 and β = 0.1.Type curves are shown in Fig.2.

Table 3 .
Probable values for the dimensionless coefficient of the nonlinear term.

Table 4 .
Basic data from the simulated well.
Numerical solutions comparisons for oil and gas flows.For convenience, the same formation and well production conditions were used in the numerical simulations.Pseudo-pressure (see curve 1 ) and pressure transients (see curves 3 ) for gas well production, and pressure transients for oil well production (see curves 2 ) were simulated.Pressure depletion of the oil well is greater than that of a gas well under the same conditions.For oil well production, propagating time of the pressure wave to the closed external boundary (t eo ) is approximately 0.6 h.For gas well production, propagating time to the boundary (t eg ) is about 4.3 h.

Table 5 .
Numerical comparisons between oil and gas well production.: t is production time; ψ w is gas pseudo-pressure of wellbore, MPa 2 (mPa s) −1 ; ψ is gas pseudo-pressure drop of wellbore, MPa 2 (mPa s) −1 ; p w is wellbore pressure; p is wellbore pressure drop. Explanations

Table 6 .
Formation and well parameters of the example well.

Table 7 .
Main interpretation parameters of the example well.C s is wellbore storage coefficient; S is skin factor; k is formation permeability; r e is distance of constant pressure boundary from wellbore; β is dimensionless coefficient of nonlinear term; DV is differential value; RDV is relative differential value.1.Effects of nonlinearity upon pressure transients are obvious and nonlinear models more accurately portrait the flow processes of fluid in porous media.2.Locations of pressure and nonlinear model derivative curves for liquid flow are lower than those derived from linear models.3.Differences between nonlinear and linear model pressure transients increase with time and the nonlinear coefficient.4.Influences of nonlinearity are greater for formations with low permeability, heavy oil or thin thickness.5.Nonlinear transient flow behaviour of gas is different from that of oil because of dissimilar fluid properties. is liquid compressibility, MPa −1 ; ρ 0 , ϕ 0 , p 0 are reference values, which are usually used under standard conditions.Substitute Eq. (A6) into Eq.(A4):Ct= C ρ + C f .(A14)whereC t is total compressibility of rock and liquid, MPa −1 .Substitute Eqs.(A8)-(A10) and Eq.(A13) into Eq.(A.3): Explanations: www.nonlin-processes-geophys.net/20/311/2013/ Nonlin.Processes Geophys., 20, 311-327, 2013 Equation (A19) is the governing differential equation containing the quadratic gradient term in Cartesian coordinates.