1 Introduction
Multiphase flow is widely encountered in the fields of hydrology and environmental engineering sciences in the context of soil and groundwater pollutions of non-aqueous phase liquids (NAPL), such as chlorinated solvents; they constitute a large and serious environmental problem (Cohen and Mercer, 1993). In addition, in soil and groundwaters, they are subject to natural attenuation. Identification of pollution sources is difficult due to the fact that organic pollutants can rapidly migrate down to the bottom of the aquifer and/or along paths different from the water (Bano et al., 2009; Benremita and Schäfer, 2003; Bohy et al., 2004, 2006; Dridi et al., 2009). Among the different attenuation processes involved, such as dispersion, adsorption, volatilisation, chemical and biological destruction, biodegradation is often considered most important (Nex et al., 2006; Sinke and Le Hecho, 1999). Taking into account all of these mechanisms to quantify and characterise groundwater pollution by NAPLs, multicomponent transport equations require a fast and accurate resolution of three-phase flow primary and secondary variables.
Two different strategies are currently used to model multiphase flow in porous media (Helmig, 1997). Using the first approach, in the case of a NAPL migrating in an unsaturated zone, the water-oil-gas three-phase system can be modelled by three pressure equations obtained by introducing the Darcy velocities of each fluid phase into the individual mass balance equations. This system of pressure equations can be transformed into a pressure-saturation form by using relationships between the wetting phase saturation and the capillary pressure, which is defined as the pressure difference across the interface between a non-wetting and a wetting phase. The pressure of the second and third fluid phases can then be removed by expressing this pressure in terms of the saturation and pressure of the other phase. This latter approach is useful when phase disappearance occurs and saturation becomes zero. The system of partial differential equations describing a three-phase flow is highly non-linear due to the nature of the relative permeability and capillary pressure functions needed to close the system.
Alternative forms of these governing flow equations have been investigated to develop better computational algorithms. This has led to the fractional flow approach, which originates from the petroleum industry. In this approach, the total fluid flow describes the individual phases as a fraction of the total flow (Binning and Celia, 1999). Through the fractional flow formulation, the immiscible displacement of oil, gas and water can usually be expressed in terms of three coupled equations, namely a mean pressure equation (Nayagum et al., 2001, 2004) or global pressure equation (Antoncev and Monahov, 1978; Chavent and Jaffré, 1986) and two saturation equations.
In this article, a fractional three-phase flow formulation is chosen using a global pressure approach. Numerical simulation of three-phase compressible flow in a porous medium usually requires the knowledge of the corresponding three-phase data. Experimental values for relative three-phase permeabilities and capillary pressures are usually known only for three two-phase sides of a water-oil-gas triangular diagram T. Three-phase relative permeabilities are generally derived from these sets of two-phase data by interpolation formulas (Baker, 1988; Stone, 1970, 1973). The compressible flow equations reformulated in terms of the global pressure and saturations of the gas and water phases are derived under the assumption that the three-phase data satisfy the total differential (TD) condition (Amaziane and Jurak, 2008; Chavent and Jaffré, 1986). This condition allows the two gradients of two-phase capillary pressure to be rewritten as a single gradient of a mathematical function called the global capillary function.
Comparisons with classical approaches such as pressure-pressure or pressure-saturation formulations have shown the computational efficiency of the global pressure-saturation formulation, as it reduces the coupling between pressure and saturation equations (Chen, 2005; Chen and Ewing, 1997). The aim of this article is to compare two different parameterisation approaches used to assess the reformulated three-phase data needed to simulate three-phase flow.
To assess the three-phase data, a step-by-step algorithm based on the approximated global pressure formulation was originally proposed by Chavent and Salzano (1985), and Chavent and Jaffré (1986) for compressible three-phase flow. It is derived under the TD condition so that each of the three-phase volume factors are evaluated at the global pressure level instead of their own phase pressure. The two corresponding degrees of freedom of this original interpolation are fractional flow functions f_{j} for each fluid phase j and the global mobility d. Resulting TD three-phase data are ultimately determined with the resolution of two constrained optimization problems.
Recently, a new class of TD-interpolations algorithm was introduced by Chavent (2009), Chavent et al. (2008), and di Chiara Roupert et al. (2010), which allows the use of an exact global pressure formulation. The phase volume factors are evaluated at their own phase pressure. The two degrees of freedom of this new class of TD interpolation include the global capillary function ${P}_{c}^{g}$ and the global mobility function d, which are computed by a C^{1} and C^{0} finite element interpolations, respectively. Corresponding procedures are designed in such a way that the resulting three-phase data satisfy important mathematical and physical properties, including:
- • monotonicity, regularity and bound constraints of the resulting three-phase data;
- • the TD condition introduced in the mathematical model to simplify the pressure equation in the three-phase flow model.
In this article, we analyse the implementation of both approaches and assess the numerical results obtained in the case of two-phase data for a typical homogeneous porous medium. In the first section, the exact global pressure approach is employed using the conventions introduced by Chavent (2009) and Chavent et al. (2008). This leads to a simplified water and gas saturation equation and an expression of the pressure equation under Darcy's law. Then, we briefly discuss the implementation of both interpolation approaches to obtain three-phase TD data, which coincide with the given two-phase data set on the boundaries of the ternary diagram. Finally, in the fourth section we show numerical results, including the TD-interpolated water fractional flows, global mobilities and the global capillary function obtained by both approaches for an incompressible test case. Finally, we highlight practical aspects concerning the use of these two approaches in the context of a flow simulator.
2 Governing equations for three-phase flow in porous media
Let us recall the equations associated with the reformulated three-phase flow under the TD condition. The following notations are used for three fluid phases: j = 1 indicates water, j = 2 indicates oil, j = 3 and indicates gas.
2.1 Global pressure equation
We introduce the volume factor B_{j} evaluated at a pressure level p:
$${B}_{j}(p)={\rho}_{j}(p)/{\rho}_{j}({p}_{0})\text{,}$$ | (1) |
$${j}_{j}=-{d}_{j}({P}_{j})k{r}_{j}(S)K.(\nabla {P}_{j}-{\rho}_{j}({P}_{j})g\nabla Z)\text{,}$$ | (2) |
$$\varphi \frac{\partial}{\partial t}(\sum _{j=1}^{3}{B}_{j}{S}_{j})+\nabla .(q)=\sum _{j=1}^{3}{Q}_{j}\text{,}$$ | (3) |
$$P={P}_{2}+{P}_{c}^{g}(S,P)\text{,}$$ | (4) |
$$\left\{\begin{array}{l}for\text{\hspace{0.28em}}all\text{\hspace{0.28em}}{S}_{1}(x,t),{S}_{3}(x,t)\text{\hspace{0.28em}}and\text{\hspace{0.28em}}P(x,t),\hfill \\ \nabla {P}_{c}^{g}(S,P)={f}_{1}(S,P)\nabla {P}_{c}^{12}+{f}_{3}(S,P)\nabla {P}_{c}^{32}({S}_{3})+\frac{\partial {P}_{c}^{g}}{\partial P}(S,P)\nabla P,\hfill \end{array}\right.$$ | (5) |
$$q=-dK.\left[\omega \nabla P-\overline{\rho}g\nabla Z\right]\text{,}$$ | (6) |
$$\left\{\begin{array}{c}\hfill {f}_{j}(s,p)=k{r}_{j}(s){d}_{j}(p)/d(s,p),j=\mathrm{1,2,3}\hfill \\ \hfill d(s,p)=k{r}_{1}(s){d}_{1}(p)+k{r}_{2}(s){d}_{2}(p)+k{r}_{3}(s){d}_{3}(p)\hfill \\ \hfill \overline{\rho}(s,p)=\sum _{j=1}^{3}{\rho}_{j}(p){f}_{j}(s,p)\hfill \\ \hfill \sum _{j=1}^{3}{f}_{j}(s,p)=1.\hfill \end{array}\right.$$ | (7) |
The original global pressure formulation employs the global pressure level p to evaluate the volume factors B_{j}. In the exact global pressure formulation, phase pressure levels p_{j} are used. Based on the capillary pressure conventions, relative phase mobilities are written as , ${d}_{1}={d}_{1}({p}_{2}+{P}_{c}^{12}({s}_{1}))$, d_{2} = d_{2}(p_{2}), and ${d}_{3}={d}_{3}({p}_{2}+{P}_{c}^{32}({s}_{3}))$; a similar framework is used for ρ_{1}, ρ_{2}, and ρ_{3}.
2.2 Water and gas saturation equations
The corresponding water and gas saturation equation satisfying the TD condition can be written for j = 1,3:
$$\begin{array}{l}\varphi \frac{\partial}{\partial t}({B}_{j}{S}_{j})+\nabla .({q}_{j})={Q}_{j}\hfill \\ \text{where}\text{\hspace{0.28em}}{q}_{j}={\omega}^{-1}{f}_{j}q-d{f}_{j}K.\left[\nabla {P}_{c}^{jg}+g\left[\Delta {\rho}_{j}+\overline{\rho}\frac{1-\omega}{\omega}\right]\nabla Z\right],\hfill \end{array}$$ | (8) |
One can see that some secondary variables, such as fractional flows f_{j}, global mobility d, and the global capillary function ${P}_{c}^{g}$, must be quantified to satisfy the TD condition (Chavent, 2009) and thus to solve the three-phase reformulated flow. In the next section, we give a short description of the two approaches used to interpolate TD three-phase data using the water-oil-gas diagram T.
3 Description of two TD-interpolation approaches and the data required for both methods
3.1 TD interpolation by optimization (approach A)
We first describe the data required for the original global pressure formulation introduced by Chavent and Jaffré (1986). For each pair of fluids belonging to the water-oil side (12) and the gas-oil side (32) (see Fig. 1), we require relative permeabilities $k{r}_{i}^{ij}$, $k{r}_{j}^{ij}$ capillary pressure curves ${P}_{c}^{12}$, ${P}_{c}^{32}$, and phase mobilities d_{j} for each fluid phase j = 1,2,3. A convenient continuation of the two-phase oil fractional flow curves is chosen inside T. The three-phase oil fractional flow f_{2} is chosen to be equal to the one obtained by Stone's model (Stone, 1970). Then, factional flows and are obtained using the TD equivalent condition (Chavent and Jaffré, 1986).
This TD equivalent condition yields two simple linear equations, which can be solved on the ternary diagram using (7) and the three-phase oil fractional flow f_{2}. Therefore, we are able to build the global capillary function ${P}_{c}^{g}$ (Chavent, 2009; di Chiara Roupert et al., 2010). The global pressure varies between water pressure and gas pressure, then ${P}_{c}^{g}(\mathrm{1,0},P)=0$ at the water vertex. The TD condition (see (5)) can be re-written for all s in T in an equivalent form:
$$\left\{\begin{array}{l}\frac{\partial {P}_{c}^{g}}{\partial {s}_{1}}(s,p)={f}_{1}(s,p-{P}_{c}^{g})\frac{d{P}_{c}^{12}}{d{s}_{1}}({s}_{1}),\hfill \\ \frac{\partial {P}_{c}^{g}}{\partial {s}_{3}}(s,p)={f}_{3}(s,p-{P}_{c}^{g})\frac{d{P}_{c}^{32}}{d{s}_{3}}({s}_{3}).\hfill \end{array}\right.$$ | (9) |
Therefore, the global capillary function ${P}_{c}^{g}$ must satisfy a differential equation along any given smooth curve C in T (see Fig. 1). From the previous step, the TD fractional flow functions allow the construction of the global capillary function ${P}_{c}^{g}$ using the two two-phase capillary pressure gradients and using a convenient path (see Fig. 1) along the water-oil side (curve C_{1}) and gas-oil side (curve C_{3}):
$${P}_{c}^{g}(s,p)={\int}_{1}^{{s}_{1}}{f}_{1}(\sigma ,0,p)\frac{d{P}_{c}^{12}}{d{s}_{1}}(\sigma )d\sigma +{\int}_{0}^{{s}_{3}}{f}_{3}({s}_{1},\sigma ,p)\frac{d{P}_{c}^{32}}{d{s}_{3}}(\sigma )d\sigma \text{.}$$ | (10) |
Following Chavent and Salzano (1985) and Chavent and Jaffré (1986), the next step determines d on the third side, that is, the water-gas side of T. The global mobility of the water-gas side d^{13}must be determined carefully. Otherwise, corresponding relative permeabilities $k{r}_{1}^{13}$ for water and $k{r}_{3}^{13}$ for gas computed from (7) may range out of the interval; they also may decrease in directions in which one would expect them to increase. Therefore, a constrained optimization problem is built up on the relative permeabilities using previous TD fractional flows for water and for gas to obtain the corresponding TD relative permeabilities. Once the ternary diagram T is closed using the initial values on the water-oil boundary T^{12}, the gas-oil boundary T^{32} and optimised TD data on T^{13}, the global mobility is assessed for T. Similarly, two additional constrained optimization problems are solved to compute TD three-phase permeabilities and, therefore, global mobility; see equation (7). At the end of optimization procedure, the corresponding TD three-phase permeabilities shall have a Laplacian value that is not too large but ranges between 0 to 1 and strictly increases towards the phase-related vertex of T (Jégou, 1997; Chavent and Salzano, 1985; Chavent and Jaffré, 1986). Bound constraints and regularity behaviour are set for the objective functions of the TD data using a Quasi-Newton algorithm (Byrd et al., 1994). The resulting problem is weakly parameterized, a damped or trust region Gauss-Newton method like Levenbeg-Marquardt, or a “Dog Leg” method could also be implemented to solve sequentially the optimization issue. Details of the construction of objective functions, their gradients and linear-bound constraints are given in di Chiara Roupert (2009).
At the end of the procedure, global mobility d(s,p), fractional flows f_{j}(s,p) and global capillary function ${P}_{c}^{g}(s,p)$ are known and, thus, are used to solve the global pressure equation (4), the global volumetric flow rate (6) and the saturation equations (8).
3.2 New TD interpolation by finite elements (approach B)
In the exact global pressure formulation, d, f_{1}, f_{2}, f_{3}, and $\overline{\rho}$ are functions expressed by the oil pressure level p_{2}; see equation (4). Initial two-phase flow data are set on water-oil side T^{12}, gas-oil side T^{32} and water-gas side T^{13}, gas-oil side and water-gas side (see Fig. 2). The reader may refer to di Chiara Roupert et al. (2010) for details of the algorithm we use here.
We solve first the initial value problems on the three two-phase sides of the ternary diagram T. The three ODEs are built up from tangential derivatives of ${P}_{c}^{g}$ over the three sides of the ternary diagram. Their expressions are deduced from relations (9) and solved with an adaptative predictor-corrector scheme (di Chiara Roupert, 2009; Quarteroni et al., 2007). Therefore, the global capillary function is calculated on the three two-phase boundaries.
To satisfy the TD condition on the boundaries, the global capillary function must be slightly modified on the gas-oil side to provide an exact matching value at the gas vertex (Chavent, 2009; di Chiara Roupert, 2009; di Chiara Roupert et al., 2010).
The three two-phase global mobilities d on boundaries are then computed (see Eq. (7)) with P replaced by ${P}_{2}=P-{P}_{c}^{g}(s,p)$. The global mobility is then obtained from previous Dirichlet boundary condition by solving a harmonic equation on T by means of C^{0} piecewise linear finite elements.
In the last step, the global capillary pressure is extended to the entire ternary diagram T. It is necessary to complete the Dirichlet conditions known from ODE calculations with Neumann conditions along the same boundaries (see Fig. 2). Thus the problem amounts to enquire on each node the function ${P}_{c}^{g}$ and its normal derivative on each mid-side node belonging to the border (Chavent, 2009; di Chiara Roupert et al., 2010). Normal derivatives are computed as linear combinations of derivatives $\partial {P}_{c}^{g}/\partial {s}_{1}$ and $\partial {P}_{c}^{g}/\partial {s}_{3}$ from (9). A bi-harmonic equation over is needed to provide a continuous global capillary pressure field. It is solved using a C^{1} Hsieh-Clough-Tocher finite element method (Bernadou and Hassan, 1980; Clough and Tocher, 1965), which ensures that the fractional flows derived from ${P}_{c}^{g}$ are continuous over the entire ternary diagram. At the end of the C^{1} finite element interpolation, ${P}_{c}^{g}$, $\partial {P}_{c}^{g}/\partial {s}_{1}$ and $\partial {P}_{c}^{g}/\partial {s}_{3}$ are provided. For each fluid phase, the corresponding fractional flow is then derived for j = 1,3 using:
$${f}_{j}(s,p)=(\partial {P}_{c}^{g}/\partial {s}_{j}(s,p))/(d{P}_{c}^{j2}/d{s}_{j}({s}_{j}))\text{,}$$ | (11) |
4 A comparison of numerically derived three-phase TD data
In this section, we compare the water fractional flow f_{1}, the global mobility d and the global capillary pressure ${P}_{c}^{g}$ obtained from approach A and approach B for an incompressible three-phase flow (ω=1). Input data were chosen for a typical fine sand; we considered the van Genuchten model function (Van Genuchten, 1980) coupled with Parker's normalisation method (Parker et al., 1987) to express the capillary pressure-saturation relationship in the two-phase systems for water-oil ${P}_{c}^{12}$ and gas-oil ${P}_{c}^{32}$. Specifically, we used both the Van Genuchten parameters (α and n) in the two-phase gas-water system and the corresponding surface tensions of water-oil and gas-oil. The perfect gas law is used at a temperature of T = 293 K. All parameters used in the case study are given in Table 1. Here, the oil phase is made of trichloroethylene (TCE). In terms of wettability, water is the wetting phase in the water-oil system and oil is the wetting phase in the gas-oil system. Consequently, water is the wetting phase in the water-gas system. The resulting TD water fractional flows are shown in Figs. 3 and 4. Differences occur on the water-gas side T^{13} because optimization does not take into consideration any relative permeabilities or initial capillary pressures. In Fig. 4, numerical results obtained from the bi-harmonic problem are given; ${P}_{c}^{g}$, $\partial {P}_{c}^{g}/\partial {s}_{1}$ and $\partial {P}_{c}^{g}/\partial {s}_{3}$ and are computed using the three two-phase data given on the sides as Dirichlet and Neuman boundary conditions.
Données initiales utilisées pour les trois systèmes diphasiques.
Symbol | Value | Description | Unit |
α | 0.0001 | α Van Genuchten parameter, fine sand | [Pa^{−1}] |
m | 0.93 | m Van Genuchten parameter, fine sand | [–] |
μ _{1} | 1.2 | Water dynamic viscosity | [×10^{−3}Pa.s] |
μ _{2} | 0.55 | Oil dynamic viscosity | [×10^{−3}Pa.s] |
μ _{3} | 0.015 | Gas dynamic viscosity | [×10^{−3}Pa.s] |
σ_{13}/σ_{32} | 2.5 | Gas-oil dimensionless surface tension | [–] |
σ_{13}/σ_{12} | 2.1 | Water-oil dimensionless surface tension | [–] |
B _{ j } | 1 | Volume factor for j = 1, 2, 3 | [–] |
p | 1 | Global pressure level | [bar] |
TD fractional flows computed from (11) exactly match the given initial models of two-phase capillary pressures and relative permeabilities (see Fig. 2). To compare both approaches, we calculated relative differences at the point of the ternary diagram at which s_{1} = 45% and s_{3} = 40%. At this point, the TD fractional flow obtained with approach A is about 25% higher than the one evaluated with approach B. In this case, the water volumetric flow vector (8) is higher in approach A than B. Figs. 5–7 show the global capillary function obtained from the optimization approach and the C^{1} finite element approach, respectively. Resulting values vary in the same range. One can see that the initial value ${P}_{c}^{g}(\mathrm{1,0},p)=0$ is satisfied at the maximum water saturation. Both functions grow very rapidely near the gas vertex due to two-phase capillary gradients. Each computed function ${P}_{c}^{g}=P-{P}_{2}$ shows a regular behaviour near the maximum oil saturation. It shows that for both methods, global pressure is a more appropriate unknown than oil pressure, as p_{2} − p_{j} is singular near the jth vertex. Figs. 8 and 9 show global mobility as obtained from the optimization approach and the C^{0} finite element approach. In Fig. 8, the resulting TD global mobility is expected to be regular and smooth over T. Due to the water-gas optimised value found on T^{13}, the TD global mobility field is quite different in the water-gas region compared to the gas-oil region. To discuss differences between approach A and B, relative differences are calculated at s_{1} = 25%, s_{1} = 65%. At this point, the TD global mobility as predicted by approach B is about 50% higher than the value calculated by approach A. Hence, in that particular case, the TD global mobility as obtained by the new TD interpolation scheme (approach B) allows the global volumetric flow vector (6) to achieve higher values than that provided by approach A. In the vicinity of the oil vertex, TD global mobility is the same as that obtained using Stone's model, which is not really surprising as it has been chosen as the oil fractional flow model f_{2} in approach A (see Section 3.1). Here, mobilities which are built in approach A are based on an approximated formulation of the global pressure while results of approach B are obtained with the exact formulation: the volume factors B_{j} are measured at pressure level p_{j} and not at the global pressure level p. If the two approaches are both TD compatible, the global mobility obtained in approach A (see Figs. 8 and 9) will tend to underestimate the total flux for s_{2} < 40% compared to approach B. This situation represents the majority of cases encountered in real pollution situations with a plume at residual oil saturation (s_{2} <<< 40%). Furthermore, we present specific results for ${P}_{c}^{g}$ data on T^{13} (see Fig. 7, left). In the case of approach B, the resulting water-gas capillary pressure curve exactly fits the Van Genuchten model, whereas at low water saturations, the capillary pressure values of approach A are several times higher than the physical data. This difference is a consequence of the sharp contrast of the global mobility towards the gas vertex. This difference is due to the TD relative permeabilities optimized on the water-gas side (approach A). In approach A, the total flow will thus be underestimated for three-phase flows where the gas saturation is high.
5 Conclusion
In this article, we analysed two three-phase parameterisation approaches to assess the TD secondary variables required for efficient multiphase flow simulators using a global pressure formulation. Usually, experimental three-phase data are rarely investigated inside the three-phase water-oil-gas diagram (Bano et al., 2009; Oak et al., 1990). From a practical point of view, these two TD interpolations have given similar CPU times within the order of a second. They represent robust and innovative approaches for efficiently modelling three-phase flow equations. Both methods satisfy the TD condition and allow the three-phase compressible flow equation to be rewritten in a more suitable fractional form; see equations (6) and (8). However, by construction, the resulting two-phase water-gas data, such as the capillary pressure curve, differ between the original and the new TD interpolation. When flow takes place only in a two-phase water-gas system, TD three-phase data obtained by the optimization (approach A) might be significantly different from that derived using the finite element interpolation (approach B). It also may differ from the data from methods using classical capillary curves, such as the Van Genuchten model, as is often employed in the field of environmental sciences. As shown in Fig. 7, the finite element approach seems to be more suitable for explicitly taking into account two-phase capillary pressure and relative permeability data in a water-gas system. Furthermore, optimization parameters may be difficult to select. They depend on the governing parameters of the simulated cases (eg. viscosities, pressure levels, etc.) and more specifically initial values for the two relative permeability models used on the water-gas side optimization. Moreover, numerical tests (see approach A) showed that the global mobility field heavily depends on dynamic viscosity of the gas phase. Therefore, a finite element interpolation approach combined with a flow simulator seems to be better suited to generate three-phase flow variables than the optimization approach. In other words, the finite element approach appears to be a robust and efficient solution for three-phase flow simulation. At the reservoir scale, we can expect these differences to have little impact on the flow behaviour, as it is knwon that the shock saturation mainly dominates two phase flows. Note that in the TD approach, Chavent and Jaffré (1986) showed that the hyperbolicity condition is not fully verified particularly for large oil saturation values in the first Stone model. This condition should be satisfied by any three-phase data set if one wants small capillary pressure problems to be well posed for any values of volume factors. Besides, Stone's models are not TD compatible and therefore do not allow reformulation of the flow in a simplified way (see (6). Further research will be focussed on more systematic numerical flow comparisons with other three-phase correlations.
Acknowledgements
The authors would like to thank Professor Frédéric Delay from the University of Poitiers and the four other anonymous referees for their valuable comments and suggestions which helped improve the manuscipt significantly. Financial support for this research from the Programme National de Recherche en Hydrologie (INSU-CNRS) and the Agence De l’Environnement et de la Maîtrise de l’Energie (ADEME) is gratefully acknowledged. We would also like to thank ADEME and Burgéap SA for financing the PhD of R. di Chiara Roupert.