Published online 24 January 2007
Published in Vadose Zone J 6:93-104 (2007)
DOI: 10.2136/vzj2006.0024
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
An Improved Semi-Analytical Solution for Verification of Numerical Models of Two-Phase Flow in Porous Media
Radek Fu
íka,
Ji
í Miky
kaa,
Michal Benea,* and
Tissa H. Illangasekareb
a Department of Mathematics, Faculty of Nuclear Science and Physical Engineering, Czech Technical University in Prague, Trojanova 13, 120 00 Prague, Czech Republic
b Center for Experimental Study of Subsurface Environmental Processes (CESEP), Division of Environmental Science and Engineering, Colorado School of Mines, 1400 Illinois Street, Golden CO 80401, CO
* Corresponding author (benes{at}km1.fjfi.cvut.cz)
Received 13 February 2006.
 |
ABSTRACT
|
|---|
A closed-form solution for one-dimensional two-phase flow through a homogeneous porous medium is presented that is applicable to water flow in the vadose zone and flow of nonaqueous phase fluids. The solution is a significant improvement to the one originally presented by McWhorter and Sunada, allowing the analysis of wetting phase entry saturations ranging from residual to full. Our aims are to provide a detailed analysis of how the solution to the governing partial differential equation of two-phase flow can be obtained from a functional integral equation arising from the analytical treatment of the problems and to present an improved algorithm for the implementation of this solution. The integral functional equation is obtained by imposing a set of assumptions for the boundary conditions. The proposed method can be used to obtain solutions that incorporate a wide range of saturation values at the entry point. The semi-analytical solution will be useful in the verification of vadose zone flow and multi-phase flow codes designed to simulate more complex two-phase flow problems in porous media where capillary effects must be included.
|
INTRODUCTION
|
|---|
COMPLEX multi-dimensional numerical models of multi-phase flow through porous media such as those described by Helmig (1997), Mikyka et al. (2004), or Mikyka and Illangasekare (2005) require verification to assure that the governing equations are solved correctly and the codes do not contain programming errors. This step of code verification is a necessary step in modeling protocols used in practice (e.g., Anderson & Woessner (2002)). Code simulations are compared to closed-form analytical solutions of the governing equations to estimate numerical errors and other inaccuracies of numerical schemes. Two well-known one-dimensional solutions of the two-phase flow equations include the Buckley-Leverett solution of flow without capillary effects (e.g., described by Helmig (1997) and LeVeque (2002); or see references in McWhorter and Sunada (1990)), and the exact integral solution derived by McWhorter and Sunada (1990) with subsequent discussions by Chen et al. (1992), McWhorter and Sunada (1992), and Fucík et al. (2005), which includes both advective and capillary effects. In this paper, we discuss the exact integral equation for the wetting-phase saturation obtained by McWhorter and Sunada (1990). This equation must be numerically integrated to yield the saturation distribution along the length of the soil column, and a value for entry saturation is needed as an input boundary condition. The solution to the problem as presented by McWhorter and Sunada (1990) has limitations in those situations where the entry wetting-phase saturations are high.
As numerical models are designed to simulate conditions that include high entry wetting saturations (e.g., wetting front propagation, water flooding for enhanced recovery), analytical models used for code verification should have the capability to simulate this flow condition. We present an improvement of the technique, that allows the exact integral solution to be reliably obtained under conditions where the McWhorter and Sunada (1990) approach fails to converge. Our approach provides insight into the solution behavior and explains the limitations of the previously known method of resolution of the integral equation. This generalized approach is applicable to unsaturated zone (water - air) or saturated zone (water - NAPL) models when both phases are assumed to be incompressible. We perform a series of qualitative and quantitative computations that show our algorithm agrees with previously obtained results while demonstrating the improved performance.
|
TWO-PHASE FLOW MODEL
|
|---|
In this section we introduce basic notation and set up the governing equations.
Transport equation with capillarity
We consider a one-dimensional problem describing flow of two incompressible and immiscible fluids through a porous medium where the non-wetting phase (indexed n) is horizontally displaced by the wetting fluid (water, indexed w) (therefore neglecting the influence of gravity). Darcy's law, when written for each of the fluid phases, has the following form:
 | [1] |
where q
, 
, and p are the flux, mobility and pressure of the phase , respectively, where we use 
{w,n}. The -phase mobility is defined as
 | [2] |
where k is the permeability and µ is the dynamic viscosity of phase (Bastian, 1999). The total flux qt is defined as the sum of the fluxes of each of the phases (qt = qw + qn) The capillary relation, pc = pn – pw, with a given function pc = pc(Sw) of the effective wetting phase saturation Sw, links the wetting and the nonwetting balance equations. The effective saturation of the phase is defined by
 | [3] |
where Swr and Snr are the residual wetting and non-wetting phase saturations, respectively, and s is the saturation of phase . The effective saturation is always between 0 and 1, which simplifies the description of the dependent variable by the definition Sw + Sn = 1 (Helmig, 1997).
Introducing the wetting and nonwetting phase fractional flow functions
 | [4] |
and diffusivity functions
 | [5] |
 | [6] |
we obtain the expression for the -phase flux as
 | [7] |
The mass-balance equation has the following form (the fluid mass density is assumed constant):
 | [8] |
where 
is the porosity.
The two-phase flow equation is obtained by substituting (7) into the mass-balance equation (8) to yield
 | [9] |
which corresponds to equation (2) in McWhorter and Sunada (1990). Substituting Sw = 1 – Sn, equation (9) becomes
 | [10] |
Equations (9) and (10) are equivalent and they can be used in the formulation of either a wetting phase or a nonwetting phase displacement problem. A general form of the two-phase flow equation is given as
 | [11] |
where we obtain equations (9) or (10) using respective substitutions for the functions f, D and S. For the one-dimensional displacement problem, the initial and boundary saturations (at x = 0 and x 
+
) must be defined.
McWhorter and Sunada (1990) presented the closed-form analytical solution for equation (11) for both one-dimensional and radial displacement. The radial displacement problem presented in McWhorter and Sunada (1990) is not discussed in this paper because a different type of the integral equation arises in that case.
We will discuss conditions under which the flow equation can be solved analytically to provide a simple one-dimensional benchmark solution for verification of more complex two-phase flow codes.
Transport equation without capillarity
The last term in equation (11) vanishes when the capillary effects represented by the term dpc(S)/dS in the diffusivity function D(S) are neglected, resulting in the Buckley-Leverett equation for two-phase flow (Helmig, 1997)
 | [12] |
The first-order hyperbolic equation (12) represents a limiting case for equation (11) when D(S) 0. The boundary and initial conditions are defined as
 | [69] |
The analytical solution to equation (11) is based on the method of characteristics and is given by
 | [13] |
The function f(S) has an inflection point, so that the solution is implicitly given by equation (13) for S0 
S St (inverted saturation profile), where St is the Welge tangent saturation (or post-shock value; see LeVeque (2002)) that is determined from the relation
 | [14] |
Capillary and relative-permeability model functions
Denoting the intrinsic permeability of the medium by k, the relative permeability for the wetting phase is defined by krw = kw/k and the relative permeability for the nonwetting phase by krn = kn/k. The functions krw, krn and the capillary-pressure expression are used in the following formulations.
The Brooks-Corey model (Brooks and Corey, 1964) relating capillary pressure pc to saturation is given by
 | [15] |
where and P0 are parameters characterising the soil and phase properties; P0 is called the entry pressure.
Application of the Burdine (1953) formulation to the Brooks-Corey model results in relative permeability functions for the wetting and non-wetting phases in the form
 | [16] |
The van Genuchten (1980) capillary-pressure pc expression is given as
 | [17] |
where the parameters m and n are often related by m = 1–1/n.
Application of the Mualem (1976) relative-permeability functions to the van Genuchten model results in
 | [18] |
|
QUASI-ANALYTICAL SOLUTION
|
|---|
Usefulness of a benchmark solution depends on its relative ease of use. We therefore consider the possibility of improving a closed-form solution to equation (11) based on the approach originally presented by McWhorter and Sunada (1990). In this section, The closed-form solution of McWhorter and Sunada (1990) is presented in this section to provide a basis for improvement. An enhancement that enables a wider range of entry saturations to be considered than the McWhorter and Sunada (1990) approach is presented in Section 4.
Problem formulation
A quasi-stationary solution of (11) under a particular set of conditions is presented. We assume that for all x (0, +) and t (0, +)
 | [19] |
 | [20] |
 | [21] |
with S0 > Si. If S0 < Si, we must use the other formulation (i.e., (10) instead of (9)) or introduce a fractional flow function, Fnw, as in McWhorter and Sunada (1990). An advantage of our approach is that we use the same code to compute the wetting as well as the non-wetting phase displacement problem simply by defining respective functions f and D appropriately.
The displacing phase (indexed ) is introduced to the column at x = 0 with volumetric flux given by
 | [22] |
where A > 0. The function g(t) must have the form g(t) = t–1/2, as will be shown in Section 3.4. The displaced phase flux at the inlet (x = 0) and the outlet (x +) are unknown. The boundary at x + is semi-permeable, characterized by a scalar coefficient R [0,1], where R = 0 implies that the boundary is impermeable and R = 1 implies no resistance to the flow at the boundary (unidirectional flow).
It follows from (7) and from the assumption of incompressibility of both phases that
 | [23] |
Therefore, qt is spatially uniform but may vary with time, i.e., for all x 0 and t 0 we get
 | [24] |
The total flux achieves its maximum value, qt(t) = Ag(t), at the outlet when R = 1. On the other hand, the total flux vanishes when R = 0. This represents bidirectional displacement where the displaced fluid is draining only at the inlet (x = 0).
McWhorter and Sunada (1990) considered only the limiting cases of R = 0 and R = 1, but we note that the approach is valid for R [0,1]. The displacing phase is thus injected in the counter-current flow direction of the total flux qt.
By combining equations (7), (22), and (24), we obtain the relationship
 | [25] |
Basic assumptions
We assume that the solution exists in the form S = S(), where
 | [26] |
This substitution is possible, provided the basic assumption
 | [27] |
This assumption allows the dependence S = S() to be inverted so that = (S).
Assuming that S as a function of is sufficiently smooth, partial differentiation of (26) yields
 | [28] |
and
 | [29] |
where g'(t) and '(S) stand for the derivatives dg(t)/dt and d(S)/dS, respectively.
Expression for Function F
We define the fractional flow function, F = F(t,x), as
 | [30] |
and introduce the normalized fractional flow function, 
(S), where
 | [31] |
Combining equations (29), (30) and (31) allows for the redefinition of F in terms of S,
 | [32] |
Stationary differential equation
Introduction of expression (32) for F to equation (11) yields
 | [33] |
Substituting equation (26), (28), and (29) into (33), we obtain the equation
 | [34] |
where F'(S) stands for dF(S)/dS.
Only a function g of the form
 | [35] |
allows the removal of the explicit time dependence of the terms in equation (34) because the term g'(t)/g3(t) equals C1. The value of C1 is arbitrary as long as it is negative. As the value of A depends on C1, it is therefore possible to choose for instance C1 = –1/2.
Differentiating (34) with respect to S and substituting equation (32) for F(S) yields the second-order ordinary differential equation
 | [36] |
where F''(S) stands for d2F(S)/dS2.
The boundary conditions for the ordinary differential equation (36) are F(S0) = 1 and F(Si) = 0, which follow from (19) and (21), respectively. Note that matching the initial condition (21) to the condition F(Si) = 0 is possible only if g(0) = +. This implies that the only possible form of the input flux function g(t) is g(t) = t–1/2, which in turn implies that C2 = 0 by (35).
Moreover, the boundary condition defined in (19) gives us F'(S0) = 0. However, the problem is not overdetermined because the condition F(Si) = 0 will be used to establish the relationship between A and S0.
Solution of the transport equation
Once the function F(S) is known, we can derive the inverted form of (11) from the relation
 | [37] |
which is in a form similar to the Buckley-Leverett analytical solution (13), given by
 | [38] |
This expression is valid for all values of S [Si,S0] because the function dF(S)/dS can be inverted as a consequence of the basic assumption expressed in (27).
In order to demonstrate the relationship between the Buckley-Leverett and the McWhorter-Sunada exact solutions, we define the Buckley-Leverett fractional flow function FBL as follows:
 | [39] |
where St is the Welge tangent saturation given by (14). It is obvious that the function FBL does not satisfy the basic assumption (27) due to the relationship (37) and the linear part of FBL. However, the solutions (13) and (38) are formally the same when FBL is substituted for F.
|
INTEGRAL SOLUTION
|
|---|
Derivation
The equation (36) cannot be solved until the relationship between A and S0 is determined. Following McWhorter and Sunada (1990), we integrate (36) twice and include F'(S0) = 0 and F(S0) = 1 to obtain
 | [40] |
The condition F(Si) = 0 allows for the establishment of the relationship between A and S0 as follows
 | [41] |
The integral equation (40) can be rewritten by means of (41) into the form
 | [42] |
Differentiating this integral equation, we obtain the function F'(S)
 | [43] |
The magnitude of the diffusion term D(S) does not influence the function F because multiplicative constants in the term D(S) can be cancelled in (42) as well as in (43). It affects only the value of A in (41).
Iteration scheme
In agreement with McWhorter and Sunada (1990), the unknown function F(S) is computed from the integral equation (42) by iteration. The iterative process is as follows:
 | [44] |
As in McWhorter and Sunada (1990), we suggest using F0 
1 as a first guess. The function Fk is considered to be the solution of (40) when successive iterations are sufficiently small in a norm. In our case, we use the L norm and terminate the iterative process when
 | [45] |
The integrals in (44) are evaluated numerically, therefore the exact solution is often referred to as quasi-analytical solution. The iterative process is rapid and convergent for all values of S0 in case of the bidirectional flow (R = 0). However, serious difficulties occur when S0 and R are close to one as the following first-iteration analysis demonstrates.
Test problems
Table 1 gives the parameter values that are used to demonstrate our approach. Water is the wetting phase in our computational experiments, while various realistic or theoretical non-wetting liquids are used. The term NAPL stands for non-aqueous phase liquid and DNAPL is denser-than-water NAPL.
The first setup consists of the use of Brooks-Corey model functions and artificially selected values of the soil parameters (see Helmig (1997)). In Setup 1, we choose µn = 0.020. Note that efficiency of our approach increases with decreasing µw/µn.
Setups 2 and 3 contain the parameters of laboratory test soils used in our ongoing research (sand #30 as in Turner (2004), p.43) and a test NAPL Soltrol 220.
First iteration
The iterative scheme presented by McWhorter and Sunada (1990) exhibits unsatisfactory behavior when the denominator in the integrand D(v)/(F(v) – (v)) in (42) approaches zero. This happens when both S0 and R are close to 1. We offer an analytical justification of this phenomenon. It can be shown that F(S) > (S) for all S (Si,S0). Hence, this relationship must stand for all approximations Fk of the function F.
The first iteration of the function F is obtained by substituting F0 1 into the right hand side of integral equation (44). Using Brooks-Corey model functions (16) and (15), the first iteration F1 for the Sw formulation is expressed analytically as follows:
 | [46] |
The second iteration of the function F cannot be computed by substituting F1 into the right hand side of integral equation (44) for certain values of S0, because the function F1 intersects the function and the integrand D(v)/(F(v)–(v)) becomes unbounded. This is illustrated in Figure 1
, where we set Si = 0, R = 1 and S0 {0.5,0.7,0.9,1}.
Equation (46) implies that the first iteration of F is only dependent on . We studied the behavior of F1 with respect to for common values of , finding that does not affect the formation of the instability of the iterative process in any remarkable way.
Although the values µw and µn do not influence F1, they have an important impact on the instability formation of the iterative process by the function , as shown in Figure 1. As an illustration we study the non-wetting phase displacement problem with R = 1 and Si = 0, i.e. fw. Whenever the function intersects the first iteration F1, a singularity in the integrand D(v)/(F(v)–(v)) in (42) occurs.
The viscosity ratio, M = µw/µn, is the key parameter that affects the stability of the iterative process because it shifts the inflexion point of the function towards S0 or Si (see Figure 1). The singularity may occur at any saturation in the interval (Si,S0), not just in the vicinity of S0. The other parameter that influences the formation of the instability after the first iteration is the initial saturation, Si, which appears in both the function and F1.
In displacement problems involving NAPLs that are less viscous than water, as is demonstrated in Tables 2, 3 and 4, the original iterative process fails for values of S0 near 1. In order to demonstrate limits of the functionality of the original iterative scheme, we introduce the critical value denoted by S0* that represents the lowest value of S0 for which the original iterative scheme (44) fails after the first iteration. We determine S0* experimentally by bisectioning an interval [S0(1),S0(2)], where S0(1) and S0(2) corresponds to S0 for which the iterative scheme (44) is stable and fails, respectively. In the next step, we test if the scheme (44) works for S0(3) = (S0(1) + S0(2))/2 and then we shift the left (S0(1) 
S0(3)) or the right (S0(2) S0(3)) boundary of the interval, so that the scheme (44) works for S0 = S0(1) and fails for S0 = S0(2). Iterations are terminated when the length of the interval [S0(1), S0(2)] is below a prescribed tolerance.
Results given in Tables 2, 3, and 4 suggest that the instability issue of the original process is not peripheral for highly viscous non-wetting fluids. For example, the original iterative process fails for values of S0 greater than 0.82 in the case of the test NAPL Soltrol 220 (Setup 2, Table 1), which is more viscous than water, µn = 0.0035 kg m–1 s–1, so that M = 0.286.
Based on Figure 1, the original iterative process will fail for the values of S0 S0*.
Modified integral equation
We propose the following modified method to avoid unstable behaviour of the numerical iterative process. Denoting the principal part of the integrand in (42) as G = D/(F – ), we can rewrite equation (42) as
 | [47] |
which allows us to deduce two types of iterative schemes; method A, given by the scheme
 | [48] |
and method B, given by the scheme
 | [49] |
We suggest using G0 = D/(1 – ) as the initial guess, which is equivalent to the case where F0 1, as proposed by McWhorter and Sunada (1990).
The integrals in (48) and (49) are evaluated numerically, taking advantage of the form of the integrand as follows. Let {Gj}j=0M be an equidistant discretization of the function G in the interval [Si,S0], defined as Gj = G(Si + j h), where h = (S0 – Si)/M. The numerical solution of the integral equations (48) and (49) requires computation of the integral
 | [50] |
We suggest introducing partial numerical integrals aj and bj, given as
| [51] |
.
Linear interpolation of {Gj}j=0M in the interval [Si,S0] allows the value of aj and bj to be expressed as
 | [52] |
and
 | [53] |
The integral (50) in the modified iterative schemes (48) and (49) is approximated by il for discrete values of saturation (S = Si + l h) by
 | [54] |
Since both F(Si) = 0 and (Si) = 0 by definition, it follows from G = D/(F – ) that the value of G(Si) is undefined (note that D(Si) > 0 if Si > 0). The value Gk0 = G(Si) = 0 is used in the scheme for all k.
Application of the discretization Dl = D(Si + l h) and l = (Si + l h) and using the expression (54) in the method A (48) yields
 | [55] |
Analogously, the method B (49) is given by the scheme
 | [56] |
The presented form of the iterative scheme benefits from the type of the integral equations (42), (48) and (49). In all iterations, only M numbers aj and bj need to be evaluated.
Values of the functions Fk are computed from
 | [57] |
as G(S) > 0 for all S (Si,S0). It is better to determine the first derivative F' based on the expression (43) in the form
 | [58] |
than using the numerical differentiation since the terms aj and I0 are already evaluated.
Behavior of the modified iterative scheme
In this section we focus on the unidirectional case with R = 1 and will illustrate the functionality of the modified method. We observe a monotone growth of successive estimates of G in all computations, i.e. Gk 
Gk+1 in [Si,S0], and fast convergence for all cases where the original iterative method succeeds (Table 5).
View this table:
[in this window]
[in a new window]
|
Table 5. Number of iterations required to obtain the function F and the value of A with accuracy  A = 10–15. In some situations with Si > 0, the modified iterative method B fails randomly.
|
|
The modified method converged for cases where the original method failed. In the modified method, successive estimates of G decreased in the L norm, but the number of iterations needed to reach a required precision of G increases considerably as both S0 and R approach one. Although there are negligible variations in successive estimates Fk in this situation, the value of A converges slowly, as shown in Figure 2
. Successive differences of the function F estimates decrease very rapidly in the beginning of the iterative process, while the value of A increases considerably. Values of F between 10–20 to 10–35 are necessary to obtain successive differences of the approximations of A below A = 10–15, which reaches the common computer round-off error.
We suggest using the difference between successive approximations of A as the stopping criterion for the iterative process. This is represented formally as
 | [59] |
We use the test models with highly viscous NAPL to demonstrate robustness of our modified iterative scheme in situations where the original iterative scheme fails even after the first iteration.
We found that the method B (49) fails randomly due to numerical division by zero, when Si > 0, because finite-precision evaluation of the fraction in (49) is indistinguishable from 1 for S very close to Si. If Si = 0, however, the process is stable because the diffusivity term D(0) = 0 lets G(S) vanish in the vicinity of Si. It is obvious that the value of G(Si) is undefined for all Si [0,S0], since by definition F(Si) = (Si) = 0. We suggest excluding the value of G(Si) from the discretization of the function G in the numerical computation because F(Si) = 0. Table 5 shows that the number of iterations required to reach machine precision of successive estimates of A increases as S0 1 for both method A and B. This is due to the extremely small difference between the function and F in the neighborhood of one, as noted by McWhorter and Sunada (1990). Results given in Table 6 and Figure 3
demonstrate how the function F approaches the Buckley-Leverett-based function FBL introduced in (39). Moreover, convergence also takes place for the first and second derivatives of F, i.e. F' FBL' and F'' FBL'' as S0 1.
The number of iterations increases as S0 1, because the integrals in the iterative scheme determined numerically become inaccurate as limited by the precision of the computer. More importantly, the limit function FBL does not obey the basic assumption that S is a strictly monotone function of and its second derivative FBL'' is discontinuous. Numerical experiments showed that the function F'', given as
 | [60] |
is bounded by FBL'' (see Figure 3). The convergence of F to FBL as S0 1 can be studied only through numerical experiments since analytical techniques are not available.
For S0 close to 1, a large number of iterations is needed to achieve convergence of A (see Figure 4
). Above a certain value of S0, the modified iterative process will not converge due to loss of numerical accuracy. However, estimates of the function F and its first and second derivative may converge even though A will not, since the fraction in (44)
 | [61] |
suppresses any effect of changing A on the function F.
View larger version (21K):
[in this window]
[in a new window]
|
Figure 4. Evolution of A in the modified iterative process, method A; test Setup 1, Si = 0. As S0 approaches 1, the iterative process requires higher number of iterations to reach convergence of A. In the very proximity of S0 = 1, the value of A is far from convergence even after 108 iterations. The situation for the method B is analogous.
|
|
The lower subfigures of Figure 3 indicate that the function G approaches the function FBL'' multiplied by a constant involving A2 (see (60)). Since FBL'' (S) = 0 for all S in (Si,St), this possible limit function of G as S0 1 does not solve the modified iterative schemes. This is due to D(S) 
0 in the interval (Si,St), since zero values of G are not admissible in the modified integral equation (47).
Consequently, the integral equation (42) cannot be solved numerically for values of S0 and R when they are too close to 1.
Limiting value of A
The convergence of A is an important part of the computational scheme, especially as it depends on S0. The iterative process may need a large number of iterations for A to converge if S0 and R are close to one, while the estimates of the function F vary negligibly. Therefore, we pursue the discussion of McWhorter and Sunada (1990), (1992) and Chen et al. (1992) concerning the limit
 | [62] |
In this section, we consider only the unidirectional displacement case when R = 1. McWhorter and Sunada (1990) claimed that the limit (62) is infinite as a consequence of F close to S0. However, this was questioned by Chen et al. (1992), claiming that the limit is always finite since the integrand
 | [63] |
is bounded as S0 1. In the reply to this comment, McWhorter and Sunada (1992) confirmed that the limit (62) is always finite because the integrand (63) is bounded for v = S0.
On the other hand, our work shows that the value of A increases without bounds as S0 approaches 1, as demonstrated in Figure 4. We extend our observations related to F'' approaching FBL'' as S0 1 as follows.
The term '(S) can be evaluated by combining (32) and (36) to yield
 | [64] |
We substitute this expression of '(S) into (29) to obtain
 | [65] |
The total flux condition (25) can be written in the terms of S0 only, as follows:
 | [70] |
This equation can be further simplified by employing the Sw formulation into
 | [66] |
and thus one can state
 | [67] |
The limit (67) is infinite for both the Brooks-Corey and van Genuchten models. That is
 | [68] |
which agrees with McWhorter and Sunada (1990). Note that the limit (68) is also infinity for the Sn formulation. This result implies that G must be unbounded at some value of S. It can be seen in Figure 3 that G grows dramatically as S0 1 in the region of St (the cusp at the front of the Buckley-Leverett shock). Since a cusp has an undefined second derivative, convergence to a cusp implies that the solution is unbounded in the vicinity of the cusp.
|
SOLUTION OVERVIEW
|
|---|
In this section, we demonstrate how the modified iterative methods using (48) and (49) can delineate the relationship between the McWhorter and Sunada and Buckley-Leverett analytical solutions. We perform computations for Setup 1 with Si = 0 with various values of R and S0. In order to compare the McWhorter and Sunada exact solution (37) with the Buckley-Leverett analytical solution (13), we use the value of A corresponding to the McWhorter and Sunada exact solution for R = 1.
Figure 5
shows how the cases of bi-directional displacement (R = 0, diffusive term only in (11)), partially unidirectional displacement (R = 0.8, both advective and diffusive terms in (11)), and unidirectional displacement (R = 1, both advective and diffusive terms in (11)) are related to the Buckley-Leverett solution of the advection equation (12). As S0 approaches 1, the diffusive term in (11) has less effect on the solution. Table 7 displays values of A for various combinations of R and S0.
View larger version (25K):
[in this window]
[in a new window]
|
Figure 5. McWhorter exact solutions (the method A) and Buckley-Leverett analytical solutions for various S0; test Setup 1, Si = 0. As S0 1, the unidirectional displacement solution (R=1) approaches the Buckley-Leverett solution, while the head of the bi-directional displacement solution (R=0) advances negligibly.
|
|
The modified iterative process allows solutions for strongly advective terms in (11), whereas the original iterative process fails in situations where the diffusive term is still significant. Since the critical value S0* for Setup 1 with Si = 0 and R = 1 is S0* = 0.69, solutions with R = 1 shown in Figure 5, except the case S0 = 0.6, are only obtainable by our modified iterative method.
|
CONCLUSIONS
|
|---|
The article is devoted to a detailed discussion of the benchmark solution described by McWhorter and Sunada (1990). We propose a reliable procedure for solving the implicit functional relationship that is the result of the analytical treatment of the advection-diffusion equation. This algorithm extends the use of the semi-analytical approach to a wider range of entry saturations than the original algorithm proposed by McWhorter and Sunada (1990). The use of our algorithm is limited by the round-off errors of the numerical computations and the number of iterations required for solution.
From our analysis, it follows that the original iterative method proposed by McWhorter and Sunada (1990) can be used to obtain solutions of the unidirectional displacement problem (R=1) only in a restricted interval of the entry saturations S0. The restricted interval can be determined by examining the first iteration. Our modified iterative method removes this restriction and offers a solution for larger range of entry saturations.
Method A (equation (48)) can be used to compute the solution for any admissible parameters except the values of S0 and R extremely close to 1 while method B (equation (49)) randomly fails if Si > 0. Therefore, the iterative method described by method A (equation (48)) can be used exclusively for use of the McWhorter and Sunada quasi-analytical solution.
The comparison of the McWhorter-Sunada fractional flow function F with the Buckley-Leverett fractional flow function, FBL, allows us to determine the limit of A as S0 1 and therefore to confirm the statement given by McWhorter and Sunada (1990), in contrast to the contentions of Chen et al. (1992) and McWhorter and Sunada (1992).
The practical value of our results is that they contribute to a detailed analysis of the analytical benchmark solution often useful for verification of more complex numerical models and in providing a tool for comparison under conditions of high wetting-phase saturations. Such a code verification was conducted by Mikyka and Illangasekare (2005) where this improved solution was used.
|
ACKNOWLEDGMENTS
|
|---|
The authors were partly supported by the project "Applied Mathematics in Technical and Physical Sciences" MSM 6840770010, by the project "Environmental modelling" KONTAKT ME878, and by the project "Jindrich-Necas Center for Mathematical Modelling" LC06052, all of the Ministry of Education of the Czech Republic, and by the National Science Foundation through the award 0222286 (CMG RESEARCH: "Numerical and Experimental Validation of Stochastic Upscaling for Subsurface Contamination Problems Involving Multi-phase Volatile Chlorinated Solvents").
|
REFERENCES
|
|---|
- Anderson, M.P., and W.W. Woessner. 2002. Applied Groundwater Modeling, Simulation of Flow and Advective Transport. San Diego: Elsevier.
- Bastian, P. 1999. Numerical Computation of Multiphase Flows in Porous Media. Habilitation. Christian-Albrechts-Universität Kiel.
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrology Paper 3:27.
- Burdine, N.T. 1953. Relative permeability calculations from pore-size distribution data. Tech. rept. Petroleum Transaction, AIME.
- Chen, Z.-X., G.S. Bodvarsson, and P.A. Witherspoon. 1992. Comment on ‘Exact Integral Solutions for Two-Phase Flow’. Water Resour. Res. 28:1477–1478.[CrossRef]
- Fuík, R., J. Mikyka, and T.H. Illangasekare. 2005. Evaluation of saturation-dependent flux on two-phase flow using generalized semi-analytic solution. Pages 25–37 of: Proceedings of Czech-Japanese Seminar in Applied Mathematics 2004, editors M. Bene, J. Mikyka, T. Oberhuber. Czech Technical University in Prague. ISBN 80-01-03181-0.
- Helmig, R. 1997. Multiphase Flow and Transport Processes in the Subsurface : A Contribution to the Modeling of Hydrosystems. Berlin: Springer Verlag.
- LeVeque, R.J. 2002. Methods for Hyperbolic Problems. Cambridge, New York, Melbourne, Madrid, Cape Town: Cambridge University Press.
- McWhorter, D.B., and D.K. Sunada. 1990. Exact Integral Solutions for Two-Phase Flow. Water Resour. Res. 26(3):399–413.
- McWhorter, D.B., and D.K. Sunada. 1992. Reply to "Comment on ‘Exact integral solutions for two-phase flow’" by Z.-X. Chen, G.S. Bodvarsson, and P.A. Witherspoon. Water Resour. Res. 28(5):1479–479.[CrossRef]
- Mikyka, J., and T.H. Illangasekare. 2005. Application of a Multiphase FlowModel for Simulations of NAPL Behavior at Inclined Material Interfaces. Pages 117–127 of: Proceedings of Czech-Japanese Seminar in Applied Mathematics 2004, editors M. Bene, J. Mikyka, T. Oberhuber. Czech Technical University in Prague. ISBN 80-01-03181-0.
- Mikyka, J., M. Bene, A. Turner, and T.H. Illangasekare. 2004. Development and Validation of a Multiphase Flow Model for Applications in NAPL Behaviour in Highly Heterogeneous Aquifer Formations. Pages 215–218 of: Proceedings on FEM MODFLOW and More, editors K. Ková, Z. Hrkal and J.Bruthans. Charles University in Prague.
- Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.[CrossRef]
- Turner, A.D. 2004. Behavior of dense non-aqueous phase liquids at soil interfaces of heterogeneous formations: Experimental methods and physical model testing. Numerical Computation of Multiphase Flows in Porous Media. Master's thesis. Colorado School of Mines, Golden, Colorado.
- Van Genuchten, M.T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Web of Science]
This article has been cited by other articles:
|
|
|
|
 
R. Fucik, J. Mikyska, M. Benes, and T. H. Illangasekare
Semianalytical Solution for Two-Phase Flow in Porous Media with a Discontinuity
Vadose Zone J.,
August 13, 2008;
7(3):
1001 - 1007.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
