Well Productivity Index for Compressible Fluids and Gases

In this paper we discuss the notion of the diffusive capacity for the generalized Forchheimer flow of fluid through porous media. The diffusive capacity is an integral characteristic of the flow motivated by the engineering notion of the productivity index (PI), Dake 1983, Raghavan 1993, Christopher et al. 2014. The PI characterizes the well capacity with respect to drainage area of the well and in general is time dependent. We study its time dynamics for two types of fluids: slightly compressible and strongly compressible fluid (ideal gas). In case of the slightly compressible fluid the PI stabilizes in time to the specific value, determined by the so-called pseudo steady state solution, Aulisa et al. 2009, 2011, 2012. Here we generalize our results from Aulisa et al. 2012 on long term dynamics of the PI in case of arbitrary order of the nonlinearity of the flow. In this paper we study the mathematical model of the PI for compressible gas flow for the first time. In contrast to slightly compressible fluid this functional mathematically speaking is not time-invariant. At the same time it stays"almost"constant for a long period of time, but then it rapidly blows up as time approaches the certain critical value. This value depends on the initial data (initial reserves) of the reservoir. The"greater"are the initial reserves, the larger is this critical value. We present numerical and analytical results for the time asymptotic of the PI and its stability with respect to the initial data. Using comparison theorems for porous media equation from V\'azquez 2007 we obtain estimates between the PI's for the original gas flow and auxiliary flow with a distributed source. The latter one generates the time independent PI, and can be calculated using formula similar to one in case of slightly compressible fluid.


(Communicated by Josef Málek)
Abstract. We discuss the notion of the well productivity index (PI) for the generalized Forchheimer flow of fluid through porous media. The PI characterizes the well capacity with respect to drainage area of the well and in general is time dependent. In case of the slightly compressible fluid the PI stabilizes in time to the specific value, determined by the so-called pseudo steady state solution, [5,3,4]. Here we generalize our results from [4] in case of arbitrary order of the nonlinearity of the flow. In case of the compressible gas flow the mathematical model of the PI is studied for the first time. In contrast to slightly compressible fluid the PI stays "almost" constant for a long period of time, but then it blows up as time approaches the certain critical value. This value depends on the initial data (initial reserves) of the reservoir. The "greater" are the initial reserves, the larger is this critical value. We present numerical and theoretical results for the time asymptotic of the PI and its stability with respect to the initial data.
1. Historical remarks and review of the results. The classical equation describing the fluid flow in porous media is the Darcy's law, stating the linear relation between the pressure gradient ∇p and the velocity u. Darcy himself observed in [11] that the area of applicability of linear relation is very limited. When, for instance, the fluid has high velocity or in the presence of fractures in the media, the nonlinear models are necessary to capture the properties of the flow.
One of the widely-used nonlinear models is the Forchheimer equation in the form g(|u|)u = −∇p, where g(s) is a polynomial, [7,20]. Originally Forchheimer in his work [13] proposed three particular equations to match the experimental data: two term, three term and power laws, with g(s) being up to second order degree polynomial. To embrace the recent findings on the nonlinearity of the fluid flow (see [18,24]) and simplify the mathematical handling, it is convenient to consider the generalization of classical Forchheimer equations to the case where g(s) is the generalized polynomial with non-negative coefficients and, possibly, non-integer powers (see [2]). We call this family of equations the g-Forchheimer equations. The g-Forchheimer equation combined with the equation of state and the conservation of mass results in a single degenerate parabolic equation for the pressure p(x, t) only. Note, that in case of natural reservoirs the inertial and viscous terms in Brinkman-Forchheimer equation (e.g. [23]) have very small impact on engineering parameters such as the productivity index, and these terms are usually neglected, [20,10].
In this paper we discuss two types of fluids: slightly compressible fluid, characterized by the equation of state ρ(p) ∼ exp(γp) with very small compressibility constant (γ ∼ 10 −8 ), and strongly compressible fluid (in particular, ideal gas) characterized by the equation of state ρ(p) ∼ p, [20]. While in general porosity depends on pressure, see [21,26,8], we only consider it to be a function of spatial variable x. In case of slightly compressible fluid we studied different properties of g-Forchheimer equations in [2,4,14]. In this paper we extend the results of our work [4] on asymptotic behavior of the pressure function to the case when the degree of the g-polynomial is arbitrary. In case of ideal gas we discuss both numerical and analytical results for the time dynamics of the solution of the corresponding parabolic equation.
Keeping in mind that the applications of our findings are in geophysics and in particular in reservoir engineering, we restrict out studies to fluid flow in a reservoir with bounded domain U . The reservoir is bounded by the exterior impermeable boundary Γ e and the interior well-boundary Γ i . Γ i is subject to various boundary conditions depending on the flow regime. We introduce the capacity type functional to study the asymptotic behavior of the fluid flow in the reservoir with respect to time t. This functional is motivated by the well Productivity Index (PI) and is often used by reservoir engineers to measure well capacity, see [10,25,15]. Productivity index is the total amount of fluid per unit pressure drawdown (the difference between reservoir and well pressures) that can be extracted by the well from a reservoir (see [25]) and as such describes the ability of the reservoir to deliver fluids to the well. It is defined as J(t) = Q(t)/P DD(t), where Q(t) is the total flux through the boundary Γ i , and the pressure drawdown P DD(t) is equal to the difference between averages of the pressure in the domain and on the boundary Γ i . Note, that in our previous articles [3,4,5] instead of PI we used the notion of the diffusive capacity that has more general physical meaning for other evolutional problems. In current article we use the term PI that is identical to the diffusive capacity and has clear engineering interpretation.
In general the PI is time dependent for both slightly and strongly compressible fluid flows, [10,25]. However its time dynamics differs greatly in these two cases.
For slightly compressible fluids and specific regimes of production the PI stabilizes in time to a constant value. This value can be determined using the solution of a particular boundary value problem. Namely, the time-invariant PI is associated with the pressure distribution p s (x, t) = − Qs |U | t + W (x). Here the initial pressure distribution W (x) is the solution of the specific steady state BVP, Q s is a constant flux through the well-boundary, and |U | is the volume of the reservoir domain. Such pressure is called the pseudo-steady state (PSS) pressure and satisfies the split boundary condition on the well − Qs |U | t + ϕ(x) for some known function ϕ(x).
For arbitrary initial data one can not expect the PSS pressure distribution and constant PI. However, as it appears from the engineering practice (see [10,25]) the constant PSS PI serves as an attractor for the transient PI. We proved this fact in [3,4] under a number of conditions. In this paper we will generalize our previous results. We consider two types of boundary value problems.
In the IBVP-I we impose total flux Q(t) on the well-boundary and consider the trace of the pressure function to be split as γ(t) + ψ(x, t), with Γi ψ ds = 0. In this case γ(t) can be considered as the average of the trace function and ψ as the deviation of the actual trace from its average. Note that we impose conditions only on the function ψ, while γ is unknown. For t → ∞ the boundary data {ψ(x, t), Q(t)} are assumed to be localized in the neighborhood of the time independent {ϕ(x), Q s }.
In the IBVP-II no assumptions on the well-boundary flux are made. Instead, we specify the Dirichlet condition γ(t) + ψ(x, t) on the well boundary. In this case both γ and ψ are known functions. For t → ∞ the boundary data {ψ(x, t), γ(t)} are assumed to be localized in the neighborhood of functions {ϕ(x), − Qs |U | t}. The main result for both IBVP-I and II is that the time dependent PI asymptotically stabilizes in the neighborhood of the PI for the PSS regime associated with the pair ϕ(x) and Q s (see Theorems 4.11 and 5.4). The corresponding value for the steady state PI can be calculated just by solving an auxiliary Dirichlet time independent BVP.
In [4] we obtained this result under the degree condition, stating that the degree of g-polynomial deg(g) ≤ 4 n−2 , where n is the dimension of space. Mathematically this condition arises from the theory of Sobolev spaces and insures the continuous embedding W 1,2−a (U ) ⊂ L 2 (U ). While for n = 2 this constraint holds for any degree deg(g), for n = 3 it holds only for deg(g) ≤ 4. In this paper we study the asymptotic behavior of the PI without any constraint on deg(g). In this case the classical Poincaré-Sobolev inequality does not hold, and we use weighted inequality for the mixed term |∇p| 2−a |p| α−2 , α = 2. We obtain estimates for the bounds of the L α -norm for both the difference between transient and PSS solutions, and the difference between transient and PSS solution time derivatives.
In Sec. 6 we discuss the concept of the PI for an ideal gas flow and some results on its time dynamics. In this case the productivity index is defined, [6,9], as J g (t) = Q(t)/P DD(t), where P DD(t) is the difference between the average of p 2 in the domain U and on the boundary Γ i . In contrast to the case of slightly compressible flow, we show numerically that until certain critical time T crit the transient PI remains almost constant and then as time approaches T crit it blows up. This result obtained on actual field data corresponds with the engineering observations. Time T crit depends on the initial reserves/initial data: the "greater" the initial reserves are the greater T crit is. Similar to our approach in case of slightly compressible fluid, we use the auxiliary pressure p 0 (x, t), resulting in time independent PI to investigate the behavior of time dependent PI. The p 0 (x, t) is the solution of the equation with positive function on the RHS. This function can be considered as fluid injection inside the reservoir. When gas reserves are considerably larger than pressure drawdown on the well, this source term is negligible, and the PI for p 0 (x, t) is almost identical with the general time dependent PI.
The paper is organized as follows.
• In Sec. 2 we discuss the various aspects of nonlinear Forchheimer equations and associated parabolic PDE. We give the definition of the PI, (10), on the solution of this parabolic equation.
• The IBVP-I for total flux and IBVP-II for Dirichlet boundary conditions 6 we introduce the concept of PI (153) for an ideal gas flow, based on the so-called pseudo-pressure function. We show numerically that the gas PI stays the same until it sharply blows up when time approaches the critical value T crit dependent on the parameter characterizing the initial reserves. This fact corresponds well with the field data. We associate this constant PI with the auxiliary pressure in (154) and prove the estimate for the difference between the actual and auxiliary pressure in terms of PI, Sec. 6.1. Finally, in Theorem 6.10 we obtain the stability of the PI with respect to the parameter characterizing the initial gas reserves.
2. Problem statement and preliminary properties. Consider a fluid in a porous medium occupying a bounded domain U ⊂ R n . Let x ∈ R n , n ≥ 2, be a spatial variable and t ∈ R be a time variable. Let u(x, t) ∈ R n be the fluid velocity and p(x, t) ∈ R be the pressure.
We consider a generalized Forchheimer equation where g : R + → R + . In particular we consider function g to be the generalized polynomial with non-negative coefficients. Namely g(s) = a 0 s α0 + a 1 s α1 + · · · + a k s α k = a 0 + k j=1 a j s αj , with k ≥ 0, the real exponents satisfy α 0 = 0 < α 1 < α 2 < . . . < α k , and the coefficients a 0 , a 1 , . . . , a k > 0. The largest exponent α k is the degree of g and is denoted by deg(g). One can notice that the equation (1) with g(s) defined as in (2) includes the linear Darcy's equation and all classical forms of Forchheimer equations [13] (for details see our previous works [2,4]).
In case of slightly compressible fluid we consider the case when the Degree Condition in formula (2.14), [4], is not satisfied, namely In this case there is no continuous embedding W 1,2−a (U ) ⊂ L 2 (U ) and corresponding Poincaré inequality does not hold. Clearly, if n = 3 condition (3) will hold for the nonlinearities α k > 4.
From (1) one can obtain the non-linear Darcy equation explicitly solved for the velocity u: where Along with (1), which is considered as a momentum equation, the dynamical system is subject to the continuity equation where ρ(x, t) ∈ R + is density of the fluid. The fluid is considered to be slightly compressible, subject to the equation of state where γ ∼ 10 −8 is the compressibility constant, see [20]. Combining (1), (6) and (7) we obtain the pressure equation (see [2] for details): Since in natural reservoirs γ ∼ 10 −8 , in engineering practice the second term is always neglected in comparison to the first one, [25,10]. We thus arrive to the truncated degenerate parabolic equation It is not difficult to show that for small γ the difference between the solutions of full and truncated equations is also small. By scaling the time variable we will assume further on that γ = 1. In case of the flow of ideal gas the equation (1) for n = 2 will take the form (see, for example, [22]) where β is the Forchheimer coefficient. As in case of slightly compressible fluid, the above equation can be solved for u: Combining the equation above with the equation of state ρ(p) = M p (without loss of generality we take M = 1) and continuity equation (6) we obtain (for the details see Sec. 6) We will study equations (8) and (9) in the open domain U ⊂ R n with the C 2 boundary ∂U = Γ = Γ e ∪ Γ i . The Γ e is considered to be external impermeable boundary of the reservoir with u · N = 0, where N is the outward normal vector to Γ e . The Γ i is the internal boundary of the well with the total flux or Dirichlet boundary conditions. To study the long term dynamics of the g-Forchheimer flow in the domain U we introduce special capacity type functional. We call it the Productivity Index (PI) (see [10,25,28]) and define as where Q(t) is the total flux through the well-boundary Γ i and P DD(t) is called the pressure drawdown in the domain U . In case of slightly compressible fluid these quantities are defined as are the pressure averages in the domain and on the well-boundary correspondingly. For gas filtration in porous media these quantities are defined as, see [6,9]: We end this section with the recollection of some properties of K(ξ) that will be used in this study (see [2,4]). If g(·) is defined as in (2), then the inverse function K(ξ), ξ ∈ [0, ∞) in (5) is decreasing and satisfies Moreover the function K(|y|)y, y ∈ R n , is monotone, i.e. for any two functions where Following [2,4] we define the functional Function H(ξ) can be compared with K(ξ) (see [2]) as follows: Combining (19) with (15) there exist constants C 1 and C 2 such that 3. Two types of boundary conditions and pseudo-steady state solution.
In Sections 3, 4 and 5 we consider the case of slightly compressible fluid. Let U ⊂ R n be an open set with the C 2 boundary ∂U = Γ = Γ e ∪ Γ i , see Fig. 1. The Γ e is considered to be external impermeable boundary of the reservoir and the Γ i is the internal boundary of the reservoir (well surface). With respect to boundary Γ i two different problems will be considered: IBVP-I with imposed total flux (generalized Neumann condition) and IBVP-II with imposed Dirichlet boundary condition.
Throughout this paper, we consider solutions with sufficient regularities both in x and t variables such that our calculations can be performed legitimately. Though the results are applicable to to weak solutions with enough regularities up to the boundary, we consider classical solution, namely p ∈ C 2,1 x,t (Γ). Since the IBVP-I (21)-(24) lacks uniqueness due to conditions (22) and (23), following [4] we further restrict the boundary data. In particular, let the trace of p(x, t) on well-boundary Γ i be of the form with Γi ψ(x, t) ds = 0.
Note, that function γ(t) is not specified and will be determined by the flux Q(t).
Such restriction of the boundary trace ensures the uniqueness of the solution of IBVP (21)-(24) (see Remark 3.2 in [4]). Transient boundary data ψ(x, t) and Q(t) will be compared with the timeindependent boundary data ϕ(x) and Q s . As in [4] we consider functions Ψ(x, t) and Φ(x) to be defined in the domain U with the traces on Γ i being ψ(x, t) and ϕ(x) correspondingly. The results are obtained under the constraints on the parameters of difference of boundary data II. IBVP-II for Dirichlet boundary condition. The function p(x, t) is a solution of the IBVP-II if it satisfies (21), (23), (24) with the condition on the well boundary for known function ψ(x, t). Transient boundary data ψ(x, t) and γ(t) will be compared with the boundary data ϕ(x) and −At, A = Q s /|U |, where Q s is constant flux on Γ i . We again consider functions Ψ(x, t) and Φ(x) to be the W 1 2 (U ) extensions of ψ(x, t) and ϕ(x) on the domain U . The results are obtained under the constraints on the differences For the existence and regularity theory of degenerate parabolic equations, see e.g. [16,19,14] and references therein.

PSS solution.
For the time-independent boundary data ϕ(x), Q s we define the pseudo steady state (PSS) solution of IBVP for equation (21) for a given function ϕ(x) and a constant total flux Q s on the boundary Γ i . Function W (x) is called the basic profile corresponding to the flux Q s and is defined as a solution of BVP The existence and uniqueness of a solution for the given pair ϕ(x), Q s is proved in [3]. It is not difficult to see that in case of the PSS solution when p(x, t) = p s (x, t) the functional J(t) in (10) is time independent and We assume that ϕ(x) and Q s are such that the denominator of (34) is not equal zero. Obviously PSS solution is the only solution for given ϕ(x), Q s which obey both types of boundary conditions. Therefore PSS will serve us a reference solution for comparison for both IBVP-I and II.
The proof uses the following Lemma which reduces the question of the convergence of J(t) to the convergence of the gradient of the solution in the appropriate norm and of the total fluxes to a corresponding constant values. Lemma 3.2. As t → ∞ the difference between transient and PSS PI's Proof. From (10) and (34) it follows is the difference of pressure drawdowns for transient and PSS regimes. The pressure drawdown for the PSS solution is time independent, and is not equal zero. Hence The difference in pressure drawdowns ∆ p (t) can be written as Applying Trace theorem and Poincaré inequality, as U z dx = 0, one has The condition (35) follows.
The results for the solution of IBVP-I will be obtained under the conditions on the deviation of the boundary data (27) in Sec. 4, and for the solution of IBVP-II in terms of deviations (29) in Sec. 5.

4.1.
Assumptions on boundary data. In order to formulate the assumptions in the Theorem 3.1 we first introduce the following notations. Suppose (40) The certain way the conditions A1-A5 will be formulated is dictated by the fact, that while the main result will be obtained under all the assumptions below, some auxiliary estimates require less restrictive constraints. As a note for each assumption we state the resulting bounds on parameters F 1 , F 2 , A 1 , A 2 and A 3 .
Along with the solution of IBVP-I p(x, t) and the PSS solution p s (x, t) of IBVP-I with boundary data Φ(x), Q s , we will use the following shifts of the solutions: It follows that q satisfies the BVP Similarly q s satisfies the BVP Note that According to Lemma 3.2 we need the convergence ∆ Q (t) → 0 and ∇(p − p s ) L 2−a → 0 as t → ∞. The affinity of Q(t) and Q s will be imposed as a condition on the boundary data (see Assumptions A4). The estimates necessary to prove the convergence of ∇(p−p s ) L 2−a are outlined in the following Lemma. The references to the corresponding results further in the paper are given. 2.
Proof. According to the inequality (5.20) in proof of Theorem 5.6 in [4] we have: Here C Φ = C Φ (p, W ) is as in (17) and is bounded by virtue of (57) and (59). It is then clear that the RHS of (60) converges to 0 under conditions (55)-(59).

4.2.
Bounds for the solutions. As before let q(x, t) be the solution of BVP (50), (51) and q s (x, t) be the solution of PSS BVP (52), (53). Let z = q − q s . We will prove that under certain condition U |z| α dx is bounded at time infinity.
In order to obtain the suitable differential inequality we will use the following result (see Lemma We first prove the following lemma to estimate the boundary data. Proof. We have z| Γi = (q − q s ) Γi = ∆ γ (t) + F 1 (t). Then by Trace theorem Then

WELL PRODUCTIVITY INDEX FOR COMPRESSIBLE FLUIDS AND GASES 13
Applying Young's inequality with ε we get Proof. Subtracting equations (50) and (52) from each other, multiplying on |z| α−1 sign(z) and integrating over U one has From boundary conditions (51) and (53) it follows that z| Γi = ∆ γ (t) + F 1 (t), and one has 1 α We split the first integral in RHS of (66) in four separate integrals and estimate them one-by-one. Namely, where (notice that ∇p s = ∇W ) Since ∇z = ∇(p − W ) − ∇(∆ Ψ ) and α−2 α < α−a α = 1 γ0 , similar to (4.5) in [14] we have The integral I pw can be estimated similar to (4.6)-(4.9) in [14]: The first integral in (69) can be estimated similar to (4.8) The second integral in (69) is estimated similar to (4.7) Similar the third integral in (69) can be estimated as Since 1 2−a + 1−a 2−a = 1, from Hölder and Young's inequalities it follows Combining the estimates above in (69) we have the estimate for the integral I pw Finally, for the integral I wp we have Applying Hölder's and Young's inequalities for the first integral in the RHS of (75) we get The second integral in the RHS of (75) can be estimated similar to (72) Then in (75) the integral I wp can be estimated as Substituting (68), (74) and (76) in (67) we get the estimate for the first integral in the RHS of (66) The second integral in the RHS of (66) can be estimated as Then similar to the estimate for I pw and (72) we get . Finally combining (77), (79), (80) and Lemma 4.3 in (66) we get after choosing small enough ε Applying (61) to the first term in the RHS of (81) we get .
In order to prove the estimate on U |z| α dx we recall the Lemma A.1 from [14]. The following result will follow from Lemma 4.4.

4.3.
Bounds for the gradient of solutions. We will obtain the bounds for the integral U H(|∇p|) dx, where H(|∇p|) is defined in (18). This result is necessary for our estimates of the difference in time derivative of fully transient and PSS solutions.
Consequently, if Q(t) and Ψ(x, t) satisfy conditions A2 then Proof. In formula (4.46) in [4] we have for z = q − q s : where Notice that Subtracting (88) from (86) we get d dt The term U z 2 dx can be estimated using (82) as α > 2. Then in (89) we have Applying Gronwall's inequality one has for t ≥ 0 According to formula (2.35) in [14] lim sup Then using estimates (92) and (82) For the first integral in the RHS of (94) we have According to Lemma 4.4 in [4] with F 2 (t) as in (38) we have Choosing ε 1 and ε 2 to be small enough and using (94) in (93) we obtain (84). Since F 2 (t) is uniformly bounded under assumptions A2, the estimate (85) follows.

Estimate of time derivatives.
We will obtain the estimate on the difference of time derivative of fully transient and PSS solutions.
The quantity in the first integral in the RHS of (100) is equal to From (13) it follows that Then in (100) we have To estimate the second term in the RHS of (101) we apply Hölder inequality and Young's inequality with ε = (1 + a)/(1 − a). Since K(|∇p|) is bounded (see (14)) applying once more Hölder inequality with powers β/(β − 2) and β/2 we get Similar for the third term in the RHS of (101) one has Then in (101) we have Since the first term in the RHS of (104) is negative, it can be neglected. We then have where The result follows from the Gronwall's inequality under assumptions A3-β.
In order to prove that lim sup t→∞ U |p t + A| 2 dx = 0 we will need the Weighted Poincaré inequality (see Lemma 2.5., [14]): Lemma 4.9. Let ξ = ξ(x) ≥ 0 and function u(x) be defined on U and is vanishing on the boundary ∂U . Assume α ≥ 2. Given two numbers θ and θ 1 such that θ > 2 (2 − a) * and max 1, where Proof. In (5.15) of [4] we have where the constant C depends on U H(|∇p|) dx and is bounded under assumptions A2.
Taking in above ε 2 → 0 we complete proof of the Theorem 4.10.
Combining it with Lemma 3.2 we get the main result for the solutions of IBVP-I: Theorem 4.11. If boundary data Q(t) and Ψ(x, t) satisfy assumptions A5 then 5. IBVP-II: Asymptotic convergence of the productivity index. Let p(x, t) be the solution of IBVP-II (21), (23), (24) with Dirichlet boundary condition (28) with boundary data Ψ(x, t). Let J(t) be the corresponding PI. Let p s (x, t) = −At + W (x), W (x)| Γi = ϕ(x), be the PSS solution of this IBVP defined by the boundary data −At, A = Q s /|U |, and Φ(x). Let J P SS be the corresponding PI.

Assumptions on the boundary data.
Let for any β ≥ 2.

EUGENIO AULISA, LIDIA BLOSHANSKAYA AND AKIF IBRAGIMOV
The results in this section are obtained under the following conditions on the boundary data (121) According to Lemma 3.2 in order to prove J(t) → J P SS as t → 0 we need to prove the following estimates:

5.2.
Estimates on the solution of IBVP-II. Along with the solutions p and p s we will use the shifts p(x, t) = p(x, t) − Ψ(x, t) and p s (x, t) = p s (x, t) − Φ(x) + At. These shifts obviously satisfy Lemma 5.1. Suppose boundary data Ψ(x, t) satisfies assumptions D1, then Proof. Similar to (100) we have Similar to (102) -(105) we have The result follows from the Gronwall's inequality under assumptions D1.
Further we will need the following result obtained in [14]: under assumptions D2.1 (122) (see [14], Th. 4.5) Lemma 5.2. Suppose boundary data Ψ(x, t) satisfies assumptions D1 and D2. Then Proof. First, notice that Differentiating both sides of (128) in t, multiplying by p t and integrating over U we have due to zero boundary condition Similar to estimate (5.10) in [4] we get since ∇ p = ∇p − ∇Ψ Applying the Weighted Poincaré inequality (106) with u = p t , p t | Γi = 0 and ξ = ∇p with powers α = 2 and θ = 2 2−a to the first integral on the RHS of (137) we get: Here Under assumptions D2.2 the RHS of (140) converges to 0 as t → ∞ and the result of Lemma follows from (135).

24
EUGENIO AULISA, LIDIA BLOSHANSKAYA AND AKIF IBRAGIMOV Further we will use another result from [14]: under assumptions D3.1 (124) (see [14], Th. 4.3) Theorem 5.3. Suppose the boundary data Ψ(x, t) and Φ(x) satisfies assumptions D1 − D3. Then Proof. From (16) we get (133) and (59) Since p − p s | Γi = 0, integration by parts in the RHS of (143) gives The last inequality holds in view of (133) and (59). Thus 6. Productivity index concept for the flow of ideal gas. In this section we discuss the concept of the PI for an ideal gas flow in the porous media for a wellreservoir system. We will repeat some discussion from Sec. 2 and provide more details. The equation of state for an ideal gas takes the form (see [1], [20]) where ρ is the density of the fluid and M is a proportionality constant. Without loss of generality we assume M = 1. As a momentum equation we consider the second order Forchheimer equation (1) which takes the form (see, for example, [22]) where β is the Forchheimer coefficient. As in case of slightly compressible fluid, the above equation can be solved for u and rewritten in terms of a nonlinear permeability function K 2 (·): Combining (145) and (147) together with the continuity equation (6) yields the parabolic equation for pressure As before the domain U models the reservoir with boundary split in two parts Γ i and Γ e . Γ i is the well-boundary while Γ e is an exterior non-permeable boundary of the reservoir. We impose Dirichlet Pseudo Steady State (PSS) boundary condition on Γ i zero mass flux boundary condition on Γ e ρu · N | Γe = 0 (150) and the initial condition for φ 0 (x) ≥ 0.
Remark 6.1. Constant B is a positive parameter (generally large) characterizing the initial reserves of gas in the reservoir domain U . The constant A is associated with the amount of gas extracted at the well-bore Γ i . Thus the quantity B − At quantifies the gas reserves in the reservoir at the moment t.
For the flow of an ideal gas in porous media, engineers define the Productivity Index as, see [6,9]: where N is an outward normal to Γ i . Then productivity index for the well-reservoir system is defined by As we saw earlier in case of slightly compressible flow, under the boundary conditions (149) and (150) the value of the PI, defined as in (10) is stabilizing to a constant value (34), which is defined by the PSS solution and depends only on the parameter A. Obviously, in case of the gas flow described by (148) this feature is not valid anymore and the PI is not stabilizing in time. Nevertheless the approach developed in previous sections is applicable to the gas flow as well. Namely we will introduce special auxiliary pressure function characterized by the time independent PI. We will then study the relation between the general time dependent PI and this time independent one.
Let Q 0 = A |U | be a given constant mass-rate of gas production and consider the auxiliary pressure distribution given by where W (x) is the solution of the following BVP Integrating (155) over the domain U , using integration by parts and boundary condition (157) yields to the integral flux condition From (154) we can also derive that Combining this property, definition (154) and BVP (155)-(157) it can be easily verified that the auxiliary pressure p 0 satisfies the following equations where From (158) and (159), it also follows the total mass flux condition for p 0 Using this last result and the explicit representation of p 0 (x, t) in formula (153) leads to the following Proposition. Proposition 6.3. The Productivity Index for an ideal gas flow defined on p 0 (x, t) is time independent and is given by Numerical computations performed for several basic reservoir geometries, show that if the initial data in the system (147)-(150) is given by p(x, 0) = p 0 (x, 0), then the corresponding productivity indices J g [p 0 ] and J g [p](t) are almost identical for a long time, see Fig. 2, as long as the quantity or, equivalently, as long as In case of compressible flow, we have to fix a critical time

28
EUGENIO AULISA, LIDIA BLOSHANSKAYA AND AKIF IBRAGIMOV and for t > T crit the negative boundary data (149) leads to the violation of ellipticity. This behavior is qualitatively justified by the fact that as long as (166) holds, function f 0 (x, t) in the RHS of (160) is negligible and p(x, t) and p 0 (x, t) behave similarly. On the other hand, as (B − At) → 0 then f 0 (x, t) approaches the constant value A, and the two solutions diverge from each other. Then two productivity indices J g [p 0 ] and J g [p](t) also diverge from each other. This phenomenon is observed on the actual field data and has a clear practical explanation. Notice that the denominator in the formula for J g [p 0 ] (165) is a measure of the pressure drawdown needed to maintain constant production Q 0 . As long as the gas reserves are considerably larger than the the pressure drawdown then the distributed source term f 0 (equivalent to reservoir fluid injection) needed to maintain constant production Q 0 is negligible. Otherwise when the gas reserves are comparable in magnitude with the pressure drawdown, then a possible way to maintain constant production rate is by resupplying the reservoir by fluid injections.
In the remaining part of this section we will theoretically investigate the difference between the functions p(x, t) (the actual solution of the problem) and p 0 (x, t) (the solution of the auxiliary problem) depending on the key parameter T crit = B/A. Unfortunately for the Forchheimer case we were not yet able to obtain the appropriate estimates for the differences between p 0 (x, t) and p(x, t). We will report mathematically rigorous result only for the case of compressible Darcy flow. We will show that for the fixed T 0 for all times t ∈ [0, T 0 ] when p(x, t) > 0 on U × [0, T 0 ] the productivity indices J g [p] and J g [p 0 ] are becoming closer to each other with the increasing parameter B. The obtained results make us believe this comparison can be proved in the general Forchheimer case as well.
We consider the case of positive solutions. For that assume that our time space domain belongs to the time layer 0 ≤ t ≤ T 0 < T crit . Let D = U × (0, T 0 ]. The following useful inequality follows directly from maximum principle Proof. Inequality (169) follows from the maximum principle (see, for example, [12]), since the nonlinear equation (148) is uniformly parabolic in the domain D, and the boundary function in (149) is decreasing with time.
6.1. Analytical comparison of the solutions for case of the Darcy flow. First the following maximum principle follows from the results in [27].
The following comparison principle holds: if From the above Lemma follows Proposition 6.7. Under the conditions of Lemma 6.6 then the following comparison holds Proof. Since f 0 (x, t) > 0 for all (x, t) (see (163)) and according to Lemma 6.5 p 2 0 (x, t) − p 2 (x, t) > 0 for all (x, t), in the RHS of (177) we have By Poincaré inequality we have Finally, estimating the RHS of (177) using (184), neglecting the first term in LHS, and using (185) for the second one, we get: Choosing appropriate ε we get (183).
6.2. Stability of PI with respect to initial and boundary data. In this section we will show that for the given time interval 0 < t ≤ T 0 < T crit the difference between the PI's for the actual and the auxiliary problems becomes small as the parameter B for the initial reserves becomes large. As before we consider only the linear Darcy case with F = 0 in (146). Let p(x, t) be the classical solution of IBVP (148)-(150) with the corresponding productivity index J g [p](t) defined by (153). Let p 0 (x, t) be the auxiliary pressure given by (154) with the corresponding productivity index J g [p 0 ] defined by (165). In linear case the original pressure p(x, t) inherits the following properties of the auxiliary pressure p 0 (x, t) (for the details see [27]): for all t, 0 < t ≤ T 0 < T crit |∇p(x, t)| + ∆p(x, t) ≤ C < ∞; (186) |p t (x, t)| ≤ C < ∞.