Abstract
1D water oil displacement in porous media is usually described by the Buckley-Leverett equation or the Rapoport-Leas equation when capillary diffusion is included. The rectilinear geometry is not representative for near well oil displacement problems. It is therefore of interest to describe the radially symmetric Buckley-Leverett or Rapoport-Leas equation in cylindrical geometry (radial Buckley-Leverett problem). We can show that under appropriate conditions, one can apply a similarity transformation \((r,t) \rightarrow \eta = r^{2}/(2t) \) that reduces the PDE in radial geometry to an ODE, even when capillary diffusion is included (as opposed to the situation in the rectilinear geometry (Yortsos, Y.C. (Phys. Fluids 30(10),2928–2935 1987)). We consider two cases (1) where the capillary diffusion is independent of the saturation and (2) where the capillary diffusion is dependent on the saturation. It turns out that the solution with a constant capillary diffusion coefficient is fundamentally different from the solution with saturation-dependent capillary diffusion. Our analytical approach allows us to observe the following conspicuous difference in the behavior of the dispersed front, where we obtain a smoothly dispersed front in the constant diffusion case and a power-law behavior around the front for a saturation-dependent capillary diffusion. We compare the numerical solution of the initial value problem for the case of saturation-dependent capillary diffusion obtained with a finite element software package to a partially analytical solution of the problem in terms of the similarity variable η.
Article PDF
Similar content being viewed by others
Notes
The derivative approaches zero because \(S_{w}^{in}\) is a continuously decreasing function of ξ, approaching a constant as \(\xi \rightarrow -\infty \) (see the last equation in (134))
The derivative approaches zero because \(S_{w}^{in}\) is a continuously increasing function of ξ, approaching a constant as \(\xi \rightarrow \infty \) (see Eq. 74)
References
Chavent, G., Jaffré, J.: Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media, vol. 17. Elsevier, New York (1986)
Dake, L.P.: Fundamentals of reservoir engineering. 1978 (1983)
Deng, L., King, M.J., et al.: Capillary corrections to buckley-leverett flow. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers (2015)
Evje, S., Karlsen, K.H.: Viscous splitting approximation of mixed hyperbolic-parabolic convection-diffusion equations. Numer. Math. 83(1), 107–137 (1999)
Karlsen, H.K., Brusdal, K., Dahle, H.K., Evje, S., Lie, K.-A.: The corrected operator splitting approach applied to a nonlinear advection-diffusion problem. Comput. Methods Appl. Mechan. Eng. 167(3-4), 239–260 (1998)
Krushkov, S.N.: First-order quasilinear equations in several independent variables. Math. Sb. 123, 228–255 (1970)
LeVeque, R.J.: Numerical methods for conservation laws, vol. 3. Springer, New York (1992)
Ling, K.: Fractional flow in radial flow systems: a study for peripheral waterflood. J. Petrol. Explor. Product. Technol. 6(3), 441–450 (2016)
Martin, B.: Imperial college lectures. In: Petroleum engineering, the-volume 2: reservoir engineering, vol. 2. World Scientific Publishing Company, Singapore (2017)
Oleinik, O.A.: Discontinuous solutions of non-linear differential equations. Uspekhi Matematicheskikh Nauk 12(3), 3–73 (1957)
Rapoport, L.A., Leas, W.J., et al.: Properties of linear waterfloods. J. Petrol. Tech. 5(05), 139–148 (1953)
Smoller, J.: Shock waves and reaction—diffusion equations, vol. 258. Springer Science & Business Media, New York (2012)
Van Dyke, M.: Perturbation methods in fluid mechanics, vol. 136. Academic Press, New York (1964)
Henry, J., et al.: Welge a simplified method for computing oil recovery by gas or water drive. J. Petrol. Tech. 4(04), 91–98 (1952)
Yortsos, Y.C.: Stability of displacement processes in porous media in radial flow geometries. Phys. Fluids 30(10), 2928–2935 (1987)
Yortsos, Y.C., Fokas, A.S., et al.: An analytical solution for linear waterflood including the effects of capillary pressure. Soc. Pet. Eng. J. 23(01), 115–124 (1983)
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix: 1 Rankine-Hugoniot condition
Mass conservation across the shock specifies the shock velocity in terms of the shock saturation; this condition is called the Rankine-Hugoniot condition. In the case of radial flow, the shock velocity vs is replaced by ηs = rvs, where ηs can be interpreted as (2D) radial flux up to a factor of 2π.
1.1 1.1 The Rankine-Hugoniot condition in the rectilinear (1D) case
We will first give a short summary of the 1D problem. In the 1D (incompressible) problem with pc = 0, we have mass conservation:
and Darcy
with boundary conditions and initial condition:
where uinj is the constant injection velocity. Adding both (99) and using the boundary conditions (101), we find:
which means that our problems (99)–(100) reduces to:
Setting \(\eta =\frac {x}{t}\) we find a rarefaction wave (similarity solution):
followed by a shock at η = ηs to the constant state Sw = Swc. The self similar variable η can be interpreted as a velocity, i.e., saturations above the shock saturation \({S_{w}^{s}}\) travel with speed \(u=\frac {u_{inj}}{\varphi }\frac {df_{w}}{dS_{w}}\). The velocity of the shock vs can be expressed in terms of the saturations in front of the shock and behind the shock using mass conservation as follows. The shock is located at position x0 at time t0 and moves to the right to the position x0 + Δx at time t0 + Δt. This means that we have at t = t0 + Δt:
slightly in front of the shock and
behind the shock (Sw is approximately constant in [x0,x0 + Δx) if Δx is small enough). Furthermore, we have the (1D) volumetric velocities in front of/behind the shock:
We use the following mass balance in the region [x0,x0 + Δx]:
where the net inflow is given by the LHS and the mass accumulation is given on the RHS. Taking the limit \({\Delta } t \rightarrow 0\) yields the shock velocity:
Notice that Eq. 109 only specifies vs (or ηs) in terms of the saturations \(S_{w}^{-}\) and \(S_{w}^{+}\); these saturation values around the shock are given by another physical mechanism called the entropy condition, described in Section B.
1.2 1.2 The Rankine-Hugoniot condition in the radial case
The derivation of the Rankine-Hugoniot in the radial case goes along the same lines as the 1D case. In the radial (incompressible) problem, we have:
and Darcy
with boundary condition and initial condition
Adding both (110) and using the boundary conditions (112), we find:
which means that our problems (110)–(111) reduce to
(we keep the product rutot throughout the section, even though it equals 1, to facilitate comparison with the 1D-case). Setting \(\eta =\frac {r^{2}}{2t}\), we find a rarefaction wave:
followed by a shock at η = ηs to the constant state Sw = Swc. The self similar variable η can be interpreted as 2D-volumetric flux (η = ru), i.e., saturations above \({S_{w}^{s}}\) travel with “speed” \(ru=\frac {ru_{tot}}{\varphi }\frac {df_{w}}{d S_{w}}\).
The shock is located at position r0 at time t0 and at position r0 + Δr at time t0 + Δt. This means that we have
slightly in front of the shock and
behind the shock (Sw is approximately constant in [r0,r0 + Δr) if Δr is small enough). Furthermore, we have the (2D) volumetric velocities in front of/behind the shock:
The mass balance (net inflow equals mass accumulation) in the annulus between r = r0 and r = r0 + Δr yields:
Taking the limit \({\Delta } t \rightarrow 0\) yields the shock velocity:
Notice that Eq. 120 only specifies vs (or ηs = rvs) in terms of the saturations \(S_{w}^{-}\) and \(S_{w}^{+}\); these saturation values are given by the entropy condition in Section 2.
Appendix 2: The entropy condition
The Rankine-Hugoniot condition does not specify the shock saturation. Behind the shock, the velocity (in the 1D case) is given as:
This means that we have at the shock saturation \({S_{w}^{s}}=S_{w}^{-}\)
If \(v_{s}< \frac {u_{tot}}{\varphi }\left .\frac {df_{w}}{S_{w}}\right \vert _{{S_{w}^{s}}}\) higher saturations would overtake the shock and the solution would become multivalued. According to the entropy condition (140) derived in Section 2.1 we also have:
which means that we have:
i.e.,
Due to the specific functional form of fw there is usually one solution that shocks to the base state, which means that \(S_{w}^{+} = S_{wc}\), i.e.,
This means that Eq. 126 has only one solution \({S_{w}^{s}}\), which corresponds to the solution found using Welge’s tangent construction.
In the radial case, a similar argument holds; we need:
otherwise higher saturations would have a higher 2D volumetric velocity and would overtake the shock, leading to a multivalued solution. Furthermore, we will show in Section 2.2 (see Eq. 143) that:
and we can conclude that the shock saturation is given by:
just as in the 1D case.
1.1 2.1 The entropy condition in the 1D case
The idea of the entropy condition is that the physical solution of problems (99)–(101) is given by the solution of equation:
with boundary condition and initial condition
in the limit \(\epsilon \rightarrow 0\).
We can solve problems (130)–(131) using the method of matched asymptotic expansions around the shock. The outer solution is given by a rarefaction wave, possibly followed by a constant state, i.e., approaching the shock we find:
where we use the moving coordinate η = x − vst. We solve the inner problem in the comoving frame and we set:
to obtain the inner problem:
Integrating once from \(-\infty \) to ξ yields:
where we used that \(\frac {dS_{w}^{in}}{d\xi }|_{{S_{w}^{s}}} \rightarrow 0\) as \(\xi \rightarrow -\infty \)Footnote 1. Solving for \(\frac {dS_{w}^{in}}{d\xi }\) yields:
and using a Taylor approximation for fw around \({S_{w}^{s}}\) we obtain:
i.e.,
where the inequality holds because \(S_{w}^{in}\) is a decreasing function of ξ; \(\frac {dS_{w}^{in}}{d\xi }\) is determined by \(S_{w}^{in}\) (see Eq. 136), so if \(\frac {dS_{w}^{in}}{d\xi }\) would be positive for some value of \(S_{w}^{in}\), the saturation will never drop below this value of \(S_{w}^{in}\) and reach Swc. We have \(S_{w}^{in}-{S_{w}^{s}} \leq 0\) which means that we need:
which leads to the following inequality for vs
used in the previous section.
1.2 2.2 The entropy condition in the radial case
We use the method of matched asymptotic expansions (see Section 3.3) to derive (128). Integrating (81) from \(\infty \) to ξ yields:
where we used that \(\left .\frac {dS_{w}^{in}}{d\xi }\right \vert _{{S_{w}^{s}}} \rightarrow 0\) as \(\xi \rightarrow \infty \).Footnote 2 Using a Taylor approximation for fw around \({S_{w}^{s}}\) and solving for \( \frac {dS_{w}}{d\xi }\) we find:
Notice that the LHS of Eq. 142 is positive (the saturation is an increasing function of ξ due to the choice of ξ), which means that we find:
Appendix 3: Mathematical details of Section 3.2.3
In this section, we will discuss briefly how Eq. 60 is derived using Eqs. 58 and 59.
For Sw ≈ Swc, only the first and last terms of Eq. 42 are nonzero. The second term vanishes due to the incompressibility assumption; for the third term, we have:
because rutot = 1 in the incompressible case (see Eq. 54). Furthermore, we have λo ≈ 0 ⇒ fw ≈ 1, which means that the third term indeed vanishes.
This means that we have to balance the first and fourth terms to get a nontrivial solution.
For the first term \(-\varphi \eta \frac {dS_{w}}{d\eta }\), we find:
and for the fourth term \(\frac {d}{d\eta } \left (2\eta D_{cap}\frac {dS_{w}}{d\eta } \right )\), we find three contributions:
where
and
Note that the contributions of Eqs. 148 and 149 are of the same order and that the contribution of Eq. 147 is of higher order and can be neglected. Balancing powers of ηf − η of Eqs. 145, 148, and 149, we indeed find a nontrivial solution, provided the balance:
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Meulenbroek, B., Gargar, N.K. & Bruining, H. An engineering approach to study the effect of saturation-dependent capillary diffusion on radial Buckley-Leverett flow. Comput Geosci 25, 637–653 (2021). https://doi.org/10.1007/s10596-020-09993-y
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10596-020-09993-y