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

Abstract

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.

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 ($k_s$), and the water endpoint relative permeability ($a_w$). The skin zone permeability is the rock permeability in the zone around the well which was stimulated or damaged during the wellbore drilling operation, while the water endpoint permeability is a measure of how easy water can flow in a specific porous media when there is immobile oil present. In the pursuance of modeling the wellbore pressure response during a water injection test, the Rapport-Leas equation [1], a nonlinear pseudo-parabolic convection-dispersion equation, is used to determine the water saturation distribution in the reservoir as a function of time by assuming a one-dimensional homogeneous medium containing incompressible fluids. Water saturation ($S_w$) is the fraction of water in a given pore space, and it is expressed in water volume by pore volume. In [2, 3, 4], it has been shown how to obtain the wellbore pressure response for the case when capillary 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 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 [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].

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

Figure 2.

Relationship between the water saturation ($S_w$), solid blue curve, the dimensionless total flow rate profile ($q_D$), and dotted red curve, during injection.

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.

2.1. Saturation

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

$$\frac{\partial S_w}{\partial t}+\frac{\theta q_t}{2\mathit{\text{πrhϕ}}}\frac{\partial F_w\left(S_w\right)}{\partial r}=0,$$E1

where throughout we assume that porosity ($\varphi$) is homogeneous; $q_t$ is the total liquid rate in RB/D; $\theta$ represents in general a unit conversion factor where in the oil field units used here, $\theta$ =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\left(r\right)h{\lambda }_p\left(S_w\right)}{\alpha }\left(r\frac{\partial p_p}{\partial r}\right),\text{for}p=o,w.$$E2

For field units used throughout, $\alpha$ = 141.2. $p_p$, is the phase p pressure. The ${\lambda }_p$ is the phase p mobility, given by the ratio of the phase permeability ($k_{\mathit{rw}}$ or $k_{\mathit{ro}}$), which are functions of the water saturation, by the phase viscosity (${\mu }_w$ or ${\mu }_o$). To find the water fractional flow ($F_w$), we can subtract Eq. (2) for water from Eq. (2) for oil, to get

$$\frac{\alpha q_o}{k\left(r\right)h{\lambda }_o\left(S_w\right)}-\frac{\alpha q_w}{k\left(r\right)h{\lambda }_w\left(S_w\right)}=-\left(r\frac{\partial p_o}{\partial r}-r\frac{\partial p_w}{\partial r}\right).$$E3

Rearranging Eq. (3), substituting the capillary pressure p_c given by the difference of the oil pressure ($p_o$) and the water pressure ($p_w$), and dividing the resulting equation by the total flow rate, $q_t$, yield

$$F_o-\frac{F_w{\lambda }_o\left(S_w\right)}{{\lambda }_w\left(S_w\right)}=-\frac{k\left(r\right)h{\lambda }_o\left(S_w\right)}{\alpha q_t}\left(r\frac{\partial p_c\left(S_w\right)}{\partial r}\right),$$E4

where $F_o$ and $F_w$ are the oil and water fractional flow given by $q_o$ and $q_w$, respectively. We assume throughout that water is the wetting phase. Finally, substituting $F_o=1-F_w$ in Eq. (4) and solving for $F_w$, we have the following expression for the water fractional flow including capillary pressure effects

$$F_w\left(S_w\right)=\frac{1+\frac{k\left(r\right)h{\lambda }_o\left(S_w\right)}{\alpha q_t}\left(r\frac{\partial p_c}{\partial r}\right)}{1+\frac{{\lambda }_o\left(S_w\right)}{{\lambda }_w\left(S_w\right)}}=\frac{1}{1+\frac{{\lambda }_o\left(S_w\right)}{{\lambda }_w\left(S_w\right)}}+\frac{\frac{k\left(r\right){\mathit{hk}}_{\mathit{ro}}}{\alpha q_t{\mu }_o}\left(r\frac{\partial p_c}{\partial r}\right)}{1+\frac{{\lambda }_o\left(S_w\right)}{{\lambda }_w\left(S_w\right)}}=f_w+\mathrm{ϵ}\mathit{rk}\left(r\right)f_w\left(S_w\right)k_{\mathit{ro}}\left(S_w\right)\frac{\partial p_c\left(S_w\right)}{\partial r},$$E5

where $f_w$ is the water mobility ratio (Figure 3), i.e., the ratio of water mobility and the total mobility (${\lambda }_t$), given by

$$f_w\left(S_w\right)=\frac{1}{1+\frac{{\lambda }_o\left(S_w\right)}{{\lambda }_w\left(S_w\right)}}=\frac{{\lambda }_w\left(S_w\right)}{{\lambda }_o\left(S_w\right)+{\lambda }_w\left(S_w\right)}\frac{{\lambda }_w\left(S_w\right)}{{\lambda }_o\left(S_w\right)+{\lambda }_w\left(S_w\right)},$$E6

which usually assumes an S-shape. $\mathrm{ϵ}$ is the perturbation parameter, defined by

$$\mathrm{ϵ}=\frac{h}{\alpha q_t{\mu }_o}$$E7

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

$$k\left(r\right)=\left\{\begin{array}{ll}{k}_{\mathrm{s}},& r_w\le r<r_{\text{skin}}\\ k,& r\ge r_{\text{skin}},\end{array}\right.$$E8

where $r_w$ is the wellbore radius, $r_{\text{skin}}$ is the radius of the damaged zone, and $k_s$ is the permeability in the skin zone. Grouping all the parameters that are the function of water saturation, we can define (Figure 4)

$$\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}\left(S_w\right)=-f_w\left(S_w\right)k_{\mathit{ro}}\left(S_w\right)\frac{{\mathit{dp}}_c}{{\mathit{dS}}_w}\left(S_w\right),$$E9

and rewrite Eq. (10) as

$$F_w\left(S_w\right)=f_w\left(S_w\right)-\mathrm{ϵ}\mathit{rk}\left(r\right)\frac{d\mathrm{\Psi }\left(S_w\right)}{{\mathit{dS}}_w}\frac{\partial S_w}{\partial r}.$$E10

Figure 3.

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

Figure 4.

$\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}$ versus water saturation for a partially water wet reservoir (b). For a strong water, we reservoir, $\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}\to \infty$ at $S_w=S_{\mathit{iw}}$.

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

$$p_c\left(S_w\right)=p_t{\left(\mathit{sv}+\frac{S_w-S_{\mathit{wi}}}{1-S_{\mathit{wi}}-S_{\mathit{\text{or}}}}\right)}^{-\frac{1}{\lambda }},$$E11

to represent capillary pressure. Here, $S_{\mathit{iw}}$ is the immobile water saturation and $S_{\mathit{\text{or}}}$ is the residual oil saturation. $\lambda$, where $0.4\le \lambda \le 4.0$, is a measure of the pore size distribution (the greater the $\lambda$ value, the more uniform is the pore size distribution), and $p_t$ is the threshold pressure. The threshold pressure is a measure of the maximum pore size [9], i.e., the minimum capillary pressure at which a continuous nonwetting phase exists in the imbibition case and a continuous wetting phase exists in the drainage case [10]. The greater is the maximum pore size, the smaller is the pressure threshold. According to [11], the extrapolation of the capillary pressure curve obtained from experimental data to $S_w=1$ yields the correct threshold value. In practice, we introduce a small variable, sv, to limit the maximum value of $p_c$ to a finite value. We can relate the relative permeabilities and the capillary pressure through $\lambda$ by using the [12] model for the water phase (wetting phase)

$$k_{\mathit{rw}}=a_w{\left(\frac{S_w-S_{\mathit{iw}}}{1-S_{\mathit{iw}}-S_{\mathit{\text{or}}}}\right)}^{\frac{2+\lambda }{\lambda }},$$E12

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

$$k_{\mathit{ro}}={\left(1-\frac{S_w-S_{\mathit{iw}}}{1-S_{\mathit{iw}}-S_{\mathit{\text{or}}}}\right)}^2\left(1-{\left(\frac{S_w-S_{\mathit{iw}}}{1-S_{\mathit{iw}}-S_{\mathit{\text{or}}}}\right)}^{\frac{2+\lambda }{\lambda }}\right).$$E13

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

$$C=\frac{\theta q_t}{\mathit{\pi h\varphi }},$$E14

yields

$$\frac{\partial S_w}{\partial t}+\frac{C}{2r}\frac{\partial f_w}{\partial r}-\mathrm{ϵ}\frac{C}{2r}\frac{\partial }{\partial r}\left(\mathit{rk}\left(r\right)\frac{\mathrm{\partial \Psi }}{\partial r}\right)=0,$$E15

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,$$E16

where $f_w$ is considered to be an S-shaped function along this chapter. During injection, for a partially water wet reservoir, the capillary pressure dispersive effect will be non-negligible only in a small region around the water front (hypodispersion phenomenon) [14, 15] where the capillary pressure derivative and the saturation gradient are significant (Figure 5). The capillary pressure smears the water front during injection balancing the self-sharpening tendency of the shock. [16] have developed an exact analytical solution for linear waterflood including the effects of capillary pressure, but their solution is limited to a particular functional form to represent relative permeabilities and capillary pressure curves and does not consider radial flow, which makes their solution very restrictive. As done by [17, 18] for Cartesian coordinates and by [19] for streamlines and streamtubes, the perturbation caused by the capillary pressure effects can be modeled as a shock layer (water front) which moves with the same speed as the shock wave. By applying the method of matched asymptotic expansions [20, 21], we can 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

$$S_w\left(r,0\right)=S_{\mathit{iw}},$$E17

$$F_w\left(r_w,t\right)=1,$$E18

$$\underset{r\to \infty }{\text{lim}}{S}_w=S_{\mathit{iw}}.$$E19

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 $\mathrm{ϵ}\to 0$, i.e., when the capillary pressure effects are null. The combination of the self-sharpening tendency of the shock ($S_{\mathit{wf}}>S_{\mathit{iw}}$) with the dispersive effect from the capillary pressure balance against each other leads to the shock layer [22]. Note: in order to guarantee pressure continuity, we have assumed that the capillary pressure gradient is zero at the wellbore, which means that $F_w=f_w=1$ is the wellbore, so $S_w\left(r_w,t\right)=1-S_{\mathit{\text{or}}}$. This boundary condition will be used for both the Buckley-Leverett and the Rapoport-Leas solutions. Ref. [23] presented the idea of using the method of matched asymptotic equations to solve the Rapoport-Leas equation, while [18, 24, 19] showed how the mass balance could be used to present a closed solution for the saturation distribution. Ref. [24] derived an approximate solution for the [1] in Cartesian coordinates for both water and oil injections into a core considering end effects by also applying the method of matched asymptotic expansions. The method of asymptotic expansions uses the inner and outer saturation solutions combined with a matching function in order to obtain a composite solution which avoids abruptly switching from the outer to the inner solution or vice versa. The inner and outer solutions are each capable of representing the real solution in two distinct regions—the “inner region” and the “outer region” of the boundary layer (Figure 6). Similarly, we approximate the saturation solution of Rapoport-Leas equation by forming a composite solution given by the combination of three saturations: $S_w^{\mathit{BL}}$, the solution obtained when the capillary pressure effects are neglected; $S_w^{\mathit{SL}}$, the saturation distribution in the shock layer obtained by magnifying the dispersion effects in the saturation governing equation; and $S_w^{\mathit{SH}}$, the shock wave represented by a Heaviside function:

$$S_w\left(r,t\right)\simeq S_w^{\mathit{BL}}\left(r,t\right)+S_w^{\mathit{SL}}\left(r,t\right)-S_w^{\mathit{SH}}\left(r,t\right),$$E20

where $\mathit{BL}$ stands for Buckley-Leverett, $\mathit{SL}$ for shock layer, and $\mathit{SH}$ for shock function.

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.

Figure 6.

True saturation distribution in the reservoir ($S_w$) compared with the outer solution ($S_w^{\mathit{BL}}$) and the inner solution ($S_w^{\mathit{SL}}$). The dashed square shows the saturation transition zone between the outer and the inner solution where none of these two are capable of approximate $S_w$.

2.1.1. Outer solution ($S_w^{\mathit{BL}}$)

The outer solution, $S_w^{\mathit{BL}}$, is obtained by letting $\mathrm{ϵ}\to 0$ in Eq. (15):

$$S_w^{\mathit{BL}}\left(r,t\right)=\underset{\mathrm{ϵ}\to 0,\left(r,t\right)\text{fixed}}{\text{lim}}{S}_w\left(r,t,\mathrm{ϵ}\right).$$E21

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^{\mathit{BL}}\left(r,0\right)=\left\{\begin{array}{ll}1-S_{\mathit{\text{or}}},& \text{for}r\le r_w\\ S_{\mathit{iw}},& \text{for}r>r_w,\end{array}\right.$$E22

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

$$S_w^{\mathit{BL}}\left(r,t\right)=\left\{\begin{array}{ll}1-S_{\mathit{\text{or}}},& r^2⩽r_w^2\\ {\left(\frac{{\mathit{df}}_w}{{\mathit{dS}}_w}\right)}^{-1}\left(\frac{1}{C}\frac{\left(r^2-r_w^2\right)}{t}\right),& r_w^2<r^2⩽r_w^2+\mathit{Dt}\\ S_{\mathit{iw}},& r^2>\mathit{Dt}+r_w^2,\end{array}\right.$$E23

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

$$D=C\frac{\left[f_w\left(S_{\mathit{wf}}\right)-f_w\left(S_{\mathit{iw}}\right)\right]}{\left[S_{\mathit{wf}}-S_{\mathit{iw}}\right]},$$E24

where $S_{\mathit{wf}}$ and $S_{\mathit{iw}}$ are the shock saturations. In this case, in order to satisfy the conservation of mass, the shock speed should correspond to the slope of a tangent line to the water fractional flow curve, i.e.,

$$\frac{f_w\left(S_{\mathit{wf}}\right)-f_w\left(S_{\mathit{iw}}\right)}{S_{\mathit{wf}}-S_{\mathit{iw}}}=\frac{{\mathit{df}}_w\left(S_{\mathit{wf}}\right)}{{\mathit{dS}}_w}.$$E25

The details of this solution can be found in [5]. Figure 7 shows the shock jump slope tangent to the fractional flow curve at $S_w=S_{\mathit{wf}}$ and the saturation distribution in the reservoir at a time $t$. The rarefaction wave family spans from $1-S_{\mathit{\text{or}}}$ to $S_{\mathit{wf}}$ from $r_w$ to $r=25$ ft, the water front position, i.e., the shock front position, r_s. Ahead of the water front position, there is an immobile water. Figure 8 compares this solution, the outer solution, with the true solution; there is the convection-dispersion saturation profile. Here, we call the true solution the solution obtained from a numerical simulator.

Figure 7.

The shock jump slope tangent (blue curve) to the S-shaped fractional flow curve at $S_w=S_{\mathit{wf}}$ (a) and the saturation profile in the reservoir at a time t (b). The rarefaction waves family spans from $1-S_{\mathit{\text{or}}}$ 1 to $S_{\mathit{wf}}$ and from $r_w$ to $r$ = 25 ft, the water front position, i.e., the shock front position, $r_{f,\mathit{inj}}$. Ahead of the water front position, there is immobile water.

Figure 8.

Saturation distribution during the injection period for the outer solution ($S_w^{\mathit{BL}}$), without capillary pressure, and for the true solution ($S_w$), with capillary pressure. Both profiles agree in the region far from the water front, the region outside the dashed square.

2.1.2. Inner solution ($S_w^{\mathit{SL}}$)

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 $S_w^{\mathit{SL}}$, we magnify the shock layer by using a stretching traveling wave coordinate. Similarly, as defined in [24, 18, 19]

$$w=w\left(r,t\right)=\frac{r^2-r_s^2\left(t\right)}{\mathrm{ϵ}},$$E26

where r_s is the shock front position:

$$r_s^2\left(t\right)=r_w^2+\mathit{CDt},$$E27

w is zero at $r=r_s$ and goes to $\pm \infty$ as $\mathrm{ϵ}\to 0$. We rewrite Eq. (15) in terms of moving coordinates, $\left(r,t\right)\to \left(w,\tau \right)$, where $\tau =\tau \left(t\right)=t$. Using Eq. (26) in the transformed equation and multiplying the resulting equation by $\mathrm{ϵ}$ yield

$$\frac{\mathrm{ϵ}}{C}\frac{\partial S_w}{\partial \tau }-D\frac{\partial S_w}{\partial w}+\frac{\partial f_w}{\partial w}-\frac{\partial }{\partial w}\left(2\left(\mathrm{ϵ}w+r_s^2\left(\tau \right)\right)k\left(\mathrm{ϵ}w+r_s^2\left(\tau \right)\right)\frac{\mathrm{\partial \Psi }}{\partial w}\right)=0.$$E28

Figure 9.

Saturation distribution during the injection period for the inner solution ($S_w^{\mathit{SL}}$) and for the true solution ($S_w$), with capillary pressure. Both profiles agree in the region around the water front, i.e., in the shock boundary layer which we have defined as the inner region.

The inner solution is obtained by letting $\mathrm{ϵ}\to 0$ in Eq. (28):

$$S_w^{\mathit{SL}}\left(w,\tau \right)=\underset{\mathrm{ϵ}\to 0,\left(w,\tau \right)\text{fixed}}{\text{lim}}{S}_w\left(\mathrm{ϵ}w+r_s^2\left(\tau \right),\tau ,\mathrm{ϵ}\right),$$E29

as presented in [26]. Therefore, neglecting the terms of order $\mathrm{ϵ}$ in Eq. (28), we have

$$-D\frac{\partial S_w^{\mathit{SL}}}{\partial w}+\frac{\partial f_w}{\partial w}-\frac{\partial }{\partial w}\left(2r_s^2\left(\tau \right)k\left(r_s^2\left(\tau \right)\right)\frac{\mathrm{\partial \Psi }}{\partial w}\right)=0.$$E30

Note that here we are treating the permeability k as function of the shock position radius, r_s, only, by assuming that in the limit of the inner solution, $\mathrm{ϵ}\left(r\to r_s\left(\tau \right)\right)$. 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

$$-{\mathit{DS}}_w^{\mathit{SL}}+f_w\left(S_w^{\mathit{SL}}\right)-2r_s^2\left(\tau \right)k\left(r_s^2\left(\tau \right)\right)\frac{d\mathrm{\Psi }\left(S_w^{\mathit{SL}}\right)}{{\mathit{dS}}_w}\frac{\partial S_w^{\mathit{SL}}\left(w,\tau \right)}{\partial w}=a\left(\tau \right),$$E31

where $a\left(\tau \right)$ 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

$$\left\{\begin{array}{ll}w\to \infty :& S_w^{\mathit{SL}}=S_{\mathit{iw}},\frac{\partial S_w^{\mathit{SL}}}{\partial w}=0,\\ w\to -\infty :& S_w^{\mathit{SL}}=S_{\mathit{wf}},\frac{\partial S_w^{\mathit{SL}}}{\partial w}=0,\end{array}\right.$$E32

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\left(\tau \right)=-{\mathit{DS}}_{\mathit{iw}},$$E33

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

$$a\left(\tau \right)=-{\mathit{DS}}_{\mathit{wf}}+f_w\left(S_{\mathit{wf}}\right);$$E34

implying that $-D\left(S\_\left\{\mathit{iw}\right\}-S\_\left\{\mathit{wf}\right\}\right)-f\_w\left(S\_\left\{\mathit{wf}\right\}\right)=0$, which it is indeed correct from the definition of D in Eq. (24). As we can see from Eqs. (33) and (34), $a\left(\tau \right)$ 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\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)$ yield

$$2r_s^2\left(\tau \right)k\left(r_s^2\left(\tau \right)\right)\frac{\frac{d\mathrm{\Psi }\left(S_w^{\mathit{SL}}\right)}{{\mathit{dS}}_w}}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}\frac{\partial S_w^{\mathit{SL}}\left(w,\tau \right)}{\partial w}=1.$$E35

Integrating Eq. (35) from $w_{\mathit{\text{well}}}=w\left(r_w,\tau \right)$ to any $w$ at any time $\tau$ gives us the relationship between any $S_w^{\mathit{SL}}$ and $w$:

$$2r_s^{2}k\left(r_s^2\right){\int }_{S_w^{\mathit{SL}}\left(w_{\mathit{\text{well}}}\right)}^{S_w^{\mathit{SL}}}\frac{\frac{d\mathrm{\Psi }\left(S_w^{\mathit{SL}}\right)}{{\mathit{dS}}_w}}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}{\mathit{dS}}_w^{\mathit{SL}}={\int }_{w_{\mathit{\text{well}}}}^w\mathit{dw}.$$E36

At $S_w^{\mathit{SL}}=S_{\mathit{wf}}$, the integral in the left side of Eq. (36) diverges as the integrand denominator goes to 0. This behavior is consistent with our boundary condition assumptions for $S_w^{\mathit{SL}}$ (Eq. (32)). At $S_w^{\mathit{SL}}=S_{\mathit{iw}}$, the integral in the left side of Eq. (36) does converge when the integrand numerator also goes to zero (Figure 4), a behavior which is consistent when trying to model a hypodispersion phenomenon (for a partially water wet reservoir). In another reservoir wettability scenario, e.g., a strong water wet rock, the capillary pressure would not be bounded at Siw, and the integral in the left side of Eq. (36) would diverge. Note: we still do not know the value of $S_w^{\mathit{SL}}$ 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 ($S_w^{\mathit{SH}}$)

The matching saturation $S_w^{\mathit{SH}}$ is defined using the matching principle by applying Prandtl’s technique [26]:

$$\underset{r^2\to r_w^2+\mathit{CDt}}{\text{lim}}{S}_w^{\mathit{BL}}\left(r,t\right)=\underset{w\to \pm \infty }{\text{lim}}{S}_w^{\mathit{SL}}\left(w,t\right);$$E37

and in the injection case is given by

$$S_w^{\mathit{SH}}\left(r,t\right)=\left\{\begin{array}{ll}{S}_{\mathit{iw}},& r^2\ge r_s^2\left(t\right)=r_w^2+\mathit{CDt},\\ S_{\mathit{wf}},& r_w^2\le r_s^2\left(t\right).\end{array}\right.$$E38

which is plotted in Figure 10 against the outer and inner solutions. As we were searching for, $S_w^{\mathit{SH}}$ 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.

Figure 10.

The matching saturation function ($S_w^{\mathit{SH}}$) compared with the outer solution ($S_w^{\mathit{BL}}$) (a) and with the inner solution ($S_w^{\mathit{SH}}$) (b). $S_w^{\mathit{SH}}$ 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). $S_w^{\mathit{SH}}$ matches with the region inner solution in the outer region, the region far from the water front.

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 ($S_w^{\mathit{BL}}$) solution and the composite solution ($S_w$) must obey material balance, the following equality

$$q_{t}t={\int }_{r_w^2}^{\infty }\left(S_w\left(r,t\right)-S_{\mathit{iw}}\right)\pi {\mathit{hdr}}^2={\int }_{r_w^2}^{\infty }\left(S_w^{\mathit{BL}}\left(r,t\right)-S_{\mathit{iw}}\right)\pi {\mathit{hdr}}^2.$$E39

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

$${\int }_{r_w^2}^{\infty }\left(S_w^{\mathit{BL}}+S_w^{\mathit{SL}}-S_w^{\mathit{SH}}-S_{\mathit{iw}}\right)\pi {\mathit{hdr}}^2={\int }_{r_w^2}^{\infty }\left(S_w^{\mathit{BL}}\left(r,t\right)-S_{\mathit{iw}}\right)\pi {\mathit{hdr}}^2,$$E40

which, upon simplification, gives

$${\int }_{r_w^2}^{\infty }\left(S_w^{\mathit{SL}}-S_w^{\mathit{SH}}\right){\mathit{dr}}^2=0.$$E41

Rearranging Eq. (41) using Eq. (38) for $S_w^{\mathit{SH}}$ gives

$${\int }_{r_w^2}^{\infty }{S}_w^{\mathit{SL}}{\mathit{dr}}^2={\int }_{r_w^2}^{r_s^2}{S}_{\mathit{wf}}{\mathit{dr}}^2+{\int }_{r_s^2}^{\infty }{S}_{\mathit{iw}}{\mathit{dr}}^2=S_{\mathit{wf}}\left(r_s^2-r_w^2\right)+\left(S_{\mathit{iw}}{\int }_{r_w^2}^{\infty }d{r}^2-S_{\mathit{iw}}{\int }_{r_w^2}^{r_s^2}d{r}^2\right).$$E42

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

$${\int }_{r_w^2}^{\infty }{S}_w^{\mathit{SL}}{\mathit{dr}}^2=\left(S_{\mathit{wf}}-S_{\mathit{iw}}\right)\mathit{CDt}+S_{\mathit{iw}}{\int }_{r_w^2}^{\infty }d{r}^2.$$E43

Transforming Eq. (43) from $\left(r,t\right)\to \left(w,\tau \right)$ and using Eq. (26), Eq. (43) becomes

$$\mathrm{ϵ}{\int }_{-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}}^{\infty }{S}_w^{\mathit{SL}}\left(w\right)\mathit{dw}=\left(S_{\mathit{wf}}-S_{\mathit{iw}}\right)\mathit{CD\tau }+\mathrm{ϵ}{S}_{\mathit{iw}}{\int }_{-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}}^{\infty }\mathit{dw}.$$E44

From Eq. (35),

$$\mathit{dw}=2r_s^{2}k\left(r_s^2\right)\frac{\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w}{\mathit{dS}}_w^{\mathit{SL}}.$$E45

Substituting Eq. (45) in Eq. (44) and solving the resulting equation divided by $\mathrm{ϵ}{S}_{\mathit{iw}}$ for

$${\int }_{-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}}^{\infty }\mathit{dw}=\frac{2r_s^{2}k\left(r_s^2\right)}{S_{\mathit{iw}}}{\int }_{S_w^{\mathit{SL}}\left(-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}\right)}^{S_{\mathit{iw}}}{S}_w^{\mathit{SL}}\frac{\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}\left(S_w^{\mathit{SL}}\right)}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}{\mathit{dS}}_w^{\mathit{SL}}-\frac{\left(S_{\mathit{wf}}-S_{\mathit{iw}}\right)\mathit{CD\tau }}{\mathrm{ϵ}{S}_{\mathit{iw}}}.$$E46

Setting $S_w=S_{\mathit{iw}}$ in the upper limits of the integrals of Eq. (36) and exchanging the two sides of the equation yield

$${\int }_{-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}}^{\infty }\mathit{dw}=2r_s^{2}k\left(r_s^2\right){\int }_{S_w^{\mathit{SL}}\left(-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}\right)}^{S_{\mathit{iw}}}\frac{\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}\left(S_w^{\mathit{SL}}\right)}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}{\mathit{dS}}_w^{\mathit{SL}}.$$E47

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

$$2r_s^{2}k\left(r_s^2\right){\int }_{S_w^{\mathit{SL}}\left(-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}\right)}^{S_{\mathit{iw}}}\frac{\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}\left(S_w^{\mathit{SL}}\right){\mathit{dS}}_w^{\mathit{SL}}}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}=\frac{2r_s^{2}k\left(r_s^2\right)}{S_{\mathit{iw}}}{\int }_{S_w^{\mathit{SL}}\left(-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}\right)}^{S_{\mathit{iw}}}\frac{\frac{S_w^{\mathit{SL}}d\mathrm{\Psi }}{{\mathit{dS}}_w}\left(S_w^{\mathit{SL}}\right){\mathit{dS}}_w^{\mathit{SL}}}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}-\frac{\left(S_{\mathit{wf}}-S_{\mathit{iw}}\right)\mathit{CD\tau }}{\mathrm{ϵ}{S}_{\mathit{iw}}}.$$E48

Multiplying Eq. (48) by $\mathrm{ϵ}{S}_{\mathit{iw}}$ and rearranging the resulting equation give

$$2r_s^{2}k\left(r_s^2\right)\mathrm{ϵ}{\int }_{S_w^{\mathit{SL}}\left(-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}\right)}^{S_{\mathit{iw}}}\left(S_w^{\mathit{SL}}-S_{\mathit{iw}}\right)\frac{\frac{d\mathrm{\Psi }}{{\mathit{dS}}_w}\left(S_w^{\mathit{SL}}\right)}{D\left(S_{\mathit{iw}}-S_w^{\mathit{SL}}\right)+f_w\left(S_w^{\mathit{SL}}\right)}{\mathit{dS}}_w^{\mathit{SL}}=\left(S_{\mathit{wf}}-S_{\mathit{iw}}\right)\mathit{CD\tau }.$$E49

Once the value $S_w^{\mathit{SL}}\left(-\frac{\mathit{CD\tau }}{\mathrm{ϵ}}\right)$ (i.e., the inner solution saturation in the wellbore $S_w^{\mathit{SL}}\left(w_{\mathit{\text{well}}}\right)$) is determined numerically by solving Eq. (49) using the bisection method at each time τ, Eq. (36) is used to determine the saturation profile in the stabilized zone. It is important to note that, as $S_w^{\mathit{SL}}$ should reach $S_{\mathit{wf}}$ and $S_{\mathit{iw}}$ asymptotically as $w\to \pm \infty$, here we did not have to fix a finite distance in which the traveling wave would reach its open bounds as done by [18, 19, 24]. With our approach, as shown in the validation section (Figure 12), we can obtain essentially a perfect match with the numerical solution, with a “smoother” water front, which is expected from the dispersive effect of capillary pressure, contrary to the sharp transition between the saturation at the water front foot ($r_{\mathit{wf}}$) which is the finite position at which water can be considered immobile—and the initial water saturation in the oil zone exhibited by solutions of previous authors.

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, $q_t\left(r_w,t\right)$ 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\left(r\right)\mathit{hr}}{\alpha }\left({\lambda }_t\frac{\partial p_o}{\partial r}-{\lambda }_w\frac{\partial p_c}{\partial r}\right),$$E50

where $p_w=p_o-p_c$. Eq. (50) can be solved for the oil pressure gradient by integrating it from the wellbore radius to infinite, assuming an infinite-acting reservoir. The bottom hole pressure difference from the reservoir initial pressure ($p_{\mathit{oi}}$) can then be expressed as

$$\mathrm{\Delta }{p}_{\mathit{wf}}\left(t\right)=p_{\mathit{wf}}\left(t\right)-p_{\mathit{oi}}=\underset{r_w}{\overset{\infty }{\int }}\frac{\alpha q_t\left(r,t\right)}{h{\lambda }_t\left(r,t\right)k\left(r\right)}\frac{\mathit{dr}}{r}-\underset{r_w}{\overset{\infty }{\int }}{f}_w\frac{{\mathit{dp}}_c}{{\mathit{dS}}_w}\frac{\partial S_w}{\partial r}\mathit{dr},$$E51

where it is assumed that $p_o\left(r_w,t\right)=p_w\left(r_w,t\right)$, i.e., $p_c=0$ at $r=r_w$, in order to satisfy the compatibility condition [27], i.e., to guarantee phase pressure continuity at the wellbore. Eq. (51) can be rewritten as

$$\mathrm{\Delta }{p}_{\mathit{wf}}\left(t\right)=\underset{r_w}{\overset{\infty }{\int }}\frac{\alpha q_t\left(r,t\right)}{h{\lambda }_t\left(r,t\right)k\left(r\right)}\frac{\mathit{dr}}{r}-\underset{r_w}{\overset{r_{\mathit{wf}\left(t\right)}}{\int }}{f}_w\frac{{\mathit{dp}}_c}{{\mathit{dS}}_w}\frac{\partial S_w}{\partial r}\mathit{dr},$$E52

by assuming that the second term in the right-hand side of Eq. (52) is zero from $r_{\mathit{wf}}$ to $\infty$, considering $f_w\left(S_{\mathit{iw}}\right)=0$ and $\frac{\partial S_w}{\partial r}=0$ for $r>r_{\mathit{wf}}\left(t\right)$, since the water in the region ahead of the water front foot is assumed immobile. $r_{\mathit{wf}}$ can be defined as the position at which $\left(S_w-S_{\mathit{iw}}\right)<\delta$, where $\delta$ is a very small number. For the hypodispersion phenomenon, we can find a finite $r_{\mathit{wf}}$ where $\delta \to 0$, i.e., at which $S_w=S_{\mathit{iw}}$. Using the [6] steady-state theory, which assumes that, $q_t\left(r,t\right)=q_t\left(r_w,t\right)$, for $r\le r_{\mathit{wf}}\left(t\right)$, Eq. (52) becomes

$$\mathrm{\Delta }{p}_{\mathit{wf}}\left(t\right)=\frac{\alpha q_t\left(r_w,t\right)}{h}\underset{r_w}{\overset{r_{\mathit{wf}}\left(t\right)}{\int }}\frac{1}{{\lambda }_t\left(r,\mathrm{\Delta }{t}_{\mathit{\text{prod}}}\right)k\left(r\right)}\frac{\mathit{dr}}{r}+\frac{\alpha }{h}\underset{r_{\mathit{wf}}\left(t\right)}{\overset{\infty }{\int }}\frac{q_t\left(r,t\right)}{{\lambda }_t\left(r,t\right)k\left(r\right)}\frac{\mathit{dr}}{r}-\underset{r_w}{\overset{r_{\mathit{wf}\left(t\right)}}{\int }}{f}_w\frac{{\mathit{dp}}_c}{{\mathit{dS}}_w}\frac{\partial S_w}{\partial r}\mathit{dr},$$E53

where for any practical set of values of physical properties [28] indicate that this assumption is valid. Adding and subtracting the term $\frac{\alpha }{h}\underset{r_w}{\overset{r_{\mathit{wf}}\left(t\right)}{\int }}\frac{q_t\left(r_w,t\right)}{{\widehat{\lambda }}_{o}k\left(r\right)}\frac{\mathit{dr}}{r}$, where ${\widehat{\lambda }}_o=\frac{k_{\mathit{ro}}\left(S_{\mathit{wi}}\right)}{{\mu }_o}$ is the endpoint oil mobility at $S_w=S_{\mathit{wi}}$, Eq. (53) can be rewritten as

$$\begin{array}{cc}\mathrm{\Delta }{p}_{\mathit{wf}}\left(t\right)=& \frac{\alpha }{h}\underset{r_w}{\overset{\infty }{\int }}\frac{q_t\left(r,t\right)}{{\widehat{\lambda }}_o\left(r,t\right)k\left(r\right)}\frac{\mathit{dr}}{r}+\frac{\alpha q_t\left(r_w,t\right)}{h}\underset{r_w}{\overset{r_{\mathit{wf}}\left(t\right)}{\int }}\left(\frac{1}{{\lambda }_t\left(r,t\right)}-\frac{1}{{\widehat{\lambda }}_o}\right)\frac{\mathit{dr}}{k\left(r\right)r}-\underset{r_w}{\overset{r_{\mathit{wf}\left(t\right)}}{\int }}{f}_w\frac{{\mathit{dp}}_c}{{\mathit{dS}}_w}\frac{\partial S_w}{\partial r}\mathit{dr}\\ & =\mathrm{\Delta }{\widehat{p}}_o\left(t\right)+\frac{\alpha q_t\left(r_w,t\right)}{h{\widehat{\lambda }}_o}\underset{r_w}{\overset{r_{\mathit{wf}}\left(t\right)}{\int }}\left(\frac{{\widehat{\lambda }}_o}{{\lambda }_t\left(r,t\right)}-1\right)\frac{1}{k\left(r\right)}\frac{\mathit{dr}}{r}\\ -\underset{r_w}{\overset{r_{\mathit{wf}\left(t\right)}}{\int }}{f}_w\frac{{\mathit{dp}}_c}{{\mathit{dS}}_w}\frac{\partial S_w}{\partial r}\mathit{dr}.\end{array}$$E54

$\mathrm{\Delta }{\widehat{p}}_o\left(t\right)$ 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

$$\mathrm{\Delta }{\widehat{p}}_o\left(t\right)=p_{\mathit{wf},o}\left(t\right)-p_i=\frac{\alpha q_t}{\mathit{kh}{\widehat{\lambda }}_o}\left[\frac{1}{2}\text{ln}\left(\frac{\mathit{\beta k}{\widehat{\lambda }}_{o}t}{\varphi {\widehat{c}}_{\mathit{\text{to}}}{r}_w^2}\right)+0.4045+s\right].$$E55

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

$${\widehat{c}}_{\mathit{\text{to}}}=c_o\left(1-S_{\mathit{wi}}\right)+c_w{S}_{\mathit{wi}}+c_r.$$E56

3. Validation

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}^{-7}$ 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.

Property	Value	Unit	Property	Value	Unit
$q_t$	3000	RB/DAY	$B_o$	1.003	RB/STB
$h$	60	ft	$B_w$	1.002	RB/STB
$r_w$	0.35	ft	$c_o$	$8\times {10}^{-6}$	1/psi
$r_e$	6800	ft	$c_w$	$3.02\times {10}^{-6}$	1/psi
$k$	300	ft	$c_r$	$5\times {10}^{-6}$	1/psi
$s$	0	mD	${\mu }_o$	3.0	cp
$S_{\mathit{iw}}$	0.10		${\mu }_w$	0.5	cp
$S_{\mathit{\text{or}}}$	0.25		$\lambda$	2
$p_i$	2500	psi	$p_t$	0.5	psi
$\varphi$	0.22

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

4. Conclusions

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

Acknowledgments

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.

Nomenclature

β	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)
Sw	Water saturation
Siw	Immobile water saturation
Sor	Residual oil saturation
t	Time (h)
C	Constant given by θqtπhϕ
D	Shock speed
h	Reservoir thickness (ft)
r	Radius (ft)
w	Traveling wave coordinate
α	Unit conversion factor (141.2)
ϵ	Perturbation parameter
λ	Pore size distribution index
ϕ	Porosity
θ	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)