Two phase partially miscible flow and transport modeling in porous media ; application to gas migration in a nuclear waste repository

# Two phase partially miscible flow and transport modeling in porous media ; application to gas migration in a nuclear waste repository

Alain Bourgeat Alain Bourgeat Université de Lyon, Université Lyon1, CNRS UMR 5208 Institut Camille Jordan, F - 69200 Villeurbanne Cedex, France Mladen Jurak Department of mathematics, University of Zagreb, Croatia, Farid Smaï Université de Lyon, Université Lyon1, CNRS UMR 5208 Institut Camille Jordan, F - 69200 Villeurbanne Cedex, France    Mladen Jurak Alain Bourgeat Université de Lyon, Université Lyon1, CNRS UMR 5208 Institut Camille Jordan, F - 69200 Villeurbanne Cedex, France Mladen Jurak Department of mathematics, University of Zagreb, Croatia, Farid Smaï Université de Lyon, Université Lyon1, CNRS UMR 5208 Institut Camille Jordan, F - 69200 Villeurbanne Cedex, France    Farid Smaï Alain Bourgeat Université de Lyon, Université Lyon1, CNRS UMR 5208 Institut Camille Jordan, F - 69200 Villeurbanne Cedex, France Mladen Jurak Department of mathematics, University of Zagreb, Croatia, Farid Smaï Université de Lyon, Université Lyon1, CNRS UMR 5208 Institut Camille Jordan, F - 69200 Villeurbanne Cedex, France
July 20, 2019
###### Abstract

We derive a compositional compressible two-phase, liquid and gas, flow model for numerical simulations of hydrogen migration in deep geological repository for radioactive waste. This model includes capillary effects and the gas high diffusivity. Moreover, it is written in variables (total hydrogen mass density and liquid pressure) chosen in order to be consistent with gas appearance or disappearance. We discuss the well possedness of this model and give some computational evidences of its adequacy to simulate gas generation in a water saturated repository.

###### Keywords:
Two-phase flow, porous medium, modeling, underground nuclear waste management

## 1 Introduction

The simultaneous flow of immiscible fluids in porous media occurs in a wide variety of applications. The most concentrated research in the field of multiphase flows over the past four decades has focused on unsaturated groundwater flows, and flows in underground petroleum reservoirs. Most recently, multiphase flows have generated serious interest among engineers concerned with deep geological repository for radioactive waste. There is growing awareness that the effect of hydrogen gas generation, due to anaerobic corrosion of the steel engineered barriers (carbon steel overpack and stainless steel envelope) of radioactive waste packages, can affect all the functions allocated to the canisters, waste forms, buffers, backfill, host rock. Host rock safety function may be threaten by overpressurisation leading to opening fractures of the host rock, inducing groundwater flow and transport of radionuclides.

The equations governing these flows are inherently nonlinear, and the geometries and material properties characterizing many problems in many applications can be quite irregular and contrasted. As a result, numerical simulation often offers the only viable approach to the mathematical modeling of multiphase flows. In nuclear waste management, the migration of gas through the near field environment and the host rock, involves two components, water and pure hydrogen ; and two phases ”liquid” and ”gas”. It is then not clear if conventional models, like for instance the Black-Oil model, used in petroleum or groundwater engineering are still valid for such a situation. Our ability to understand and predict underground gas migration is crucial to the design of reliable waste storages. This is a fairly new frontier in multiphase porous-media flows, and again the inherent complexity of the physics leads to governing equations for which the only practical way to produce solutions may be numerical simulation. This exposition provides an overview of the types of standard models that are used in the field of compressible multiphase multispecies flows in porous media and includes discussions of the problems coming from  gas being one of the components. Finally the paper also addresses one of the outstanding physical and mathematical problems in multiphase flow simulation: the appearance disappearance of one of the phases, leading to the degeneracy of the equations satisfied by the saturation. It has been seen recently in a benchmark organized by the French agency in charge of the Nuclear Waste management in France Bench () that none of the usual codes used in that field where able to simulate adequately the appearance or/and disappearance of one of the phases. In order to overcome this difficulty, we will discuss a formulation based on new variables which doesn’t degenerate. We will demonstrate through two test cases, the ability of this new formulation to actually cope with the appearance or/and disappearance of one phase. The scope of the paper is limited to isothermal flows in rigid porous media and will not include the possible process of ”pathway dilation” as described in Hors ().

## 2 Conceptual and mathematical model

Our goal is first to present a survey of the conventional models used for describing two-phase two-components flow in porous media. Most of these models have been designed and widely used in petroleum engineering (see for instance: Allen (), Bear:1 (), ChJa (), DPeace (), ECLIPSE ()). We consider herein a porous medium saturated with a fluid composed of 2 phases : liquid and gas. According to the application we have in mind, the fluid is a mixture of two components: water (mostly liquid) and hydrogen (, mostly gas). Water is present in the gas phase through vaporization, and hydrogen is present in the liquid phase through dissolution. The fluids are compressible and the model is assumed to be isothermal. For simplicity, we assume first the porous medium to be rigid, meaning the porosity is only a function of the space variable ; and second, we neglect the pressure induced dilation of gas pathways. Hydrogen being highly diffusive we have described in detail how the diffusion could be taken in account in the models.

### 2.1 Petrographic and fluid properties

#### 2.1.1 Fluid phases

The two phases will be denoted by indices , for liquid and for gas. Associated to each phase in the porous media are the following quantities:

• liquid and gas phase pressures: , ;

• liquid and gas phase saturations: , ;

• liquid and gas phase mass densities: , ;

• liquid and gas phase viscosities: , ;

• phase volumetric flow rates,  and ;

Then, the Darcy-Muskat law says:

 ql=−K(\bf x)krl(Sl)μl(∇pl−ρlg),qg=−K(\bf x)krg(Sg)μg(∇pg−ρgg), (1)

where is the absolute permeability tensor, and are the relative permeability functions, and is the gravity acceleration.

The above phase mass densities and viscosities are all functions of phase pressures and of phase composition. Moreover, phase saturations satisfy

 Sl+Sg=1 (2)

and the pressures are connected through a given experimental capillary pressure law:

 pc(Sg)=pg−pl. (3)

From definition (3) we should notice that is strictly increasing function of gas saturation, .

#### 2.1.2 Fluid components

Water and pure hydrogen components will be denoted by indices and . Since the liquid phase could be composed of water and dissolved hydrogen we need to introduce the water mass density in the liquid phase , and the hydrogen mass density in the liquid phase . Similarly, we introduce the water and hydrogen mass densities in the gas phase: , . Note that the upper index is always the component index, and the lower one denotes the phase. We have, then

 ρl=ρwl+ρhl,ρg=ρwg+ρhg. (4)

Since the composition of each phase is generally unknown we introduce the mass fraction of the component in the phase

 ωiα=ρiαρα;ωwα+ωhα=1,α∈{g,l}. (5)

### 2.2 Mass conservation equation

The conservation of the mass applies to each component and reads, for any arbitrary control volume :

 dmidt+Fi=Fi,i∈{w,h}. (6)

where

• is the mass of the component , in the control volume , at given instant ;

• is the rate at which the component is leaving (migrating from) the volume , at given instant ;

• is the source term (the rate at which the component is added to by the source).

The mass of each component is a sum over all phases and :

 mi=∫RΦ(Slρlωil+Sgρgωig)d\bf x,i∈{w,h}.

In each phase, the migration of a component is due to the transport by the phase velocity and to the molecular diffusion:

 Fi =∫∂R(ρlωilql+ρgωigqg+jil+jig)⋅nd\bf x;i∈{w,h};

where is the unit outer normal to . The phase flow velocities, and are given by the Darcy-Muskat law (1), and the component diffusive flux in phase is denoted , , , and will be defined in the next paragraph, by equations (12). From (6) we get the differential equations:

 Φ∂∂t(Slρlωwl+Sgρgωwg)+div(ρlωwlql+ρgωwgqg+jwl+jwg)=Fw, (7) Φ∂∂t(Slρlωhl+Sgρgωhg)+div(ρlωhlql+ρgωhgqg+jhl+jhg)=Fh. (8)

#### 2.2.1 Diffusion fluxes

From and , the water and hydrogen molar masses, using definitions (5) we define the following water and hydrogen molar concentrations in each phase :

 chα=SαρhαMh=SαραωhαMh,cwα=SαρwαMw=SαραωwαMw. (9)

Then the phase molar concentrations, , is:

 cα=chα+cwα=Sαρα(ωhαMh+ωwαMw). (10)

Usually, component diffusive flux in phase is assumed to be depending on , the component molar fraction in phase , defined from (9) and (10):

 Xhα=chαcα=ωhαωhα+(Mh/Mw)ωwα,Xwα=cwαcα=ωwαωwα+(Mw/Mh)ωhα;∑i=h,wXiα=1, (11)

for . Molar diffusive flux of component in phase is given by

 Jiα=−cαDiα∇Xiα;α∈{g,l},i∈{w,h}

and give molar diffusive flow rate of component through the unit area. Coefficients and (unit ) are Darcy scale molecular diffusion coefficients of components in phase . Mass flux of component in phase , in equations (7), (8) is then obtained by multiplying the above component molar diffusive fluxes, , by the molar mass of the component , and by the rock porosity:

 jhα=−ΦMhcαDhα∇Xhα,jwα=−ΦMwcαDwα∇Xwα. (12)
###### Remark 1

: Number of unknowns in the system is eight: ; but, up to now, we have only four equations (2), (3), (7) and (8), and therefore, four additional equations are needed to close the system.

###### Remark 2

: Note that and are not exactly molecular diffusion coefficients in phase , corresponding to molecule-molecule interactions in free space but , , , are effective diffusion coefficients, obtained from the strict molecular diffusion coefficients by a kind of averaging through the whole porous medium (see section 2.6 in Bear:1 (), or Milli () or Jin ()). Moreover, for simplicity, we have not included in diffusion of component in phase any dependancy on phase saturations; and then did not consider possible non linear effects coming from coupling between advective-diffusive transport (”dusty gas” model) and molecular streaming effects (Knutsen diffusion) (see for instance Tough ()).

###### Remark 3

In a binary system diffusive fluxes satisfy , for , and therefore we have

 MhDhα=MwDwα,α∈{g,l}. (13)

### 2.3 Phase equilibrium Black-oil model

The additional equations needed to close the system of equations (2), (3), (7) and (8) will come from the assumption that the two phases are in equilibrium; equilibrium meaning that at any time the quantity of hydrogen dissolved in the water is maximal for the given pressure, and similarly, the quantity of evaporated water is maximal for the given pressure. Composition of each phase is then uniquely determined by its phase pressure and saturation. In this flow situation we say that water and gas phases are saturated. Nevertheless it could happen that one of the phases disappears: either water can completely evaporate or hydrogen can be completely dissolved in the water. Then in these situations the composition of the remaining phase is not uniquely determined by its phase pressure, and the phase composition becomes an independent variable (instead of the saturation which is now constant, 0 or 1). This flow situation correspond to the so-called unsaturated flow. For unsaturated flow, standard practice in petroleum reservoir engineering is to introduce the following quantities:

• liquid and gas formation volume factors, , ;

• solution gas/liquid phase Ratio ;

• vapor water/gas phase Ratio .

The explanation of formation volume factors and solution component/phase ratios, as used in the oil reservoir modelling, is as follows. Considering the volume of liquid at reservoir conditions (reservoir temperature and pressure); when this volume of liquid is transported through the tubing to the surface it separates at standard conditions to a volume of liquid , and a volume of gas (coming out of the liquid, due to the pressure drop). At standard (i.e. stock tank) conditions, the liquid phase contains only oil component (here only water component) and the gas phase contains only gas component (here only hydrogen). Then applying conservation of mass,

 ΔVreslρl=ΔVstdlρstdl+ΔVstdgρstdg;

and denoting the solution gas/liquid Ratio , we may write:

 ΔVreslρl=ΔVstdl(ρstdl+Rsρstdg),

and the liquid phase mass density decomposition

 ρl=ρstdl+RsρstdgBl, (14)

where is the liquid formation volume factor, . Similarly, from and the gas formation volume factor , we get the gas phase mass density decomposition

 ρg=ρstdg+RvρstdlBg; (15)

where the vapor water/gas phase ratio . Now, in order to express the diffusion fluxes in terms of the new variables and we have to rewrite the molar concentrations defined in (9):

 chl=SlRsρstdgMhBl,cwl=SlρstdlMwBl,chg=SgρstdgMhBg,cwg=SgRvρstdlMwBg. (16)

These concentrations, defined in (16) may be slightly different from those previously defined in (9) if the gas at standard conditions is not composed only of hydrogen, and the liquid, also at standard conditions, is note made only of water. The phase molar concentrations, from (16), are then given by:

 cl=chl+cwl=SlBl(RsρstdgMh+ρstdlMw)=SlBlρstdgMh(Rs+F),
 cg=chg+cwg=SgBg(ρstdgMh+RvρstdlMw)=SgBgρstdgMh(1+FRv),

where is given by:

 F=MhρstdlMwρstdg. (17)

The component molar fractions in phase, are now defined as:

 Xhl =chlcl=RsρstdgRsρstdg+(Mh/Mw)ρstdl=RsRs+F, (18) Xwl =cwlcl=ρstdlρstdl+(Mw/Mh)Rsρstdg=FRs+F, Xhg =chgcg=ρstdgρstdg+(Mh/Mw)Rvρstdl=11+FRv, Xwg =cwgcg=RvρstdlRvρstdl+(Mw/Mh)ρstdg=FRv1+FRv.

Mass diffusive fluxes of components, defined in (12), depend on component molar fractions in phases; they have then to be rewritten. For instance, the mass diffusion flux of hydrogen in water, takes now the form:

 jhl/ρstdg =−ΦMhρstdgSlBlρstdgMh(Rs+F)Dhl∇Xhl=−ΦSlBlFRs+FDhl∇Rs.

Other diffusion fluxes could be obtained by similar calculations, leading to the following formulas:

 ϕhl=jhl/ρstdg =−ΦSlBlFRs+FDhl∇% Rs,ϕwl=jwl/ρstdl=ΦSlBl1Rs+FDwl∇Rs, (19) ϕhg=jhg/ρstdg =ΦSgBgF1+FRvDhg∇% Rv,ϕwg=jwg/ρstdl=−ΦSgBg11+FRvDwg∇Rv. (20)
###### Remark 4

Let us recall, like in Remark 3, that the diffusion coefficients satisfy , ; and, like in (11), for any phase , in (18).

The Darcy fluxes (1) can now be rewritten in the following form:

 ql=−Kkrlμl(∇pl−ρstdl+RsρstdgBlg),qg=−Kkrgμg(∇pg−ρstdg+RvρstdlBgg),

and the component fluxes, normalized by standard densities, are,

 ϕw=1Blql+RvBgqg+ϕwl+ϕwg,ϕh=RsBl% ql+1Bgqg+ϕhl+ϕhg. (21)

Finally, we can write mass conservation (7), (8) in the following form:

 Φ∂∂t(SlBl+RvSgBg)+div (ϕw)=Fw/ρstdl, (22) Φ∂∂t(SlRsBl+SgBg)+div (ϕh)=Fh/ρstdg. (23)

These equations have to be completed by equations (2) and (3). For instance, in (22) and (23), we can take saturation and one of the pressures as independent variables; for example and , and for the other terms the following functional dependencies:

 Bl(pl),Bg(pg),Rs(pl),Rv(pg),μl(pl),μg(pg),krw(Sg),krg(Sg).

According to their definition, , , and () are constants, and and are depending only on space position.

#### 2.3.1 Unsaturated flow

Equations (22) and (23), with saturation and pressure as unknowns, are valid if there is no missing phase (saturated flow); but if one of the phases is missing (unsaturated flow), equations and unknowns have to be adapted. There are two possible unsaturated cases, according to either the gas phase or the liquid phase disappears :

1. Gas phase missing (Hydrogen totally dissolved in water): then we have , . Generalized gas phase Darcy’s velocity is equal to zero since . Independent variables are now , the liquid phase pressure, and , the solution gas/liquid phase Ratio; then we must write and , for , where is equilibrium solution gas/liquid phase ratio.

2. Liquid phase missing: then, , . Generalized liquid phase Darcy’s velocity is zero since ; independent variables are now , the gas phase pressure, and the water vapor/gas phase ratio, become the new independent variables. However pressure can be kept as independent variable since could be expressed through the capillary pressure law (3). We must also write and , for , where is the equilibrium water vapor/gas phase ratio.

These above conditions can be summarized as

 Sg≥0,^Rs(pl)−Rs≥0,(^Rs(pl)−Rs)Sg=0, (24) Sl≥0,^Rv(pg)−Rv≥0,(^Rv(pg)−Rv)Sl=0. (25)
###### Remark 5

In this last section, Phase Equilibrium Black-oil model, we did not take in account a possible interplay between dissolution and capillary pressure.

### 2.4 Thermodynamical equilibrium Henry-Raoult model

Another way of closing the system of equations (2), (3), (7) and (8) is to use phase thermodynamical properties for characterizing equilibrium . We use first ideal gas law and Dalton law,

 pg=pwg+phg, (26)

where and are the vaporized water and hydrogen partial pressures in the gas phase; and

 pwg=ρwgMwRT,phg=ρhgMhRT. (27)

is the temperature, is universal gas constant and , are the water and hydrogen molar masses.

Next, we apply Henry’s and Raoult’s laws which say that, at equilibrium, the vapor pressure of a substance varies linearly with its mole fraction in solution. In Henry’s law the constant of proportionality is obtained by experiment and in Raoult’s law the constant is the pressure of the component in its pure state. Here we will assume, for simplicity, that the quantity of dissolved hydrogen in the liquid is small; then these laws reduce to the linear Henry’s law, which says that the amount of gas dissolved in a given volume of the liquid phase is directly proportional to the partial pressure of that same gas in the gas phase:

 ρhl=H(T)Mhphg, (28)

where is the Henry’s law constant, depending only on the temperature. For the liquid phase, we apply Raoult’s law which says that the water vapor pressure is equal to the vapor pressure of the pure solvent, at given temperature, multiplied by the mole fraction of the solvent. The water vapor partial pressure of the pure solvent depends only on the temperature and therefore is a constant, denoted here by , so we have from definition (9)

 pwg=^pwg(T)Xwl=^pwg(T)ρwlρwl+(Mw/Mh)ρhl. (29)

Further on, we can include in formula (29) the presence of capillary pressure, by using Kelvin’s equation (see Abbas ()) which gives:

 pwg=^pwg(T)ρwlρwl+(Mw/Mh)ρhle−Mwpc/(RTρl). (30)

If now, to equations (26)-(28) and (30) we add the relations

 ρhl+ρwl=ρl,ρhg+ρwg=ρg, (31)

and the water compressibility, defined by

 ρwl=ρstdlBl(pl); (32)

then, we have 8 equations: (26), (27), (28), (30), (31) and (32) and 10 unknowns:

 pl,pg,pwg,phg,ρhl,ρwl,ρl,ρhg,ρwg,ρg.

We may, for instance, parametrize all these 10 unknowns, by the two phase pressures and . For this we should combine the Henry law (28), the Raoult-Kelvin law (30) and (26) leading to the system of two equations for and :

 pwg =^pwg(T)ρwlρwl+(MwH(T))phge−Mw(pg−pl)/(RT(ρwl+MhH(T)phg) pg =pwg+phg.

It is easy to show that this above system of two equations has a unique solution for any ; and solving this system we obtain

 pwg=f(pg,pl),phg=g(pg,pl).

By (27) we can then write , and as functions of and , and by (28) and (32) we can finally express , and as functions of the phase pressures.

###### Remark 6

The gas phase will appear only if there is a sufficient quantity of dissolved hydrogen in the liquid phase, and this quantity is exactly the dissolved gas quantity, at equilibrium, given by Henry’s law. But when the dissolved gas (hydrogen) quantity is smaller than the quantity of hydrogen at equilibrium, then the Henry law does not apply; and is then equal to zero and could not be taken as unknown. In this situation, instead of saturation, we may take as independent variable. We notice that Henry-Raoult model based on thermodynamical equilibrium leads to a similar concepts as the ones developed for establishing Black Oil model for reservoir modeling.

###### Remark 7

In the Henry-Raoult model, if there is no vaporized water, , the criteria for non saturated flow is simple and reads

 ρhlMhH(T)

if there is a threshold pressure, i.e. . The same criterion can be used also if there is vaporized water, as long as the pressure in porous medium is much larger than vaporized water pressure.

###### Remark 8

In Henry-Raoult model, above, we were assuming for simplicity, no complete evaporation of the water.

#### 2.4.1 Comparison of Equilibrium models

We will now consider in more details a special case, where we may, from the Henry-Raoult model, get back the Black oil model. In both models we will then use diffusive fluxes developed in Section 2.3. For simplicity, we will neglect in (29) influence of capillary pressure and solved gas on water vapor partial pressure, leading to being a constant. Then we have from (28), (32) and (31)

 ρl=ρhl+ρwl=H(T)Mhphg+ρstdlBl(pl)=ρstdl+Bl(pl)H(T)Mh(pg−^pwg(T))Bl(pl), (34)

and also from (27),

 ρg =ρhg+ρwg=MhRT(pg−^pwg(T))+MwRT^pwg(T). (35)

Therefore, equations (34) and (35) can be written in a form similar to (14)-(15) by defining the gas formation volume factor , solution gas/liquid phase ratio and vapor water/gas phase ratio , from above thermodynamical relations:

 Rs=Bl(pl)H(T)Mhρstdg(pg−^pwg(T)) (36) Bg=RTρstdgMh(pg−^pwg(T)),Rv=1F^pwg(T)pg−^pwg(T), (37)

where is given by (17).

Compared to the Black-Oil model, Section 2.3, here, it is clear from (36) that solution gas/liquid phase ratio depends on and and not only on .

###### Remark 9

For the case without any water vapor, the corresponding thermodynamical model is obtained simply by taking in (36), (37) and in equations (22) and (23). In the same way, if we neglect dissolved hydrogen, the corresponding thermodynamical model is obtained simply by taking in (36), (37), and in (22) and (23).

#### 2.4.2 Model assuming water incompressibility and no water vaporization

In the Henry-Raoult model we assume now that the water is incompressible and the water vapor quantity is neglectful; i.e. the gas phase contains only hydrogen, then (see Remark 7) and in (37) leading to . Formulas (36), (37), could be rewritten as

 Bl≡1, Rs=Chpg,1Bg=Cvpg; (38)

where we have denoted

 Ch=H(T)Mhρstdg,Cv=MhRTρstdg; (39)

and equations (14)-(15) become:

 ρl(Rs)=ρstdl+Rsρstdg,ρg(pg)=Cvρstdgpg. (40)

Similarly to constant defined by (17), we use the density ratio

 G=ρstdlρstdg, (41)

and equations (22), (23) reduce to

 Φ∂Sl∂t+div(ql−1GJ)=Fw/ρstdl, (42) Φ∂∂t(SlRs+CvpgSg)+%$div$(Rsql+Cvpgqg+J)=Fh/ρstdg, (43) ql=−Kkrlμl(∇pl−(ρstdl+Rsρstdg)g),qg=−Kkrgμg(∇pg−Cvρstdgpgg), (44) J=−ΦSlFRs+FDhl∇Rs, (45)

where we have denoted and used Remark 3 to get , from which it follows

Here there is only hydrogen in phase gas, and Henry’s law assumes thermodynamical equilibrium in which quantity of dissolved hydrogen is proportional to gas phase pressure. In saturated case (where two phases are present) Henry’s law reads , and we can then work with variables and in equations (42)–(45).

When the gas phase disappears, the gas pressure drops to the liquid pressure augmented by entry pressure, , the liquid can contain any quantity of dissolved hydrogen between zero and , from Henry’s law (see (28) and definition (39)).

And then, when the gas phase is absent (one of the unsaturated cases) , we will replace saturation, as we did in the Black-Oil model for the same unsaturated case in Section 2.3.1, by a new variable (see definitions (28), (36), (39)), such that

 Rsρstdg=ρhl

is the mass density of dissolved hydrogen in the liquid phase.

Assuming at standard conditions the gas phase contains only hydrogen (hydrogen component massgas phase mass) and the liquid phase contains only water, with water incompressibility and no vaporized water, we see that, in definition (36), which is exactly the definitiion given in the Black-Oil model in Section 2.3.

The physical meaning of then stays the same either in Henry-Raoult or in Black-Oil models even in the unsaturated case. Moreover, from mass conservation law it follows that the dissolved hydrogen mass density is continuous when the gas phase vanishes and this can be expressed as an unilateral condition:

 0≤Sg≤1,0≤Rs≤Chpg,Sg(Chpg−Rs)=0.

In unsaturated region, where (that is ) we replace variable by , and equations (42)–(45) degenerate to:

 div(ql−1G%$J$)=Fw/ρstdl; (46) Φ∂Rs∂t +div(Rsql+J)=Fh/ρstdg; (47) ql =−Kλl(1)(∇pl−(ρstdl+Rsρstdg)g); (48) J =−ΦFRs+FDhl∇Rs. (49)

### 2.5 Saturated/unsaturated state, general formulation

Finally in the saturated region we used and as variables in (42)–(45) but in unsaturated region we should use other variables, and , in (46)–(49). In order to avoid the change of variables and equations in different regions as above we prefer to introduce a new variable

 X=(1−Sg)Rs+CvpgSg; (50)

in view to make equation (43) parabolic in .

This new variable is well defined both in saturated and unsaturated regions. Moreover, from (38) and (28),