Reservoir, rock, and fluid properties for simulation and analytical solution.
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 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. 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.
- method of matched asymptotics
- boundary layer approximation
- nonlinear pseudo-parabolic partial differential equation
- convection-dispersion phenomenon
- multiphase flow in porous media
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 (), the skin zone permeability (), and the water endpoint relative permeability (). 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 , 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 () 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 pressure effects are negligible, i.e., when dispersion effects are not significant. In this case, the Rapport-Leas equation reduces to the Buckley-Leverett  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 infinite-acting 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 . 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 .
2. 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 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.
The water mass balance equation, in radial coordinates, leads to the following nonlinear partial differential equation :
where throughout we assume that porosity () is homogeneous; 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, , and the radius, , are in ft; and time, t, is in hours. Let us use Darcy’s equation in radial coordinates without gravity for the oil () and water () flow rate in RB/D given by
For field units used throughout, = 141.2. , is the phase p pressure. The is the phase p mobility, given by the ratio of the phase permeability (or ), which are functions of the water saturation, by the phase viscosity (or ). To find the water fractional flow (), we can subtract Eq. (2) for water from Eq. (2) for oil, to get
Rearranging Eq. (3), substituting the capillary pressure pc given by the difference of the oil pressure () and the water pressure (), and dividing the resulting equation by the total flow rate, , yield
where and are the oil and water fractional flow given by and , respectively. We assume throughout that water is the wetting phase. Finally, substituting in Eq. (4) and solving for , we have the following expression for the water fractional flow including capillary pressure effects
where is the water mobility ratio (Figure 3), i.e., the ratio of water mobility and the total mobility (), given by
which usually assumes an S-shape. is the perturbation parameter, defined by
and the permeability is a function of radius because we consider a skin-damaged zone:
where is the wellbore radius, is the radius of the damaged zone, and is the permeability in the skin zone. Grouping all the parameters that are the function of water saturation, we can define (Figure 4)
and rewrite Eq. (10) as
For simplicity, we use the Brooks and Corey model  given by
to represent capillary pressure. Here, is the immobile water saturation and is the residual oil saturation. , where , is a measure of the pore size distribution (the greater the value, the more uniform is the pore size distribution), and is the threshold pressure. The threshold pressure is a measure of the maximum pore size , 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 . The greater is the maximum pore size, the smaller is the pressure threshold. According to , the extrapolation of the capillary pressure curve obtained from experimental data to yields the correct threshold value. In practice, we introduce a small variable, sv, to limit the maximum value of to a finite value. We can relate the relative permeabilities and the capillary pressure through by using the  model for the water phase (wetting phase)
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 , given by
where 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.  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  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 combine the solution of the Buckley-Leverett equation (Eq. (16)) with this steady traveling wave to generate an approximate solution of the Rappaport and Leas equation, i.e., the solution of the convection–dispersion saturation equation. In order to solve the  equation for the injection period, with the following initial and boundary conditions
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 , i.e., when the capillary pressure effects are null. The combination of the self-sharpening tendency of the shock () with the dispersive effect from the capillary pressure balance against each other leads to the shock layer . Note: in order to guarantee pressure continuity, we have assumed that the capillary pressure gradient is zero at the wellbore, which means that is the wellbore, so . This boundary condition will be used for both the Buckley-Leverett and the Rapoport-Leas solutions. Ref.  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.  derived an approximate solution for the  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: , the solution obtained when the capillary pressure effects are neglected; , the saturation distribution in the shock layer obtained by magnifying the dispersion effects in the saturation governing equation; and , the shock wave represented by a Heaviside function:
where stands for Buckley-Leverett, for shock layer, and for shock function.
2.1.1. Outer solution ()
The outer solution, , is obtained by letting 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 .
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 semi-shock has a constant speed, satisfying the Rankine-Hugoniot condition :
where and 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 . Figure 7 shows the shock jump slope tangent to the fractional flow curve at and the saturation distribution in the reservoir at a time . The rarefaction wave family spans from to from to ft, the water front position, i.e., the shock front position, rs. 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.
2.1.2. Inner solution ()
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 with the same speed as the shock itself (Figure 9). In order to find , we magnify the shock layer by using a stretching traveling wave coordinate. Similarly, as defined in [24, 18, 19]
where rs is the shock front position:
The inner solution is obtained by letting in Eq. (28):
Note that here we are treating the permeability k as function of the shock position radius, rs, only, by assuming that in the limit of the inner solution, . 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 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 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 , which it is indeed correct from the definition of D in Eq. (24). As we can see from Eqs. (33) and (34), is a constant and it will be called simply by a from now on. Substituting the constant (Eq. (33)) in Eq. (31) and dividing it by yield
Integrating Eq. (35) from to any at any time gives us the relationship between any and :
At , 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 (Eq. (32)). At , 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 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.
2.1.3. Matching solution ()
The matching saturation is defined using the matching principle by applying Prandtl’s technique :
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, 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 “non-correspondent” 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.
2.1.4. 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 the saturation distribution based in the mass balance. Since both the Buckley-Leverett () solution and the composite solution () must obey material balance, the following equality
which, upon simplification, gives
From Eq. (35),
Setting in the upper limits of the integrals of Eq. (36) and exchanging the two sides of the equation yield
Multiplying Eq. (48) by and rearranging the resulting equation give
Once the value (i.e., the inner solution saturation in the wellbore ) 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 should reach and asymptotically as , 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 () 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.
2.2. Wellbore pressure
As mentioned previously, after finding the saturation distribution, we can obtain the wellbore pressure by applying the pressure solutions presented by . During injection at a constant flow rate, RB/D, where at the beginning of the water injection, by integrating Darcy’s law as in [6, 2], given by
where . 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 () can then be expressed as
by assuming that the second term in the right-hand side of Eq. (52) is zero from to , considering and for , since the water in the region ahead of the water front foot is assumed immobile. can be defined as the position at which , where is a very small number. For the hypodispersion phenomenon, we can find a finite where , i.e., at which . Using the  steady-state theory, which assumes that, , for , Eq. (52) becomes
where for any practical set of values of physical properties  indicate that this assumption is valid. Adding and subtracting the term , where is the endpoint oil mobility at , Eq. (53) can be rewritten as
is the single-phase oil transient pressure drop, the known pressure drop solution that is obtained if we inject oil into an oil reservoir (injection period), whose well-known approximate solution can be approximated as
Here, is a unit conversion factor in which oil field unit is 0.0002637 and the single-phase total compressibility is
We have compared our pressure and saturation solution including capillary pressure effects with the commercial numerical simulator IMEX, using the properties shown in Table 1. Figure 12 compares the saturation distribution obtained from our analytical solution with the one obtained with IMEX, while Figure 13 shows the comparison of the wellbore pressure response from our analytical solution and IMEX during injection test. In order to be able to match saturation and pressure obtained from our solution with IMEX, we have to use a very refined grid (0.01 ft) around the wellbore in the zone invaded by water and then increase it exponentially to a very large external radius (10,000 times the wellbore radius) in order to reproduce an infinite acting reservoir. In addition, we have to start with very short time steps, day. Figure 14 presents the log-log diagnostic plots of injection. We can see that at early times of injection, there is a plateau (stabilization) in the wellbore pressure derivative plot, which by inspection reflects the original total mobility (endpoint oil mobility), while, at late times, we find that the derivative shows stabilized radial flow, which by inspection reflects the endpoint water mobility.
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.
This research was conducted under the auspices of TUPREP, the Tulsa University Petroleum Reservoir Exploitation Projects, and it was prepared with financial support from the Coordination for the Improvement of Higher Education Personnel (CAPES) within the Brazilian Ministry of Education.
Conflict of interest
The authors declare that there is no conflict of interest.
|β||Unit conversion factor (0.0002637)|
|Δpo||Single-phase oil pressure drop (psi)|
|λo||Oil mobility (1/cp)|
|λt||Total mobility (1/cp)|
|λw||Water mobility (1/cp)|
|μo||Oil viscosity (cp)|
|μw||Water viscosity (cp)|
|aw||Water endpoint relative permeability|
|Fo||Oil fractional flow|
|Fw||Water fractional flow|
|fw||Water mobility ratio|
|k||Absolute permeability (mD)|
|ks||Skin permeability (mD)|
|kro||Oil relative permeability (mD)|
|krw||Water relative permeability (mD)|
|pc||Capillary pressure (psi)|
|pi||Reservoir initial pressure (psi)|
|po||Oil pressure (psi)|
|pt||Pressure threshold (psi)|
|pw||Water pressure (psi)|
|qo||Oil flow rate (RB/D)|
|qt||Total liquid rate (RB/D)|
|rs||Shock front position (ft)|
|rskin||Skin zone radius (ft)|
|rw||Wellbore radius (ft)|
|Siw||Immobile water saturation|
|Sor||Residual oil saturation|
|C||Constant given by θqtπhϕ|
|h||Reservoir thickness (ft)|
|w||Traveling wave coordinate|
|α||Unit conversion factor (141.2)|
|λ||Pore size distribution index|
|θ||Unit conversion factor (5.6146/24)|
|co||Oil compressibility (1/psi)|
|cr||Rock compressibility (1/psi)|
|cw||Water compressibility (1/psi)|
|cto||Single-phase total compressibility (1/psi)|