#### 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.

#### 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 bywhere is the fluid constituent mass density (local ratio between the fluid constituent mass and the corresponding volume of mixture), is its velocity in the mixture, is the partial stress tensor associated with it (analogous to Cauchy stress tensor in Continuum Mechanics), is the momentum supply acting on it due to its interaction with the remaining constituents of the mixture, and is the body force. The momentum supply to the mixture must be zero; in other words: .

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 and the fluid mass density (Continuum Mechanics viewpoint) ; in other words, and represents the porous matrix porosity.

The dynamic interaction among the constituents is accounted for by the momentum supply which can be represented bywhere is the fluid viscosity (in a Continuum Mechanics approach), represents the specific permeability of the porous matrix, and is a diffusion coefficient.

The first term of 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 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 towhere is a pressure (assumed constant while the flow is unsaturated) and is the identity tensor.

At this point all the quantities are assumed to depend solely on time* t* and on position* x* and the only nonvanishing component of the fluid constituent velocity, , is denoted by* v*, so (1)-(3) give rise to the following system:

In order to write system (4) in a more convenient form, a reference pressure is defined such that

Since 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]:in which is a positive constant, represents the fluid fraction, represents the velocity, represents the pressure (for simplicity ), and the fluid density 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]:to encompass saturated and unsaturated flows within the same relationship, allowing a small supersaturation.

The main idea consists of choosing a very large constant 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 it is not possible to establish, a priori, a bound for the fluid fraction. In other words, we have that

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:in which and are positive constants and is a positive constant that can be very small. From the above equation it is easy to see thatand, so, we have an a priori limit for the supersaturation, even for very large pressures, since

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:

The Riemann Problem associated with (13) arises when, for all , the following initial data is assumed:

where are constants.

The eigenvalues for this problem are given by

Taking into account (10) we may write

Since the second derivative of with respect to , given byis always positive valued, we have that the solution (in a generalized sense) of the Riemann Problem consists of connecting the left state and the right state to an intermediate state , by means of rarefactions or shocks, according to one of the four ways listed below:

The state will be connected to the state by a 1-rarefaction if, and only if, the eigenvalue is an increasing function of the ratio between these states. The state will be connected to the state by a 2-rarefaction if, and only if, the eigenvalue is an increasing function of the ratio between these states. Otherwise, we have 1-shock instead of 1-rarefaction and 2-shock instead of 2-rarefaction.

When the states and are connected by a 1-rarefaction, the first Riemann invariant is a constant between them. On the other hand, when the states and are connected by a 2-rarefaction, the second Riemann invariant is a constant between them. The Riemann invariants and , associated with and , are obtained from the following differential equations (that arises from (13)):

After some calculations, we obtain

It is to be noticed that (10) was proposed looking for an easy representation for the invariants.

When we ensure that the solution is 1-rarefaction/2-rarefaction. In this case, the intermediate state is obtained from

This system leads us to and, hence, is obtained from one of the equations of (22).

When (21) holds we always have

When two states are connected by a discontinuity, they must satisfy the Rankine-Hugoniot jump conditions given bywhere “” denotes the jump and denotes the shock (discontinuity) speed, as well as the entropy conditions (which ensure that the states cannot be connected by a rarefaction).

The states and will be connected by a 1-shock, if satisfying the jump conditions and

The states and will be connected by a 2-shock, if satisfying the jump conditions and

Since (17) holds, the entropy conditions are satisfied, provided the Rankine-Hugoniot jump conditions hold and (for the 1-shock) or (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 Table 1.

If the states and are connected by a 1-shock while the states and are connected by a 2-shock, then the jump conditions becomein which is the 1-shock speed, while the 2-shock speed is . The above equations and the entropy conditions lead us toin which

The intermediate fluid fraction, , is the unique root of while is obtained from (17), for instance, from

Two states are connected by a 1-shock/2-shock ifand, in these cases, the following relations always hold:

When neither (21) nor (33) holds, we have 1-rarefaction/2-shock (if ) or 1-shock/2-rarefaction (if ). These cases occur when

The intermediate state, in case of 1-rarefaction/2-rarefaction, is obtained from

In case of 1-shock/2-shock, the intermediate state is obtained from

In case of 1-rarefaction/2-shock, the intermediate state is obtained from

Finally, in case of 1-shock/2-rarefaction, the intermediate state is obtained from

#### 5. An Example

In order to illustrate the features of the proposed constitutive relation (10), let us consider the following initial data:

represented in Figure 1.

All the results presented in the sequence are obtained assuming a porosity with a supersaturation of 5% (in other words with ).

Figures 2–5 represent the fluid fraction behavior as a function of the position for four time instants, considering some selected values of and whereas Figures 6–9 depict the velocity behavior for these constants with the same values of and . In both cases the following values are considered for the constants : , , , and , while the constant assumes the following values: , , and .

In Figures 2–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.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors Maria Laura Martins-Costa and Rogério Martins Saldanha da Gama gratefully acknowledge the financial support provided by Brazilian agency CNPq while the author Daniel Cunha da Silva gratefully acknowledges the financial support provided by Brazilian agency CAPES.