An engineering approach to study the effect of saturation-dependent capillary diffusion on radial Buckley-Leverett flow

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)→η=r2/(2t)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(r,t) \rightarrow \eta = r^{2}/(2t) $\end{document} 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 η.


Introduction
Water drive recovery of oil is one of the most important secondary recovery methods. The displacement in rectilinear (1D) geometry can be described by the Buckley-Leverett (BL) model. The BL model disregards capillary diffusion but uses a saturation-dependent fractional flow to model the water and oil flux functions. To obtain a unique solution, it is necessary to consider the limit of vanishing capillary diffusion ( [6,9,10,14]), leading to the entropy condition. The solution in 1D consists of a rarefaction wave, [8]). Capillary diffusion can be straightforwardly included in the cylindrical model. This model can be solved using a finite element package ( [1]), if the capillary diffusion is taken large enough. We solve it using COMSOL as detailed in Section 4. The numerical results indicate that the shock is replaced by a continuous curve that follows a power-law behavior: continuous but not differentiable at the toe where S w = S wc . A larger capillary diffusion leads to broader fronts. We also studied the problem using a similarity transformation (r, t) → r 2 /(2t = η). Contrary to the 1D case, the model equations including the capillary diffusion term can be expressed entirely in terms of η. This leads to a system of coupled ODEs for the saturation and pressure. In the incompressible case, the saturation equation decouples from the pressure equation. We obtain a second-order ODE for the saturation [15] and an ODE for the pressure with saturation-dependent coefficients. Yortsos [15] studies the stability of the saturation profile in two asymptotic cases. Yortsos [15] did not consider explicit solutions for general cases of saturation-dependent capillary diffusion in a radially symmetric Buckley-Leverett problem. In this manuscript, we carry out an extension of the work by Yortsos in the sense that we consider generic saturation-dependent capillary diffusion. A second innovation is that we consider explicit solutions for S w and more importantly for the pressure, which can be helpful when these solutions are applied for the pressure buildup. We follow the approach due to Yortsos, where we determine self-similar solutions for cases of small well radii and long-time behavior. Furthermore, also in Yortsos, we approximate the location of the boundary condition on the injection well by shifting the injection well boundary to the origin of the computational domain. This approximation is accurate moderately far from the well and for the longtime behavior. Our more general approach can be used to gain quick insight into the relation between several input and output parameters, as well as for testing of numerical simulators.
The model equations are presented in Section 2 and analyzed in Section 3 using our similarity transformation. The case with a saturation-dependent capillary diffusion is solved with a combination of analytical and numerical methods in Section 3.1. The case with small, but constant, capillary diffusion is solved in Section 3.3 using the method of matched asymptotic expansions. Numerical results are presented in Section 4 and subsequently discussed in Section 5. We end the paper with some conclusions. Although most of the theory in the Appendices is wellknown, we have added the Appendices for the convenience of the reader who needs to recap this theory.
We finally note that the current manuscript is meant to have a descriptive nature, rather than being mathematically rigorous.

Mathematical model
In this section, we formulate our model equations in dimensionless form; we use capital letters or variables with superscript ∼ for (in)dependent variables with dimensions, small letters for dimensionless (in)dependent variables, and we use the subindex c for characteristic scales. Notice that only three scales need to be chosen: length, time, and pressure: ρ c will drop out entirely because gravity is neglected.
The mass conservation equation for water reads in fulldimensional form: The term U w ∂ρ w ∂R is proportional to the square of the pressure gradient ∂p w ∂R which is small and can thus be neglected (since U w ∼ ∂p w ∂R and ∂ρ w ∂R = dρ w dp w ∂p w ∂R see [2], Ch. 5, equation (5.18)). Furthermore, we divide by ρ w and we obtain: Similarly, we find for oil (neglecting ∂ρ o ∂R ): Adding (2) and (3) yields: where U tot = U w + U o . Using isothermal compressibilities for oil and water, we have: We finally obtain: Furthermore, we use Darcy's law, i.e.: where we used the mobility This means that we find: The capillary pressure P cap = P o − P w can be used to express U w in terms of U tot , i.e.: where we used f w = w tot and the capillary diffusion coefficientD cap = − w o tot dP cap dS w . This can be substituted in Eq. 2 to obtain: The model is supplemented with appropriate initial and boundary conditions in Section 2.1; constitutive relations for the mobilities and capillary pressure are given in Section 2.2. The relevant scales for time, length, and pressure (t c , r c and p c ) are discussed in Section 2.3 and used to rewrite the model equations in dimensionless form in Section 2.4. We end with a summary of the model equations for easy reference in Section 2.5.

Initial and boundary conditions
Inflow of pure water at a given rate yields boundary conditions for S w and RU tot at the inlet R = R w , i.e.: Initially in the compressible case, we have oil at constant pressure, which yields initial conditions for S w and P w at T = 0, i.e.: S w = S wc , at T = 0 (14) The initial water pressure is some constant. Due to the gauge invariance, the specific value of the constant does not matter.
In the incompressible case, we only have an initial condition for S w ; however, the time dependence of the pressure enters through the time-dependent saturation profile, but not due to an initial pressure condition for the incompressible case. Notice that Eq. 42 requires a numerical boundary condition for S w as R → ∞; we set: as water travels a finite distance in a finite time, which means that the water saturation effectively remains at its initial value as R → ∞.

Overview of the constitutive relations
Water mobility: Oil mobility: Total mobility and fractional flows: Capillary pressure: We will use the dimensionless parameters from Table 1 in the constitutive relations.

Length scales and parameters
We must now assign the reference values, t c , r c , and p c . Our interest is in well testing and therefore we choose for the time scale of interest t c = one day = 8.64 × 10 4 s. From this, we can derive the characteristic radius r c around the well that will be affected during the well test: Note that due to this specific choice, numeric constants are absorbed in the length scale which will simplify our expression ru tot = 1 in Section 3.2. The reference pressure can be obtained independent of our choice of t c from: In summary:  Table 2.

Model equations in dimensionless form
Using R = r r c , T = t t c , U tot = u r c t c , and P w = p w p c , the conservation equation for water (11) becomes: where we set c w = p ccw and D cap =D cap Similarly, we rewrite the equation for the sum of oil and water conservation (6) as follows: We use the capillary pressure P cap = P o − P w to eliminate P o from the Darcy (9) and find: which yields the dimensionless expression for u tot : where we have usedD cap = D cap r 2 c t c and the dimensionless mobilities λ w and λ o : This means that we have according to Eqs. 17 and 18: for the (dimensionless) mobilities and according to Eq. 20: ; due to our choice of scales, this combination now equals 1, which means that the water and oil mobility reduce to: where the ratio μ o μ w arises due to our scaling procedure. The injection flux Q inj 2πH t c r 2 c from Eq. 13 also equals 1 due to our choice of scales. The only parameters left are the combinations p cb = P cb p c and c w =c w p c and c o =c o p c , which are given in Table 3.
We summarize the model for easy reference.

Similarity transformation
We express (32)-(34) in terms of η using the similarity transformation (r, t) → (r 2 /2t = η). The product (ru tot ) is used as a dependent variable, i.e., the product (ru tot ) depends on η (see Eq. 44). Physically, this product represents the flux, i.e., the flux depends on the variable η (whereas the velocity depends on (r, t)). Using this transformation, Eq. 32 yields: Similarly, (33) becomes: Finally, (34) yields: Note that p w ∼ ln η for small η, which means that the first term on the RHS of Eq. 44 yields a finite, nonzero contribution.
The boundary conditions at the injection well r = r w are mapped to η = η w (t), where: We consider a well with a radius that is much smaller than the dimensions of the outer region. Furthermore, we are interested in the intermediate time behavior. The well radius is R w = 0.05 m, then the behavior of the solution for T > 30 min, which corresponds to t = 0.02, gives We are interested in the behavior of the solution far away from the well, i.e., η ≈ 5 and therewith η w (t) 5; and therefore, we use the approximation η w (t) ≈ 0. This means that we use the following boundary conditions as an approximation: In Section 3.2.5, we will analyze this approximation for the boundary condition quantitatively and we will see that the error is less than 1%. The initial condition (at t = 0) and the boundary condition (at r → ∞) for S w are both mapped to η → ∞, i.e.: The initial condition for the pressure in the compressible case is mapped to η → ∞, i.e.: We will solve problem (42)-(49) for saturation-dependent capillary diffusion in Section 3.1 and for constant capillary diffusion in Section 3.3. We focus mainly on the incompressible case where c o = c w = 0.

Summary of the self similar (in)compressible problem
Equations (42)-(44) can be rewritten as a set of four coupled nonlinear ordinary differential equations, i.e.: where f 3 is given by Eq. 44, f 2 is given by Eq. 43 and f 1 is given by Eq. 42. Notice that Eq. 42 is converted into two first-order equations. The functions f 1 , f 2 , f 3 are given by the following expressions: We have two conditions at η = 0 (see Eq. 47) and two conditions at η → ∞ (see Eqs. 48 and 49). However, our analysis is for incompressible flows (see [2] Chapters 5 and 7); this analysis is valid for t > ϕμc eff r 2 c /k ≈ 1s for Fig. 1 Water saturation versus η. The figure shows the behavior in the region I, i.e., 0 ≤ η ≤ η , region II, i.e., η ≤ η ≤η, and region III, i.e.,η ≤ η ≤ η f . In region II, we use the numerical solution, whereas in regions I and III, we use an analytical approximation to solve the equations (see text) our parameters. Therefore, compressibility effects can be disregarded for processes beyond 1 s.

Analysis in the incompressible case (c w = c o = 0)
Setting c w = c o = 0, we have f 2 = 0, which means that: Furthermore, f 1 now only depends on S w and dS w dη = S w and not on p w . Consequently, the equation for S w is decoupled from p w now, which means that we can first solve the second-order ODE for S w and use the result to find the pressure profile p w . We do have the following problem though. Close to S w = 1 and S w = S wc , D cap becomes very small because This means that numerical integration close to S w = 1 and S w = S wc is difficult, due to the presence of D cap in the denominator of f 1 (see e.g. [4], [5] for similar problems). For this reason, we split the domain in three regions (see Fig. 1): Region I Small η and S w ≈ 1 : S w ≤ S w ≤ 1 and 0 ≤ η ≤ η Region II Intermediate η and S w :S w ≤ S w ≤ S w and η ≤ η ≤η Here, η = η f is the endpoint of the region where S w > S wc , i.e., S w (η f ) = S wc and for η ≥ η f we have S w = S wc 2 . The parameters S w andS w are numerical parameters that are chosen as close to S w = 1 and S w = S wc as possible; several consecutive values are chosen until no visible difference between the corresponding solutions of the model equations is observed. Once the parameters S w andS w are fixed, all the other parameters (η ,η, η f ) follow as part of our solution to the model equations.
We will solve our problem analytically in regions I and III and we will use numerical integration in Region II as discussed in more detail in the following Sections 3.2.1-3.2.3.

Small η and S w ≈ 1: S w ≤ S w ≤ 1
First, we choose a value S w close to one; we will investigate the influence of this (numerical) choice later. We then use a linear approximation for our solution: where the condition S w (0) = 1 fixes η in terms of S w : The value of S w is determined as follows: a value of S w is picked and adjusted such that the total mass is conserved in the entire domain.

Intermediate η, S w ≥ S w ≥S w
We integrate numerically from η = η , where we have S w = S w and S w = S w until η =η, whereη is chosen such that S w atη is close to S wc . We will investigate the influence of this (numerical) choice later.

Large η ≈ η f and S w ≈ S wc :S w ≥ S w ≥ S wc
Around the toe, η = η f , we have S w ≈ S wc and consequently the oil mobility λ o ≈ 0. This means that we have: We assume power law behavior of S w around η f where the exponent p remains to be determined. Substitution of Ansatz (59) into Eq. 42 and balancing different powers of S w − S wc , we find (see Appendix C): where K 1 and η f are constant. This means that S w connects continuously to the initial state S w = S wc at η = η f (i.e., no shock). The derivatives of S w though do blow up near the toe.
The constants K 1 and η f are determined by matching the analytical solution in region III to the numerical solution in region II at η =η.

Global mass conservation
The unknown parameter S w has to be chosen such that we select the correct saturation profile, i.e., such that mass is globally conserved. The total amount of mass injected until time T , m in equals: Using cylindrical coordinates for the triple integrals, we have: In the incompressible case, we have constant water density; furthermore, we use R = r r c and T = t t c to obtain: Using η = r 2 2t ⇒ dη = r t dr, we find: Due to our choice of scales, we have furthermore, we have S w = S wc for η ≥ η f , which means that Eq. 67 reduces to: The value of S w is adjusted until (68) is satisfied. Figure 2 shows the influence of the choice of S w on the results. For the scale used to display the figure, we observe that the solution with S w = 0.95 slightly deviates from the other three curves, which appear (almost) completely on top of each other. This means that our solution S w (η) is practically independent of our choice of S w if S w is chosen large enough, i.c., above 0.97. Note furthermore that for η < 0.11, we have S w > 0.99, which means that the error on the boundary condition due to approximating η w (t) by zero is in the worst case 1%.

Saturation-dependent capillary diffusion results
In Fig. 3, we explore the effect of different bubbling pressures P cb on the saturation profile. For small values the curves approach a steep shock-like behavior. For larger values of the bubbling pressure, which acts as a parameter for enhanced capillary diffusion, the saturation profiles are broadened. In the limit of zero capillary pressure, we find the Buckley-Leverett solution as shown in Fig. 4, which also shows that capillary diffusion broadens the front. Figure 5 shows the water pressure profile. Our solution is only valid for t ≥ 2 · 10 −2 as discussed in Section 3. This occurs frequently in the case of well testing, where the very early data are ignored (see Dake [2], chapter VII, p 159 ff).
In this period, the total pressure increase of p ≈ 14, P = pp c ≈ 8.5 bar is limited due to the practical duration of the test. η shows that small differences can be observed for S w = 0.95 and S w = 0.96, but in practice no discernible differences for S w ≥ 0.97. Therefore, we use S w = 0.97 in our computations

Constant small capillary diffusion: matched asymptotic expansions
In general, the capillary diffusion is saturation dependent, i.e.: However, in this section, we will assume a constant (small) capillary diffusion coefficient; we will take a few representative values for D cap and compute the corresponding saturation profiles. These saturation profiles are compared with the profiles (with a saturation-dependent capillary diffusion) obtained in Section 3.2.5. We take incompressible conditions for water and oil, i.e., c o = c w = 0. This means that Eq. 43 yields: Furthermore, we set D cap = ε constant; Eq. 42 yields: where η f ≈ η s (the location of the shock in the unperturbed profile). We will use the method of matched asymptotic expansions to obtain an approximation to the solution of Eq. 71. We solve the outer problem to obtain S out w (η) in Section 3.3.1 and we solve the inner problem to obtain S in The total solution of Eq. 71 is then given by: Fig. 3 Dependence of S w (η), for S w = 0.97 for different values of the bubbling pressure P cb . We use P cb : 5, 10,20, 50, and 100 times the bases case value of 4.3 ×10 4 P a. The sharpest profile corresponds to P cb = 5 times the base value. For larger values of P cb , the front gets more disperse and the whole profile gets broader Fig. 4 Comparison of the saturation profile S w (η) for S w = 0.97 to the Buckley-Leverett solution. We observe that the overall mass balance is satisfied where the matching condition states that the inner limit of the outer solution = the outer limit of the inner solution ( [13]). The matching condition Water pressure p w (t) at the inlet (R w = 0.1 m) for 0.02 ≤ t ≤ 5. The lower limit is chosen according to the arguments at the end of Section 3.2.5. The upper limit is taken as t = 5, being a reasonable limit of the duration of the test. In this period, the total pressure increase of p ≈ 14, P = p · p c ≈ 8.5 bar is limited due to the practical duration of the test determines S m w for η < η s and the matching condition determines S m w for η > η s .

Solution of the outer problem
In the outer problem, we neglect contributions proportional to ε, and thus Eq. 71 reduces to: The solution of Eq. 76 is given by the Buckley-Leverett profile: where we use that ru tot = 1; see Eq. 70. For η > η s we have S w = S wc . Furthermore, the shock saturation S s w is obtained via Welge's tangent construction: Details of the derivation of Welge's tangent construction are outlined in the Appendices 1 and 2 for easy referencing for the general audience. Mathematical derivations can be found in [12] and [7]. Using Eq. (77), we obtain η s ≈ 5.551. Later on, we need η f for which we use mass conservation. Notice, indeed that we have: due to mass conservation; see Eq. 68.

Solution of the inner problem
The shock at η = η s is spread out and replaced by a front in the inner region small. For the inner problem, we define the new variable ξ , i.e.: which yields the following inner problem (to lowest order in ε): where we used η = η s + O(ε) in the inner region. We have an additional boundary condition due to the matching condition: Furthermore, we have the trivial condition: This means that our total solution (outer solution plus inner solution minus matching saturation) becomes: because S out w = S wc = S m w (see Eq. 83) for η > η s . Integrating (81) once yields: i.e., where we set the constant K = dS in w dξ | S sc and We use separation of variables to solve (86) and find We use the matching condition (82) to determine the constant K. The integral: has to diverge due to the matching condition (82); ξ = ∞ if and only if S in w = S s w . Note that the integration interval is finite and that g(u) is continuous; the integral only blows up if the denominator equals 0, i.e., if g(u) + K = 0. Notice that we have g(S s w ) = 0; this means that we need K = 0 in order to have a divergent integral around S w = S s w , i.e.: Equation (90) is integrated numerically to obtain ξ(S in w ) and inverted to obtain S in w (ξ ). Fig. 6 Inner, outer and total solution of the saturation profile for ε = 10 −2 . Consecutively the solution over the entire η range, and a zoom in the front region are given. At η ≈ 4.8 the transition from the outer to the inner solution occurs width of the front, which is approximately given by η f − η s (see Fig. 6).

Numerical results
We consider a fully coupled, implicit numerical solution approach based on finite elements, which is solved with the mathematical module of COMSOL to solve the model equations in weak form. We consider the spatial domain 0 ≤ x ≤ L of length L = 100m, where the Dirichlet boundary condition is taken at the production side, x = L. The grid size in the numerical simulation is 0.001 m, which corresponds to a spatial resolution of 100,000. This is fine enough to resolve the salient features. The reservoir parameters are given in Tables 2  and 3.
Numerical results for the base case are shown in Fig. 9. Here, the numerical simulation confirms the mentioned analytical results. As shown in Fig. 9, for smaller capillary diffusion (lower P cb ), the front is steep, whereas for higher capillary diffusion, the front is smeared out. Fig. 7 Inner, outer, and total solution of the saturation profile for ε = 10 −3 . On the left is the solution over the entire η range; on the right is a zoom in the front region. It is difficult to observe the transition from the inner to the outer solution due to the small value of ε 5 Discussion Figure 10 compares the analytical and numerical solutions of problems (32)-(41). The solutions practically coincide over the entire domain except where the downstream solution connects to the initial condition.This is due to the fact that capillary diffusion vanishes at this point (toe). This causes power law behavior for the saturation S w , i.e., S w is a continuous function of η at the toe, but the For smaller values of the capillary diffusion coefficient, i.e., less than five times P cb , numerical problems occur at the injection point due to the vanishing capillary diffusion. Figure 11 compares the analytical solution of problems (32)-(41) with the analytical solution of problems (71)-(75), i.e., the problem with a saturation-dependent capillary diffusion to the problem with a constant diffusion Fig. 9 The saturation profile S w (η) is given for different values of the bubbling pressure P cb . We use P cb : 5, 10,20, 50, and 100 times the bases case value of 4.3 ×10 4 P a. The sharpest profile corresponds to P cb = 5 times the base value. For larger values of P cb , the front gets more disperse and the whole profile gets broader. For smaller values of P cb , COMSOL is not able to find a solution Fig. 10 Combination of Figs. 3 (saturation profile obtained by the analytical method) and 9 (saturation profile obtained using COMSOL) for P cb 10 and 100 times of the base value. We observe that both methods yield the same saturation profile coefficient of comparable magnitude. Close to the injection point, the solutions more or less coincide. Close to the toe, we observe that the problem with constant diffusion has a smooth saturation profile, whereas the solution with the saturation-dependent diffusion coefficient shows power law behavior. Note that a very small constant diffusion coefficient admits a viable solution, which can be obtained with the method of matched asymptotic expansions.

Conclusion
The equations for radial Buckley-Leverett flow can be reduced to a first-order ODE using a similarity coordinate Fig. 11 Pressure p w as function of η obtained by the analytical method and using COMSOL for P cb 10 times the base value. We observe that both methods yield the same pressure profile until η = 5. Note that for η > 5, the time t < 2 × 10 −3 at r = r w , which is outside the range of validity of the analytical solution (and outside the range of interest) η = r 2 /(2t) transformation of the model equations. The resulting ODE can be solved in a similar way as the 1D problem. The solution consists of a rarefaction wave, followed by a shock to the constant initial state. When we include capillary diffusion, we obtain the Rapoport-Leas equation, which can be reduced to a second-order ODE using the same similarity transformation. This problem can be solved with a combination of analytical and numerical techniques. In this case, we observe that the shock is replaced by a curve that shows power law behavior; the saturation profile is continuous, but the derivatives blow up at the toe. We compare the similarity solution of the second-order ODE to a numerical solution of the initial value problem over the entire domain. The solutions for the saturation profile show excellent agreement. As long as we stay away from the connate water saturation, where the capillary pressure blows up, the solutions for the water pressure also show excellent agreement. If we replace the saturation-dependent capillary diffusion by a constant diffusion coefficient, we obtain a simpler problem, which can be solved with the method of matched asymptotic expansions. In this case, the shock is replaced by a smooth profile.
Our analysis is for incompressible flows (see [2], Chapters 5 and 7). It is shown in [2] that this analysis is valid for t > ϕμc eff r 2 c /k corresponding with our parameter values to approximately 1 s. Therefore, compressibility effects can be disregarded for processes beyond 1 s. For other parameters, e.g., for gas shales, this analysis can be used including the self-similar transformation for t < ϕμc eff r 2 c /k. This results in a system of four coupled ODEs in terms of η. In general, the similarity transformation reduces the number of independent variables from two to one; this may considerably simplify the equations used to interpret the results of this type of well test problems.
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:// creativecommonshorg/licenses/by/4.0/.