Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation: The Saturation Convection-Dispersion Equation Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation: The Saturation Convection-Dispersion Equation

In this work, we apply the method of matched asymptotic expansions to solve the one- dimensional saturation convection-dispersion equation, a nonlinear pseudo-parabolic partial differential equation. This equation is one of the governing equations for two- phase flow in a porous media when including capillary pressure effects, for the specific initial and boundary conditions arising when injecting water in an infinite radial piece- wise homogeneous horizontal medium containing oil and water. The method of matched asymptotic expansions combines inner and outer expansions to construct the global solu- tion. In here, the outer expansion corresponds to the solution of the nonlinear first-order hyperbolic equation obtained when the dispersion effects driven by capillary pressure became negligible. This equation has a monotonic flux function with an inflection point, and its weak solution can be found by applying the method of characteristics. The inner expansion corresponds to the shock layer, which is modeled as a traveling wave obtained by a stretching transformation of the partial differential equation. In the transformed domain, the traveling wave solution is solved using regular perturbation theory. By combining the solution for saturation with the so-called Thompson-Reynolds steady-state theory for obtaining the pressure, one can obtain an approximate analytical solution for the wellbore pressure, which can be used as the forward solution which analyzes pressure data by pressure-transient analysis.


Introduction
In this chapter, we show how to generate a semi-analytical solution for the wellbore pressure response during a water injection test. In the petroleum industry, well testing is a common practice which consists of wellbore pressure and wellbore flow rate data acquisition in order to estimate parameters that govern flow in the porous media, i.e., the reservoir rock which stores the hydrocarbons. Well tests give an insight into the oil and gas field production potential and profitability and allow the estimation of reservoir parameters. Estimated parameters can be used to calibrate the reservoir numerical simulation model that are used to describe the fluid flow in these reservoirs and forecast their performance as well as to maximize the productivity of the wells. Injections are important tests on reservoirs containing high amount of harmful gases like carbon dioxide and sulfur dissolved in the oil, causing conventional production testing in the exploratory phase of offshore field development inviable. Multiphase flow is the norm in petroleum reservoirs, and an injection test consists of a period of water or gas injection into an oil reservoir ( Figure 1), a common technique known as waterflooding or gasflooding that is used to displace oil to a producing well. Data from an injection test can be used to estimate the reservoir rock absolute permeability (k), the skin zone permeability (k s ), and the water endpoint relative permeability (a w ). The skin zone permeability is the rock permeability in the zone around the well which was stimulated or damaged during the wellbore drilling operation, while the water endpoint permeability is a measure of how easy water can flow in a specific porous media when there is immobile oil present. In the pursuance of modeling the wellbore pressure response during a water injection test, the Rapport-Leas equation [1], a nonlinear pseudo-parabolic convection-dispersion equation, is used to determine the water saturation distribution in the reservoir as a function of time by assuming a one-dimensional homogeneous medium containing incompressible fluids. Water saturation (S w ) is the fraction of water in a given pore space, and it is expressed in water volume by pore volume. In [2][3][4], it has been shown how to obtain the wellbore pressure response for the case when capillary Figure 1. Sketch of the injection test. Reservoir is assumed to be at rest at the beginning of the test with constant pressure and immobile water saturation distribution (a). Water is injected at constant flow rate leading to pressure change that propagates from the well (b). pressure effects are negligible, i.e., when dispersion effects are not significant. In this case, the Rapport-Leas equation reduces to the Buckley-Leverett [5] equation, a nonlinear hyperbolic equation. In this work, we have extended their model to include capillary pressure. Although the wellbore pressure during injection seems to be insensitive to capillarity effects insensitive to the accuracy of the calculated saturation distribution in the reservoir, knowledge of the correct saturation profile at the end of injection represents the initial condition and hence is required to calculate the saturation distribution during subsequent tests as shut-in (falloff) and flowback (production) test which would allow the estimation of relative permeabilities and capillary pressure curves. Once the water saturation distribution is determined for each time, the corresponding pressure solution can be obtained by integrating the expression for the pressure gradient, given by Darcy's law, from the wellbore radius to infinity while assuming an infiniteacting reservoir. Because Darcy's law does not assume incompressible flow, the pressure solution is transient and does not need to assume incompressible flow even though the saturation profile is generated from a incompressible assumption. To actually evaluate this integral which represents the pressure solution, however, we must assume that the reservoir rate profile becomes constant in a region from the wellbore to a radius such that all the injected water is contained within the reservoir volume within this radius ( Figure 2); this radius increases with time [6]. The region within this radius is referred to as the steady-state region or zone. Intuitively, the assumption that this steady-state zone exists appears to be more tenuous as the total compressibility of the system increases. However, the assumption that this steady-state zone exists has shown to yield accurate semi-analytical pressure solutions for gas-condensate systems [7].

Mathematical model
The solutions presented assume infinite-acting one-dimensional radial flow from and to a fully penetrating vertical well with no gravity effects. We apply the method of matched asymptotic expansions to solve the one-dimensional saturation convection-dispersion equation, a nonlinear pseudo-parabolic partial differential equation. This equation is one of the governing equations for two-phase flow in a porous media when including capillary pressure effects, for Relationship between the water saturation (S w ), solid blue curve, the dimensionless total flow rate profile (q D ), and dotted red curve, during injection. the specific initial and boundary conditions arising when injecting water in an infinite radial piecewise homogeneous horizontal medium containing oil and water. The method of matched asymptotic expansions combines inner and outer expansions to construct the global solution. In here, the outer expansion corresponds to the solution of the nonlinear first-order hyperbolic equation obtained when the dispersion effects driven by capillary pressure became negligible. This equation has a monotonic flux function with an inflection point, and its weak solution can be found by applying the method of characteristics. The inner expansion corresponds to the shock layer, which is modeled as a traveling wave obtained by a stretching transformation of the partial differential equation. By combining the solution for saturation with the so-called Thompson-Reynolds steady-state theory, one can obtain an approximate analytical solution for the wellbore pressure, which can be used as the forward solution which analyzes pressure data by pressure-transient analysis. Let us start by finding the saturation distribution in the reservoir during injection and show how to find pressure.

Saturation
The water mass balance equation, in radial coordinates, leads to the following nonlinear partial differential equation [5]: where throughout we assume that porosity (ϕ) is homogeneous; q t is the total liquid rate in RB/D; θ represents in general a unit conversion factor where in the oil field units used here, θ =5.6146/24; the reservoir thickness, h, and the radius, r, are in ft; and time, t, is in hours. Let us use Darcy's equation in radial coordinates without gravity for the oil (o) and water (w) flow rate in RB/D given by For field units used throughout, α = 141.2. p p , is the phase p pressure. The λ p is the phase p mobility, given by the ratio of the phase permeability (k rw or k ro ), which are functions of the water saturation, by the phase viscosity (μ w or μ o ). To find the water fractional flow (F w ), we can subtract Eq. (2) for water from Eq. (2) for oil, to get Rearranging Eq. (3), substituting the capillary pressure p c given by the difference of the oil pressure (p o ) and the water pressure (p w ), and dividing the resulting equation by the total flow rate, q t , yield where F o and F w are the oil and water fractional flow given by q o and q w , respectively. We assume throughout that water is the wetting phase. Finally, substituting F o ¼ 1 À F w in Eq. (4) and solving for F w , we have the following expression for the water fractional flow including capillary pressure effects where f w is the water mobility ratio (Figure 3), i.e., the ratio of water mobility and the total mobility (λ t ), given by which usually assumes an S-shape. e is the perturbation parameter, defined by and the permeability is a function of radius because we consider a skin-damaged zone: where r w is the wellbore radius, r skin is the radius of the damaged zone, and k s is the permeability in the skin zone. Grouping all the parameters that are the function of water saturation, we can define ( Figure 4) For simplicity, we use the Brooks and Corey model [8] given by to represent capillary pressure. Here, S iw is the immobile water saturation and S or is the residual oil saturation. λ, where 0:4 ≤ λ ≤ 4:0, is a measure of the pore size distribution (the greater the λ value, the more uniform is the pore size distribution), and p t is the threshold pressure. The threshold pressure is a measure of the maximum pore size [9], i.e., the minimum capillary pressure at which a continuous nonwetting phase exists in the imbibition case and a continuous wetting phase exists in the drainage case [10]. The greater is the maximum pore size, the smaller is the pressure threshold. According to [11], the extrapolation of the capillary pressure curve obtained from experimental data to S w ¼ 1 yields the correct threshold value. In practice, we introduce a small variable, sv, to limit the maximum value of p c to a finite value. We can relate the relative permeabilities and the capillary pressure through λ by using the [12] model for the water phase (wetting phase) and the [8] model for the oil phase (nonwetting phase) [13]  . dΨ dSw versus water saturation for a partially water wet reservoir (b). For a strong water, we reservoir, dΨ dSw ! ∞ at S w ¼ S iw . Now that we have defined the fractional flow rate and its parameters, let us go back to our governing equation for saturation (Eq. (1)). Inserting Eq. (10) into Eq. (1) and defining yields which is the nonlinear "pseudo-parabolic" governing equation for saturation. If we insert same common values for the parameters in Eq. (7) to have an idea of its order of magnitude, we can see that epsilon is a very small number. This suggests that the effect of the third term in Eq. (15) may be treated as a perturbation to the first-order hyperbolic equation [5], given by where f w is considered to be an S-shaped function along this chapter. During injection, for a partially water wet reservoir, the capillary pressure dispersive effect will be non-negligible only in a small region around the water front (hypodispersion phenomenon) [14,15] where the capillary pressure derivative and the saturation gradient are significant ( Figure 5). The capillary pressure smears the water front during injection balancing the self-sharpening tendency of the shock. [16] have developed an exact analytical solution for linear waterflood including the effects of capillary pressure, but their solution is limited to a particular functional form to represent relative permeabilities and capillary pressure curves and does not consider radial flow, which makes their solution very restrictive. As done by [17,18] for Cartesian coordinates and by [19] for streamlines and streamtubes, the perturbation caused by the capillary pressure effects can be modeled as a shock layer (water front) which moves with the same speed as the shock wave. By applying the method of matched asymptotic expansions [20,21], we can we divide the domain into two regions, outer and inner regions (Figure 6), where the inner region, the region around the water front, is modeled as a shock layer which propagates with the same speed as the shock would be obtained when e ! 0, i.e., when the capillary pressure effects are null. The combination of the self-sharpening tendency of the shock (S wf > S iw ) with the dispersive effect from the capillary pressure balance against each other leads to the shock layer [22]. Note: in order to guarantee pressure continuity, we have assumed that the capillary pressure gradient is zero at the wellbore, which means that F w ¼ f w ¼ 1 is the wellbore, so S w r w ; t ð Þ¼1 À S or . This boundary condition will be used for both the Buckley-Leverett and the Rapoport-Leas solutions. Ref. [23] presented the idea of using the method of matched asymptotic equations to solve the Rapoport-Leas equation, while [18,24,19] showed how the mass balance could be used to present a closed solution for the saturation distribution. Ref. [24] derived an approximate solution for the [1] in Cartesian coordinates for both water and oil injections into a core considering end effects by also applying the method of matched asymptotic expansions. The method of asymptotic expansions uses the inner and outer saturation solutions combined with a matching function in order to obtain a composite solution which avoids abruptly switching from the outer to the inner solution or vice versa. The inner and outer solutions are each capable of representing the real solution in two distinct regions-the "inner region" and the "outer region" of the boundary layer ( Figure 6). Similarly, we approximate the saturation solution of Rapoport-Leas equation by forming a composite solution given by the combination of three saturations: S BL w , the solution obtained when the capillary pressure effects are neglected; S SL w , the saturation distribution in the shock layer obtained by magnifying the dispersion effects in the saturation governing equation; and S SH w , the shock wave represented by a Heaviside function: where BL stands for Buckley-Leverett, SL for shock layer, and SH for shock function.

Outer solution (S BL w )
The outer solution, S BL w , is obtained by letting e ! 0 in Eq. (15): That is the nonlinear hyperbolic convection equation known as the Buckley-Leverett saturation equation given by Eq. (16) which is obtained when capillary and gravity effects are neglected. The well-known unique admissible weak solution of this Riemann problem, with the following initial condition can be obtained by the application of the method of characteristics and is given by [5].
that is, by a family of rarefaction waves, a semi-shock wave, and a constant saturation zone where water is immobile. The shock jump is caused by the S-shaped form of the fractional flow curve, which leads to a gradient catastrophe and consequently a shock solution. This semishock has a constant speed, satisfying the Rankine-Hugoniot condition [25]: where S wf and S iw are the shock saturations. In this case, in order to satisfy the conservation of mass, the shock speed should correspond to the slope of a tangent line to the water fractional flow curve, i.e., The details of this solution can be found in [5]. Figure 7 shows the shock jump slope tangent to the fractional flow curve at S w ¼ S wf and the saturation distribution in the reservoir at a time t. The rarefaction wave family spans from 1 À S or to S wf from r w to r ¼ 25 ft, the water front position, i.e., the shock front position, r s . Ahead of the water front position, there is an immobile water. Figure 8 compares this solution, the outer solution, with the true solution; there is the convection-dispersion saturation profile. Here, we call the true solution the solution obtained from a numerical simulator.

Inner solution (S SL w )
As mentioned, the inner solution intends to represent the saturation profile in the "inner" region around the water front, which is a shock layer (a boundary layer) around the shock traveling Figure 7. The shock jump slope tangent (blue curve) to the S-shaped fractional flow curve at S w ¼ S wf (a) and the saturation profile in the reservoir at a time t (b). The rarefaction waves family spans from 1 À S or 1 to S wf and from r w to r = 25 ft, the water front position, i.e., the shock front position, r f , inj . Ahead of the water front position, there is immobile water. with the same speed as the shock itself ( Figure 9). In order to find S SL w , we magnify the shock layer by using a stretching traveling wave coordinate. Similarly, as defined in [24,18,19] where r s is the shock front position: w is zero at r ¼ r s and goes to AE∞ as e ! 0. We rewrite Eq. (15) in terms of moving coordinates, The inner solution is obtained by letting e ! 0 in Eq. (28): as presented in [26]. Therefore, neglecting the terms of order e in Eq. (28), we have Note that here we are treating the permeability k as function of the shock position radius, r s , only, by assuming that in the limit of the inner solution, e r ! r s τ ð Þ ð Þ . Intuitively, this assumption does not seem valid when the shock layer is crossing heterogeneity interfaces, i.e., interfaces between two different permeability zones. However, for the water injection in a field Figure 9. Saturation distribution during the injection period for the inner solution (S SL w ) and for the true solution (S w ), with capillary pressure. Both profiles agree in the region around the water front, i.e., in the shock boundary layer which we have defined as the inner region. scale, the skin zone will be crossed by the water front in a very short time, and we will only need to use the pseudo-parabolic equation (Eq. (15)) to find saturation for the end of injection period (to be used as initial condition for falloff and flowback tests, as mentioned in the introduction). Consequently, we can simplify the problem as shown above. Integrating the ordinary differential equation given by Eq. (30) with respect to w for any fixed time τ and applying the chain rule gives where a τ ð Þ is constant for the injection case, as we will show later. As mentioned, the inner solution is modeled as a traveling wave with a constant speed-the shock speed-and the boundary conditions (for the inner solution) given by as the inner solution goes asymptotically to the shock saturations. This necessity of this behavior will be clearer very soon when we compare the inner solution with the matching saturation solution. Using the first boundary condition given by Eq. (32) in Eq. (31) leads to while using the second boundary condition given by Eq. (32) yields implying that ÀD S_ iw f g À S_ wf f g ð Þ À f _w S_ wf f g ð Þ¼0, which it is indeed correct from the definition of D in Eq. (24). As we can see from Eqs. (33) and (34), a τ ð Þ is a constant and it will be called simply by a from now on. Substituting the constant a (Eq. (33)) in Eq. (31) and dividing it by D S iw À S SL Integrating Eq. (35) from w well ¼ w r w ; τ ð Þ to any w at any time τ gives us the relationship between any S SL w and w: At S SL w ¼ S wf , the integral in the left side of Eq. (36) diverges as the integrand denominator goes to 0. This behavior is consistent with our boundary condition assumptions for S SL w (Eq. (32)). At S SL w ¼ S iw , the integral in the left side of Eq. (36) does converge when the integrand numerator also goes to zero (Figure 4), a behavior which is consistent when trying to model a hypodispersion phenomenon (for a partially water wet reservoir). In another reservoir wettability scenario, e.g., a strong water wet rock, the capillary pressure would not be bounded at Siw, and the integral in the left side of Eq. (36) would diverge. Note: we still do not know the value of S SL w at w well. In order to find a closed form for this problem, mass balance can be used, but first let us present the matching saturation, since this solution will be necessary for the mass balance.

Matching solution (S SH w )
The matching saturation S SH w is defined using the matching principle by applying Prandtl's technique [26]: and in the injection case is given by which is plotted in Figure 10 against the outer and inner solutions. As we were searching for, S SH w matches with the outer solution in the inner region and with the inner solution in the outer region, being able to subtract their effect in the composite solution in their "noncorrespondent" zones. Figure 11 compares the saturation distribution during the injection period for the true solution obtained from the numerical simulator IMEX with the outer, inner, and matching saturation solutions.

Mass balance
Now that we have defined all the three saturations that are composed of the approximate solution for the convection dispersion saturation equation, let us try to find a closed form for Setting S w ¼ S iw in the upper limits of the integrals of Eq. (36) and exchanging the two sides of the equation yield As the left sides of Eqs. (46) and (47) are the same, the right sides of these two equations must be equal which gives Multiplying Eq. (48) by eS iw and rearranging the resulting equation give Once the value S SL w À CDτ e À Á (i.e., the inner solution saturation in the wellbore S SL w w well ð Þ) is determined numerically by solving Eq. (49) using the bisection method at each time τ, Eq. (36) is used to determine the saturation profile in the stabilized zone. It is important to note that, as S SL w should reach S wf and S iw asymptotically as w ! AE∞, here we did not have to fix a finite distance in which the traveling wave would reach its open bounds as done by [18,19,24]. With our approach, as shown in the validation section (Figure 12), we can obtain essentially a perfect match with the numerical solution, with a "smoother" water front, which is expected from the dispersive effect of capillary pressure, contrary to the sharp transition between the saturation at the water front foot (r wf ) which is the finite position at which water can be considered immobile-and the initial water saturation in the oil zone exhibited by solutions of previous authors.

Wellbore pressure
As mentioned previously, after finding the saturation distribution, we can obtain the wellbore pressure by applying the pressure solutions presented by [2]. During injection at a constant flow rate, q t r w ; t ð Þ RB/D, where t ¼ 0 at the beginning of the water injection, by integrating Darcy's law as in [6,2], given by where p w ¼ p o À p c . Eq. (50) can be solved for the oil pressure gradient by integrating it from the wellbore radius to infinite, assuming an infinite-acting reservoir. The bottom hole pressure difference from the reservoir initial pressure (p oi ) can then be expressed as where it is assumed that p o r w ; t ð Þ¼p w r w ; t ð Þ, i.e., p c ¼ 0 at r ¼ r w , in order to satisfy the compatibility condition [27], i.e., to guarantee phase pressure continuity at the wellbore. Eq. (51) can be rewritten as by assuming that the second term in the right-hand side of Eq. (52) is zero from r wf to ∞, considering f w S iw ð Þ ¼ 0 and ∂Sw ∂r ¼ 0 for r > r wf t ð Þ, since the water in the region ahead of the water front foot is assumed immobile. r wf can be defined as the position at which S w À S iw ð Þ< δ, where δ is a very small number. For the hypodispersion phenomenon, we can find a finite r wf where δ ! 0, i.e., at which S w ¼ S iw . Using the [6] steady-state theory, which assumes that, q t r; t ð Þ ¼ q t r w ; t ð Þ, for r ≤ r wf t ð Þ, Eq. (52) becomes   Table 1. Reservoir, rock, and fluid properties for simulation and analytical solution. Figure 13. Comparison of the wellbore pressure response from analytical solution and IMEX during water injection test with capillary pressure (a). Figure 14. The log-log diagnostic plots for wellbore pressure data (blue curve) and its derivative with respect to time (red curve) during water injection (b).

Conclusions
In this work, an accurate approximate analytical solution was constructed for wellbore pressure during water injection test in a reservoir containing oil and immobile water. Our solution was validated by comparing the bottom hole pressure calculated from the analytical model with the data obtained from a commercial numerical simulator. Our solution presented here for water injection together with the wellbore pressure and flow rate history for subsequent tests as shut-in and flowback can be used as forward model in a nonlinear regression in order to estimate relative permeabilities and capillary pressure curves in addition to the rock absolute permeability, the skin zone permeability, and the water endpoint relative permeability.