1. 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 (ks), and the water endpoint relative permeability (aw). 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 (Sw) is the fraction of water in a given pore space, and it is expressed in water volume by pore volume. In [2–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 [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

Figure 2. Relationship between the water saturation (Sw), solid blue curve, the dimensionless total flow rate profile (qD),

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

http://dx.doi.org/10.5772/intechopen.76828

119

yield accurate semi-analytical pressure solutions for gas-condensate systems [7].

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

2. Mathematical model

and dotted red curve, during injection.

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).

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation… http://dx.doi.org/10.5772/intechopen.76828 119

Figure 2. Relationship between the water saturation (Sw), solid blue curve, the dimensionless total flow rate profile (qD), and dotted red curve, during injection.

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].

## 2. Mathematical model

1. Introduction

118 Perturbation Methods with Applications in Science and Engineering

propagates from the well (b).

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 (ks), and the water endpoint relative permeability (aw). 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 (Sw) is the fraction of water in a given pore space, and it is expressed in water volume by pore volume. In [2–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 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.

#### 2.1. Saturation

The water mass balance equation, in radial coordinates, leads to the following nonlinear partial differential equation [5]:

$$\frac{\partial \mathcal{S}\_w}{\partial t} + \frac{\partial q\_t}{2\pi r h \phi} \frac{\partial F\_w(\mathcal{S}\_w)}{\partial r} = 0,\tag{1}$$

where Fo and Fw are the oil and water fractional flow given by qo and qw, respectively. We assume throughout that water is the wetting phase. Finally, substituting Fo ¼ 1 � Fw in Eq. (4) and solving for Fw, we have the following expression for the water fractional flow including

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

k rð Þhkro αqtμ<sup>o</sup>

where f <sup>w</sup> is the water mobility ratio (Figure 3), i.e., the ratio of water mobility and the total

<sup>¼</sup> <sup>λ</sup>wð Þ Sw λoð Þþ Sw λwð Þ Sw

> <sup>e</sup> <sup>¼</sup> <sup>h</sup> αqt μo

k rð Þ¼ <sup>k</sup>s, rw <sup>≤</sup> <sup>r</sup> <sup>&</sup>lt; <sup>r</sup>skin k, r ≥ rskin,

where rw is the wellbore radius, rskin is the radius of the damaged zone, and ks is the permeability in the skin zone. Grouping all the parameters that are the function of water saturation,

> dpc dSw

ð Þ¼� Sw f <sup>w</sup>ð Þ Sw kroð Þ Sw

and the permeability is a function of radius because we consider a skin-damaged zone:

which usually assumes an S-shape. e is the perturbation parameter, defined by

r ∂pc ∂r 

¼ f <sup>w</sup> þ erk rð Þf <sup>w</sup>ð Þ Sw kroð Þ Sw

http://dx.doi.org/10.5772/intechopen.76828

λwð Þ Sw λoð Þþ Sw λwð Þ Sw ∂pcð Þ Sw <sup>∂</sup><sup>r</sup> ,

, (6)

ð Þ Sw , (9)

(5)

121

(7)

(8)

<sup>1</sup> <sup>þ</sup> <sup>λ</sup><sup>o</sup> ð Þ Sw λwð Þ Sw

capillary pressure effects

mobility (λt), given by

we can define (Figure 4)

<sup>1</sup> <sup>þ</sup> k rð Þhλoð Þ Sw αqt

> <sup>1</sup> <sup>þ</sup> <sup>λ</sup>oð Þ Sw λwð Þ Sw

> > f <sup>w</sup>ð Þ¼ Sw

r ∂pc ∂r 

<sup>¼</sup> <sup>1</sup> <sup>1</sup> <sup>þ</sup> <sup>λ</sup><sup>o</sup> ð Þ Sw λwð Þ Sw þ

Figure 3. Water fractional flow curve (dark solid curve) and its derivative (dotted curve).

1 <sup>1</sup> <sup>þ</sup> <sup>λ</sup>oð Þ Sw λwð Þ Sw

dΨ dSw

Fwð Þ¼ Sw

where throughout we assume that porosity (ϕ) is homogeneous; qt 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

$$q\_p = -\frac{k(r)h\lambda\_p(S\_w)}{a} \left( r \frac{\partial p\_p}{\partial r} \right), \qquad \text{for } p = o, w. \tag{2}$$

For field units used throughout, α = 141.2. pp, is the phase p pressure. The λ<sup>p</sup> is the phase p mobility, given by the ratio of the phase permeability (krw or kro), which are functions of the water saturation, by the phase viscosity (μ<sup>w</sup> or μo). To find the water fractional flow (Fw), we can subtract Eq. (2) for water from Eq. (2) for oil, to get

$$-\frac{\alpha q\_o}{k(r)h\lambda\_o(S\_w)} - \frac{\alpha q\_w}{k(r)h\lambda\_w(S\_w)} = -\left(r\frac{\partial p\_o}{\partial r} - r\frac{\partial p\_w}{\partial r}\right). \tag{3}$$

Rearranging Eq. (3), substituting the capillary pressure pc given by the difference of the oil pressure (po) and the water pressure (pw), and dividing the resulting equation by the total flow rate, qt , yield

$$F\_o - \frac{F\_w \lambda\_o(S\_w)}{\lambda\_w(S\_w)} = -\frac{k(r)h\lambda\_o(S\_w)}{\alpha q\_t} \left( r \frac{\partial p\_c(S\_w)}{\partial r} \right),\tag{4}$$

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation… http://dx.doi.org/10.5772/intechopen.76828 121

Figure 3. Water fractional flow curve (dark solid curve) and its derivative (dotted curve).

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 reser-

The water mass balance equation, in radial coordinates, leads to the following nonlinear partial

where throughout we assume that porosity (ϕ) is homogeneous; qt 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

> r ∂pp ∂r

For field units used throughout, α = 141.2. pp, is the phase p pressure. The λ<sup>p</sup> is the phase p mobility, given by the ratio of the phase permeability (krw or kro), which are functions of the water saturation, by the phase viscosity (μ<sup>w</sup> or μo). To find the water fractional flow (Fw), we

Rearranging Eq. (3), substituting the capillary pressure pc given by the difference of the oil pressure (po) and the water pressure (pw), and dividing the resulting equation by the total flow

> ¼ � k rð Þhλoð Þ Sw αqt

¼ � r

∂po <sup>∂</sup><sup>r</sup> � <sup>r</sup>

r ∂pcð Þ Sw ∂r 

∂pw ∂r

� <sup>α</sup>qw k rð Þhλwð Þ Sw

∂Fwð Þ Sw

<sup>∂</sup><sup>r</sup> <sup>¼</sup> <sup>0</sup>, (1)

, for p ¼ o, w: (2)

: (3)

, (4)

θqt 2πrhϕ

∂Sw ∂t þ

qp ¼ � k rð Þhλpð Þ Sw α

can subtract Eq. (2) for water from Eq. (2) for oil, to get

αqo k rð Þhλoð Þ Sw

> Fo � Fwλoð Þ Sw λwð Þ Sw

voir during injection and show how to find pressure.

120 Perturbation Methods with Applications in Science and Engineering

2.1. Saturation

differential equation [5]:

rate in RB/D given by

rate, qt

, yield

where Fo and Fw are the oil and water fractional flow given by qo and qw, respectively. We assume throughout that water is the wetting phase. Finally, substituting Fo ¼ 1 � Fw in Eq. (4) and solving for Fw, we have the following expression for the water fractional flow including capillary pressure effects

$$F\_w(\mathbf{S}\_w) = \frac{1 + \frac{k(r)k\lambda\_c(\mathbf{S}\_w)}{aq\_t} \left(r \frac{\partial p}{\partial r}\right)}{1 + \frac{\lambda\_c(\mathbf{S}\_w)}{\lambda\_w(\mathbf{S}\_w)}} = \frac{1}{1 + \frac{\lambda\_d(\mathbf{S}\_w)}{\lambda\_w(\mathbf{S}\_w)}} + \frac{\frac{k(r)k\lambda\_w}{aq\_t\mu\_s} \left(r \frac{\partial p}{\partial r}\right)}{1 + \frac{\lambda\_d(\mathbf{S}\_w)}{\lambda\_w(\mathbf{S}\_w)}} = f\_w + \epsilon rk(r) f\_w(\mathbf{S}\_w) k\_{ro}(\mathbf{S}\_w) \frac{\partial p\_c(\mathbf{S}\_w)}{\partial r},\tag{5}$$

where f <sup>w</sup> is the water mobility ratio (Figure 3), i.e., the ratio of water mobility and the total mobility (λt), given by

$$f\_w(\mathcal{S}\_w) = \frac{1}{1 + \frac{\lambda\_v(\mathcal{S}\_w)}{\lambda\_w(\mathcal{S}\_w)}} = \frac{\lambda\_w(\mathcal{S}\_w)}{\lambda\_o(\mathcal{S}\_w) + \lambda\_w(\mathcal{S}\_w)} \frac{\lambda\_w(\mathcal{S}\_w)}{\lambda\_o(\mathcal{S}\_w) + \lambda\_w(\mathcal{S}\_w)}\,\tag{6}$$

which usually assumes an S-shape. e is the perturbation parameter, defined by

$$\mathbf{e} = \frac{h}{\alpha \eta\_t \mu\_o} \tag{7}$$

and the permeability is a function of radius because we consider a skin-damaged zone:

$$k(r) = \begin{cases} k\_{s\prime} & r\_w \le r < r\_{\text{skin}} \\ k\_{\prime} & r \ge r\_{\text{skin}\prime} \end{cases} \tag{8}$$

where rw is the wellbore radius, rskin is the radius of the damaged zone, and ks is the permeability in the skin zone. Grouping all the parameters that are the function of water saturation, we can define (Figure 4)

$$\frac{d\Psi}{dS\_w}(S\_w) = -f\_w(S\_w)k\_{ro}(S\_w)\frac{dp\_c}{dS\_w}(S\_w),\tag{9}$$

Figure 4. <sup>d</sup><sup>Ψ</sup> dSw versus water saturation for a partially water wet reservoir (b). For a strong water, we reservoir, <sup>d</sup><sup>Ψ</sup> dSw ! ∞ at Sw ¼ Siw.

and rewrite Eq. (10) as

$$F\_w(\mathcal{S}\_w) = f\_w(\mathcal{S}\_w) - \epsilon r k(r) \frac{d\Psi(\mathcal{S}\_w)}{d\mathcal{S}\_w} \frac{\partial \mathcal{S}\_w}{\partial r} \,. \tag{10}$$

Now that we have defined the fractional flow rate and its parameters, let us go back to our

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

<sup>C</sup> <sup>¼</sup> <sup>θ</sup>qt

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)

> C 2r ∂f <sup>w</sup>

where f <sup>w</sup> 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

Figure 5. Capillary pressure curve (a) and saturation profile at a time t during water injection (b). The dashed green lines

represent the capillary pressure derivative (a) and the saturation gradient (b) at the water front.

may be treated as a perturbation to the first-order hyperbolic equation [5], given by

∂Sw ∂t þ

rk rð Þ <sup>∂</sup><sup>Ψ</sup> ∂r 

<sup>π</sup>h<sup>ϕ</sup> , (14)

<sup>∂</sup><sup>r</sup> <sup>¼</sup> <sup>0</sup>, (16)

¼ 0, (15)

http://dx.doi.org/10.5772/intechopen.76828

123

governing equation for saturation (Eq. (1)). Inserting Eq. (10) into Eq. (1) and defining

∂Sw ∂t þ

C 2r ∂f <sup>w</sup> <sup>∂</sup><sup>r</sup> � <sup>e</sup> <sup>C</sup> 2r ∂ ∂r

yields

For simplicity, we use the Brooks and Corey model [8] given by

$$p\_c(\mathcal{S}\_w) = p\_t \left( s\upsilon + \frac{\mathcal{S}\_w - \mathcal{S}\_{wi}}{1 - \mathcal{S}\_{wi} - \mathcal{S}\_{or}} \right)^{-\frac{1}{\lambda}} \tag{11}$$

to represent capillary pressure. Here, Siw is the immobile water saturation and Sor 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 pt 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 Sw ¼ 1 yields the correct threshold value. In practice, we introduce a small variable, sv, to limit the maximum value of pc 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)

$$k\_{rw} = a\_w \left(\frac{\mathcal{S}\_w - \mathcal{S}\_{iw}}{1 - \mathcal{S}\_{iw} - \mathcal{S}\_{or}}\right)^{\frac{2+\lambda}{\lambda}},\tag{12}$$

and the [8] model for the oil phase (nonwetting phase) [13]

$$k\_{ro} = \left(1 - \frac{S\_{w} - S\_{iw}}{1 - S\_{iw} - S\_{or}}\right)^{2} \left(1 - \left(\frac{S\_{w} - S\_{iw}}{1 - S\_{iw} - S\_{or}}\right)^{\frac{2 + \lambda}{\lambda}}\right). \tag{13}$$

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

$$\mathbf{C} = \frac{\Theta \eta\_t}{\pi h \phi'} \tag{14}$$

yields

and rewrite Eq. (10) as

Figure 4. <sup>d</sup><sup>Ψ</sup>

Sw ¼ Siw.

Fwð Þ¼ Sw <sup>f</sup> <sup>w</sup>ð Þ� Sw <sup>e</sup>rk rð Þ <sup>d</sup>Ψð Þ Sw

dSw versus water saturation for a partially water wet reservoir (b). For a strong water, we reservoir, <sup>d</sup><sup>Ψ</sup>

to represent capillary pressure. Here, Siw is the immobile water saturation and Sor 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 pt 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 Sw ¼ 1 yields the correct threshold value. In practice, we introduce a small variable, sv, to limit the maximum value of pc to a finite value. We can relate the relative permeabilities and the capillary pressure through λ by using the [12]

> Sw � Siw 1 � Siw � Sor � �<sup>2</sup>þ<sup>λ</sup>

λ

<sup>1</sup> � Sw � Siw 1 � Siw � Sor � �<sup>2</sup>þ<sup>λ</sup>

!

For simplicity, we use the Brooks and Corey model [8] given by

122 Perturbation Methods with Applications in Science and Engineering

model for the water phase (wetting phase)

pcð Þ¼ Sw pt sv þ

krw ¼ aw

1 � Siw � Sor � �<sup>2</sup>

and the [8] model for the oil phase (nonwetting phase) [13]

kro <sup>¼</sup> <sup>1</sup> � Sw � Siw

dSw

Sw � Swi 1 � Swi � Sor � ��<sup>1</sup>

∂Sw

λ

<sup>∂</sup><sup>r</sup> : (10)

dSw ! ∞ at

, (11)

, (12)

: (13)

λ

$$
\frac{\partial \mathcal{S}\_w}{\partial t} + \frac{\mathcal{C}}{2r} \frac{\partial f}{\partial r} - \epsilon \frac{\mathcal{C}}{2r} \frac{\partial}{\partial r} \left( rk(r) \frac{\partial \Psi}{\partial r} \right) = 0,\tag{15}
$$

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

$$\frac{\partial S\_w}{\partial t} + \frac{C}{2r} \frac{\partial f\_w}{\partial r} = 0,\tag{16}$$

where f <sup>w</sup> 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

Figure 5. Capillary pressure curve (a) and saturation profile at a time t during water injection (b). The dashed green lines represent the capillary pressure derivative (a) and the saturation gradient (b) at the water front.

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 [1] equation for the injection period, with the following initial and boundary conditions

$$\mathcal{S}\_{w}(r,0) = \mathcal{S}\_{iw\prime} \tag{17}$$

given by the combination of three saturations: SBL

magnifying the dispersion effects in the saturation governing equation; and SSH

<sup>w</sup> ð Þþ <sup>r</sup>; <sup>t</sup> SSL

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

where BL stands for Buckley-Leverett, SL for shock layer, and SH for shock function.

<sup>w</sup> , is obtained by letting e ! 0 in Eq. (15):

<sup>w</sup> ð Þ¼ <sup>r</sup>; <sup>t</sup> lim <sup>e</sup>!<sup>0</sup>,ð Þ <sup>r</sup>;<sup>t</sup> fixed

�

can be obtained by the application of the method of characteristics and is given by [5].

C

<sup>1</sup> � Sor, r<sup>2</sup> <sup>⩽</sup>r<sup>2</sup>

<sup>r</sup><sup>2</sup> � <sup>r</sup><sup>2</sup> w � � t � �, r<sup>2</sup>

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-

<sup>D</sup> <sup>¼</sup> <sup>C</sup> <sup>f</sup> <sup>w</sup> Swf � � � <sup>f</sup> <sup>w</sup>ð Þ Siw � �

where Swf and Siw 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

> <sup>¼</sup> df <sup>w</sup> Swf � � dSw

Siw, r<sup>2</sup> <sup>&</sup>gt; Dt <sup>þ</sup> <sup>r</sup><sup>2</sup>

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

> 1 � Sor, for r ≤ rw Siw, for r > rw,

> > w

<sup>w</sup> < r<sup>2</sup> ⩽r<sup>2</sup>

<sup>w</sup> þ Dt

w,

Swf � Siw � � , (24)

: (25)

<sup>w</sup> ð Þ� <sup>r</sup>; <sup>t</sup> <sup>S</sup>SH

Swð Þ <sup>r</sup>; <sup>t</sup> <sup>≃</sup>SBL

SBL

SBL <sup>w</sup> ð Þ¼ r; 0

8 >>><

>>>:

df <sup>w</sup> dSw � ��<sup>1</sup> 1

shock has a constant speed, satisfying the Rankine-Hugoniot condition [25]:

<sup>f</sup> <sup>w</sup> Swf � � � <sup>f</sup> <sup>w</sup>ð Þ Siw Swf � Siw

pressure effects are neglected; SSL

2.1.1. Outer solution (SBL

The outer solution, SBL

initial condition

flow curve, i.e.,

wave represented by a Heaviside function:

w )

SBL <sup>w</sup> ð Þ¼ r; t <sup>w</sup> , the solution obtained when the capillary

<sup>w</sup> ð Þ r; t , (20)

http://dx.doi.org/10.5772/intechopen.76828

Swð Þ r; t; e : (21)

<sup>w</sup> , the shock

125

(22)

(23)

<sup>w</sup> , the saturation distribution in the shock layer obtained by

$$F\_w(r\_w, t) = 1,\tag{18}$$

$$\lim\_{r \to \infty} \mathcal{S}\_w = \mathcal{S}\_{iw}.\tag{19}$$

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 (Swf > Siw) 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 Fw ¼ f <sup>w</sup> ¼ 1 is the wellbore, so Swð Þ¼ rw; t 1 � Sor. 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

Figure 6. True saturation distribution in the reservoir (Sw) compared with the outer solution (SBL <sup>w</sup> ) and the inner solution (SSL <sup>w</sup> ). The dashed square shows the saturation transition zone between the outer and the inner solution where none of these two are capable of approximate Sw.

given by the combination of three saturations: SBL <sup>w</sup> , the solution obtained when the capillary pressure effects are neglected; SSL <sup>w</sup> , the saturation distribution in the shock layer obtained by magnifying the dispersion effects in the saturation governing equation; and SSH <sup>w</sup> , the shock wave represented by a Heaviside function:

$$\mathcal{S}\_{w}(r,t) \simeq \mathcal{S}\_{w}^{BL}(r,t) + \mathcal{S}\_{w}^{SL}(r,t) - \mathcal{S}\_{w}^{SH}(r,t),\tag{20}$$

where BL stands for Buckley-Leverett, SL for shock layer, and SH for shock function.

#### 2.1.1. Outer solution (SBL w )

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 [1] equation

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 (Swf > Siw) 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 Fw ¼ f <sup>w</sup> ¼ 1 is the wellbore, so Swð Þ¼ rw; t 1 � Sor. 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

Swð Þ¼ r; 0 Siw, (17)

Fwð Þ¼ rw; t 1, (18)

lim<sup>r</sup>!<sup>∞</sup> Sw <sup>¼</sup> Siw: (19)

<sup>w</sup> ) and the inner solution

for the injection period, with the following initial and boundary conditions

124 Perturbation Methods with Applications in Science and Engineering

Figure 6. True saturation distribution in the reservoir (Sw) compared with the outer solution (SBL

<sup>w</sup> ). The dashed square shows the saturation transition zone between the outer and the inner solution where none of

(SSL

these two are capable of approximate Sw.

The outer solution, SBL <sup>w</sup> , is obtained by letting e ! 0 in Eq. (15):

$$\mathcal{S}\_{w}^{BL}(r,t) = \lim\_{\mathfrak{e} \to 0\_{\ell}\left(r,t\right)} \mathcal{S}\_{w}(r,t,\mathfrak{e}).\tag{21}$$

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

$$S\_w^{BL}(r,0) = \begin{cases} 1 - S\_{or\prime} & \text{for } r \le r\_w \\ S\_{iw\_{\prime\prime}} & \text{for } r > r\_{w\prime} \end{cases} \tag{22}$$

can be obtained by the application of the method of characteristics and is given by [5].

$$S\_w^{BL}(r,t) = \begin{cases} 1 - \mathbb{S}\_{\theta\nu} & r^2 \leqslant r\_w^2\\ \left(\frac{df\_w}{dS\_w}\right)^{-1} \left(\frac{1}{C} \frac{\left(r^2 - r\_w^2\right)}{t}\right), & r\_w^2 < r^2 \leqslant r\_w^2 + Dt\\ \mathbb{S}\_{\text{div}\prime} & r^2 > Dt + r\_{w\prime}^2 \end{cases} \tag{23}$$

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]:

$$D = \mathcal{C} \frac{\left[f\_w(\mathcal{S}\_{w\circ}) - f\_w(\mathcal{S}\_{iw})\right]}{\left[\mathcal{S}\_{w\circ} - \mathcal{S}\_{iw}\right]},\tag{24}$$

where Swf and Siw 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.,

$$\frac{\int\_{w} \left(\mathcal{S}\_{w\f}\right) - f\_w(\mathcal{S}\_{iw})}{\mathcal{S}\_{w\f} - \mathcal{S}\_{iw}} = \frac{df\_w(\mathcal{S}\_{w\f})}{d\mathcal{S}\_w}.\tag{25}$$

The details of this solution can be found in [5]. Figure 7 shows the shock jump slope tangent to the fractional flow curve at Sw ¼ Swf and the saturation distribution in the reservoir at a time t. The rarefaction wave family spans from 1 � Sor to Swf from rw to r ¼ 25 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.

with the same speed as the shock itself (Figure 9). In order to find SSL

where rs is the shock front position:

the resulting equation by e yield

e C ∂Sw

we have defined as the inner region.

<sup>∂</sup><sup>τ</sup> � <sup>D</sup> <sup>∂</sup>Sw ∂w þ ∂f <sup>w</sup> <sup>∂</sup><sup>w</sup> � <sup>∂</sup> ∂w

The inner solution is obtained by letting e ! 0 in Eq. (28):

�<sup>D</sup> <sup>∂</sup>SSL w ∂w þ ∂f <sup>w</sup> <sup>∂</sup><sup>w</sup> � <sup>∂</sup> ∂w 2r 2 <sup>s</sup> ð Þ<sup>τ</sup> k r<sup>2</sup>

SSL

layer by using a stretching traveling wave coordinate. Similarly, as defined in [24, 18, 19]

w is zero at r ¼ rs and goes to �∞ as e ! 0. We rewrite Eq. (15) in terms of moving coordinates, ð Þ! r; t ð Þ w; τ , where τ ¼ τðÞ¼ t t. Using Eq. (26) in the transformed equation and multiplying

> 2 ew þ r 2 <sup>s</sup> ð Þ<sup>τ</sup> <sup>k</sup> <sup>e</sup><sup>w</sup> <sup>þ</sup> <sup>r</sup>

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, eð Þ r ! rsð Þτ . 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

with capillary pressure. Both profiles agree in the region around the water front, i.e., in the shock boundary layer which

<sup>r</sup><sup>2</sup> � <sup>r</sup><sup>2</sup> <sup>s</sup> ð Þt

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

w ¼ w rð Þ¼ ; t

r 2 <sup>s</sup> ðÞ¼ t r 2

<sup>w</sup> ð Þ¼ <sup>w</sup>; <sup>τ</sup> lim <sup>e</sup>!<sup>0</sup>,ð Þ <sup>w</sup>;<sup>τ</sup> fixed

Figure 9. Saturation distribution during the injection period for the inner solution (SSL

as presented in [26]. Therefore, neglecting the terms of order e in Eq. (28), we have

<sup>w</sup> , we magnify the shock

127

¼ 0: (28)

<sup>e</sup> , (26)

http://dx.doi.org/10.5772/intechopen.76828

<sup>w</sup> þ CDt, (27)

∂w

<sup>s</sup> ð Þ<sup>τ</sup> ; <sup>τ</sup>; <sup>e</sup> , (29)

¼ 0: (30)

<sup>w</sup> ) and for the true solution (Sw),

2 <sup>s</sup> ð Þ<sup>τ</sup> <sup>∂</sup><sup>Ψ</sup>

2

<sup>s</sup> ð Þ<sup>τ</sup> <sup>∂</sup><sup>Ψ</sup> ∂w

Sw ew þ r

#### 2.1.2. Inner solution (SSL 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 Sw ¼ Swf (a) and the saturation profile in the reservoir at a time t (b). The rarefaction waves family spans from 1 � Sor 1 to Swf and from rw to r = 25 ft, the water front position, i.e., the shock front position, rf ,inj. Ahead of the water front position, there is immobile water.

Figure 8. Saturation distribution during the injection period for the outer solution (SBL <sup>w</sup> ), without capillary pressure, and for the true solution (Sw), with capillary pressure. Both profiles agree in the region far from the water front, the region outside the dashed square.

with the same speed as the shock itself (Figure 9). In order to find SSL <sup>w</sup> , we magnify the shock layer by using a stretching traveling wave coordinate. Similarly, as defined in [24, 18, 19]

$$w = w(r, t) = \frac{r^2 - r\_s^2(t)}{\epsilon},\tag{26}$$

where rs is the shock front position:

The details of this solution can be found in [5]. Figure 7 shows the shock jump slope tangent to the fractional flow curve at Sw ¼ Swf and the saturation distribution in the reservoir at a time t. The rarefaction wave family spans from 1 � Sor to Swf from rw to r ¼ 25 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 solu-

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 Sw ¼ Swf (a) and the saturation profile in the reservoir at a time t (b). The rarefaction waves family spans from 1 � Sor 1 to Swf and from rw to r = 25 ft, the water front position, i.e., the shock front position, rf ,inj. Ahead of the water front position, there is immobile water.

for the true solution (Sw), with capillary pressure. Both profiles agree in the region far from the water front, the region

<sup>w</sup> ), without capillary pressure, and

Figure 8. Saturation distribution during the injection period for the outer solution (SBL

tion obtained from a numerical simulator.

w )

126 Perturbation Methods with Applications in Science and Engineering

2.1.2. Inner solution (SSL

outside the dashed square.

$$r\_s^2(t) = r\_w^2 + \mathsf{C}Dt,\tag{27}$$

w is zero at r ¼ rs and goes to �∞ as e ! 0. We rewrite Eq. (15) in terms of moving coordinates, ð Þ! r; t ð Þ w; τ , where τ ¼ τðÞ¼ t t. Using Eq. (26) in the transformed equation and multiplying the resulting equation by e yield

$$\frac{\partial}{\partial \mathbf{C}} \frac{\partial \mathbf{S}\_w}{\partial \tau} - D \frac{\partial \mathbf{S}\_w}{\partial w} + \frac{\partial f\_w}{\partial w} - \frac{\partial}{\partial w} \left( 2 \left( \epsilon w + r\_s^2(\tau) \right) k \left( \epsilon w + r\_s^2(\tau) \right) \frac{\partial \Psi'}{\partial w} \right) = 0. \tag{28}$$

The inner solution is obtained by letting e ! 0 in Eq. (28):

$$S\_w^{SL}(w, \tau) = \lim\_{\mathfrak{e} \to 0\_\ell(w, \mathfrak{e}) \text{ fixed}} S\_w(\mathfrak{e}w + r\_s^2(\tau), \tau, \mathfrak{e}), \tag{29}$$

as presented in [26]. Therefore, neglecting the terms of order e in Eq. (28), we have

$$-D\frac{\partial S\_w^{SL}}{\partial w} + \frac{\partial f\_w}{\partial w} - \frac{\partial}{\partial w} \left( 2r\_s^2(\tau)k(r\_s^2(\tau))\frac{\partial \Psi'}{\partial w} \right) = 0. \tag{30}$$

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, eð Þ r ! rsð Þτ . 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 (SSL <sup>w</sup> ) and for the true solution (Sw), 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

$$-DS\_w^{SL} + f\_w \left( S\_w^{SL} \right) - 2r\_s^2(\tau) k \left( r\_s^2(\tau) \right) \frac{d\Psi \left( S\_w^{SL} \right)}{dS\_w} \frac{\partial S\_w^{SL}(w, \tau)}{\partial w} = a(\tau), \tag{31}$$

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

$$\begin{cases} w \to \rightsquigarrow : \qquad S\_w^{SL} = S\_{iw\prime} \quad \frac{\partial S\_w^{SL}}{\partial w} = 0, \\ w \to -\rightsquigarrow : \quad S\_w^{SL} = S\_{w\prime\prime} \quad \frac{\partial S\_w^{SL}}{\partial w} = 0, \end{cases} \tag{32}$$

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

$$a(\pi) = -DS\_{iw\nu} \tag{33}$$

SSL

value of SSL

the mass balance.

technique [26]:

SSH

(SSH <sup>w</sup> ) (b). SSH

2.1.3. Matching solution (SSH

The matching saturation SSH

and in the injection case is given by

and matching saturation solutions.

Figure 10. The matching saturation function (SSH

2.1.4. Mass balance

<sup>w</sup> )

SSH <sup>w</sup> ð Þ¼ r; t

lim <sup>r</sup>2!r<sup>2</sup> <sup>w</sup>þCDt SBL

(

<sup>w</sup> ¼ Siw, 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

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

used, but first let us present the matching saturation, since this solution will be necessary for

<sup>w</sup> ð Þ¼ r; t lim

Siw, r<sup>2</sup> ≥ r<sup>2</sup>

which is plotted in Figure 10 against the outer and inner solutions. As we were searching for,

<sup>w</sup> 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,

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

<sup>w</sup> ) compared with the outer solution (SBL

<sup>w</sup> matches with the outer solution in the inner region and with the inner solution in the outer region, as desired.

<sup>w</sup> ≤ r<sup>2</sup> <sup>s</sup> ð Þt :

Swf , r<sup>2</sup>

<sup>w</sup>!�<sup>∞</sup> SSL

<sup>s</sup> ðÞ¼ <sup>t</sup> <sup>r</sup><sup>2</sup>

<sup>w</sup> at w well. In order to find a closed form for this problem, mass balance can be

<sup>w</sup> is defined using the matching principle by applying Prandtl's

<sup>w</sup> þ CDt,

<sup>w</sup> ð Þ w; t ; (37)

http://dx.doi.org/10.5772/intechopen.76828

<sup>w</sup> ) (a) and with the inner solution

(38)

129

while using the second boundary condition given by Eq. (32) yields

$$a(\pi) = -DS\_{w\!\!f} + f\_w(\mathbb{S}\_{w\!\!f});\tag{34}$$

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 Siw � <sup>S</sup>SL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> <sup>S</sup>SL w � � yield

$$2r\_s^2(\tau)k\left(r\_s^2(\tau)\right)\frac{\frac{d\Psi\left(S\_w^{\rm SL}\right)}{dS\_w}}{D\left(S\_{iw}-S\_w^{\rm SL}\right)+f\_w\left(S\_w^{\rm SL}\right)}\frac{\partial S\_w^{\rm SL}(w,\tau)}{\partial w} = 1.\tag{35}$$

Integrating Eq. (35) from wwell ¼ w rð Þ <sup>w</sup>; τ to any w at any time τ gives us the relationship between any SSL <sup>w</sup> and w:

$$2r\_s^2 k(r\_s^2) \int\_{S\_w^{\rm SL}(w\_{\rm mol})}^{S\_w^{\rm SL}} \frac{\frac{d\Psi\left(S\_w^{\rm SL}\right)}{dS\_w}}{D\left(\mathcal{S}\_{\rm inv} - \mathcal{S}\_w^{\rm SL}\right) + f\_w\left(\mathcal{S}\_w^{\rm SL}\right)} dS\_w^{\rm SL} = \int\_{w\_{\rm red}}^w dw. \tag{36}$$

At SSL <sup>w</sup> ¼ Swf , 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 SSL <sup>w</sup> (Eq. (32)). At SSL <sup>w</sup> ¼ Siw, 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 SSL <sup>w</sup> 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 (SSH <sup>w</sup> )

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

<sup>s</sup> ð Þ<sup>τ</sup> � � <sup>d</sup><sup>Ψ</sup> <sup>S</sup>SL

<sup>w</sup> <sup>¼</sup> Siw, <sup>∂</sup>SSL

<sup>w</sup> <sup>¼</sup> Swf , <sup>∂</sup>SSL

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

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

að Þ¼� τ DSwf þ f <sup>w</sup> Swf

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

> <sup>d</sup><sup>Ψ</sup> SSL ð Þ <sup>w</sup> dSw

Integrating Eq. (35) from wwell ¼ w rð Þ <sup>w</sup>; τ to any w at any time τ gives us the relationship

<sup>d</sup><sup>Ψ</sup> SSL ð Þ <sup>w</sup> dSw

<sup>w</sup> ¼ Swf , the integral in the left side of Eq. (36) diverges as the integrand denominator goes

w � �

> w � � dSSL

<sup>w</sup> ¼ ðw wwell

∂SSL <sup>w</sup> ð Þ w; τ

D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> SSL

D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> SSL

to 0. This behavior is consistent with our boundary condition assumptions for SSL

w � � dSw

> w <sup>∂</sup><sup>w</sup> <sup>¼</sup> <sup>0</sup>,

w <sup>∂</sup><sup>w</sup> <sup>¼</sup> <sup>0</sup>,

að Þ¼� τ DSiw, (33)

� �; (34)

<sup>∂</sup><sup>w</sup> <sup>¼</sup> <sup>1</sup>: (35)

dw: (36)

<sup>w</sup> (Eq. (32)). At

∂SSL <sup>w</sup> ð Þ w; τ

<sup>∂</sup><sup>w</sup> <sup>¼</sup> <sup>a</sup>ð Þ<sup>τ</sup> , (31)

(32)

applying the chain rule gives

by D Siw � <sup>S</sup>SL

between any SSL

At SSL

w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> <sup>S</sup>SL

w � � yield

2r 2 <sup>s</sup> ð Þ<sup>τ</sup> k r<sup>2</sup> <sup>s</sup> ð Þ<sup>τ</sup> � �

<sup>w</sup> and w:

2r 2 <sup>s</sup> k r<sup>2</sup> s � � <sup>ð</sup>SSL w SSL <sup>w</sup> ð Þ wwell

�DSSL

<sup>w</sup> <sup>þ</sup> <sup>f</sup> <sup>w</sup> <sup>S</sup>SL w � � � <sup>2</sup><sup>r</sup>

128 Perturbation Methods with Applications in Science and Engineering

boundary conditions (for the inner solution) given by

8 >><

>>:

while using the second boundary condition given by Eq. (32) yields

2 <sup>s</sup> ð Þ<sup>τ</sup> k r<sup>2</sup>

<sup>w</sup> ! <sup>∞</sup> : <sup>S</sup>SL

<sup>w</sup> ! �<sup>∞</sup> : <sup>S</sup>SL

The matching saturation SSH <sup>w</sup> is defined using the matching principle by applying Prandtl's technique [26]:

$$\lim\_{w^2 \to r\_w^2 + \mathcal{C}Dt} S\_w^{\mathcal{BL}}(r, t) = \lim\_{w \to \pm \infty} S\_w^{\mathcal{SL}}(w, t); \tag{37}$$

and in the injection case is given by

$$\mathcal{S}\_{w}^{SH}(r,t) = \begin{cases} \mathcal{S}\_{\text{inv}} & r^2 \ge r\_s^2(t) = r\_w^2 + \mathsf{CD}t, \\ \mathcal{S}\_{wf} & r\_w^2 \le r\_s^2(t). \end{cases} \tag{38}$$

which is plotted in Figure 10 against the outer and inner solutions. As we were searching for, SSH <sup>w</sup> 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.

#### 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

Figure 10. The matching saturation function (SSH <sup>w</sup> ) compared with the outer solution (SBL <sup>w</sup> ) (a) and with the inner solution (SSH <sup>w</sup> ) (b). SSH <sup>w</sup> matches with the outer solution in the inner region and with the inner solution in the outer region, as desired.

Figure 11. Saturation distribution during the injection period with (true solution) and without capillary pressure (outer solution), the traveling wave (inner solution), and the matching saturation (a). SSH <sup>w</sup> matches with the region inner solution in the outer region, the region far from the water front.

the saturation distribution based in the mass balance. Since both the Buckley-Leverett (SBL w ) solution and the composite solution (Sw) must obey material balance, the following equality

$$\eta\_t t = \int\_{r\_w^2}^{\infty} (S\_w(r, t) - S\_{iw}) \tau h dr^2 = \int\_{r\_w^2}^{\infty} (S\_w^{RL}(r, t) - S\_{iw}) \tau h dr^2. \tag{39}$$

must hold. From Eq. (20) and Eq. (39), it follows that

$$\int\_{r\_w^2}^{\infty} \left( \mathcal{S}\_w^{BL} + \mathcal{S}\_w^{SL} - \mathcal{S}\_w^{SH} - \mathcal{S}\_{iw} \right) \pi h dr^2 = \int\_{r\_w^2}^{\infty} \left( \mathcal{S}\_w^{RL}(r, t) - \mathcal{S}\_{iw} \right) \pi h dr^2,\tag{40}$$

which, upon simplification, gives

$$\int\_{r\_w^2}^{\infty} (\mathcal{S}\_w^{SL} - \mathcal{S}\_w^{SH}) dr^2 = 0. \tag{41}$$

From Eq. (35),

ð∞ �CD<sup>τ</sup> e

the equation yield

be equal which gives

2r 2 <sup>s</sup> k r<sup>2</sup> s � �e ðSiw SSL <sup>w</sup> �CD<sup>τ</sup> ð Þ <sup>e</sup>

Once the value SSL

note that, as SSL

solutions of previous authors.

2r 2 <sup>s</sup> k r<sup>2</sup> s � � <sup>ð</sup>Siw SSL <sup>w</sup> �CD<sup>τ</sup> ð Þ <sup>e</sup>

dw <sup>¼</sup> <sup>2</sup>r<sup>2</sup>

ð∞ �CD<sup>τ</sup> e

> dΨ dSw <sup>S</sup>SL w � �dSSL w

D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> <sup>S</sup>SL

<sup>w</sup> � CD<sup>τ</sup> e

dw ¼ 2r 2 <sup>s</sup> k r<sup>2</sup> s � � <sup>ð</sup>Siw SSL <sup>w</sup> �CD<sup>τ</sup> ð Þ <sup>e</sup>

<sup>s</sup> k r<sup>2</sup> s � � Siw <sup>ð</sup>Siw SSL <sup>w</sup> �CD<sup>τ</sup> ð Þ <sup>e</sup>

dw ¼ 2r 2 <sup>s</sup> k r<sup>2</sup> s � �

dΨ dSw D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup>

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

Substituting Eq. (45) in Eq. (44) and solving the resulting equation divided by eSiw for

dΨ dSw SSL w � �

D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> SSL

Setting Sw ¼ Siw in the upper limits of the integrals of Eq. (36) and exchanging the two sides of

As the left sides of Eqs. (46) and (47) are the same, the right sides of these two equations must

dΨ dSw SSL w � �

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

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 (rwf ) which is the finite position at which water can be considered immobile—and the initial water saturation in the oil zone exhibited by

� � (i.e., the inner solution saturation in the wellbore SSL

<sup>w</sup> should reach Swf and Siw asymptotically as w ! �∞, here we did not have to

D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> SSL

<sup>s</sup> k r<sup>2</sup> s � � Siw <sup>ð</sup>Siw SSL <sup>w</sup> �CD<sup>τ</sup> ð Þ <sup>e</sup>

SSL w

w � � <sup>¼</sup> <sup>2</sup>r<sup>2</sup>

Multiplying Eq. (48) by eSiw and rearranging the resulting equation give

SSL <sup>w</sup> � Siw � � dSSL

w � � dSSL

> SSL <sup>w</sup> dΨ dSw <sup>S</sup>SL w � �dSSL w

D Siw � <sup>S</sup>SL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> <sup>S</sup>SL

w � � dSSL

w � � dSSL

w

dΨ dSw SSL w � �

D Siw � SSL w � � <sup>þ</sup> <sup>f</sup> <sup>w</sup> SSL

<sup>w</sup> : (45)

http://dx.doi.org/10.5772/intechopen.76828

<sup>w</sup> � Swf � Siw � �CD<sup>τ</sup> eSiw

: (46)

131

<sup>w</sup> : (47)

� � � Swf � Siw � �CD<sup>τ</sup> eSiw

<sup>w</sup> <sup>¼</sup> Swf � Siw � �CDτ: (49)

:

(48)

<sup>w</sup> ð Þ wwell ) is

Rearranging Eq. (41) using Eq. (38) for SSH <sup>w</sup> gives

$$\int\_{r\_w^2}^{\prime\prime} S\_w^{\rm SI} dr^2 = \int\_{r\_w^2}^{r\_s^2} S\_{w\dagger} dr^2 + \int\_{r\_s^2}^{\prime\prime} S\_{iw} dr^2 = S\_{w\dagger} \left( r\_s^2 - r\_w^2 \right) + \left( S\_{iw} \int\_{r\_w^2}^{\prime\prime} dr^2 - S\_{iw} \int\_{r\_w^2}^{r\_s^2} dr^2 \right). \tag{42}$$

Using Eq. (27) in Eq. (42), it follows that

$$\int\_{r\_w^2}^{\infty} S\_w^{SL} dr^2 = \left(\mathcal{S}\_{w\!f} - \mathcal{S}\_{iw}\right) \mathcal{C} Dt + \mathcal{S}\_{iw} \int\_{r\_w^2}^{\infty} dr^2. \tag{43}$$

Transforming Eq. (43) from ð Þ! r; t ð Þ w; τ and using Eq. (26), Eq. (43) becomes

$$
\epsilon \int\_{-\frac{\triangle\pi}{\bullet}}^{\infty} S\_w^{\rm SL}(w) dw = \left(\mathcal{S}\_{w\circ} - \mathcal{S}\_{iw}\right) \text{CD}\tau + \epsilon \mathcal{S}\_{iw} \int\_{-\frac{\triangle\pi}{\bullet}}^{\infty} d\upsilon. \tag{44}
$$

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation… http://dx.doi.org/10.5772/intechopen.76828 131

From Eq. (35),

the saturation distribution based in the mass balance. Since both the Buckley-Leverett (SBL

Figure 11. Saturation distribution during the injection period with (true solution) and without capillary pressure (outer

solution and the composite solution (Sw) must obey material balance, the following equality

ð∞ r2 w SBL

> ð∞ r2 w SBL

<sup>w</sup> ð Þ� r; t Siw � �πhdr<sup>2</sup>

<sup>w</sup> ð Þ� r; t Siw � �πhdr<sup>2</sup>

� �dr<sup>2</sup> <sup>¼</sup> <sup>0</sup>: (41)

ð∞ r2 w

ð∞ r2 w dr<sup>2</sup>

> ð∞ �CD<sup>τ</sup> e

dr<sup>2</sup> � Siw

!

ðr2 s r2 w dr<sup>2</sup>

: (43)

dw: (44)

ð Þ Swð Þ� <sup>r</sup>; <sup>t</sup> Siw <sup>π</sup>hdr<sup>2</sup> <sup>¼</sup>

<sup>w</sup> � Siw � �πhdr<sup>2</sup> <sup>¼</sup>

> ð∞ r2 w SSL <sup>w</sup> � SSH w

ð∞ r2 s

<sup>w</sup> gives

2 <sup>s</sup> � r 2 w � � <sup>þ</sup> Siw

� �CDt <sup>þ</sup> Siw

� �CD<sup>τ</sup> <sup>þ</sup> <sup>e</sup>Siw

Siwdr<sup>2</sup> <sup>¼</sup> Swf <sup>r</sup>

<sup>w</sup> dr<sup>2</sup> <sup>¼</sup> Swf � Siw

Transforming Eq. (43) from ð Þ! r; t ð Þ w; τ and using Eq. (26), Eq. (43) becomes

<sup>w</sup> ð Þ w dw ¼ Swf � Siw

qt t ¼ ð∞ r2 w

in the outer region, the region far from the water front.

130 Perturbation Methods with Applications in Science and Engineering

ð∞ r2 w SBL <sup>w</sup> <sup>þ</sup> <sup>S</sup>SL

which, upon simplification, gives

ð∞ r2 w SSL <sup>w</sup> dr<sup>2</sup> <sup>¼</sup>

Rearranging Eq. (41) using Eq. (38) for SSH

ðr2 s r2 w

Using Eq. (27) in Eq. (42), it follows that

e ð∞ �CD<sup>τ</sup> e SSL

Swf dr<sup>2</sup> <sup>þ</sup>

ð∞ r2 w SSL

must hold. From Eq. (20) and Eq. (39), it follows that

<sup>w</sup> � SSH

solution), the traveling wave (inner solution), and the matching saturation (a). SSH

w )

: (39)

<sup>w</sup> matches with the region inner solution

, (40)

: (42)

$$dw = 2r\_s^2 k \left( r\_s^2 \right) \frac{\frac{d\Psi}{dS\_w}}{D \left( S\_{iw} - S\_w^{SL} \right) + f\_w} dS\_w^{SL} \,. \tag{45}$$

Substituting Eq. (45) in Eq. (44) and solving the resulting equation divided by eSiw for

$$\int\_{-\varDelta\_{\rm tr}}^{\varDelta} dw = \frac{2r\_s^2 k \left(r\_s^2\right)}{S\_{\rm inv}} \int\_{S\_w^{\rm SL}\left(-\varDelta\_{\rm tr}\right)}^{S\_{\rm inv}} S\_w^{\rm SL} \frac{d^{\rm SU}\left(S\_w^{\rm SL}\right)}{D\left(S\_{\rm inv} - S\_w^{\rm SL}\right) + f\_w\left(S\_w^{\rm SL}\right)} dS\_w^{\rm SL} - \frac{\left(S\_{\rm wf} - S\_{\rm inv}\right) \mathsf{C} \mathsf{D}\tau}{\mathsf{e} S\_{\rm inv}}.\tag{46}$$

Setting Sw ¼ Siw in the upper limits of the integrals of Eq. (36) and exchanging the two sides of the equation yield

$$\int\_{-\frac{\Box\eta}{\epsilon}}^{\prime} dw = 2r\_s^2 k(r\_s^2) \int\_{S\_w^{\rm SL}\left(-\frac{\Box\tau}{\epsilon}\right)}^{S\_{\rm inv}} \frac{\frac{d\Psi}{dS\_w} \left(S\_w^{\rm SL}\right)}{D\left(S\_{\rm inv} - S\_w^{\rm SL}\right) + f\_w\left(S\_w^{\rm SL}\right)} dS\_w^{\rm SL}.\tag{47}$$

As the left sides of Eqs. (46) and (47) are the same, the right sides of these two equations must be equal which gives

$$2\sigma\_s^2 k(r\_s^2) \int\_{S\_w^2}^{\mathbb{S}\_{\text{int}}} \frac{\frac{d\mathcal{V}}{d\mathcal{S}\_w} \left(\mathcal{S}\_w^{\text{SL}}\right) dS\_w^{\text{SL}}}{D\left(\mathcal{S}\_{\text{in}} - \mathcal{S}\_w^{\text{SL}}\right) + f\_w\left(\mathcal{S}\_w^{\text{SL}}\right)} = \frac{2\sigma\_s^2 k(r\_s^2)}{\mathcal{S}\_{\text{in}}} \int\_{S\_w^2 \left(-\frac{\varpi\_s}{\bullet}\right)}^{\mathbb{S}\_{\text{in}}} \frac{\frac{\mathcal{S}\_w^{\text{SL}} d\mathcal{V}}{d\mathcal{S}\_w} \left(\mathcal{S}\_w^{\text{SL}}\right) dS\_w^{\text{SL}}}{D\left(\mathcal{S}\_{\text{in}} - \mathcal{S}\_w^{\text{SL}}\right) + f\_w\left(\mathcal{S}\_w^{\text{SL}}\right)} - \frac{\left(\mathcal{S}\_{wf} - \mathcal{S}\_{\text{in}}\right) \mathsf{C} \mathcal{D} \mathsf{r}}{\mathfrak{C} \mathcal{S}\_{\text{in}}}.\tag{48}$$

Multiplying Eq. (48) by eSiw and rearranging the resulting equation give

$$\left(2\tau\_s^2 k(r\_s^2)\epsilon\right)\_{S\_w^{\rm SL}\left(-\frac{\rm CDr}{\bullet}\right)}^{\rm S\_{\rm in}} \left(S\_w^{\rm SL} - S\_{\rm in}\right) \frac{\frac{d\Psi}{d\mathcal{S}\_w} \left(S\_w^{\rm SL}\right)}{D\left(S\_{\rm in} - S\_w^{\rm SL}\right) + f\_w\left(S\_w^{\rm SL}\right)} dS\_w^{\rm SL} = \left(S\_{w\mathcal{f}} - S\_{\rm in}\right) \mathcal{CD}\tau. \tag{49}$$

Once the value SSL <sup>w</sup> � CD<sup>τ</sup> e � � (i.e., the inner solution saturation in the wellbore SSL <sup>w</sup> ð Þ wwell ) 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 SSL <sup>w</sup> should reach Swf and Siw asymptotically as w ! �∞, 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 (rwf ) 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.

Figure 12. Comparison of the saturation distribution from analytical solution and IMEX during the injection period with capillary pressure.

#### 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 [2]. During injection at a constant flow rate, qt ð Þ rw; t RB/D, where t ¼ 0 at the beginning of the water injection, by integrating Darcy's law as in [6, 2], given by

$$q\_t = -\frac{k(r)hr}{a} \left(\lambda\_t \frac{\partial p\_o}{\partial r} - \lambda\_w \frac{\partial p\_c}{\partial r}\right),\tag{50}$$

ΔpwfðÞ¼ t

ΔpwfðÞ¼ t

compressibility is

3. Validation

αqt ð Þ rw; t h

> α h ð ∞

solution can be approximated as

rw

qt ð Þ r; t λbo

<sup>¼</sup> <sup>Δ</sup>bpoð Þþ <sup>t</sup>

ðrwfð Þt

1 <sup>λ</sup><sup>t</sup> <sup>r</sup>;Δtprod � �k rð Þ dr r þ α h ð ∞

h Ð rwf ð Þt rw qt ð Þ rw;t b λok rð Þ dr

λbo <sup>λ</sup>tð Þ <sup>r</sup>; <sup>t</sup> � <sup>1</sup> ! <sup>1</sup>

khλb<sup>o</sup>

<sup>Δ</sup>bpoð Þ<sup>t</sup> 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

> 1 2

Here, β is a unit conversion factor in which oil field unit is 0.0002637 and the single-phase total

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, 10�<sup>7</sup> 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.

rwfð Þt

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

ðrwfð Þt

rw

where for any practical set of values of physical properties [28] indicate that this assumption is

qt ð Þ r; t λtð Þ r; t k rð Þ

1 <sup>λ</sup>tð Þ <sup>r</sup>; <sup>t</sup> � <sup>1</sup>

k rð Þ dr r �

ln <sup>β</sup>kλbot <sup>ϕ</sup>bctor<sup>2</sup> w

!

" #

<sup>b</sup>cto <sup>¼</sup> coð Þþ <sup>1</sup> � Swi cwSwi <sup>þ</sup> cr: (56)

dr r �

<sup>r</sup> , where <sup>λ</sup>b<sup>o</sup> <sup>¼</sup> kro ð Þ Swi

λbo

ðrwf tð Þ

rw f w dpc dSw

þ 0:4045 þ s

! dr

ðrwf tð Þ

http://dx.doi.org/10.5772/intechopen.76828

∂Sw ∂r

is the endpoint oil

∂Sw ∂r dr

dr: (54)

: (55)

dr, (53)

133

rw f w dpc dSw

μo

k rð Þr � ðrwf tð Þ

rw f w dpc dSw

∂Sw ∂r

rw

valid. Adding and subtracting the term <sup>α</sup>

mobility at Sw ¼ Swi, Eq. (53) can be rewritten as

ð Þ <sup>r</sup>; <sup>t</sup> k rð Þ dr

<sup>Δ</sup>bpoðÞ¼ <sup>t</sup> pwf , <sup>o</sup>ð Þ� <sup>t</sup> pi <sup>¼</sup> <sup>α</sup>qt

αqt ð Þ rw; t hλb<sup>o</sup>

r þ αqt ð Þ rw; t h

ðrwfð Þt

rw

where pw ¼ po � pc. 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 (poi) can then be expressed as

$$
\Delta p\_{w\not\!f}(t) = p\_{w\not\!f}(t) - p\_{o\not\!i} = \int\_{r\_w}^{\sim} \frac{aq\_t(r,t)}{h\lambda\_t(r,t)k(r)} \frac{dr}{r} - \int\_{r\_w}^{\sim} f\_w \frac{dp\_c}{dS\_w} \frac{\partial S\_w}{\partial r} dr,\tag{51}
$$

where it is assumed that poð Þ¼ rw; t pwð Þ rw; t , i.e., pc ¼ 0 at r ¼ rw, in order to satisfy the compatibility condition [27], i.e., to guarantee phase pressure continuity at the wellbore. Eq. (51) can be rewritten as

$$
\Delta p\_{wf}(t) = \int\_{r\_w}^{\infty} \frac{\alpha q\_t(r, t)}{h \lambda\_t(r, t) k(r)} \frac{dr}{r} - \int\_{r\_w}^{r\_{wf(t)}} f\_w \frac{dp\_c}{dS\_w} \frac{\partial S\_w}{\partial r} dr,\tag{52}
$$

by assuming that the second term in the right-hand side of Eq. (52) is zero from rwf to ∞, considering <sup>f</sup> <sup>w</sup>ð Þ¼ Siw 0 and <sup>∂</sup>Sw <sup>∂</sup><sup>r</sup> ¼ 0 for r > rwfð Þt , since the water in the region ahead of the water front foot is assumed immobile. rwf can be defined as the position at which ð Þ Sw � Siw < δ, where δ is a very small number. For the hypodispersion phenomenon, we can find a finite rwf where δ ! 0, i.e., at which Sw ¼ Siw. Using the [6] steady-state theory, which assumes that, qt ð Þ¼ r; t qt ð Þ rw; t , for r ≤ rwfð Þt , Eq. (52) becomes

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation… http://dx.doi.org/10.5772/intechopen.76828 133

$$\Delta p\_{w\f}(t) = \frac{a q\_t(r\_w, t)}{h} \int\_{r\_w}^{r\_{w\f}(t)} \frac{1}{\lambda\_t(r, \Delta t\_{prod}) k(r)} \frac{dr}{r} + \frac{a}{h} \int\_{r\_{w\f}(t)}^{r} \frac{q\_t(r, t)}{\lambda\_t(r, t) k(r)} \frac{dr}{r} - \int\_{r\_w}^{r\_{w\f}(t)} f\_w \frac{dp\_c}{dS\_w} \frac{\partial \mathbf{S}\_w}{\partial r} dr,\tag{53}$$

where for any practical set of values of physical properties [28] indicate that this assumption is valid. Adding and subtracting the term <sup>α</sup> h Ð rwf ð Þt rw qt ð Þ rw;t b λok rð Þ dr <sup>r</sup> , where <sup>λ</sup>b<sup>o</sup> <sup>¼</sup> kro ð Þ Swi μo is the endpoint oil mobility at Sw ¼ Swi, Eq. (53) can be rewritten as

$$\Delta p\_{w\uparrow}(t) = \left. \frac{a}{\hbar} \right] \frac{\overset{\text{or}}{\underset{r\_w}{\rightleftharpoons}} q\_t(r,t)k(r)}{\widehat{\lambda}\_o} \left( r, t \middle| k(r) \right) \frac{dr}{r} + \left. \frac{aq\_t(r\_w, t)}{\hbar} \right\rbrack \left( \frac{1}{\lambda\_t(r,t)} - \frac{1}{\widehat{\lambda}\_o} \right) \frac{dr}{k(r)r} - \left. \int\_{r\_w} f\_w \frac{dp\_c}{dS\_w} \frac{\partial \mathcal{S}\_w}{\partial r} dr \right|$$

$$= \Delta \widehat{p}\_o(t) + \left. \frac{aq\_t(r\_w, t)}{\hbar \widehat{\lambda}\_o} \right\rbrack\_{r\_w}^{r\_w(t)} \left( \frac{\widehat{\lambda}\_o}{\lambda\_t(r,t)} - 1 \right) \frac{1}{k(r)} \frac{dr}{r} - \left. \int\_{r\_w} f\_w \frac{dp\_c}{dS\_w} \frac{\partial \mathcal{S}\_w}{\partial r} dr \right. \tag{54}$$

<sup>Δ</sup>bpoð Þ<sup>t</sup> 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

$$
\Delta\widehat{p}\_o(t) = p\_{wf,o}(t) - p\_i = \frac{\alpha q\_t}{k \hbar \widehat{\lambda}\_o} \left[ \frac{1}{2} \ln \left( \frac{\beta k \widehat{\lambda}\_o t}{\phi \widehat{c}\_{to} r\_w^2} \right) + 0.4045 + s \right]. \tag{55}
$$

Here, β is a unit conversion factor in which oil field unit is 0.0002637 and the single-phase total compressibility is

$$
\widehat{\mathcal{L}}\_{tb} = \mathfrak{c}\_o (1 - \mathcal{S}\_{wi}) + \mathfrak{c}\_w \mathcal{S}\_{wi} + \mathfrak{c}\_r. \tag{56}
$$

## 3. Validation

2.2. Wellbore pressure

Darcy's law as in [6, 2], given by

132 Perturbation Methods with Applications in Science and Engineering

Eq. (51) can be rewritten as

considering <sup>f</sup> <sup>w</sup>ð Þ¼ Siw 0 and <sup>∂</sup>Sw

ð Þ¼ r; t qt

assumes that, qt

flow rate, qt

capillary 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

Figure 12. Comparison of the saturation distribution from analytical solution and IMEX during the injection period with

λt ∂po <sup>∂</sup><sup>r</sup> � <sup>λ</sup><sup>w</sup>

where pw ¼ po � pc. 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

> αqt ð Þ r; t hλtð Þ r; t k rð Þ

where it is assumed that poð Þ¼ rw; t pwð Þ rw; t , i.e., pc ¼ 0 at r ¼ rw, in order to satisfy the compatibility condition [27], i.e., to guarantee phase pressure continuity at the wellbore.

> dr r �

by assuming that the second term in the right-hand side of Eq. (52) is zero from rwf to ∞,

water front foot is assumed immobile. rwf can be defined as the position at which ð Þ Sw � Siw < δ, where δ is a very small number. For the hypodispersion phenomenon, we can find a finite rwf where δ ! 0, i.e., at which Sw ¼ Siw. Using the [6] steady-state theory, which

ð Þ rw; t , for r ≤ rwfð Þt , Eq. (52) becomes

ðrwf tð Þ

rw f w dpc dSw

ð ∞

rw

qt ¼ � k rð Þhr α

difference from the reservoir initial pressure (poi) can then be expressed as

ΔpwfðÞ¼ t pwfðÞ� t poi ¼

ΔpwfðÞ¼ t

ð ∞

αqt ð Þ r; t hλtð Þ r; t k rð Þ

rw

ð Þ rw; t RB/D, where t ¼ 0 at the beginning of the water injection, by integrating

� �

∂pc ∂r

dr r � ð ∞

rw f w dpc dSw

> ∂Sw ∂r

<sup>∂</sup><sup>r</sup> ¼ 0 for r > rwfð Þt , since the water in the region ahead of the

∂Sw ∂r

, (50)

dr, (51)

dr, (52)

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, 10�<sup>7</sup> 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.


4. Conclusions

Acknowledgements

Ministry of Education.

Conflict of interest

Nomenclature

λ<sup>o</sup> Oil mobility (1/cp) λ<sup>t</sup> Total mobility (1/cp) λ<sup>w</sup> Water mobility (1/cp)

μ<sup>o</sup> Oil viscosity (cp)

μ<sup>w</sup> Water viscosity (cp)

Fo Oil fractional flow

Fw Water fractional flow

f <sup>w</sup> Water mobility ratio

The authors declare that there is no conflict of interest.

β Unit conversion factor (0.0002637)

Δpo Single-phase oil pressure drop (psi)

aw Water endpoint relative permeability

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.

Application of the Method of Matched Asymptotic Expansions to Solve a Nonlinear Pseudo-Parabolic Equation…

http://dx.doi.org/10.5772/intechopen.76828

135

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

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).
