# A New Constitutive Relation for Describing the Flow through Porous Media with Unsaturated–Saturated Transition

Authors: Maria Laura Martins-Costa; Daniel Cunha da Silva; Michele Cunha da Silva; Rogério Martins Saldanha da Gama
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The present work proposes an adequate constitutive relation to treat the process of filling with a fluid porous medium. This relation assures the problem to remain hyperbolic when the porous medium is saturated by the fluid. When the fluid fraction reaches porosity, the proposed constitutive relation increases the porous matrix resistance to more fluid inlet due to an important feature: it is a continuous and differentiable function with first derivative being also an increasing function. This allows assuring that the fluid fraction may exceed the porosity only by a very small value, making the constitutive relation realistic. Some examples compare this new constitutive relation with previous ones, highlighting its advantages.


---

## Body

## 1. Introduction

Along the time, an effort to adequately describe the filling up of a porous medium by a fluid has given rise to many works, especially when the transition from unsaturated to saturated flow is concerned. For instance, in Saldanha da Gama el al. [1], a physically realistic restriction which never allows the fluid fraction to be greater than the porosity, has been proposed. This constitutive relation, actually an upper bound, assures the problem to remain hyperbolic while the flow is unsaturated. However, when saturation is reached the problem loses its hyperbolicity, so numerical techniques such as Glimm’s Scheme cannot be applied to approximate it. (Actually, the work [1] is a generalization of [2].) A subsequent work [3] proposes a constitutive relation allowing a small supersaturation. Its advantage over the previous work [1] is that the hyperbolicity of the problem is never lost, even when saturation is reached. Although these two mentioned previous works adequately describe the transition from unsaturated flow to saturated one and the latter assures that the hyperbolicity feature is always maintained, in Martins-Costa el al. [3] the first derivative of the constitutive relation is a discontinuous function. Its second part is actually a straight line, so that it is not possible to increase the porous matrix resistance to the incoming of more fluid when the fluid fraction reaches the porosity.The new constitutive relation proposed in this work to describe the filling up of a porous matrix by a fluid is a continuous and differentiable function and its first derivative is an increasing function. Besides keeping the good feature of assuring the problem hyperbolicity regardless the fluid fraction value, this relation is more realistic than the previous one [3]. Its advantage is that it allows increasing the porous medium resistance to more fluid inlet when the fluid fraction reaches the porosity, preventing the fluid fraction from being significantly greater than the porosity. In other words, this relation allows controlling the (very small) quantity of fluid fraction that may exceed the porosity.Flows through unsaturated porous media are mainly caused by a force depending on the saturation gradient. Distinct models for transport in porous media are compared by Alazmy and Vafai [4]. Most problems are treated by means of a volume averaging technique (quantities as concentration and velocity components are averaged in volume and a Continuum Mechanics approach is employed), like phase change in porous media [5]. When the porous media deformation is accounted for, mainly in soil applications, the seminal work of Biot [6] and his followers on consolidation (actually a time-dependent volume reduction of the soil skeleton caused by a load, generating distinct velocities of the soil particles and of the liquid in the pores) must be accounted for. Another approach to deal with porous media deformation is to consider a hierarchy of length scales, the so-called multiscale models, in which an upscaling of the microscopic model is used to give rise to the macroscopic model. An interesting example is the hybrid theory developed by Hassanizadeh and Gray [7] to account for interfacial effects in two-phase flow in porous media.A distinct approach is used in this work: A Mixture Theory approach. Its basic assumption is that the mixture is composed of superposed continuous constituents, each one occupying the whole volume of the mixture. Also, each constituent is endowed with a kinematics and must satisfy the balance laws and the mixture as a whole must satisfy the balance laws and the Second Law of Thermodynamics [8, 9]. In this case the partially saturated porous media are regarded as a three constituents’ mixture, namely, the liquid constituent (representing the fluid), the gas constituent (with very low density, to account for the mixture compressibility), and the solid constituent (representing the porous matrix).In Martins-Costa and Saldanha da Gama [10], the partially saturated porous matrix one-dimensional description gives rise to a nonlinear system of hyperbolic equations. A porous slab with distinct conditions is considered to account for most one-dimensional cases of interest. After the solution of the associated Riemann Problem, the system is approximated by Glimm’s Scheme and an operator splitting technique. In reality, the solution of the associated Riemann Problem plays an important role in hyperbolic systems. It is worth noting that the constitutive relation proposed in the present work allows the Riemann invariants to be represented by explicit closed form formulae.Nonisothermal fluid flows through unsaturated porous media have been considered by Saldanha da Gama and Martins-Costa [11], in which the approximation of the hyperbolic part (representing the continuity and momentum equations for the fluid constituent) is used as input for the elliptic one (representing the energy equations for both the solid and the fluid constituents). The hydrodynamic problem is solved by combining an operator splitting technique and Glimm’s Scheme, which, in turn, uses the solution of a previously defined number of associated Riemann problems. Glimm’s method is also employed to describe a pollutant transport in the atmosphere (considering also a contact shock which is associated with the second eigenvalue) by Martins-Costa and Saldanha da Gama [12].The complete solution of the Riemann Problem associated with a compressible flow of the generalized Chaplygin gas was presented by Pang et al. [13], accounting for a delta shock wave. Abbassi and Namah [14] describe the oil-water interface of the displacement of oil by water injection in a porous medium (representing an oil reservoir). An auxiliary Riemann Problem is solved in order to obtain the oil-water speed explicit expression. It is worth noting that Abbassi and Namah [14] employ a Mixture Theory approach. The Dam-Break Problem has a clear analogy with the classical Riemann Problem, since dam-break flows present an initial discontinuity on both sides of the dam wall caused by two district states. Gupta and Singh [15] solved the associate Riemann Problem in order to implement Glimm’s method to approximate one-dimensional shallow water equations of the Dam-Break problem. In the sequence, they use an operator splitting technique to extend to a two-dimensional case. Martins et al. [16] solved the Local Inertial Equations in flood hydrodynamic modeling obtained by neglecting the convective terms in the Saint-Venant equations. According to the authors neglecting these terms do not impair the problem accuracy. They compute Riemann invariants to use the Method of Characteristics, derived for quasilinear differential equations.Essentially, the present work proposes a very convenient constitutive relation to describe the filling up of a porous matrix by a fluid, considering a Mixture Theory approach. This relation not only assures that the hyperbolicity of the system is maintained even when saturation is reached, like in Martins-Costa et al. [3], but also it consists of a continuous and differentiable function with an increasing function as its first derivative. This important feature allows controlling the (very low) quantity of fluid fraction that may surpass the porosity. The complete solution of the associated Riemann Problem is presented, and some meaningful comparisons are made with the constitutive relations employed in previous works [1, 3] to show the advantage of this new relation.


## 2. Mechanical Model

In a Mixture Theory context, the balance equations should be satisfied by each constituent (of the mixture) as well as by the mixture [8, 9]. The unsaturated porous matrix is modeled as chemically nonreacting continuous mixture of a rigid solid constituent (and at rest), a liquid constituent (denoted as fluid constituent), and a third gas constituent with negligible inertia to account for the compressibility of the mixture as a whole, coexisting superposed. Since the solid constituent is rigid and the third constituent has negligible inertia, it suffices to solve mass and momentum balance equations for the fluid constituent only, which are given by(1)∂ρF∂t+∇·ρFvF=0ρF∂vF∂t+∇vFvF=∇·TF+mF+ρFbFwhere ρF is the fluid constituent mass density (local ratio between the fluid constituent mass and the corresponding volume of mixture), vF is its velocity in the mixture, TF is the partial stress tensor associated with it (analogous to Cauchy stress tensor in Continuum Mechanics), mF is the momentum supply acting on it due to its interaction with the remaining constituents of the mixture, and bF is the body force. The momentum supply to the mixture must be zero; in other words: ∑i=1nmi=0.Assuming the partial stress tensor symmetrical, the balance of angular momentum is automatically satisfied.The fluid fraction is the ratio between the fluid constituent mass densityρF and the fluid mass density (Continuum Mechanics viewpoint) ρf; in other words, ϕ=ρF/ρf and ε represents the porous matrix porosity.The dynamic interaction among the constituents is accounted for by the momentum supplymF which can be represented by(2)mF=-μfKϕ2vF-μfDK∇ϕwhere μf is the fluid viscosity (in a Continuum Mechanics approach), K represents the specific permeability of the porous matrix, and D is a diffusion coefficient.The first term ofmF accounts for the drag forces effect on the fluid constituent (the well-known Darcian term in the literature—see [17]) while the second term is an attempt to model the capillary forces arising from the nonuniform fluid distribution inside the porous matrix. In some very particular conditions, it would allow recovering Fick’s law (see [10]). It is important to note the absence of interfacial properties when a Mixture Theory viewpoint is employed, although Bowen [18] has accounted for capillary pressures using a Mixture Theory approach and interaction stresses in this approach were considered by Williams [19] and Sampaio and Williams [20], being analogous to the interfacial stresses of the multiscale theories of Hassanizadeh and Gray [7].The partial stress tensor is analogous to Cauchy tensor in Continuum Mechanics. Its constitutive equation is based on Williams [19] assumption of TF being proportional to the pressure acting on the fluid constituent and to the gradient of fluid constituent velocity and on the dominance of normal fluid stresses proposed by Allen [21], leading to(3)TF=-p-ϕIwhere p- is a pressure (assumed constant while the flow is unsaturated) and I is the identity tensor.At this point all the quantities are assumed to depend solely on timet and on positionx and the only nonvanishing component of the fluid constituent velocity, vF, is denoted byv, so (1)-(3) give rise to the following system:(4)∂ϕ∂t+∂∂xϕv=0ρfϕ∂v∂t+ϕv∂v∂x=-∂∂xϕp--ρf2μfDK∂ϕ∂x-μfKϕ2vIn order to write system (4) in a more convenient form, a reference pressure p0 is defined such that(5)p0∂ϕ∂x=∂∂xϕp-+ρf2μfDK∂ϕ∂xSince this article aims at adequately treating the transition from unsaturated flows to saturated ones, the unsaturated problem must be a hyperbolic one. This is a fundamental step to allow using methodologies such as Glimm method. It is worth noting that the Darcian terms are small when compared with the inertial and concentration gradients terms, a physically reasonable hypothesis for highly permeable porous media and/or very small viscosity fluids. In order to allow the direct application of Glimm method, the Darcian term is neglected in this work.It is important to note that the Darcian term may be considered using Glimm method associated with an operator splitting technique, as previously done in Martins-Costa and Saldanha da Gama [22], for instance.


## 3. Constitutive Relation

The one-dimensional flow through an unsaturated homogeneous porous medium may be represented by the following set of equations [1, 3]:(6)∂ϕ∂t+∂∂xϕv=0∂∂tϕv+∂∂xp+ϕv2=0p=c2ϕ,c>0in which c is a positive constant, ϕ represents the fluid fraction, v represents the velocity, p represents the pressure (for simplicity p0=p), and the fluid density ρf has been considered unitary in some unit systems.When the saturation is reached, we have that the fluid fraction becomes equal to the porosityε (ϕ=ε) and the pressure becomes independent of ϕ, since the speeds becomes unbounded.In order to surpass this shortcoming, the following constitutive relationship was proposed in Martins-Costa el al. [3]:(7)p=c2ϕ,forϕ<εc2ε+c□2ϕ-ε,forϕ≥ε,c□>c>0to encompass saturated and unsaturated flows within the same relationship, allowing a small supersaturation.The main idea consists of choosing a very large constantc□ and allowing the fluid fraction to be slightly greater than the porosity. Nevertheless, the model proposed in [3] fails when the pressure becomes very large, since there is not a bound for the fluid fraction and, during a simulation, we can have unrealistic values for ϕ.In fact, taking into account the fact that(8)ϕ=p-c2εc□2+εit is not possible to establish, a priori, a bound for the fluid fraction. In other words, we have that (9)limp→∞ϕ=+∞The present work proposes a relationship between the pressure and the fluid fraction that allows establishing, a priori, an upper bound for the fluid fraction, thus allowing a reliable and always physically consistent simulation. In addition, the constitutive relation to be proposed here gives rise to a user-friendly Riemann Problem in which the Riemann invariants can be represented by an explicit closed form formula.The proposed constitutive relation is the following:(10)p=c12ϕ+c22ϕε--ϕ-2c1c2ε-ln⁡ε--ϕε-,ε-=ε+δ,δ>0in which c1 and c2 are positive constants and δ is a positive constant that can be very small. From the above equation it is easy to see that(11)p>c22ϕε--ϕ⇒pc22+pε->ϕ⇒ϕ<ε+δalwaysand, so, we have an a priori limit for the supersaturation, even for very large pressures, since(12)limp→∞ϕ=ε+δ<+∞The flow is considered saturated in the regions in whichε≤ϕ<ε+δ. When ϕ>ε, as the fluid fraction becomes larger, the speeds of propagation increase rapidly, as well as the pressure, providing a realistic response and preserving the hyperbolic behavior of the system of partial differential equations.


## 4. The Associated Riemann Problem

Let us consider now the following problem:(13)∂ϕ∂t+∂∂xϕv=0∂∂tϕv+∂∂xp+ϕv2=0p=c12ϕ+c22ϕε--ϕ-2c1c2ε-ln⁡ε--ϕε-The Riemann Problem associated with (13) arises when, for all x∈(-∞,∞), the following initial data is assumed:(14)ϕ=ϕLfort=0,-∞<x<0ϕRfort=0,0<x<∞v=vLfort=0,-∞<x<0vRfort=0,0<x<∞whereϕL,ϕR,vL,andvR are constants.The eigenvalues for this problem are given by(15)λ1=λ1^ϕ,v=v-p\andλ2=λ2^ϕ,v=v+p\Taking into account (10) we may write (16)λ1=λ1^ϕ,v=v-c1+c2ε-ε--ϕλ2=λ2^ϕ,v=v+c1+c2ε-ε--ϕSince the second derivative ofp with respect to ϕ, given by(17)p\\=2c1+c2ε-ε--ϕc2ε-ε--ϕ2,is always positive valued, we have that the solution (in a generalized sense) of the Riemann Problem consists of connecting the left state ϕL,vL and the right state ϕR,vR to an intermediate state ϕ∗,v∗, by means of rarefactions or shocks, according to one of the four ways listed below:(18)ϕL,vL→1-rarefactionϕ∗,v∗→2-rarefactionϕR,vRϕL,vL→1-rarefactionϕ∗,v∗→2-shockϕR,vRϕL,vL→1-shockϕ∗,v∗→2-rarefactionϕR,vRϕL,vL→1-shockϕ∗,v∗→2-shockϕR,vRThe stateϕL,vL will be connected to the state ϕ∗,v∗ by a 1-rarefaction if, and only if, the eigenvalue λ1 is an increasing function of the ratio x/t between these states. The state ϕR,vR will be connected to the state ϕ∗,v∗ by a 2-rarefaction if, and only if, the eigenvalue λ2 is an increasing function of the ratio x/t between these states. Otherwise, we have 1-shock instead of 1-rarefaction and 2-shock instead of 2-rarefaction.When the statesϕL,vL and ϕ∗,v∗ are connected by a 1-rarefaction, the first Riemann invariant is a constant between them. On the other hand, when the states ϕ∗,v∗and ϕR,vR are connected by a 2-rarefaction, the second Riemann invariant is a constant between them. The Riemann invariants R1 and R2, associated with λ1 and λ2, are obtained from the following differential equations (that arises from (13)):(19)-λ1dϕ+dϕv=0and-λ2dϕ+dϕv=0After some calculations, we obtain(20)R1=c1ln⁡ϕ+c21ε-ln⁡ϕε--ϕ+v=constant→firstRiemanninvariantR2=-c1ln⁡ϕ-c21ε-ln⁡ϕε--ϕ+v=constant→secondRiemanninvariantIt is to be noticed that (10) was proposed looking for an easy representation for the invariants.When(21)vR-vL≥c1ln⁡ϕLϕR+c21εln⁡ϕLε-ϕRϕRε-ϕLwe ensure that the solution is 1-rarefaction/2-rarefaction. In this case, the intermediate state is obtained from (22)c1ln⁡ϕ∗+c21εln⁡ϕ∗ε-ϕ∗+v∗=c1ln⁡ϕL+c21εln⁡ϕLε-ϕL+vL-c1ln⁡ϕ∗-c21εln⁡ϕ∗ε-ϕ∗+v∗=-c1ln⁡ϕR-c21εln⁡ϕRε-ϕR+vRThis system leads us to(23)v∗=12c1ln⁡ϕL-ln⁡ϕR+c21εln⁡ϕLε-ϕL-ln⁡ϕRε-ϕR+vR+vLand, hence, ϕ∗ is obtained from one of the equations of (22).When (21) holds we always have(24)ϕL≥ϕ∗≤ϕRvL≤v∗≤vRWhen two states are connected by a discontinuity, they must satisfy the Rankine-Hugoniot jump conditions given by(25)ϕvϕ=ϕv2+pϕv=swhere “[]” denotes the jump and s denotes the shock (discontinuity) speed, as well as the entropy conditions (which ensure that the states cannot be connected by a rarefaction).The states(ϕL,vL) and (ϕ∗,v∗) will be connected by a 1-shock, if satisfying the jump conditions and (26)s1<λ1^ϕL,vLandλ1^ϕ∗,v∗<s1<λ2^ϕ∗,v∗,for  1-shockThe states(ϕ∗,v∗) and (ϕR,vR) will be connected by a 2-shock, if satisfying the jump conditions and (27)s2>λ2^ϕR,vRandλ2^ϕ∗,v∗>s2>λ1^ϕ∗,v∗,for  2-shockSince (17) holds, the entropy conditions are satisfied, provided the Rankine-Hugoniot jump conditions hold and ϕL<ϕ∗ (for the 1-shock) or ϕ∗>ϕR (for the 2-shock).The above results allow us to conclude that there are four possible solutions: 1-rarefaction/2-rarefaction, 1-shock/2-shock, 1-rarefaction/2-shock, and 1-shock/2-rarefaction, as resumed in Table1.Table 1
The four possible solutions and its consequences.
1-rarefaction/2-rarefaction→→ϕL≥ϕ∗≤ϕRvL<v∗<vR 1-rarefaction /2-shock→→ϕL>ϕ∗>ϕRvL≤v∗≥vR 1-shock /2-rarefaction→→ϕL<ϕ∗<ϕRvL≥v∗≤vR 1-shock/2-shock→→ϕL≤ϕ∗≥ϕRvL>v∗>vRIf the states(ϕL,vL) and (ϕ∗,v∗) are connected by a 1-shock while the states (ϕR,vR) and (ϕ∗,v∗) are connected by a 2-shock, then the jump conditions become(28)ϕ∗v∗-ϕLvLϕ∗-ϕL=ϕ∗v∗2-ϕLvL2+p∗-pLϕ∗v∗-ϕLvL=s1ϕRvR-ϕ∗v∗ϕR-ϕ∗=ϕRvR2-ϕ∗v∗2+pR-p∗ϕRvR-ϕ∗v∗=s2in which s1 is the 1-shock speed, while the 2-shock speed is s2. The above equations and the entropy conditions lead us to(29)ϕ∗=ϕL-ϕ∗-ϕLp∗-pLϕ∗-ϕL1ϕ∗ϕLϕ∗=ϕR-ϕR-ϕ∗pR-p∗ϕR-ϕ∗1ϕRϕ∗in which(30)pL=c12ϕL+c22ϕLε--ϕL-2c1c2ε-ln⁡ε--ϕLε-pR=c12ϕR+c22ϕRε--ϕR-2c1c2ε-ln⁡ε--ϕRε-p∗=c12ϕ∗+c22ϕ∗ε--ϕ∗-2c1c2ε-ln⁡ε--ϕ∗ε-The intermediate fluid fraction,ϕ∗, is the unique root of (31)vR-vL+ϕ∗ϕR-ϕRϕ∗pR-p∗ϕR-ϕ∗+ϕ∗ϕL-ϕLϕ∗p∗-pLϕ∗-ϕL=0while v∗ is obtained from (17), for instance, from (32)v∗=12vL+vR-ϕ∗-ϕLp∗-pLϕ∗-ϕL1ϕ∗ϕL-ϕR-ϕ∗pR-p∗ϕR-ϕ∗1ϕRϕ∗Two states are connected by a 1-shock/2-shock if(33)vL-vR≥ϕR-ϕLpR-pLϕR-ϕL1ϕRϕLand, in these cases, the following relations always hold:(34)ϕL≤ϕ∗≥ϕRvL≥v∗≥vRWhen neither (21) nor (33) holds, we have 1-rarefaction/2-shock (if ϕL>ϕR) or 1-shock/2-rarefaction (if ϕR>ϕL). These cases occur when (35)ϕR-ϕLpR-pLϕR-ϕL1ϕRϕL>vL-vR>-c1ln⁡ϕLϕR+c21εln⁡ϕLε-ϕRϕRε-ϕLThe intermediate state, in case of 1-rarefaction/2-rarefaction, is obtained from(36)c1ln⁡ϕ∗+c21ε-ln⁡ϕ∗ε--ϕ∗+v∗=c1ln⁡ϕL+c21ε-ln⁡ϕLε--ϕL+vL-c1ln⁡ϕ∗-c21ε-ln⁡ϕ∗ε--ϕ∗+v∗=-c1ln⁡ϕR-c21ε-ln⁡ϕRε--ϕR+vRIn case of 1-shock/2-shock, the intermediate state is obtained from(37)v∗=vL-ϕ∗-ϕLp∗-pLϕ∗-ϕL1ϕ∗ϕLv∗=vR-ϕR-ϕ∗pR-p∗ϕR-ϕ∗1ϕRϕ∗In case of 1-rarefaction/2-shock, the intermediate state is obtained from(38)c1ln⁡ϕ∗+c21ε-ln⁡ϕ∗ε--ϕ∗+v∗=c1ln⁡ϕL+c21ε-ln⁡ϕLε--ϕL+vLv∗=vR-ϕR-ϕ∗pR-p∗ϕR-ϕ∗1ϕRϕ∗Finally, in case of 1-shock/2-rarefaction, the intermediate state is obtained from(39)v∗=vL-ϕ∗-ϕLp∗-pLϕ∗-ϕL1ϕ∗ϕL-c1ln⁡ϕ∗-c21ε-ln⁡ϕ∗ε--ϕ∗+v∗=-c1ln⁡ϕR-c21ε-ln⁡ϕRε--ϕR+vR


## 5. An Example

In order to illustrate the features of the proposed constitutive relation (10), let us consider the following initial data: (40)ϕ,v=0.85,10.0,for-∞<x<0ϕ,v=0.60,0.0,for0<x<∞represented in Figure1.Figure 1
Initial data.All the results presented in the sequence are obtained assuming a porosityε=0.9 with a supersaturation of 5% (in other words with δ=0.045).Figures2–5 represent the fluid fraction behavior as a function of the position for four time instants, considering some selected values of c1 and c2 whereas Figures 6–9 depict the velocity behavior for these constants with the same values of c1 and c2. In both cases the following values are considered for the constants c1: c1=1.0, c1=10.0, c1=20.0, and c1=30.0, while the constant c2 assumes the following values: c2=0.5, c2=1.0, and c2=2.0.Figure 2
Fluid fraction behavior alongx-axis for selected values of c1 and c2 at t=0.25.Figure 3
Fluid fraction behavior alongx-axis for selected values of c1 and c2 at t=0.50.Figure 4
Fluid fraction behavior alongx-axis for selected values of c1 and c2 at t=0.75.Figure 5
Fluid fraction behavior alongx-axis for selected values of c1 and c2 at t=1.00.Figure 6
Velocity behavior alongx-axis for selected values of c1 and c2 at t=0.25.Figure 7
Velocity behavior alongx-axis for selected values of c1 and c2 at t=0.50.Figure 8
Velocity behavior alongx-axis for selected values of c1 and c2 at t=0.75.Figure 9
Velocity behavior alongx-axis for selected values of c1 and c2 at t=1.00.In Figures2–5 there are two dashed lines: the upper one representing ε- and the lower one representing ε. It is worthwhile to observe that relationship (10) avoids any value of the fluid fraction greater than ε-. Such feature holds for any choice of δ, provided it is positive valued.Relationships like (7), proposed by Martins-Costa el al. [3], as well as the relationship presented in (6) (see [1]) do not ensure an upper bound for the fluid fraction and may lead us to unrealistic results.


## 6. Final Remarks

The constitutive relation for the pressure proposed in this work, besides giving rise to explicit formulae for the Riemann invariants, provides an upper bound for the fluid fraction. So, any simulation carried out based on this relationship is free from presenting inconsistent physical results. Hence, we have a reliable tool for describing the transition saturated-unsaturated (and vice versa) in rigid porous media.



---

*Source: 2347809-2019-06-27.xml*