Approximate Semi-analytical Solution for Injection–Falloff–Production Well Test: An Analytical Tool for the In Situ Estimation of Relative Permeability Curves

Abstract

An injection–falloff–production test (IFPT) was originally proposed in Chen et al. (in: SPE conference paper, 2006. doi:10.2118/103271-MS, SPE Reserv Eval Eng 11(1):95–107, 2008) as a well test for the in situ estimation of two-phase relative permeability curves to be used for simulating multiphase flows in porous media. Hence, we develop an approximate semi-analytical solution for the two-phase saturation distribution in an oil–water system during the flowback period of an IFPT according to the mathematical theory of waves. In fact, we show that the weak solution we construct for the saturation equation for the flowback period satisfies the Oleinik entropy condition and hence is unique. In addition, we allow the governing relative permeabilities during the flowback period to be different from the relative permeabilities during injection. Using the saturation solution with the steady-state pressure theory of Thompson and Reynolds, we obtain a solution for the wellbore pressure during the flowback period. By comparing results from our solution with those from a commercial numerical simulator, we show that our approximate semi-analytical solution yields accurate saturation profiles and bottom hole pressures history. The use of very small time steps and a highly refined radial grid is necessary to generate a good solution from a reservoir simulator. The approximate analytical pressure solution developed is used as a forward model to match pressure and water flow rate data from an IFPT in order to estimate reservoir rock absolute permeability and skin factor in conjunction with in situ imbibition and drainage water–oil relative permeabilities.

Appendices

Appendix 1: Derivation of the Shock Position During Flowback Period

In order to obtain the shock position expression for the flowback period, Eq. 15 can be rearranged as

$$\begin{aligned} \int \limits _{r_\mathrm{s}^2}^{r_{S_\mathrm{wf}}^2} S_\mathrm{w1} \text {d}r^2 - \int \limits _{r_\mathrm{IP}^2}^{r_{S_\mathrm{wf}}^2} S_\mathrm{w2} \text {d}r^2+ \int \limits _{r_\mathrm{IP}^2}^{r_\mathrm{s}^2} S_\mathrm{w3} \text {d}r^2 = 0, \end{aligned}$$

(32)

where we can evaluate each integral in Eq. 32 by the following integration by parts formula

$$\begin{aligned} \int \limits _{a}^{b} S_{\mathrm{w}j} \text {d}r^2 = S_{\mathrm{w}j}(r,t)r^2|_{a}^b-\int \limits _{S_{\mathrm{w}j}(a)}^{S_{\mathrm{w}j}(b)} r^2 \text {d}S_{\mathrm{w}j},\quad j=1,2,3. \end{aligned}$$

(33)

Substituting the expression for r given by Eq. 12 into Eq. 33 gives

$$\begin{aligned} \int \limits _{a}^{b} S_{\mathrm{w}j} \text {d}r^2 = S_{\mathrm{w}j}(r,t)r^2|_{a}^b-\int \limits _{S_{\mathrm{w}j}(a)}^{S_\mathrm{iw}(b)} \left( r_\mathrm{inj}^2 -\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_{\mathrm{w}j}}\right) \text {d}S_{\mathrm{w}j}. \end{aligned}$$

(34)

Equation 34 applies at for any time t greater than the time at which the falloff period ends. We neglect any changes in the saturation profile during the falloff period. The saturation radius position at the end of the injection period ($r_\mathrm{inj}$) is defined as

$$\begin{aligned} r_\mathrm{inj}^2 \equiv r_\mathrm{inj}^2 (S_\mathrm{w1},t_\mathrm{inj}) = r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{w1})}{\text {d}S_\mathrm{w1}} \end{aligned}$$

(35)

for curve $S_1$ and as

$$\begin{aligned} r_\mathrm{inj}^2 (S_{w2,3},t_\mathrm{inj}) \equiv r_{\mathrm{f,inj}}^2 = r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_{w2,3}} \end{aligned}$$

(36)

for curves $S_\mathrm{w2}$ and $S_\mathrm{w3}$, where $r_{\mathrm{f,inj}}$ is the radius of the water front at the end of injection period. By using Eqs. 35 and 36 in Eq. 34 with $j = 1$ and inserting $S_\mathrm{w1}(r_\mathrm{s},t) = S_\mathrm{w}^+$ and $S_\mathrm{w3}(r_\mathrm{s},t) = S_\mathrm{w}^-$ in the resulting equation, we obtain

$$\begin{aligned} I\equiv & {} \int \limits _{r_\mathrm{s}^2}^{r_{S_\mathrm{wf}}^2} S_\mathrm{w1} \text {d}r^2 = S_\mathrm{wf}r_{\mathrm{f,inj}}^2 - S_\mathrm{w}^+ r_\mathrm{s}^2 \nonumber \\&- \int \limits _{S_\mathrm{w}^+}^{S_\mathrm{wf}} \left( r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{w1}) }{\text {d}S_\mathrm{w1}}- \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{w1}) }{\text {d}S_\mathrm{w1}}\right) \text {d}S_\mathrm{w1}, \end{aligned}$$

(37)

and by integrating the last term in Eq. 37, it follows that

$$\begin{aligned} I= & {} S_\mathrm{wf}r_{\mathrm{f,inj}}^2 - S_\mathrm{w}^+ r_\mathrm{s}^2 - r_\mathrm{w}^2 (S_\mathrm{wf}-S_\mathrm{w}^+) \nonumber \\&-\frac{\theta (q_\mathrm{inj}t_\mathrm{inj} - q_{t,\mathrm{prod}}\varDelta t_\mathrm{prod})}{\pi h \phi }\left( f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right) . \end{aligned}$$

(38)

For the second term in Eq. 32, we use Eq. 34 with $i = 2$ to obtain

$$\begin{aligned} II&\equiv - \int \limits _{r_\mathrm{IP}^2}^{r_{S_\mathrm{w}^2}} S_\mathrm{w2} \text {d}r^2 = S_\mathrm{wIP}r_\mathrm{IP}^2 - S_\mathrm{wf}r_\mathrm{w}^2 \nonumber \\&\quad - \int \limits _{S_\mathrm{wf}}^{S_\mathrm{wIP}} \bigg ( r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj} }{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}} - \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi \phi h}\frac{\text {d}f_\mathrm{w}(S_\mathrm{w2})}{\text {d}S_\mathrm{w2}} \bigg ) \text {d}S_\mathrm{w2} \end{aligned}$$

(39)

and after integrating, it follows that

$$\begin{aligned} II&= S_\mathrm{wIP} r_\mathrm{IP}^2 - S_\mathrm{wf}r_{\mathrm{f,inj}}^2 - \left( r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj} }{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S}\right) (S_\mathrm{wIP}-S_\mathrm{wf}) \nonumber \\&\quad + \frac{\theta q_{t,\mathrm{prod}}\varDelta t_\mathrm{prod}}{\pi \phi h}(f_\mathrm{w}(S_\mathrm{wIP})-f_\mathrm{w}(S_\mathrm{wf})). \end{aligned}$$

(40)

For the third term of Eq. 32, we use Eq. 34 with $i = 3$ to obtain

$$\begin{aligned} III\equiv & {} \int \limits _{r_\mathrm{IP}^2}^{r_\mathrm{s}^2} S_\mathrm{w3} \text {d}r^2 =r_\mathrm{s}^2 S_\mathrm{w}^- - r_\mathrm{IP}^2 S_\mathrm{wIP} \nonumber \\&- \int \limits _{S_\mathrm{wIP}}^{S_\mathrm{w}^-} \left( r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi \phi h}\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf}) }{\text {d}S_\mathrm{w}} - \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{w3})}{\text {d}S_\mathrm{w3}} \right) \text {d}S_\mathrm{w3} \end{aligned}$$

(41)

which after integrating yields

$$\begin{aligned} III= & {} r_\mathrm{s}^2 S_\mathrm{w}^- - r_\mathrm{IP}^2 S_\mathrm{wIP} - \bigg (r_\mathrm{w}^2+ \frac{\theta q_\mathrm{inj} t_\mathrm{inj}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S}\bigg )\left( S_\mathrm{w}^--S_\mathrm{wIP}\right) \nonumber \\&+ \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\left( f_\mathrm{w}\left( S_\mathrm{w}^-\right) -f_\mathrm{w}(S_\mathrm{wIP})\right) . \end{aligned}$$

(42)

From Eq. 32, $I+II+III = 0$, which is equivalent to

$$\begin{aligned}&r_\mathrm{s}^2\left( S_\mathrm{w}^- - S_\mathrm{w}^+\right) - r_\mathrm{w}^2 \left( S_\mathrm{w}^- - S_\mathrm{w}^+\right) + \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\left( f_\mathrm{w}\left( S_\mathrm{w}^-\right) -f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right) \nonumber \\&\quad - \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\left( f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S}\left( S_\mathrm{wf}-S_\mathrm{w}^-\right) \right) =0 \end{aligned}$$

(43)

and finally, solving Eq. 43 for $r_\mathrm{s}$, gives

$$\begin{aligned}&r_\mathrm{s}\left( \varDelta t_\mathrm{prod},S_\mathrm{w}^-,S_\mathrm{w}^+\right) = \Big (r_\mathrm{w}^2-\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{[f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}(S^+)]}{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\nonumber \\&\quad +\frac{ \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\left( f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\left( S_\mathrm{wf}-S_\mathrm{w}^-\right) \frac{\text {d}f_\mathrm{w}}{\text {d}S}(S_\mathrm{wf})\right) }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\Big )^{\frac{1}{2}}. \end{aligned}$$

(44)

Appendix 2: Does the Shock Speed Correspond to the Slope of a Tangent Line to the Water Fractional Flow Curve?

By analogy to the injection Riemann problem, one might think that during production period, the shock speed would also be tangent to the water fractional flow curve, i.e.,

$$\begin{aligned} \frac{\text {d}f_\mathrm{w}\left( S_\mathrm{w}^-\right) }{\text {d}S_\mathrm{w}} = \frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right] }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }. \end{aligned}$$

(45)

Equation 45 was assumed to hold in Chen et al. (2006), but we show that, in general, Eq. 45 is not valid. Substituting $r_\mathrm{s}^2$ from Eq. 21 into Eq. 17, after squaring both sides of Eq. 17, yields

$$\begin{aligned} r_\mathrm{inj}^2(S_\mathrm{wf},t_\mathrm{inj})&- \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w} \left( S_\mathrm{w}^-\right) }{\text {d}S_\mathrm{w}}= r_\mathrm{w}^2-\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}(S^+)\right] }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\nonumber \\&+\frac{ \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi \phi h}\left( f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}\left( S_\mathrm{wf}-S_\mathrm{w}^-\right) \right) }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }, \end{aligned}$$

(46)

which is equivalent to

$$\begin{aligned}&r_\mathrm{w}^2 + \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\frac{\text {d}f_\mathrm{w} (S_\mathrm{wf})}{\text {d}S_\mathrm{w}} - \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w} \left( S_\mathrm{w}^-\right) }{\text {d}S_\mathrm{w}} \nonumber \\&\quad = r_\mathrm{w}^2-\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right] }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\nonumber \\&\qquad +\frac{ \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\left( f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}\left( S_\mathrm{wf}-S^-\right) \right) }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }. \end{aligned}$$

(47)

Equation 47 can be rearranged to obtain

$$\begin{aligned}&- \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^-\right) =-\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi \phi h}\frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right] }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\nonumber \\&\quad +\frac{ \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }\left( f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}\left( S_\mathrm{wf}-S^-\right) \right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^- - S_\mathrm{w}^+\right) }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }, \end{aligned}$$

(48)

or, equivalently,

$$\begin{aligned}&- \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^-\right) =-\frac{\theta q_{t,\mathrm{prod}} \varDelta t}{\pi h \phi }\frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right] }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\nonumber \\&\quad +\frac{ \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }(f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}S_\mathrm{wf}+\frac{\text {d}f_\mathrm{w} (S_\mathrm{wf})}{\text {d}S_\mathrm{w}} S_\mathrm{w}^+}{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }. \end{aligned}$$

(49)

By simply grouping the terms that are multiplied by $\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}$, we obtain

$$\begin{aligned}&- \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^-\right) =-\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-\right) - f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right] }{\left[ S^- - S^+\right] }\nonumber \\&\quad +\frac{ \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi }(f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) -\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}\left( S_\mathrm{wf}-S_\mathrm{w}^+\right) }{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }, \end{aligned}$$

(50)

which is equivalent to

$$\begin{aligned}&- \frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}(S^-) =-\frac{\theta q_{t,\mathrm{prod}} \varDelta t_\mathrm{prod}}{\pi h \phi }\frac{[f_\mathrm{w}(S^-) - f_\mathrm{w}(S^+)]}{\left[ S_\mathrm{w}^- - S_\mathrm{w}^+\right] }\nonumber \\&\quad + \frac{\theta q_\mathrm{inj}t_\mathrm{inj}}{\pi h \phi [S_\mathrm{w}^- - S_\mathrm{w}^+]}\Bigg (f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}(S^+)-\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}[S_\mathrm{wf} - S^+]\Bigg ). \end{aligned}$$

(51)

From Eq. 51, we can see that Eq. 45 is true if and only if

$$\begin{aligned} \frac{\left[ f_\mathrm{w}(S_\mathrm{wf})-f_\mathrm{w}\left( S_\mathrm{w}^+\right) \right] }{\left[ S_\mathrm{wf} - S_\mathrm{w}^+\right] }=\frac{\text {d}f_\mathrm{w}(S_\mathrm{wf})}{\text {d}S_\mathrm{w}}, \end{aligned}$$

(52)

which is not possible for the S-shaped fractional flow curve. Consequently, Eq. 45 does not hold for the flowback solution, i.e., Eq. 45 is not consistent with Eqs. 17, 19 and 21, meaning that the shock speed is not the slope of an tangent line to the water fractional flow curve. Distinct from injection, in the production (flowback) period, we have an unstable shock, resultantly from the non-uniform saturation distribution at the beginning of the production period. The speeds of the saturations behind the oil front are higher than the speed of the shock itself, which has higher speed than the saturations ahead of it. Consequently, the shock is caught up to by the waves behind it and also catches up to the waves ahead of it, as explained in Sect. 2.1.

Appendix 3: Comparison of Area Equality and Front Tracking Shock Fitting Methods

Using the saturation profile at the end of the water injection as the initial condition, the saturation profile in the reservoir at different production times can be obtained with the shock tracking method which has as its basis the integration of Eq. 13. In the front tracking approach, Euler’s method can be applied to find the shock path numerically, with initial conditions given by

$$\begin{aligned} IC: S_\mathrm{w}^-(r,0) = S_\mathrm{iw}, \quad S_\mathrm{w}^+(r,0) = S_\mathrm{wf}, \quad r_\mathrm{s}(0) = r_{\mathrm{f,inj}} \end{aligned}$$

(53)

assuming that shock speed is constant throughout each time step ($\varDelta t$),

$$\begin{aligned} \varDelta r^2 (t^n) = - \frac{\theta q_{t,\mathrm{prod}} \varDelta t}{\pi h \phi } \frac{\left[ f_\mathrm{w}\left( S_\mathrm{w}^-(t^{n-1})\right) - f_\mathrm{w}\left( S_\mathrm{w}^+(t^{n-1})\right) \right] }{\left[ S_\mathrm{w}^-(t^{n-1}) - S_\mathrm{w}^+(t^{n-1})\right] }, \end{aligned}$$

(54)

where $\varDelta r^2 (t^n)$ is the radius squared step at the current time,

$$\begin{aligned} t^n =t^{n-1} + \varDelta t, \end{aligned}$$

(55)

yielding to the updated shock position squared,

$$\begin{aligned} r_\mathrm{s}^2 (t^n) = r_\mathrm{s}^2 (t^{n-1}) + \varDelta r^2 (t^n) . \end{aligned}$$

(56)

Equation 54 is obtained from the discretization of Eq. 13 using a finite difference scheme backward in time. Given $r_\mathrm{s}^2$, $S_\mathrm{w}^-$ and $S_\mathrm{w}^+$ can be obtained by solving

$$\begin{aligned} r_\mathrm{inj}^2\left( S_\mathrm{w}^+,t_\mathrm{inj}\right) - r_\mathrm{s}^2 ( t^n)= \frac{\theta q_{t,\mathrm{prod}} t^n}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^+ (t^n)\right) \end{aligned}$$

(57)

and

$$\begin{aligned} r_\mathrm{inj}^2(S_\mathrm{wf},t_\mathrm{inj}) - r_\mathrm{s}^2 ( t^n) = \frac{\theta q_{t,\mathrm{prod}} t^n}{\pi h \phi }\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^- (t^n)\right) , \end{aligned}$$

(58)

for $\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^+ (t^n)\right)$ and $\frac{\text {d}f_\mathrm{w}}{\text {d}S_\mathrm{w}}\left( S_\mathrm{w}^- (t^n)\right)$ and then $S_\mathrm{w}^+$ and $S_\mathrm{w}^-$ can be computed from the water fractional flow derivative curve.

The disadvantage of this method is that it is necessary to use a small time step in order to obtain an accurate shock path because the method assumes $\varDelta _t$ is sufficiently small so the shock speed does not change within the time step. The best approach (less expensive computationally) to use depends on the amount of data to be matched in the well testing application. If successive data points in time are closely spaced so that very small time steps must be used to predict data at observations, the front tracking procedure would be more efficient than the area equality method.