Modeling a Geothermal Field: A Not-Trivial Starter Kit

Luca Meacci

Instituto de Ciências Matemáticas e de Computação, Universidade de São Paulo, Campus de São Carlos, São Carlos, Brazil.

DOI: 10.4236/oalib.1106836                

Abstract

In this paper, we discuss the assumptions, the balances, and the constitutive relationships in order to provide a set of tools for the mathematical modeling of a geothermal system. In particular, we present a model for pressure and saturation supposing that: 1) the geothermal fluid flows in a porous medium, 2) it is composed of pure water, 3) the simultaneous presence of the gaseous (vapor) and liquid phases occurs, and 4) the effects of capillarity action can be introduced.

Keywords

Geothermal System, Porous Medium, Pure Water Phases, Constitutive Relationships, Capillarity, Partial Differential Equations

1. Introduction

Geothermal energy is the energy in the form of heat contained within the Earth that can be extracted using manpower [1]. The geothermal energy exploitation is a topic of great interest and has potential for the future. It is estimated that the expected geothermal target for the year 2050 is 140 GWe, i.e., it would be possible to produce from geothermal up 8.3% of total world electricity production, serving 17% of world population [2]. Furthermore, “40 countries (located mostly in Africa, Central/South America, Pacific) can be 100% geothermal powered” and “the overall CO₂ saving from geothermal electricity can be in around 1000 million tons per year” [3].

A challenge with which we strive is trying to achieve this goal with description, understanding, and prediction of geothermal system behavior. Such a system is been defined schematically as “convecting water in the upper crust of the Earth, which, in a confined space, transfers heat from a heat source to a heat sink, usually the free surface” [4]. The processes involved in a geothermal reservoir are complex. Since the 1970s developments of mathematical models of these dynamics have been used for simulation purposes [5] [6] [7]. During the years, a lot of knowledge was produced on this topic [8] [9] [10], and the continuous publication of recent developments is a demonstration not only of the importance of this field of study but also of its unchanged interest [11] [12].

Nevertheless, this amount of information and concepts presented in literature may leave, especially at first glance, the researcher with an unclear perception of the main underlying and characterizing elements. The purpose of this work is to provide a short presentation of a general mathematical model that allows in a simple (but not simplistic) way to identify which balances and assumptions can be considered to model a geothermal system. The aim is therefore to introduce, albeit under simplified hypotheses (fluid composed only of pure water and isotropic system), a not trivial starter kit for modeling a geothermal field. To facilitate the reading of the work, we display in Table 1 the symbols of the physical quantities and the notations used in the text. We discuss all the elements concerning the problem, the assumption of a porous medium, that schematizes the fractured permeable rocks of the reservoir, as well as the consequences of the presence of more fluid phases. We develop and describe the mass and moment conservation equations for each phase contained in the geothermal system. We argue for the constitutive relationships and assumptions, in absence and presence of capillarity, in order to finally reduce the balance equations to a nonlinear Partial Differential Equations (PDEs) problem with pressure and saturation unknowns and to achieve the closure of the model.

2. Basics Ingredients and Methods

The basic method for mathematical model development and physics modeling of a geothermal reservoir is represented by the following equation balance equations. An elementary balance equation for the generic quantity M can be written in the following form [13]:

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_V}M\text{d}V=-{\displaystyle {\int }_{\Gamma }}F\cdot n\text{d}\Gamma +{\displaystyle {\int }_V}q\text{d}V\mathrm{.}$ (1)

The integration takes place on an arbitrary V sub-domain enclosed by a $\Gamma$ surface. $F$ denotes the flow and q the source and absorption terms. Finally, is the normal versus outer to surface. The local form of Equation (1) is given by the arbitrary of the V volume, using the divergence theorem. The result can be obtained as:

$\frac{\partial }{\partial t}M+\nabla \cdot F=q\mathrm{.}$ (2)

Table 1. Symbols and notations.

However, the fluid is not free in our case, but it is under dynamics conditions within the rock fractures that form the reservoir of the geothermal system. Therefore, in the following system, a non-deformable porous medium is introduced. Let us introduce how we can describe the balances in a porous medium which are filled with pure water in gaseous and liquid state.

A porous medium is defined as any material in which a solid matrix and an empty space are present; or rather, a part where the fluid can penetrate [14] [15]. Therefore, the empty space can be occupied by both liquid and vapor H₂O. The solid matrix is distributed throughout the porous medium (see Figure 1). This fact implies that for sufficiently large volumes of domain, arbitrarily taken within the porous medium, we will always find the presence of solid matter. We refer to such a volume as an Arbitrary Elementary Volume (AEV).

In the porous medium, the H₂O flows through a complex network of pores structure which constitutes the space of the material. This flow is therefore delimited by the microscopic interface between solid and fluid. This configuration can be essentially analyzed in two different ways:

• Microscopic Level: Conceptually, the H₂O flowing in a porous medium could be studied at a microscopic level, that is by focusing attention on what happens in the single pore. In this case, the fluid is modeled as a continuous means subjected to specific boundary conditions. These conditions are defined on the interface that delimits the domain under consideration. However, this type of approach is generally theoretical because of the complex geometry of the interface. Besides, the mathematical model solutions from measurements at these levels can hardly be checked.

Figure 1. Example of an arbitrary domain V. The presence of rock is labelled with the letter r while vapor and liquid water are labelled with the letters v and l, respectively. The Figure is modified from [16].

• Macroscopic Level: To overcome the difficulties mentioned above, another level of description is needed, namely the macroscopic one. On this scale, the quantities involved can be easily measured. To obtain a description of the flow at this level, the so-called continuous approach is adopted. According to this approach, the real porous medium is replaced by a fictitious model in which each phase (liquid and gaseous) is seen as a continuous medium that fills the entire AEV. Let us proceed with the concept of macroscopic quantity in the porous medium. Suppose, for instance, we aim to define $\rho$ quantity (the quantity may be the volume fraction occupied by grains, the mass of liquid water, etc.) of the medium at the point $x\in {ℝ}^3$ (if the domain is considered three-dimensional). The procedure is given as follows: We consider an ideal spherical ball centered in the point $x$ and we measure the average value of $\rho$ on such a sphere varying the radius length. It is expected to obtain an oscillating behavior of the measurement values; precisely, due to the random nature of the medium (being that as more or less solid parts can be included in the sphere). Nevertheless, by observing the graph of the measurements (a diagram of which is shown in Figure 2), we can notice an interval of the radius of the ball for which the $\rho$ value stabilizes. We therefore proceed in choosing the volume with such radius as the Representative Elementary Volume (REV). This means that, the phenomenon is studied within such a REV, and each quantity of the medium in a point $x$ must be considered as average quantity in the REV centered in $x$ [17].

The theory exposed in this paper refers to the macroscopic level, and the following quantities are considered as average quantities. In Section 3, we describe the balance equations which are used in continuum mechanics for these average quantities of the medium.

Figure 2. Example of a graph of the measures of the quantity $\rho$ performed on a small spherical ball of variable radius, i.e., for an increasing volume V.

3. Balance Equations

Let us first consider the mass and momentum balance equations for a system in which the geothermal fluid is composed of pure liquid and gaseous H₂O. The mass of H₂O presented in the volume unit of the porous medium is given by:

$M_{{\text{H}}_{\text{2}}\text{O}}={\displaystyle \underset{\beta }{\sum }}{S}_{\beta }{\rho }_{\beta }\varphi ,$ (3)

where

• $\beta$ is the index representing the phase ( $\beta =v$ for vapour and $\beta =l$ for liquid water);

• $\varphi$ is the porosity of the medium, that is, the fraction of volume not occupied by the solid matrix;

• $S_{\beta }$ is the saturation of phase $\beta$, that is the volume fraction occupied inside the pores by the water of the corresponding state. The saturation variables comply with the constraint

$S_v+S_l=1;$

• ${\rho }_{\beta }$ is the density of phase $\beta$, that is the mass of $\beta$ phase per volume unit occupied by the $\beta$ phase. For example, ${\rho }_l=1\text{gr}/{\text{cm}}^{\text{3}}$.

The mass flow $F$ is the sum of all the phases

$F={\displaystyle \underset{\beta }{\sum }}{F}_{\beta },$ (4)

whereas each $F_{\beta }$ flow of the phase $\beta$ is given by

$F_{\beta }={\rho }_{\beta }{S}_{\beta }{u}_{\beta }\mathrm{,}$ (5)

by the product between the density and the volumetric flux $u_{\beta }$ of the fluid in the $\beta$ phase.

The fluids move in response to forces in which the pressure forces, gravity force, and resistance forces are most important due to the viscosity and friction within the grains of the porous medium. The fundamental law of the cause of fluid fluxes in porous media was established by Henri Darcy yet in 1856 [18] based on a series of water flow experiments in a horizontal sand-filled tube. He realized that the specific flow rate of the tube depended on the pressure variation of the tube length unit which is called a pressure gradient. By generalizing Darcy’s experiments, the specific volumetric flow rate is obtained as $u$ proportional to the pressure gradient according to the law

$u=-\frac{K}{\mu }\nabla P,$ (6)

being P the pressure of H₂O. The sign minus is needed because the fluid flows from more high-pressure areas to low-pressure areas. In Equation (6), $K$ is called permeability tensor and measures the ease of the fluid to move inside the porous medium (depending on the porosity and the microscopic geometric structure). Finally, $\mu$ denotes the viscosity, i.e., the resistance of a fluid to a change in shape. We can also consider the contribution of gravitational force [10]; the full expression becomes:

$u=-\frac{K}{\mu }\left(\nabla P-\rho g\right)\mathrm{.}$ (7)

We remark that $u$ represents the volume of fluid that crosses a unit area section per unit time. The volumetric flow $u$ is also called Darcy’s velocity; but, the volumetric flow is not the average velocity with which fluid particles are moving. This quantity is known as pore velocity $v$ and is defined according to the following equation:

$u=\varphi v\mathrm{.}$ (8)

This definition implies that in a medium with low porosity, the velocity $v$ may be greater as compared to the volumetric flow $u$.

Let us now proceed to specify the meaning of pressures of the H₂O and its phases. As specified, P denotes the pressure of H₂O while with $P_v$ and $P_l$ the pressures of H₂O in vaporous state and liquid phase, respectively. By way of definition it follows that

$\begin{array}{l}{S}_l=1\iff P_l=P,\text{and}{P}_v=0,\\ S_l=0\iff P_l=0,\text{and}{P}_v=P.\end{array}$

If $0<S_l<1$ holds, $P_l$ and $P_v$ have to be specified in a constitutive manner. In that case, the liquid and vapor are considered as a two-component mixture with P representing mixture pressure in its totality. In particular, the partial stresses of the components in the fluid mixture (see [13] ) are given by

${\tilde{T}}_l=S_l{P}_{l}I\text{and}{\tilde{T}}_v=S_v{P}_{v}I,$

while the total stress $T$ exerted by the mixture fluid is

$T={\tilde{T}}_l+{\tilde{T}}_v=S_l{P}_{l}I+S_v{P}_{v}I,$ (9)

where $T=PI$.

If, as in a geothermal reservoir occurs, the H₂O is present in several phases then it is proved that each phase flows moved by its own pressure gradient and by its own weight [10]. The effective permeability of the medium relative to the single-phase is reduced with the comparison where the phase occupies all available space. This effect is taken into account by introducing the reduction permeability factors or the relative permeabilities $k_{r\beta }$ ( $\beta =v,l$ ).

The Darcy’s expression for the single $\beta$ phase becomes

$u_{\beta }=-K\frac{k_{r\beta }}{{\mu }_{\beta }}\left(\nabla P_{\beta }-{\rho }_{\beta }g\right),$ (10)

with, resuming:

• $K$ the absolute permeability tensor of the medium ( $\left[K\right]=\left[L^2\right]={\text{m}}^2$ );

• $k_{r\beta }$ the effective permeability of the single phase $\beta$, $\left[k_{r\beta }\right]=[-],0\le k_{r\beta }\le 1$. It is a reduction factor of the permeability because only a part of the space available in the pores is occupied by the $\beta$ phase. It is indeed true that

$k_{r\beta }\left(S_{\beta }=1\right)=1\text{and}{k}_{r\beta }\left(S_{\beta }=0\right)=0,$

while, for intermediate saturation values, these functions are defined experimentally. It also applies

${\displaystyle \underset{\beta }{\sum }}{k}_{r\beta }=1;$

• ${\mu }_{\beta }$ the viscosity of the $\beta$ phase, usually $\left[{\mu }_{\beta }\right]=\text{Pa}\cdot \text{s}$.

From Equation (4) and Equation (5) the total mass flow is obtained as $F$ of H₂O, given by

$F=-{\displaystyle \underset{\beta }{\sum }}K\frac{k_{r\beta }}{{\mu }_{\beta }}{\rho }_{\beta }\left(\nabla P_{\beta }-{\rho }_{\beta }g\right).$ (11)

It can be stated that the local mass conservation equation [17] [19] leads to the following form:

$\frac{\partial }{\partial t}\left\{\varphi \left[S_l{\rho }_l+\left(1-S_l\right){\rho }_v\right]\right\}+\nabla \cdot \left({\rho }_l{S}_l{u}_l+{\rho }_v{S}_v{u}_v\right)=0,$ (12)

or as an alternative

$\frac{\partial }{\partial t}\left\{\varphi \left[S_l{\rho }_l+\left(1-S_l\right){\rho }_v\right]\right\}+\nabla \cdot \left({\rho }_l{S}_l\varphi v_l+{\rho }_v{S}_v\varphi v_v\right)=0,$ (13)

where any type of wells or sources is not considered. There are three unknown variables in this equation, which are $S_l$, $P_l$ and $P_v$. Note that the vapor density ${\rho }_v$, albeit is variable, does not introduce an unknown variable because it depends on the vapor pressure $P_v$. We can assume to apply the ideal gas law

$\frac{P_v}{{\rho }_v}=rT,$ (14)

where T is temperature in Kelvin degrees (K) and r is the specific gas constant, in H₂O case equal to $r=4.6\times {10}^2\text{J}/\text{kg}\cdot \text{K}$.

We provide three remarks concerning Equation (12), according to three different cases that can occur.

Remark 1. First, if the whole water is in the liquid state ( $S_l=1$ and $S_g=0$ ), then the Equation (12) is reduced to

$\nabla \cdot u_l=\mathrm{0,}\Rightarrow \nabla \cdot \left(\varphi \frac{K}{{\mu }_l}\left(\nabla P_l-{\rho }_{l}g\right)\right)=0\mathrm{,}$

that is an elliptic equation in the unknown $P_l$ : the classic equation of porous media. Resolving it and specifying appropriate boundary conditions, the equation provides the pressure field in the domain, and using this Equation (10) we can evaluate the water flow.

Remark 2. Let us suppose now, in this second case, that $S_l=0$ while $S_v=1$. We consider a system in which all the water is in the gaseous state. We obtain

$\frac{\partial }{\partial t}\left(\varphi {\rho }_v\right)+\nabla \cdot \left({\rho }_v{u}_v\right)=0,\Rightarrow \frac{\partial }{\partial t}\left(\varphi {\rho }_v\right)-\nabla \cdot \left({\rho }_v\varphi \frac{K}{{\mu }_v}\left(\nabla P_v-{\rho }_{v}g\right)\right)=0,$

with ${\rho }_v$ and $P_v$ linked by Equation (14). A not-uniform parabolic equation in the unknown $P_v$ is formulated here. If it is assumed that $T\left(x\right)$ is known, the solution of this equation provides the vapor pressure field in the domain.

Remark 3. The third case occurs when the water is present in both the gaseous and liquid phase. According to this configuration, we must reduce the number of unknowns, which for now are: $S_l$ (or $S_v$ ), $P_l$, and $P_v$. Through constitutive equations, we can introduce a mutual relationship between the two pressures and a second relationship that binds one of them to the saturation degree.

In Sections 4 and 5, we present therefore a constitutive model depending on whether the capillarity effects are or are not neglected, respectively.

4. Constitutive Model in the Absence of Capillarity

The state of the H₂O is determined via the Clapeyron diagram [20]. The Clapeyron equation describes the variation of pressure with temperature along the equilibrium curve between two phases ( $\alpha$ and $\beta$ )

$\frac{\text{d}P}{\text{d}T}=\frac{\Lambda }{T\left({\mathcal{V}}_{\beta }^{mol}-{\mathcal{V}}_{\alpha }^{mol}\right)}\mathrm{,}$ (15)

being $\Lambda$ the latent heat (per unit of mass) of transition from one phase to another, and ${\mathcal{V}}^{mol}$ is the molar volume of the two phases. An approximate solution of this equation in the case of the liquid-vapor transition of pure H₂O is given after fitting of the experimental data [21]. We obtain a function of the saturated vapor pressure $P^{\mathrm{<mo>\; *\; </mo>}}$ which is entitled as Clapeyron pressure, given by the expression

$P^{*}\left(T\right)=961.7\exp \left\{17.35\frac{T-273}{T}\right\}\text{Pa}.$ (16)

If, at a fixed temperature, the pressure P is such that $P>P^{*}$, then the H₂O is in the liquid phase and we have $P=P_l$. Conversely, if $P<P^{*}$ then the H₂O is in the gaseous phase (vapor), and the result is $P=P_v$. Unlike again, if $P=P^{*}$ then there is the coexistence of the liquid and gaseous phase. In such a situation, the pressure of water (now viewed as a mixture of two phases) exactly corresponds to $P^{\mathrm{<mo>\; *\; </mo>}}$. Under this condition, $P_l$ and $P_v$ must be specified. By treating the fluid as a liquid-vapour mixture, from Equation (9) we have

$\left(1-S_l\right)P_v+S_l{P}_l=P^{*}.$ (17)

Therefore, if we neglect the capillary effects, the simplest constitutive assumption, encountered in literature [17], is

$P_v=P^{*},$ (18)

$P_l=P^{*}.$ (19)

In this case, the variable that describes the system behavior is therefore the saturation of one of these phases. For instance, we can choose $S_l$ ; the density of the phases follows accordingly. The liquid phase is calculated with constant density while the vapor density is given by the following Equation (14) in this way,

${\rho }_v={\rho }_v^{*}\left(T\right)\text{with}{\rho }_v^{*}\left(T\right)=\frac{P^{*}\left(T\right)}{rT}.$

Therefore, the only state variable is precisely the saturation of one of the phases.

We aim to generalize these results, by introducing the following constitutive relationship for $S_l$ in terms of P. We describe $S_l$ with the graph ${\hat{S}}_l\left(P\right)$ on the pressure, i.e.

$S_l\in {\hat{S}}_l\left(P\right)\mathrm{,}$ (20)

being

${\hat{S}}_l\left(P\right)=\{\begin{array}{ll}1\hfill & \text{if}P>P^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\hfill \\ \left(\mathrm{0,1}\right)\hfill & \text{if}P=P^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\hfill \\ 0\hfill & \text{if}P>P^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\hfill \end{array}$ (21)

as shown in Figure 3. Concerning the pressures $P_l$ and $P_v$, we obtain

$P_v=\left(1-\mathcal{H}\left(P-P^{\mathrm{<mo>\; *\; </mo>}}\right)\right)P\mathrm{,}$ (22)

$P_l=\mathcal{H}\left(P-P^{\mathrm{<mo>\; *\; </mo>}}\right)P\mathrm{,}$ (23)

where $\mathcal{H}$ is the Heaviside function

$\mathcal{H}\left(x\right)=\{\begin{array}{ll}1\hfill & \text{if}x\ge 0\mathrm{,}\hfill \\ 0\hfill & \text{if}x<0\mathrm{.}\hfill \end{array}$ (24)

Summing up, if capillary effects are neglected, the H₂O flow in the porous medium is described by the Equation (12), read as an equation in P, with $S_l$ given by Equation (21), and with $P_v$ and $P_l$ given by Equations (22) and (23), respectively.

5. Constitutive Model in the Presence of Capillarity

When the liquid and gaseous phases coexist and phenomena of capillarity are unable to neglect, a notion of capillary pressure [22] [23] can be introduced as;

$P_c=P^w-P^{nw},$ (25)

Figure 3. Graph of ${\hat{S}}_l\left(P\right)$.

where $P_w$ is the wetting fluid pressure which tends to adhere with the pour walls while $P^{nw}$ is the fluid pressure which tends to minimize the surface contacting the porous medium. In the case of the H₂O in the form of vapor and liquid in the soil, the wetting phase is the liquid phase. Therefore

$P_c=P_l-P_v,$ (26)

with $P_c$ described by the following so-called Laplace formula

$P_c=\gamma \langle \kappa \rangle \cos \alpha ,$ (27)

with

• $\gamma =\gamma \left(T\right)$ the surface tension;

• $\langle \kappa \rangle$ the average curvature of the liquid-vapour interface throughout the REV;

• $\cos \alpha$ the contact angle between the liquid and the pore wall.

Let us now consider an isotherm line in the Clapeyron plane, according to Equation (15). At the phase transition (assuming as always, the equilibrium between the phases), we observe that

${\mathcal{V}}_v^{mol}\text{d}{P}_v={\mathcal{V}}_l^{mol}\text{d}{P}_l\mathrm{,}$ (28)

being ${\mathcal{V}}_{\beta }^{mol}$, $\beta =l,v$ the molar H₂O volume of the respective phases. By differentiating Equation (26) we have

$\text{d}{P}_c=\text{d}{P}_l-\text{d}{P}_v=\left(\frac{{\mathcal{V}}_v^{mol}}{{\mathcal{V}}_l^{mol}}-1\right)\text{d}{P}_g=\frac{{\mathcal{V}}_v^{mol}}{{\mathcal{V}}_l^{mol}}\left(1-\frac{{\mathcal{V}}_l^{mol}}{{\mathcal{V}}_v^{mol}}\right)\text{d}{P}_v.$

Now, if we are “sufficiently” far from the critical point, it is legitimate to consider ${\mathcal{V}}_l^{mol}\ll {\mathcal{V}}_v^{mol}$, and then we can state that

$\text{d}{P}_c=\frac{{\mathcal{V}}_v^{mol}}{{\mathcal{V}}_l^{mol}}\text{d}{P}_v\mathrm{.}$

Let us now assume that vapor is an ideal gas; it can be considered the equation of state as,

$P_v{\mathcal{V}}_g^{mol}=\mathcal{R}T\mathrm{,}$ (29)

and we get accordingly

$\text{d}{P}_c=\frac{\mathcal{R}T}{{\mathcal{V}}_l^{mol}}\frac{\text{d}{P}_v}{P_v}\mathrm{.}$ (30)

Now we integrate Equation (30) remarking that $P_c=0$ corresponds to $P_v=P^{*}$,

${\displaystyle \underset{0}{\overset{P_c}{\int }}}\text{d}{P}_c=\frac{\mathcal{R}T}{{\mathcal{V}}_l^{mol}}{\displaystyle \underset{P^{\mathrm{<mo>\; *\; </mo>}}}{\overset{P_v}{\int }}}\frac{\text{d}{P}_v}{{P^{\prime }}_v}\mathrm{,}$

in order to arrive at

$P_v=P^{\mathrm{<mo>\; *\; </mo>}}exp\left\{\frac{{\mathcal{V}}_l^{mol}}{\mathcal{R}T}{P}_c\right\}\mathrm{,}$ (31)

that is well-known Kelvin equation. In particular, being

${\mathcal{M}}_{{\text{H}}_{\text{2}}\text{O}}^{mol}={\rho }_l{\mathcal{V}}_l^{mol}\mathrm{,}$

the Equation (31) can be written as

$P_v=P^{*}\exp \left\{\frac{{\rho }_l}{rT}{P}_c\right\}.$ (32)

The model so defined has an additional variable: $P_c$. It is necessary to introduce a constitutive law that links capillary pressure to some other variable. Following classical literature [14], it can be assumed that $P_c$ is expressible in terms of the liquid saturation $S_l$. In this line, various functional forms can be determined (generally based on the best fit of experimental data) from the literature. In most cases (or rather, the most used in a hydro-geological context) Van Genuchten and Corey expressions [24] cab be adopted. For now, we just write

$P_c=P_c\left(S_l\right)\mathrm{,}\text{con}\{\begin{array}{ll}{P}_c\left(S_l\right)<0,\hfill & \text{se}0<S_l<1,\hfill \\ P_c\left(S_l\right)=0,\hfill & \text{se}{S}_l=1,\hfill \end{array}$ (33)

recalling that it is required to specify $P_c\left(S_l\right)$ for $0<S_l<1$.

An important estimation of capillary pressure is obtained from Equation (27). Using the experimental tabulations which reproduce this relationship and the typical values of $\langle \kappa \rangle$ and $\gamma$, it is estimated that the $P_c$ inside a capillary and can be obtained as a comparison parameter which confutes or not the approximation of capillary absence. For example, the Young-Laplace formula [25] for one capillary can be used

$P_c=\frac{2\gamma \cos \alpha }{\eta },$ (34)

where $\eta$ is the radius of the capillary section, to evaluate some reference values. In Table 2, we estimate the capillary pressure, approximating $\cos \alpha$ to the unit, changing the temperature in the reservoirs, and increasing the capillary radius. The data on the surface tension of the H₂O is taken from [26].

Thus, in the case of capillary not negligible, the state of the system is described by only two variables: P and $S_l$, with the liquid phase pressure $P_l$ linked to P and to $S_l$ from this relationship

Table 2. Capillary pressure estimates (expressed in Pa) according to the Young-Laplace Equation (34) for some values of temperature and radius of the capillary and approximating $\cos \alpha \approx 1$.

$P_l=\{\begin{array}{ll}P\hfill & \text{if}P>P^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\hfill \\ P_c\left(S_l\right)+P^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)exp\left\{\frac{{\rho }_l}{rT}{P}_c\left(S_l\right)\right\}\hfill & \text{if}P=P^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\hfill \end{array}$ (35)

being $P_c\left(S_l\right)$ given by Equation (33).

Considering the capillarity action, the mass and momentum balance problem thus arrives at the following mathematical model:

$\{\begin{array}{l}\frac{\partial }{\partial t}\left\{\varphi \left[S_l{\rho }_l+\left(1-S_l\right){\rho }_v\right]\right\}+\nabla \cdot \left({\rho }_l{S}_l{u}_l+{\rho }_v\left(1-S_l\right)u_v\right)=0,\\ u_l=-K\frac{k_{rl}}{{\mu }_l}\left(\nabla P_l-{\rho }_{l}g\right),\\ u_v=-K\frac{k_{rv}}{{\mu }_v}\left(\nabla P_v-{\rho }_{v}g\right),\\ P_v=P^{*}\exp \left\{\frac{{\rho }_l}{rT}{P}_c(S_l)\right\},\\ P_l=P_c+P_v,\\ P_c=P_c\left(S_l\right),\\ S_l+S_v=1,\end{array}$ (36)

or, expressed in function of the pore velocity

$\{\begin{array}{l}\frac{\partial }{\partial t}\left\{\varphi \left[S_l{\rho }_l+\left(1-S_l\right){\rho }_v\right]\right\}+\nabla \cdot \varphi \left({\rho }_l{S}_l{v}_l+{\rho }_v\left(1-S_l\right)v_v\right)=0,\\ v_l=-K\frac{k_{rl}}{\varphi {\mu }_l}\left(\nabla P_l-{\rho }_{l}g\right),\\ v_v=-K\frac{k_{rv}}{\varphi {\mu }_v}\left(\nabla P_v-{\rho }_{v}g\right),\\ P_v=P^{*}\exp \left\{\frac{{\rho }_l}{rT}{P}_c\left(S_l\right)\right\},\\ P_l=P_c+P_v,\\ P_c=P_c\left(S_l\right),\\ S_l+S_v=1.\end{array}$ (37)

The model must be closed with the relative permeability curves specification

$k_{rl}=k_{rl}\left(S_l\right)\mathrm{,}$

$k_{rv}=k_{rv}\left(S_v\right)=k_{rv}\left(1-S_l\right),$

while the permeability dependence of the structure is given by the definition of the permeability tensor K.

6. Conclusions

The main purpose of the paper was not limited to reach a generic model of a geothermal system that we presented as a result. The significance of this work is recognized in the approach used in building the model. We think the aim of this work lies in how we can obtain such model. The paper was written in dissertation style. It treats how to work with physical quantities in porous media, how to consider the consequences of the various hypotheses, and how to insert the constitutive relationships. Inter alia, we gave attention to the clarity of the steps and the reproducibility of the deductions.

The main result of our findings is to propose a generalised mathematical model for pressure and saturation of a two-component H₂O mixture flow. This finding is consistent with the sub-problems in which the fluid consists only of a single phase. In the case of the presence of pure water in the liquid state only, the problem is reduced to the classic elliptic PDEs problem of a liquid fluid in a porous medium. When the fluid is completely gaseous, we obtain a parabolic PDEs problem that can be solved if the system temperature is known. Now, in the case of simultaneous coexistence of the two phases, the addition of a constitutive model is necessary. Under the hypothesis of neglecting the effects due to capillarity, we adopted the application of the Clapeyron pressure and the law of perfect gases, to obtain a constitutive law for liquid saturation that depends on pressure, and using the Heaviside function. Finally, assuming to consider the effects of capillarity, we took advantage of the fact that capillary pressure can be expressed as a function of saturation. Therefore, we have expressed a relationship between the pressure of the liquid phase and the saturation, with the introduction of the measurement estimates of the capillary pressure, according to the Young-Laplace formula.

The various equations and the described methods can be used as a set of tools, a not-trivial starter kit, to understand and operate with models for geothermal basins, even by a readership who does not have strictly physical and mathematical skills. For these reasons, we believe this work can contribute to potential future directions [12], specifically to the topic of current interest, represented by the numerical simulation of the evolution of the geothermal system [27].

Acknowledgements

The author would like to express his gratitude to Professor Angiolo Farina (Università Degli Studi di Firenze, Italy) for having introduced him to the study of porous media and for the very important contribution to the realization of this work. The author also acknowledges with appreciation the foundation CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) of the Ministry of Education of the Federal Republic of Brazil for receiving economic support within the funding program (Financial code PROEX-9740044/D) of academic excellence PROEX (Programa de Excelência Acadêmica).

Conflicts of Interest

The author declares no conflicts of interest.

References

[1]	Dickson, M.H. and Fanelli, M. (2013) Geothermal Energy: Utilization and Technology. Routledge, London. https://doi.org/10.4324/9781315065786
[2]	Bertani, R. (2003) What Is Geothermal Potential. IGA News, 53, 1-3.
[3]	Bertani, R. (2016) Geothermal Power Generation in the World 2010-2014 Update Report. Geothermics, 60, 31-43. https://doi.org/10.1016/j.geothermics.2015.11.003
[4]	Hochstein, M.P. (1990) Classification and Assessment of Geothermal Resources. Small Geothermal Resources: A Guide to Development and Utilization. UNITAR, New York, 31-57.
[5]	Coats, K., et al. (1977) Geothermal Reservoir Modelling. SPE Annual Fall Technical Conference and Exhibition, Denver, 9-12 October, SPE-6892-MS. https://doi.org/10.2118/6892-MS
[6]	Donaldson, I.G. and Sorey, M.L. (1979) The Best Uses of Numerical Simulators. In: Proceedings of the 5th Workshop on Geothermal Reservoir Engineering, Stanford University, Palo Alto, Vol. 241.
[7]	Mercer Jr., J.W. and Pinder, G.F. (1973) Galerkin Finite-Element Simulation of a Geothermal Reservoir. Geothermics, 2, 81-89. https://doi.org/10.1016/0375-6505(73)90011-4
[8]	Bodvarsson, G.S., Pruess, K. and Lippmann, M.J. (1985) Modeling of Geothermal Systems. SPE 1985 California Regional Meeting, Bakersfield, California, 27-29 March 1985. University of California: Lawrence Berkeley National Laboratory. https://escholarship.org/uc/item/8k00x73q
[9]	Ganguly, S. and Kumar, M.M. (2012) Geothermal Reservoirs—A Brief Review. Journal of the Geological Society of India, 79, 589-602. https://doi.org/10.1007/s12594-012-0098-8
[10]	Pruess, K. (2002) Mathematical Modeling of Fluid Flow and Heat Transfer in Geothermal Systems: An Introduction in Five Lectures. Orkustofnun, Iceland. https://orkustofnun.is/gogn/unu-gtp-report/UNU-GTP-2002-03.pdf
[11]	Dobson, P., Asanuma, H., Huenges, E., Poletto, F., Reinsch, T. and Sanjuan, B. (2017) Supercritical Geothermal Systems—A Review of Past Studies and Ongoing Research Activities. 42nd Workshop on Geothermal Reservoir Engineering, Stanford, 13-15 February 2017.
[12]	Meacci, L., Farina, A. and Primicerio, M. (2019) Fluid Mechanics—A Free Boundary Model for the Evolution of a Geothermal System. Rendiconti Lincei-Matematica e Applicazioni, 30, 125-137. https://doi.org/10.4171/RLM/838
[13]	Rajagopal, K.R. and Tao, L. (1995) Mechanics of Mixtures, Vol. 35. World Scientific, Hackensack, NJ.
[14]	Bear, J. and Verruijt, A. (1987) Modeling of Groundwater Flow and Pollution. In: Theory and Applications of Transport in Porous Media, Reidel, Dordrecht. https://doi.org/10.1007/978-94-009-3379-8
[15]	Coussy, O. (1995) Mechanics of Porous Continua. Wiley, Hoboken.
[16]	Faust, C.R. and Mercer, J.W. (1977) Theoretical Analysis of Fluid Flow and Energy Transport in Hydrothermal Systems. Vol. 77, US Department of the Interior, Geological Survey. https://doi.org/10.3133/ofr7760
[17]	Faust, C.R. and Mercer, J.W. (1979) Geothermal Reservoir Simulation: 1. Mathematical Models for Liquid- and Vapor Dominated Hydrothermal Systems. Water Resources Research, 15, 23-30. https://doi.org/10.1029/WR015i001p00023
[18]	Darcy, H.P.G. (1856) Les Fontaines publiques de la ville de Dijon. Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d’eau, etc. V. Dalamont.
[19]	Farina, A., Bodin, J., Clopeau, T., Fasano, A., Meacci, L. and Mikelić, A. (2014) Isothermal Water Flows in Low Porosity Porous Media in Presence of Vapor-Liquid Phase Change. Nonlinear Analysis: Real World Applications, 15, 306-325. https://doi.org/10.1016/j.nonrwa.2011.11.021
[20]	Zemansky, M.W. (1968) Heat and Thermodynamics: An Intermediate Textbook. McGraw-Hill Book Company, New York.
[21]	Schmidt, E. (1968) VDI-Wasserdampftafeln/VDI Steam Tables/Tables VDI de la vapeurd’eau/Tablas VDI de vapor de agua. Springer-Verlag Berlin Heidelberg. https://doi.org/10.1007/978-3-662-24776-1
[22]	Melrose, J.C. (1974) Role of Capillary Forces in Detennining Microscopic Displacement Efficiency for Oil Recovery by Waterflooding. Journal of Canadian Petroleum Technology, 13, PETSOC-74-04-05. https://doi.org/10.2118/74-04-05
[23]	Vavra, C.L., Kaldi, J.G. and Sneider, R.M. (1992) Geological Applications of Capillary Pressure: A Review. AAPG Bulletin, 76, 840-850. https://doi.org/10.1306/BDFF88F8-1718-11D7-8645000102C1865D
[24]	Pruess, K., Oldenburg, C.M. and Moridis, G. (1999) Tough2 User’s Guide Version 2. Technical Report. Lawrence Berkeley National Lab. (LBNL), Berkeley, CA. https://doi.org/10.2172/751729
[25]	Pellicer, J., Garcia-Morales, V. and Hernandez, M.J. (2000) On the Demonstration of the Young-Laplace Equation in Introductory Physics Courses. Physics Education, 35, 126. https://doi.org/10.1088/0031-9120/35/2/309
[26]	Hsieh, C.-H. (1980) Vapor Pressure Lowering in Porous Media. PhD Thesis, Stanford University, Palo Alto.
[27]	Meacci, L., Farina, A., Rosso, F., Borsi, I., Ceseri, M. and Speranza, A. (2009) A Simplified Model for the Evolution of a Geothermal Field. COMSOL Conference, Milan. https://cn.comsol.com/paper/a-semplified-model-for-the-evolution-of-a-geothermal-field-7158