Thermodynamic Modeling of Mutual Solubilities in Gas-Laden Brines Systems Containing CO₂, CH₄, N₂, O₂, H₂, H₂O, NaCl, CaCl₂, and KCl: Application to Degassing in Geothermal Processes

Abstract

With the growing interest in geothermal energy as a renewable and sustainable energy source, nowadays engineers and researchers are facing technological and environmental challenges during geothermal wells’ operation or energy recovery improvement by optimizing surface installations. One of the major problems encountered is the degassing of geothermal brines which are often loaded with dissolved gases, resulting in technical problems (scale formation, corrosion, reduced process efficiency, etc.) and environmental problems through the possible emission of greenhouse gases (CO₂, CH₄ and water vapor) into the atmosphere. In this work, a method to predict, from readily available information such as temperature and GLR, the bubble point pressure of geothermal fluids as well as the GHG emission rate depending on the surface conditions is presented. This method is based on an extended version of the Soreide and Whitson model with new parameters optimized on the solubility data of several gases (CO₂, CH₄, N₂, O₂ and H₂) in brine (NaCl + CaCl₂ + KCl). The developed approach has been successfully used for the prediction of water content of different gases and their solubilities in different types of brines over a wide temperature and pressure range, and has been applied for the prediction of bubble point pressure and GHG emissions by comparing the results with available industrial data of geothermal power plants including the Upper Rhine Graben sites.

1. Introduction

The development of renewable energy is one of the key strategies for contributing to the decrease of greenhouse gas (GHG) emissions [1]. Geothermal energy, used for heat and electric power generation [2], appears to be a sustainable option to facilitate this energy transition. Depending on the depth of the reservoir containing the geothermal fluid, two applications are distinguished: shallow geothermal and deep geothermal [3,4]. In some deep geothermal reservoirs, ground permeability is too low to properly extract the fluid from its well. Thus, the geothermal resource heat is retrieved with an Enhanced Geothermal System (EGS) [3], which boosts rock porosity to extract the proper amount of fluid for heat recovery [5].

Currently, EGS projects are ongoing in Australia, France, Germany, the United Kingdom and the United States. However, their development stages are different from one country to another. While some of them are at plant construction and well drilling phases, such as those in Australia and the United Kingdom, others are commercial-scale plants, such as those in French territory [5,6,7].

France is a leading country in EGS development, notably thanks to the Soultz-sous-Fôrets project, which started in 1987 [5]. The Soultz-sous-Fôrets site, at the Upper Rhine Graben, produces electricity (12 GWh/year) from geothermal fluids [8], and it was the first commercial-scale EGS power plant in the world [5]. Deep geothermal on the Upper Rhine Graben, with the EGS technology, is also carried out at Rittershoffen, where heat is produced (180 GWh/year) for industrial needs [6,9]. EGS projects for heat and electricity production are currently under development on the region, at Illkirch-Graffenstaden and Eckbolsheim [6,8].

Usually, geothermal fluids are pumped from their wells to the surface, where they pass through heat exchangers to recover their energy. Nevertheless, geothermal fluids are often hot brines containing dissolved gases, and when they are sent to the surface, their pressure decreases and the dissolved gases evolve from their liquid phase [2]. The pressure under which this degassing occurs is the bubble point pressure [10,11].

On the one hand, the presence of a gas phase in the geothermal fluid sent to the surface reduces the process efficiency. Firstly, heat transfer efficiency in the heat exchangers declines if the fluid entering these devices is a gas-liquid mixture rather than a single-phase liquid. Moreover, depending on the geothermal brine composition, degassing leads to changes on its pH. Then, salt scaling (e.g., calcite-CaCO₃) could arise in the heat exchangers [9] or in the production well [11,12], impairing the process efficiency. On the other hand, carbon dioxide (CO₂) and methane (CH₄), both GHGs, are often found in geothermal fluids. Thus, their degassing can also cause undesired GHG emissions [13,14,15], when the geothermal resource is at the surface.

Hence, it is necessary to determine the bubble point depth to decide where a pump should be placed in the well to increase the fluid pressure [10] and minimize the amount of gas sent to the surface. The sites of Rittershoffen and Soultz-sous-Fôrets are pumped-well geothermal plants [6,9]. In Rittershoffen, the production well contains a pump set at a depth of 460 m [6], whereas the Soultz-sous-Fôrets production well presents a pump immersed at 350 m [9]. The bubble point pressure of these sites geothermal fluids is estimated to be about 2.0 MPa [9].

In some cases, despite the efforts to avoid its degassing, the geothermal resource arrives flashed at the wellhead of plants such as those in Rittershoffen and Soultz-sous-Fôrets. Thus, a vapor-liquid separator ensures that only the liquid fraction of the geothermal fluid is sent to the heat exchangers, while the gas phase, containing the GHGs, is vented to the atmosphere.

Most of the published data on gas emissions from geothermal sites focus on carbon dioxide emissions, which may be interesting information, as this gas is the most abundant in geothermal fluids [13]. These data concern only geothermal power plants, whose average carbon dioxide emission rate is 122 gCO₂·kWh⁻¹ [16]. However, this average emission rate is lower than those typically observed in fossil fuel power plants: 450 gCO₂·kWh⁻¹ in natural gas plants; 675 gCO₂·kWh⁻¹ in oil plants; and 900 gCO₂·kWh⁻¹ in coal plants [13].

Moreover, the configuration of the plant can influence its GHG emissions. The EGS sites in the Upper Rhine Graben employ process configurations where the geothermal fluid extracted from the well is reinjected in the ground [5,6]. Gas emissions would occur only at the plant liquid-vapor separator. These type of sites should emit low GHG quantities [6,13,17]. However, in some extreme cases, such as those of geothermal power plants in Turkey, where the geothermal fluid is also reinjected in its reservoir, carbon dioxide emissions (1050 gCO₂·kWh⁻¹) [18] can exceed those of some coal power plants [19]. This stems from the fact that the geothermal reservoirs exploited in this country are rich in carbonate-dominated rocks [20].

To the best of our knowledge, there is no available data on GHG emissions arising from the degassing of geothermal fluids at the Upper Rhine Graben sites during the process of energy production. However, as this environmentally harmful phenomenon can actually occur, it would be interesting to assess it at these geothermal sites.

To estimate GHG emissions and optimize the process, knowing the breakout pressure of geothermal fluid is essential. Thermodynamic characterization of the liquid-vapor equilibrium of these fluids, which depends on downhole composition data, is necessary to predict their boiling points. Each formation brine is unique and its composition is intimately linked to the geochemistry of the place where it is recovered on the planet [21]. This huge diversity of fluid composition participates in the complexity of this kind of study. Moreover, it could be difficult to determine the brine’s composition in a reservoir due to technical difficulties in setting up downhole sampling, related to temperature and pressure conditions. Therefore, although some articles describe techniques to perform this type of sampling [22,23,24,25], downhole compositions are usually obtained with thermodynamics models from physical and chemical data related to the fluid recovered at wellhead [12,26,27,28,29]. Among these data, one of the main parameters used to deduce downhole composition is the Gas-Liquid Ratio (GLR), which can be obtained by flashing a sample of the geothermal fluid at atmospheric pressure and wellhead temperature. Then, the gas and liquid volumes issued from this flash operation are measured, and the GLR can be expressed in gas volume units at standard or normal conditions per liquid volume units at flash conditions [30,31]. Finally, the gas phase composition is analyzed by chromatography [32], and these analysis results are associated to GLR data to ensure a further characterization of the geothermal fluid composition.

As the bubble point pressure is used to decide on the position of a submerged pump on the production well, its miscalculation leads to problems for geothermal plant design. Thus, any tool predicting geothermal fluid bubble point pressure should be accurate enough to avoid these complications. In order to ensure the calculations’ accuracy, the model parameters must be fitted from experimental data obtained for thermodynamic systems representatives of the natural geothermal resources of interest. Experimental thermodynamic characterization of gas-brine systems is therefore essential to obtain the breakout pressure. In case of the Upper Rhine Graben, the geothermal brine composition is well-known [32]. This fluid is a very salty formation water (TDS: Total Dissolved Species around 90 g·L⁻¹) which is registered as Na-Ca-Cl type and also contains dissolved gas. The off-gas obtained from the natural brines is composed of 86.2% carbon dioxide, with smaller quantities of nitrogen (N₂, 9.85%) and methane (2.28%). In the literature, some studies have obtained experimental solubility data of gas mixtures in pure water: CO₂ + CH₄ + H₂O [33,34,35,36] or CO₂ + N₂ + H₂O [37,38]. Nevertheless, the studied gas mixtures were not in accordance with the Upper Rhine Graben representative gas mixture. In addition, these investigations were carried out in a salt-free aqueous solution. By considering only the main gas, carbon dioxide, a recent work on the CO₂ + H₂O + Na + Ca + K + Cl system has shown that the minor species of the brine did not impact the carbon dioxide solubility [39]. The Upper Rhine Graben water could therefore be synthetized by water containing only sodium chloride (NaCl) and calcium chloride (CaCl₂).

From this information, the aim of this paper is to use a simple modelling approach to determine the bubble point pressure and GHG emissions in the Upper Rhine Graben and other geothermal sites worldwide. To perform this modeling work, the well-known Soreide and Whitson model with new optimized parameters will be used to estimate the downhole composition from the GLR, temperature, salinity (NaCl, CaCl₂ and KCl) and dry composition of the degassed fluid. A comparison with available industrial data of geothermal power plants will also be made.

2. Thermodynamic Modeling

The idea of the approach proposed in this section is to link the information easily available to geothermal operators, such as GLR, temperature and salinity, with the thermodynamics of phase equilibria in order to estimate equilibrium properties that are very important for the design, operation, and monitoring of geothermal wells. These properties mainly concern the prediction of the bubble-point pressure of gas-laden brines and the estimation of greenhouse gas emissions, and more specifically the amount of carbon dioxide and methane emitted in the process.

2.1. From GLR to the Gas-Laden Brine Composition

To calculate the bubble-point pressure, it is necessary to know the composition of all molecular species dissolved in the geothermal fluid. To do so, two types of approaches are possible: (1) an iterative approach to find the liquid composition by minimizing an objective function based on the differences between the calculated properties and the input information (GLR and gas composition, for instance). However, this can lead to optimization problems by converging on the local minimum instead of the global minimum. Hence, it requires a robust global optimization algorithm; (2) a simpler non-iterative approach based on a material balance by including the phases and compounds existing in the system. This work is based on the latter approach.

First, to link the GLR to the gas composition in geothermal fluids (brines), a material balance must be established on a flash system at a fixed temperature and pressure. Under certain conditions of temperature and pressure, a mixture of overall (monophasic) composition $z_i$ and molar flow rate $F$ may split into several phases. This is called a flash or a phase split equilibrium. The most common case is the two-phase liquid-vapor separation (Figure 1), where a liquid phase of flow rate L and composition $x_i$ and a vapor phase of flow rate V and composition $y_i$ are formed. The different compounds $i$ are distributed by number of moles in both phases according to their volatilities.

The global and partial material balances are written as follows:

$$F=V+L$$

(1)

$$F{z}_i=V{y}_i+L{x}_i=V{y}_i+\left(F-V\right)x_i$$

(2)

By defining the molar Gas-to-Liquid Ratio ${\left(GLR\right)}_{mol}=\frac{V\left(T,P\right)}{L\left(T,P\right)}$ and the partition coefficient $K_i=\frac{y_i}{x_i}$, we obtain from Equation (1):

$$z_i=\frac{y_i\left({\left(GLR\right)}_{mol}+\frac{1}{K_i}\right)}{{\left(GLR\right)}_{mol}+1}$$

(3)

From Equation (3), we can easily calculate the composition $z_i$ in the single-phase geothermal brine if we have the vapor composition $y_i$, the GLR and the $K_i$ partition coefficient. In Equation (3), the composition $y_i$ is the wet steam composition. Since in most cases one only knows the dry steam composition, a calculation must be performed to obtain the wet steam composition by including the water content $y_w$ estimated by Equation (6) (or by an equation of state, see below) and by renormalizing the composition on all molecular species (water and gas). In Equation (4), $f_i^L$ and $f_i^V$ are the fugacities in the liquid (L) and vapor (V) phases, $a_w$ is the activity of water in brine calculated by an activity coefficient model (e.g., Pitzer model) and $P_w^{sat}$ is the water vapor pressure.

$$f_i^L\left(T,P,n_i^L\right)=f_i^V\left(T,P,n_i^V\right)$$

(4)

The determination of the GLR is performed at low pressure (typically atmospheric), and under these conditions, with the Poynting factor and the fugacity coefficients very close to 1, one can consider:

$$f_{w\left(=water\right)}^L\left(T,P,n_w^L\right)=a_w\times P_w^{sat}\mathrm{and}{f}_w^V\left(T,P,n_w^V\right)=y_w\times P$$

(5)

therefore:

$$y_w=a_w\times \frac{P_w^{sat}\left(T\right)}{P}$$

(6)

Concerning the GLR, most of the data found in the literature are in the form of a volume ratio (Nm³/m³), whereas the GLR in the above material balance equations is molar (mol/mol), which therefore requires the knowledge of the densities of the liquid and the vapor phases. The following Equation (7) can be used to convert the GLR from volume ratio to molar ratio.

$${\left(GLR\right)}_{mol}={\left(GLR\right)}_{vol}\frac{{\rho }_{Vap}\times M{w}_{Liq}}{{\rho }_{Liq}\times M{w}_{Vap}}$$

(7)

In this work, the density of the liquid phase (brine) ${\rho }_{Liq}$ was calculated by the model proposed by Al ghafri et al., 2012 [40] and that of the gas phase ${\rho }_{Vap}$ by the multiparametric equation of state GERG-2008 (Kunz and Wagner 2012 [41]). Concerning the GERG-2008 model, its Fortran code for calculating the thermophysical properties of gas mixtures is freely released by NIST on github [42].

Regarding the partition coefficient $K_i$, correlations such as those of Giggenbach 1980 [43] and D’Amore and Truesdell 1988 [44] can be used, but with prudence, since they only consider the temperature and not the effect of pressure and salinity. A better option is the one we use in this work which is to calculate the partition coefficient by the equation of state that will be used for the calculation of phase equilibrium and that will be presented below.

2.2. Equation of State (EoS)

Because of their simplicity and versatility, the most widely used thermodynamic models for determining phase diagrams are the cubic equations of state. However, among their weaknesses is the poor representation of the phase equilibria of mixtures containing water, especially the Type III phase diagrams according to the Van Konynenburg and Scott classification [45]. This very particular molecule (water) has a molecular structure and thermophysical properties very different from those of hydrocarbons, which completely changes the phase equilibria, and even more so in the case of the existence of dissolved salts in water. The need for a good representation of the phase equilibria of mixtures of gases with water or brine therefore requires modifications to the cubic equations of state, either by modifying some parameters (binary interaction coefficient Kij, cohesion parameter a, etc.) as in the Soreide and Whitson model [46] or by adding other terms to the equation of state to take into account the presence of water or saline water (CPA type models [47,48]).

The Soreide and Whitson (SW) EoS [46] was selected in this work since it is widely used for gas-brine systems and is implemented in several process modeling and reservoir simulator tools (Simulis Thermodynamics [49], MUFITS [50], Eclipse 300 [51], IHRRS [52], etc.). It is an update of the well-known Peng Robinson equation of state (Equations (8)–(13) by adding a NaCl salinity $m_{NaCl}$ (mol/kgw) dependence to some model parameters. The alpha function of water (Equation (14)) and the binary gas/water interaction parameters specific to each phase (aqueous $k_{ij}^{AQ}$ and non-aqueous $k_{ij}^{NA}$, see Equations (15)–(18)) are concerned. It makes the model “asymmetric” and adds more flexibility for the simultaneous representation of liquid and vapor phases by different parameters.

$$P\left(T,v\right)=\frac{RT}{v-b}-\frac{a\left(T\right)}{v\left(v+b\right)+b\left(v-b\right)}$$

(8)

where P is the pressure, T the temperature, R the ideal gas constant and $v$ the molar volume. The energy parameter $a$ and the molar co-volume $b$ for a pure compound i are given by:

$$a_i\left(T\right)=0.457235529\frac{R^2{T}_{c,i}^2}{P_{c,i}}\times \alpha \left(T\right)$$

(9)

$$\alpha \left(T\right)={\left[1+\left(0.37464+1.54226{\omega }_i-0.26992{{\omega }_i}^2\right)\left(1-\sqrt{\frac{T}{T_{c,i}}}\right)\right]}^2$$

(10)

$$a_m\left(T\right)={\displaystyle \sum }_{i=1}^{n_c}{\displaystyle \sum }_{i=1}^{n_c}{x}_i{x}_j\sqrt{a_i\left(T\right)a_j\left(T\right)}\left(1-k_{ij}\right)$$

(11)

$$b_i=0.0777960739\frac{R{T}_{c,i}}{P_{c,i}}$$

(12)

$$b_m={\displaystyle \sum }_i{x}_i{b}_i$$

(13)

where $P_{c,i}$, $T_{c,i}$ and ${\omega }_i$ are respectively the critical pressure, the critical temperature and the acentric factor of compound i.

$${\alpha }_w\left(T,m_s\right)={\left[1+0.4530\left[1-\frac{T}{T_{c,w}}\left(1-0.0103m_{NaCl}^{1.1}\right)\right]+0.0034\left({\left(\frac{T}{T_{c,w}}\right)}^{-3}-1\right)\right]}^2$$

(14)

$$a_m^{AQ}\left(T\right)={\displaystyle \sum }_{i=1}^{n_c}{\displaystyle \sum }_{i=1}^{n_c}{x}_i{x}_j\sqrt{a_i\left(T\right)a_j\left(T\right)}\left(1-k_{ij}^{AQ}\right)$$

(15)

$$a_m^{NA}\left(T\right)={\displaystyle \sum }_{i=1}^{n_c}{\displaystyle \sum }_{i=1}^{n_c}{y}_i{y}_j\sqrt{a_i\left(T\right)a_j\left(T\right)}\left(1-k_{ij}^{NA}\right)$$

(16)

The binary interaction parameters $k_{ij}^{AQ}$ and $k_{ij}^{NA}$ are respectively adjusted on the gas solubility data and water content data of gas + water + NaCl systems. The $k_{ij}^{NA}$ parameter is generally constant or slightly dependent on temperature while the $k_{ij}^{AQ}$ paramater depends on temperature and NaCl molality $m_{NaCl}$. The binary gas/gas interaction parameters are those of the Peng-Robinson model and are widely available in the literature (Whitson and Brule 2000 [53], Lasala 2016 [54]. The gas/water parameters are listed in the “Results and Discussion” section.

2.3. Calculation Procedure

The procedure for calculating the bubble point pressure from the GLR, temperature and dry gas composition using an equation of state is shown in Figure 2. In the first step the dry gas composition is converted to the “true” composition (including water) and the GLR is converted to a molar ratio. Then, to estimate the composition of the geothermal fluid, the partition coefficients $K_i$ are generated by a classical PT-Flash calculation (by the Rachford-Rice method) [55,56,57] using the same temperature and pressure used to determine the GLR. However, the partition coefficients are very weakly dependent on the global composition z, so values of $z_i$ estimated close to those of the real case or simply values equal to 1/nc (relative to the number of compounds nc) can be used. Once the partition coefficients and the composition in the vapor phase $y_i$ are determined, the composition of the downhole fluid can be determined by Equation (3). Hence, the calculation of the bubble point could be started, knowing the required set of input data (temperature and liquid composition $x_i=z_i$). It consists of an iterative search for the bubble point pressure, making sure each time that the vapor composition is stable, the sum of the gas molar fraction is equal to 1. The iterative part can be performed by the Newton method [56] or by substitutions [58] (presented in Figure 2), which is simpler and sufficient. With this method, the bubble point pressure can be estimated at any point in the reservoir by knowing the temperature gradient, which will allow the determination of at which depth level the bubbling takes place. With this method, bubble point pressure can be estimated at any point at the surface or in the reservoir by knowing the temperature or its gradient.

The GHG emission rate is directly related to the vaporized fraction $\beta =\frac{V}{F}$ which can be determined by a PT-Flash calculation at known temperature and pressure. Therefore, the GHG emission rate depends on the temperature and pressure conditions at the surface or in the tank and on the difference between the actual pressure and the bubble point pressure. In other words, the lower the pressure is below the bubble point pressure, the higher the vaporized fraction and the quantity of gas emitted. Generally, the flow rates F of geothermal fluids are given in terms of volume (e.g., m³/h), hence there is a need to know its density in order to be able to calculate the amount of gas emitted. The density of the liquid and vapor phases are given directly by the equation of state. This represents an advantage of using the phi-phi approach (use of an equation of state to describe both liquid and vapor phases) instead of the gamma-phi approach. The volumetric flow rate must be converted to a molar flow rate F using the molar density given by the equation of state at T, P and $z_i$ (Equation (3)), then multiplied by the vaporized fraction $\beta$ to obtain the degassed quantity.

3. Results and Discussion

3.1. Mutual Solubilities of Gas-Water and Gas-Brine (Single or Mixed Salts) Systems

a.

Water content in gas rich-phase

The asymmetry of the binary interaction parameters used in the Soreide and Whitson EoS permits one to overcome the lack of “flexibility” of the symmetrical (phi-phi) approach concerning the distinct representation of the liquid and vapor phases by using specific parameters for each phase. In addition to the improvement of the calculation of the liquid phase, the other advantage of this adaptation is a very good representation of the vapor phase in particular the water content in the gas-rich phase. In Figure 3, experimental literature data of water content in several gases of interest (CO₂, CH₄, N₂, O₂, and H₂) are compared with the SW model predictions using the optimized interaction parameters proposed in this work (see Section 3.1 d). As shown in the figure, the model accurately predicts the water content in these gases over a wide range of pressure and temperature, covering reservoir and surface conditions. The average absolute deviations between the predictions and the experimental literature data vary between 2 and 4% for the mentioned gases.

b.

CO₂ solubility in different brines (NaCl, CaCl₂, or KCl)

CO₂ is generally the major species among the gases dissolved in geothermal waters, and this is justified by its relatively high solubility compared to other gases. Therefore, in this section we focus on the characterization of its solubility under thermodynamic conditions in geothermal reservoirs. As previously mentioned, the model takes into account the salinity by including the molality of NaCl in the model parameters (alpha function of water and kij in the aqueous phase). Since sodium and chloride are the main species in natural brine, NaCl is the only salt considered in the original version of SW model. Considering NaCl as the only salt is a good first assumption to simplify the physico-chemical calculations. However, in some geological reservoirs, other salts may be dissolved in the brine such as CaCl₂, KCl, etc. These salts can exist in significant quantities and exhibit salting-out effects more or less important than NaCl for the same concentrations. Hence, it is important to consider these species for a more “realistic” representation of formation water. In this work, the binary interaction parameters used to model the CO₂ + H₂O + NaCl system are those proposed by Chabab et al., 2019 [47] which improve very significantly the representation of the CO₂ solubility in particular in very concentrated and saturated NaCl-brines compared to the original parameters.

Regarding the modeling of CO₂ + H₂O + CaCl₂ and CO₂ + H₂O + KCl systems, the following assumptions were considered to extend the SW model to these salt solutions:

-

The original alpha function of water was kept as a function of the NaCl molality and was not modified for the other salts for the following reasons: (1) To propose a salt equivalence technique considering NaCl as the reference salt (see hereafter); (2) the alpha function proposed by SW is obtained by fitting on vapor pressure data of the H₂O + NaCl system and the difference between the vapor pressures of the brines of different salts does not have significant effects on the representation of the phase equilibria especially at the thermodynamic conditions of interest (high pressure).

-

The coefficients of the CO₂-H₂O binary interaction parameter expression in the aqueous phase $k_{C{O}_2,H_{2}O}^{AQ}$ (Equation (17)) were readjusted on the CO₂ solubility data in the mentioned brines (H₂O + CaCl₂ and H₂O + KCl). The readjusted Kij coefficients are listed in

Table 1. Coefficients of the $k_{C{O}_2,H_{2}O}^{AQ}$ correlation (Equation (17)).

	NaCl	CaCl2	KCl
a	0.43575155	0.478902106	0.44269716
b	−5.766906744 × 10−2	−7.35093213 × 10−2	−5.746192569 × 10−2
c	8.26464849 × 10−3	2.05379759 × 10−2	6.655117 × 10−3
d	1.29539193 × 10−3	2.04457138 × 10−3	−1.318885435 × 10−3
e	−1.6698848 × 10−3	−3.22673398 × 10−3	6.161936936 × 10−4
f	−0.47866096	−0.50876635	−0.488832169

.

-

The CO₂-H₂O binary interaction parameter in the non-aqueous phase $k_{C{O}_2,H_{2}O}^{NA}$ (Equation (18)) was not modified for the other salts because: (1) there are no (or few) water content data of CO₂ in equilibrium with CaCl₂ and KCl brines; (2) the salting-out effect on water content is much less important than on gas solubility. Thus, for the salts CaCl₂ and KCl, the kij used are: Equation (17) for the aqueous phase with the coefficients listed in Table 1 and Equation (18) for the non-aqueous phase.

$$k_{C{O}_2,H_{2}O}^{AQ}=\left(\frac{T}{T_{c,C{O}_2}}\right)\left[a+b\left(\frac{T}{T_{c,C{O}_2}}\right)+c\left(\frac{T}{T_{c,C{O}_2}}\right)m_{salt}\right]+m_{salt}^2\left[d+e\left(\frac{T}{T_{c,C{O}_2}}\right)\right]+f$$

(17)

$$k_{C{O}_2,H_{2}O}^{NA}=0.68208385571\times {10}^{-3}T-2.066623464504\times {10}^{-2}$$

(18)

In Figure 4, the experimental literature data of CO₂ solubilities in the different salts (NaCl, CaCl₂ and KCl) are compared to the calculations generated by the SW model using the proposed optimized parameters. As shown in the figure, the model captures well the salting-out effect of the different salts on the solubility of CO₂ under the different conditions of temperature, pressure and molality of salts. The model reproduces the solubility in NaCl and CaCl₂ solutions perfectly, but a little less well in KCl solution. This could be due to several factors (data consistency, unsuitable kij expression, optimization problem (local minimum), etc.), but is probably due to the consistency of the literature data, knowing that there is a mismatch between the data of Hou et al., 2013 [65] and those of Zhao et al., 2015 [66], and that the model reproduces very well the data of Hou et al. (see

Table 2. Deviation between experimental data and SW model prediction of CO₂ solubility in brine.

	Data Reference	AAD% (SW Model) a
CO2 + H2O + CaCl2	Zhao et al., 2015 [66]	2.24
	Messabeb et al., 2017 [74]	4.13
	Lara Cruz et al., 2021 [75]	3.36
CO2 + H2O + KCl	Zhao et al., 2015 [66]	6.97
	Hou et al., 2013 [65]	2.97

^a: using the optimized parameters proposed in this work (Table 1).

).

In order to make the SW model capable of predicting the solubility of CO₂ in real brines (mixed salt systems), correlations (Duhring lines) between the molality of the salts (CaCl₂ and KCl) and their NaCl equivalents were developed in this work. This was undertaken by an iterative (Newton-Raphson) approach of finding equal CO₂ solubilities in the different salts considering NaCl as reference salt (Markham and Kobe 1941 [76], Ratnakar et al., 2020 [77]). The data resulting from the search for equal solubilities have been correlated in the following equations (Equations (19) and (20)) to make the conversion from salt molality to NaCl equivalent molality straightforward. The correspondence between salts and NaCl is independent of pressure but depends on temperature. Figure 5 shows the Duhring lines using the salt→NaCl equivalence correlations developed to estimate the NaCl equivalent of CaCl₂ and KCl at equal CO₂ solubilities for two isotherms (323 and 373 K).

$$\begin{array}{c}{m}_{NaCl}=m_{CaC{l}_2}\left(6.4824099-3.808876\times {10}^{-3}\times T\right)-923.14\times \frac{m_{CaC{l}_2}}{T}\\ -\frac{5.0858999\times m_{CaC{l}_2}}{11\times m_{CaC{l}_2}+0.0262706\times T-6.739781}\end{array}$$

(19)

$$\begin{array}{ll}{m}_{NaCl}& =0.0095690161\times T\times m_{KCl}+0.0004886086\times T\times {m_{KCl}}^2-1.1835076\times m_{KCl}\\ & -0.16724994\times {m_{KCl}}^2-1.170219186\times {10}^{-5}\times m_{KCl}\times T^2\\ & -21.4278766\times {0.00467461197}^{m_{KCl}}\times {m_{KCl}}^3\end{array}$$

(20)

In Figure 6, the solubility of CO₂ in CaCl₂ and KCl solutions calculated by the SW model using the Kij corresponding to each salt (Table 1) were compared to the calculations using their NaCl equivalent (with the help of Equations (19) and (20)) and the NaCl coefficients for the binary interaction parameter (Table 1). The results show a perfect adequacy between the calculations by the “apparent” molality and those by the “real” Kij corresponding to each salt.

c.

CO₂ solubility in mixed-salt brine (NaCl + CaCl₂ + H₂O)

Thanks to the salt→NaCl equivalence correlations proposed in this work, the SW model can be applied to the study of gas solubility in mixed-salt systems, which is more realistic than considering NaCl as the only dissolved salt in groundwater. The idea of this approach is to convert with the help of these correlations the molalities of the other salts (CaCl₂ and KCl) into equivalent NaCl molality, add them to the NaCl molality, and this total is considered in the calculations of the water alpha function and the binary interaction parameters. To test this approach, the prediction results for the solubility of CO₂ in a salt mixture (CaCl₂ + NaCl) using NaCl equivalence were compared with predictions without including this molality correction. This comparison is illustrated in Figure 7 taking as reference the experimental data recently measured by Poulain et al., 2019 [39] and Lara Cruz et al., 2021 [75], and this corresponds to five isotherms between 323 and 453 K at fixed salt molality (1.2 mol/kgw NaCl + 0.2 mol/kgw CaCl₂) and for pressures up to 40 MPa.

As shown in Figure 7, the inclusion of the salt equivalence correlation improves the prediction of CO₂ solubility in mixed-salt solution (NaCl + CaCl₂) and the enhancement becomes clearer and more important with increasing pressure (Figure 7a). Finally, these equivalence relations based on the salting-out effects of the different salts improves the prediction in multi-salt systems and especially when other salts exist in non-negligible quantities in the studied environments.

d.

CH₄, N₂, O₂, and H₂ solubility in water and NaCl brine

Besides the CO₂ presence, other gases such as CH₄, N₂, O₂ and H₂ are generally dissolved in geothermal brines in varying amounts and must be considered as they significantly change the bubble pressure of gas-laden brines. In their original paper, Soreide and Whitson have also proposed binary interaction parameters for CH₄ + H₂O and N₂ + H₂O. However, they have not proposed parameters for the H₂ + H₂O and O₂ + H₂O binaries. The parameters used for these two binaries are those recently proposed in the works of Chabab et al. [78,79]. On the other hand, the coefficients of the Kij correlations proposed by Soreide and Whitson for CH₄ + H₂O and N₂ + H₂O have been readjusted in this work. The expressions of the Kij and their coefficients are presented in

Table 3. Gas-H₂O interaction parameters in aqueous and non-aqueous phases for the SW EoS.

	Gas-H2O Interaction Parameters in Aqueous and Non-Aqueous Phases
CH4 $T_{c,C{H}_4}=190.56\mathrm{K}$ $P_{c,C{H}_4}=45.992\mathrm{bar}$ ${\omega }_{C{H}_4}=0.01142$	$\begin{array}{cc}\hfill k_{C{H}_4,H_{2}O}^{AQ}& =-1.625685\left(1+8.590105\times {10}^{-21}{m}_{NaCl}\right)\hfill \\ & +1.114873\frac{T}{T_{c,C{H}_4}}\left(1+1.812763\times {10}^{-3}{m}_{NaCl}\right)\hfill \\ & -0.169968{\left(\frac{T}{T_{c,C{H}_4}}\right)}^2\left(1-4.198569\times {10}^{-2}{m}_{NaCl}\right)\hfill \end{array}$ $k_{C{H}_4,H_{2}O}^{NA}=0.494435$
O2 $T_{c,O_2}=154.581\mathrm{K}$ $P_{c,O_2}=50.43\mathrm{bar}$ ${\omega }_{O_2}=0.0222$	$\begin{array}{cc}\hfill k_{O_2,H_{2}O}^{AQ}& =-1.167744\left(1+3.361921\times {10}^{-2}{m_{NaCl}}^{0.8}\right)\hfill \\ & +0.466607\frac{T}{T_{c,O_2}}\left(1+8.457306\times {10}^{-2}{m_{NaCl}}^{0.8}\right)\hfill \end{array}$ $k_{O_2,H_{2}O}^{NA}=0.581650$
N2 $T_{c,N_2}=126.19\mathrm{K}$ $P_{c,N_2}=33.958\mathrm{bar}$ ${\omega }_{N_2}=0.0372$	$\begin{array}{cc}\hfill k_{N_2,H_{2}O}^{AQ}& =-1.709096\left(1+1.792130\times {10}^{-2}{m_{NaCl}}^{0.8}\right)\hfill \\ & +0.450487\frac{T}{T_{c,N_2}}\left(1+0.066426{m_{NaCl}}^{0.8}\right)\hfill \end{array}$ $k_{N_2,H_{2}O}^{NA}=0.385438$
H2 $T_{c,H_2}=33.145\mathrm{K}$ $P_{c,H_2}=12.964\mathrm{bar}$ ${\omega }_{H_2}=-0.219$	$\begin{array}{cc}\hfill k_{H_2,H_{2}O}^{AQ}& =-2.34\left(1+3.88\times {10}^{-3}{m_{NaCl}}^{0.443}\right)+0.166\frac{T}{T_{c,H_2}}\left(1+0.049{m_{NaCl}}^{0.799}\right)\hfill \\ & +-12.69\exp \left(-0.474\frac{T}{T_{c,H_2}}\right)\hfill \end{array}$ $k_{H_2,H_{2}O}^{NA}=-0.3776+0.08385\frac{T}{T_{c,H_2}}$

. It should be noted that the calculations by the SW model are very sensitive to the critical data (critical temperature and pressure and acentric factor) used, so it is important to use the appropriate critical data provided with the binary interaction parameters.

In Figure 8, the prediction of the solubility of CH₄, N₂, O₂, and H₂ in water and NaCl brine is compared with the literature’s experimental data. As shown in this figure, the predictions are in very good agreement with the experimental data. Indeed, the use of the Søreide-Whitson model with the proposed optimized parameters allows the estimation of the solubility of these gases in brine over a wide range of temperature, pressure and salinity and with an average absolute deviation of about 3%.

The difference in the studied gas solubilities can be explained by several factors such as the size of the molecules [81]. Indeed, large molecules need more energy to form a cavity in the aqueous phase. Moreover, this effect is important at high pressure. As such, it can explain the strong increase in the solubility of H₂ at elevated pressure (see Figure 8), which can even exceed the solubility of other gases (O₂ and CH₄). The contribution of hydrogen bonds (solvation) to solubility is very important and explains the high dissolution of CO₂ (which is a Lewis acid) [82] compared to other gases. Another factor to explain the variation in solubility could be the polarizability, but it is less dominant (in water) compared to hydrogen bonding and coulombic interactions [83]. Finally, the notion of solubility is influenced by several parameters which can also have counterbalances of these different dependencies. Hence, the comprehension of the difference between gas solubility could be difficult to explain directly.

e.

Mixed-gas cosolubility: CO₂ + CH₄ + H₂O and CO₂ + N₂ + H₂O systems

For mixed gas cosolubility modeling, thanks to the mixing rules, the model is easily applicable to multicomponent mixtures when binary interaction parameters are available. In this work, we are interested in gas mixtures including mainly CO₂ (70–99%) with other gases in more or less important fractions such as CH₄ and N₂. Concerning the CO₂ + N₂ + H₂O ternary mixture, the only solubility data available in the literature (Liu et al., 2012 [84]) are only at low temperature (up to 318 K) which have been studied for other applications such as the transport of flue gases. In Figure 9, the cosolubility data of the mixture of 90%CO₂ + 10%N₂ in water at 318 K are compared to the predictions of the SW model. In general, the predictions are good especially for N₂ but less good for CO₂ since the SW model underestimates the solubility of CO₂ at low temperatures (see Chabab et al., 2019 [47]). This is explained by the fact that the model parameters were fitted to data at higher temperatures and which are more representative of those found in geological reservoirs.

Concerning the CO₂ + CH₄ + H₂O ternary mixture, there are some experimental data available in the literature [33,34] and especially in the temperature and pressure range representative of deep geological formations. In Figure 10, the cosolubility data of 72%CO₂ + 28%CH₄ in water at 376 K are compared to the predictions of the SW model. The predictions are very good for both CO₂ and CH₄ over a very wide pressure range (up to 50 MPa). It should be noted that these results were obtained without any additional parameters and without the inclusion of these data in the adjustment of the binary interaction parameters. This illustrates the advantage of the “predictive” aspect of the equations of state as opposed to the empirical correlations/models.

Finally, knowing the temperature, pressure, dry gas composition and salt concentrations, the phase equilibria of the gas(s)/brine(s) systems can be predicted with good accuracy with the SW model using the interaction parameters and salt equivalence correlations proposed in this work. Moreover, any other types of phase equilibrium calculations are applicable to this model such as the search for the bubble point of gas-laden geothermal brines, which is the subject of the next section.

3.2. Application: Geothermal Fields Operations

In this part, the approach proposed in Section 2 is applied to characterize the degassing of geothermal fluids. The bubble pressure of the gas-laden brines and the GHG emission rate are estimated from information available to the geothermal operators such as the gas-to-liquid ratio (GLR), the temperature (T) in the well and the dry gas compositions of the degassed compounds.

3.2.1. Estimation of the Bubble-Point Pressure from the GLR

In order to validate the developed approach, this input information (GLR, T, salinity and gas composition) and their corresponding bubble-point pressures must be known. However, in the literature, these data are not widely available and generally incomplete. Two data sources [31,85] providing some of this necessary information have been considered for comparison and validation of the proposed predictive approach. The data from the Paris area reported by Ungemach [85] contain more information but the salinity of the brine samples is not given, so we made the calculations at different equivalent NaCl molalities. In

Table 4. Comparison of predicted bubble pressures (Pb) with those obtained from PVT data of downhole samples [31,85].

Data Source	Depth (m bgl)	CO2 (% mol)	CH4 (% mol)	N2 (% mol)	O2 (% mol)	GLR (m3/m3)	T Exp (K)	Pb Exp (MPa)	NaCl Molality (mol/kgw)	Pb Calc (MPa)	AD% Pb
Vandenberghe et al. (2001) [31]	1646	0.9367	0.0340	0.0271	0.0022	1.62	347.05	2.654	2.5	2.470	6.95
	1646	0.9225	0.0321	0.0443	0.0011	1.11	347.05	1.482	2.5	2.074	39.93
	1740	0.8311	0.0627	0.1015	0.0047	0.96	347.05	2.654	2.5	3.158	18.95
	1740	0.9243	0.0354	0.0323	0.0079	0.71	347.05	1.276	2.5	1.331	4.38
	1646	0.8795	0.0443	0.0751	0.0011	1.22	347.05	2.896	2.5	3.107	7.28
Ungemach, P. (2001) [85]	160	0.5241	0.2113	0.2646	-	0.23	335.35	0.790	0.5	0.737	6.77
	1818	0.5353	0.2006	0.2641	-	0.25	337.35	0.820	0.5	0.767	6.41
	1818	0.5245	0.1937	0.2818	-	0.26	337.35	0.840	0.5	0.820	2.41
	160	0.5307	0.1994	0.2699	-	0.28	333.35	0.830	0.5	0.890	7.24

and Figure 11 and Figure 12, the bubble pressure data reported from PVT analyses of samples taken from the geothermal wells are listed, illustrated and compared to those predicted by the model. Considering a NaCl molality of 0.5 mol/kgw, in Figure 11, the model reproduces well the Ungemach data with an AAD% of 5.7% for GLR between 0.23 and 0.28 and temperatures between 333 and 337 K knowing that the uncertainties of these data are not reported either. For samples taken at the same depth (see Table 4), in contrast to the Ungemach data (Paris) [85], the Vandenberghe et al. (North Belgium) [31] data appear to have very large uncertainty (poor repeatability) with large differences in gas composition, GLR and inhomogeneity of salt content in the analyzed brine samples (from the same site) and consequently wide variation in bubble pressure. For this reason, larger deviations between the model predictions and the Vandenberghe data are observed. However, the model captures qualitatively very well the variation of the bubble point as a function of GLR and gas composition and gives coherent values with three out of five points among these data (see Table 4 and Figure 12).

3.2.2. Estimation of Bubble-Point Pressure and GHG Emission Rate in the Upper Rhine Graben Geothermal Site

Based on recently published data (Mouchot et al., 2018 [86]) from a few geothermal wells operated in the Upper Rhine Graben, the presented approach was used to predict the bubble point pressure and GHG emission rates in these sites and more precisely in the production wells GPK2 (Soultz-sous-Forêts) and GRT2 (Rittershoffen). These data are listed in

Table 5. Prediction of bubble pressures (Pb) and GHG emission rates in the Upper Rhine Graben geothermal sites.

Geothermal Site	Gas Composition at Production Well a			GLR (Nm3/m3) at 0 °C and 1 atm a	T (K)	Pb (MPa) Exp at T a	NaCl Molality (mol/kgw) b	Pb Calc (MPa) c	P (MPa)	Geothermal Brine Flow Rate (kg/s) a	Degassed Amount (kg/h) c
	CO2 (% mol) a	CH4 (% mol) a	N2 (% mol) a
Soultz-sous-Forêts	0.91	0.02	0.07	1.03	423.15	$\approx 2.3$	1.426	2.296	≥2.296	30	0
	0.91	0.02	0.07	1.03	423.15	N.A.	1.426	2.296	2.100	30	7.7763
	0.91	0.02	0.07	1.03	423.15	N.A.	1.426	2.296	1.900	30	18.6904
	0.91	0.02	0.07	1.03	423.15	N.A.	1.426	2.296	1.5	30	62.6332
Rittershoffen	0.91	0.02	0.07	1.2	443.15	$\approx 2.5$	1.3	2.59	≥2.59	72.5	0
	0.91	0.02	0.07	1.2	443.15	N.A.	1.3	2.59	2.100	72.5	84.3599
	0.91	0.02	0.07	1.2	443.15	N.A.	1.3	2.59	1.900	72.5	154.0723
	0.91	0.02	0.07	1.2	443.15	N.A.	1.3	2.59	1.5	72.5	425.6754

^a Data from Mouchot et al. (2018) [86]. ^b Equivalent salinity in NaCl (mNaCl = 1.2 and mCaCl₂ = 0.2 mol/kgw) calculated by the correlation developed in this work (Equation (19)). ^c Predicted using the developed model.

along with the predicted bubble pressure values and flow rates of the degassed fluids at different temperatures and pressures which are further detailed in

Table 6. Prediction of GHG emission rates in the Upper Rhine Graben geothermal sites at different pressures and temperatures: estimated amounts and compositions of emitted gases.

Geothermal Site	T (K)	P (MPa)	$\mathit{\beta }=\frac{\mathit{V}}{\mathit{F}}$ (Molar)	Geother-mal Brine Flow Rate $\mathit{F}$ (kg/s) a	Geother-mal Brine Flow Rate $\mathit{F}$ (mol/s) b	${\mathit{y}}_{{\mathit{H}}_2\mathit{O}}$b	${\mathit{y}}_{\mathit{C}{\mathit{O}}_2}$b	${\mathit{y}}_{\mathit{C}{\mathit{H}}_4}$b	${\mathit{y}}_{{\mathit{N}}_2}$b	Degassed Amount $\mathit{V}$ (mol/s) b	Degassed Amount $\mathit{V}$ (kg/s) b	GHG Emitted (g-eqCO2/kWh) b
Soultz-sous-Forêts (12 GWhe/y)	423.15	2.10	$4.01\times {10}^{-5}$	30	1662.561	0.23081	0.45380	0.04583	0.26956	0.0666364	0.0021601	6.71
	423.15	1.90	$9.53\times {10}^{-5}$	30	1662.561	0.25370	0.48601	0.04174	0.21855	0.1585203	0.0051918	15.86
	423.15	1.50	$3.16\times {10}^{-4}$	30	1662.561	0.31739	0.53799	0.02783	0.11679	0.5254182	0.0173981	48.06
Rittershoffen (180 GWhth/y)	443.15	2.10	$1.89\times {10}^{-4}$	72.5	4017.522	0.38136	0.43862	0.02994	0.15008	0.7593836	0.0234333	4.16
	443.15	1.90	$3.47\times {10}^{-4}$	72.5	4017.522	0.41889	0.44965	0.02408	0.10738	1.3927135	0.04279786	7.18
	443.15	1.50	$1.00\times {10}^{-3}$	72.5	4017.522	0.52317	0.41666	0.01255	0.04762	4.0359467	0.11824318	16.51

^a Data from Mouchot et al. (2018) [86]. ^b Predicted using the developed model.

. According to Mouchot et al., at the geothermal fluid temperature of 150–170 °C, a pressure of 2.3–2.5 MPa is maintained at the surface to avoid degassing and related problems such as corrosion and scale formation. This pressure, determined just by feedback, represents a pressure around the bubble point to keep the geothermal fluid in a liquid state. At a NaCl molality of 1.3–1.4 mol/kgw (obtained by conversion using the equivalence correlation proposed in this work), the model predicts bubble pressures of 2.296 and 2.59 MPa for both sites, which are very close to those reported (2.3 and 2.5 MPa) and used as lower limits in surface facilities. As a result of this excellent bubble pressure prediction, the method can also be considered for estimating the GHG rates that can be emitted and which are directly related to the vaporized fraction.

The estimates of GHG emissions expressed in grams of CO₂ equivalent per kWh of energy produced are listed in Table 6. They have been determined from the energy production of these facilities and the quantities of degassed fluids, notably CO₂ and CH₄, considering that the global warming potential of CH₄ is 25 times higher than that of CO₂. As shown in Table 6, by decreasing the pressure and/or increasing the temperature, the degassing and the GHG emissions become important. Compared to the average of 122 gCO₂/kWh (see Introduction), the gas emissions are low, and may be zero if the pressure is maintained above the bubble pressure. Even at high depressurization, the emissions at these two sites remain low compared to the average emissions from geothermal sites and even lower than those from fossil fuel-based plants. This conclusion is in line with that of Lu et al. [5] and Pratiwi et al. [6] on the fact that the Rhine Graben sites are environmentally friendly and indeed offer low-carbon geothermal energy.

4. Conclusions

In the context of geothermal energy exploitation, the knowledge of the bubble pressure of geothermal fluids is very important to avoid technical and environmental problems as well as for the design, optimization and monitoring of surface facilities. In this work, a simple and efficient procedure for the prediction of phase equilibria (bubble point and flash calculation) of gas/brine systems during geothermal fluid exploitation is proposed. Being the main factors involved in the degassing problem, water content in the gas-rich phase and gas solubility have been successfully modeled in different brines at conditions covering widely the temperature, pressure and salinity range of geothermal applications. The model reproduces accurately the data reported by industrial operators and can therefore be used to estimate the equilibrium properties of any gas-laden brine from some readily available information such as temperature, GLR and dry gas composition. Finally, to fill the great lack of equilibrium data (solubility, GLR and bubble pressure) of gas/brine systems, an experimental pilot is under development, which will allow us to generate data that can be compared and used to improve the model presented in this work.