# Modelling reinjection of two-phase non-condensable gases and water in geothermal wells

Vlasios Leontidis, Pouriya H. Niknam, Ismail Durgut, Lorenzo Talluri, Giampaolo Manfrida, Daniele Fiaschi, Serhat Akin, Martin Gainville

► To cite this version:

Vlasios Leontidis, Pouriya H. Niknam, Ismail Durgut, Lorenzo Talluri, Giampaolo Manfrida, et al. Modelling reinjection of two-phase non-condensable gases and water in geothermal wells. Applied Thermal Engineering, 2023, 223, pp.120018. 10.1016/j.applthermaleng.2023.120018. hal-03980303

HAL Id: hal-03980303

https://ifp.hal.science/hal-03980303

Submitted on 9 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Nomenclature

Symbols

<table><tr><td>A</td><td>cross-section (m²)</td></tr><tr><td>d</td><td>internal diameter (m)</td></tr><tr><td>C</td><td>liquid phase concentration (mol/kg)</td></tr><tr><td>f</td><td>friction factor</td></tr><tr><td>g</td><td>gravitational acceleration (m/s²)</td></tr><tr><td>h</td><td>convective heat transfer coefficient (W·K⁻¹·m⁻²)</td></tr><tr><td>h̃</td><td>specific enthalpy (J/kg)</td></tr><tr><td>H</td><td>Henry's law constant (mol·m⁻³·Pa⁻¹)</td></tr><tr><td>m</td><td>mass flow rate (kg/s or t/h)</td></tr><tr><td>I.D.</td><td>internal diameter (mm)</td></tr><tr><td>IP</td><td>injectivity index (kg·s⁻¹·bar⁻¹)</td></tr><tr><td>MD</td><td>measured depth (m)</td></tr><tr><td>O.D.</td><td>external diameter (mm)</td></tr><tr><td>P</td><td>pressure (bar or Pa)</td></tr><tr><td>q</td><td>rate of heat transfer from the well fluid to the formation (W/m)</td></tr><tr><td>r</td><td>radius (m)</td></tr><tr><td>R</td><td>universal gas constant (m³·Pa·K⁻¹·mol⁻¹)</td></tr><tr><td>Rθ</td><td>thermal resistance (K·m/W)</td></tr><tr><td>Re</td><td>Reynolds number</td></tr><tr><td>t</td><td>time (s)</td></tr><tr><td>T</td><td>temperature (°C)</td></tr><tr><td>TVD</td><td>true vertical depth (m)</td></tr><tr><td>u</td><td>velocity (m/s)</td></tr><tr><td>U</td><td>overall heat transfer coefficient (W·m⁻²·K⁻¹)</td></tr><tr><td>x</td><td>mole fraction (%mol)</td></tr><tr><td>X</td><td>monomer fraction</td></tr><tr><td>y</td><td>compositional concentration in the gas phase (%wt)</td></tr></table>

<table><tr><td>z</td><td>distance in the vertical direction (m)</td></tr></table>

Greek symbols

<table><tr><td>α</td><td>void fraction</td></tr><tr><td>β</td><td>radial distribution function</td></tr><tr><td>Δ</td><td>difference</td></tr><tr><td>δ</td><td>error (%)</td></tr><tr><td>ε</td><td>wall roughness (m)</td></tr><tr><td>θ</td><td>inclination angle of the well</td></tr><tr><td>λ</td><td>thermal conductivity (W·m⁻¹·K⁻¹)</td></tr><tr><td>μ</td><td>viscosity (Pa·s)</td></tr><tr><td>ρ</td><td>density (kg/m³)</td></tr><tr><td>σ</td><td>surface tension (N/m)</td></tr><tr><td>ν</td><td>molar volume (m³/mol)</td></tr><tr><td>Φ</td><td>parameter</td></tr></table>

Subscripts

<table><tr><td>av</td><td>average</td></tr><tr><td>bh</td><td>bottom-hole</td></tr><tr><td>calc</td><td>calculated</td></tr><tr><td>G</td><td>gas</td></tr><tr><td>f</td><td>friction</td></tr><tr><td>k</td><td>phase (gas or liquid)</td></tr><tr><td>L</td><td>liquid</td></tr><tr><td>m</td><td>mixture or solid material</td></tr><tr><td>res</td><td>reservoir</td></tr><tr><td>s</td><td>superficial</td></tr><tr><td>TP</td><td>two-phase</td></tr><tr><td>surf</td><td>surface</td></tr><tr><td>tot</td><td>total</td></tr><tr><td>wh</td><td>wellhead</td></tr></table>

would anyway reach the surface as natural degassing [2]. Consequently, new methods for reducing the environmental impact of geothermal exploitation resources must be found and applied, especially in sites with high NCGs content on the resource fluid.

Reinjection of the NCGs back to the reservoir, together with the extracted liquid, is the first step to confine and dispose them (dissolution in deep geothermal resources and possible mineralization through suitable geochemistry interaction). It is also a key step to operate the geothermal system as a closed loop, maintaining the integrity of the resource and ensuring completely renewable energy compliance. Injection of fluids through a well to underground geological formations is a solution that has already been applied at an industrial scale [3]. Applications include the disposal of liquid by-products and industrial waste, displacing oil from porous rocks toward production wells (Enhanced Oil Recovery, EOR), supporting the reservoir pressure and stimulating reservoirs in geothermal applications. Moreover, underground gas injection technology is nowadays attractive for the seasonal storage of natural gas and hydrogen and is a key step for the long-sequestration of CO₂ by mineralization (Carbon Capture and Storage, CCS). In the above applications, the injected fluid is under single-phase conditions (liquid, gas or supercritical) either pure or with impurities or chemical additives (surfactants or nanoparticles for EOR). Target reservoirs can be mature Oil & Gas (O&G) and vapour fields, depleted oil and gas reservoirs, deep saline aquifers, and salt caverns, depending on the application.

It has become clear over the years that the injection process must be carefully designed to avoid severe issues, such as premature thermal breakthrough, groundwater contamination, and leakage of the reinjected fluid to the surface [4–6]. The co-injection of CO₂ (or NCGs) and water (or brine) is a technique that has been promoted to avoid these issues and enhance mineralization (under suitable conditions) [7]. The

co-injection may take place under either single- or two-phase flow conditions and the two phases (gas and water) can be mixed either on the surface before the injection well [8–16] or directly in the injection well at a certain depth [17–21]. The single-phase co-injection requires that the gas is completely dissolved in the liquid stream. For CO₂ and water, this means that the quantity of gas that can be dissolved in normal operating conditions is relatively low. For example, around 22 g of CO₂ can be dissolved in 1 kg of pure water at 50 bar and 90 °C, or higher at lower temperature and/or higher pressure [22]. In almost all geothermal exploitations, the amount of CO₂ in the resource fluid is much higher. Thus, part of the gas will continue to be vented if single-phase co-injection methods are applied.

To our knowledge, the only reinjection scenarios of CO₂/H₂S with a liquid brine stream that have reached industrial stage and demonstrate the possibility of CCS through mineralization, were those developed within the Carbfix Projects (1 and 2). These relied on the co-injection of water and CO₂ as a single-phase (liquid with pre-dissolved gas at high pressure). The original Carbfix [17] approach aimed to co-inject water and soluble gases into the subsurface in two separate streams at the surface level. Gas was released as fine bubbles into the water at depth and was completely dissolved into the geothermal brine stream before it entered the porous aquifer rocks. Water was then pumped into the annulus space, and gas was injected into a central tubing. At a later stage, in order to reduce costs and streamline the original Carbfix approach, the CO₂ and H₂S dominated gas mixture was directly captured from the power plant exhaust gas stream and dissolved into geothermal water in a scrubbing tower [23]. The pressurised gas-charged water was transported to the injection well, where it was injected together with additional effluent water into the subsurface. The two streams were mixed at a certain depth to prevent any contact between the carbon steel of the casing and the gas-charged water, which was acidic and corrosive,
and to avoid any undesirable gas release from the injected fluid [24]. The main reasons that other reinjection projects were not developed further, with the exception of some CCS projects, were the technical difficulties and uncertainties, as well as the high installation and operational costs that made industries reluctant. However, with the anticipated imposed regulations to reduce gas emissions, more and more industries are expected to turn their attention to reinjection strategies.

Well flow simulations are important for understanding fluid behaviour and modelling crucial parameters in the conception and design phases of a project. Available two-phase flow models have limitations with regards to the types of fluid, the diameter of the pipe and the orientation of the flow used in their development. Since most of them were designed for either the O&G or the nuclear industry, careful consideration should be given when adopting these for geothermal applications. Several commercial and numerical codes exist today for modelling the downward multiphase flow in injection wells. OLGA, PIPESIM, LedaFlow®, and PROSPER are among the software packages that are widely used for steady-state and dynamic calculations. Some of them have already been applied in the modelling of production and injection geothermal wells, such as OLGA [19] and LedaFlow® [25]. On the other hand, there are plenty of geothermal wellbore numerical tools, such as GWELL [26], UTWELL [27], as well as in-house codes [28]. These simulators model under steady-state or transient conditions single-phase (water or steam) or multiphase (mainly water and steam) flows. All two-phase simulators are based on well-known and established correlations for the two-phase flow to predict the flow regime, the frictional pressure drop, and the void fraction. These correlations are generally based on the separated (slip) flow model, which assumes that the two phases flow separately with different velocities and share a definite interface, or the drift-flux model, which assumes one phase is dispersed in the other continuous phase and requires the determination of the distribution parameter and the drift velocity as variables to calculate void fraction, or empirical equations [29].

The presence of dissolved salt and NCGs in the case of geothermal fluids strongly impacts the prediction of the flow behaviour. For example, the solubility of CO₂ in water reduces in the presence of dissolved NaCl [22] or other NCGs, such as H₂S [30]. Thus, an important aspect of the accuracy of any numerical tool dealing with multicomponent mixtures and multiphase flows is the prediction of the phase composition and the physical properties based on the pressure and the temperature. For that reason, most tools have integrated thermodynamic packages, which perform calculations based on Equation-of-State (EoS), black-oil (B-O) models, PVT (pressure, molar volume and temperature) tables, or simplified equations. Those based on EoS are expected to be the most precise for calculations in thermodynamic equilibrium. Combining them with solubility kinetics may increase prediction accuracy in the case of multiphase flows.

This work focuses on the simultaneous injection of NCGs and water in the same well predominantly under two-phase flow conditions. Initially, the characteristics of the two participating geothermal sites and a description of the applied injection scenarios are presented. This is followed by a description of the four applied numerical tools and their configurations. Finally, the results of the injection techniques are presented, evaluating the different injection scenarios, and comparing the results of the different tools. This study examines the injection of either a single-phase mixture, due to the complete dissolution of the gas in the water, or a two-phase mixture. In addition, preliminary simulations with water-only injection were performed for both geothermal sites.

The project's objective was twofold. Firstly, to introduce a reinjection technology (referred to as deep mixing) of NCGs together with water in the same well under two-phase flow conditions, which can be applied in several geothermal sites all over the world to reduce their environmental impact. Secondly, to validate the concept by highlighting its advantages and identifying its application limits. To do so, numerical tests with three commercially available simulators and a newly developed in-house wellbore simulator were performed for two specific sites

with different geological conditions. Simulations of injecting liquid water were used as benchmark cases to validate the in-house code.

## 2. Case studies

### 2.1. Geothermal sites

The first geothermal case was the Kizildere (K) field in the Denizli and Aydin provinces of Western Turkey at the Western extreme of the Büyük Menderes Graben, where several GPPs and wells already exist [31]. An existing injection well was exploited as a candidate for reinjection of NCGs dissolved in a water stream, targeting a high-temperature and high-pressure reservoir. Currently, and within the framework of the GECO H2020 project, injections in the same well of water and CO₂ under single-phase flow conditions (complete dissolution of CO₂) are ongoing.

The second site was the Castelnuovo (C) power plant, which is still at the proposal phase and has been designed to exploit a dry-steam geothermal resource in the Montecastelli Pisano area, located in the Eastern part of Anqua-Pomarance, within the vapour dominant field of Larderello [32]. The injection targeted a high-temperature and low-pressure reservoir, and data were taken from the case study described by Vaccaro et al. [33].

### 2.2. Injection scenarios

Two reinjecting techniques were considered. First, the mixing of the two phases on the surface and the direct injection of the resulting mixture from the wellhead (surface mixing, Inj1). Second, the injection of the two phases separately but simultaneously in the annulus and the central tube of a dual-pipe (coaxial) well and mixing of them at some depth below the surface (deep mixing, Inj2). In the latter case, the gas could be injected either in the annulus part or the central tubing, and water into the other compartment. Simulations were performed under single-phase flow conditions, injecting water-only (Inj0), and results obtained using the commercial software were used to benchmark the in-house code.

For the Kizildere case study, two well configurations were investigated. The first configuration explored water-only (K-Inj0) or H₂O-CO₂/ NCG (K-Inj1) mixture injections through the casing (without an internal tubing), while the second explored the injection of water through the annulus and CO₂/NCG through the internal tubing and mixing them at a depth in the well (K-Inj2), similar to the Carbfix solution.

Surface mixing and injection from the surface of a two-phase mixture in the case of Castelnuovo (C-Inj1) was not considered for physical reasons. In fact, due to the reservoir conditions (pressure and injectivity), the water column inside the well was located several hundred meters below the surface and the upper part of the well was in "vacuum" conditions. Even after considering the nominal water-only injection flow rate, the water interface could only rise a few hundred meters. Directly injecting a two-phase mixture into this location of the well would result in severe conditions (flashing of the injecting fluid, gas accumulation, and pressure build-up) preventing downward fluid direction and further injection from the wellhead.

To ensure that the gas would be driven towards the reservoir a new injection technology and deep mixing strategy was proposed [18,34]. Gas and liquid phases were injected separately from the wellhead in the annulus part and the central tubing, respectively, a solution alternative to the Carbfix technology. One or multiple one-way (inversed gas-lift) valves were located along the well connecting the two compartments. At least one valve was placed at a depth that was inside the initial liquid column. In the multiple valves scenario, the injection starts by opening the deepest valve, located inside the liquid, causing the fluid column in the tubing to rise. When the column rose above the level of the next valve, by a minimum predefined threshold (fluid height), that valve opened and the deeper one closed. The same sequence was applied when
<table><caption>**Table 1**
Summary of modelling scenarios and numerical tools applied.</caption><thead><tr><th>Site</th><th>Injection method</th><th>Tool</th><th>Hydrodynamic approach</th><th>Thermodynamic approach</th><th>Case label</th></tr></thead><tbody><tr><td rowspan="4">Kizildere (K), Turkey</td><td>Water-only (Inj0)</td><td>PIPESIM<br>GWellFM</td><td>Hagedorn &amp; Brown<br>Custom</td><td>Compositional<br>Black-oil<br>CPA</td><td>K-Inj0</td></tr><tr><td>Surface mixing (Inj1)</td><td>PIPESIM<br>OLGA<br>GWellFM</td><td>Hagedorn &amp; Brown<br>Hagedorn &amp; Brown<br>Custom</td><td>Compositional<br>Black-oil<br>Black-oil<br>CPA</td><td>K-Inj1</td></tr><tr><td>Deep mixing (Inj2)</td><td>PIPESIM<br>GWellFM</td><td>Hagedorn &amp; Brown<br>Custom</td><td>Black-oil<br>CPA</td><td>K-Inj2</td></tr><tr><td>Water-only (Inj0)</td><td>UniSim®<br>GWellFM</td><td>Beggs &amp; Brill<br>Custom</td><td>Sour-PR<br>CPA</td><td>C-Inj0</td></tr><tr><td rowspan="2">Castelnuovo (C), Italy</td><td rowspan="2">Deep mixing (Inj2)</td><td>One valve</td><td>UniSim®<br>GWellFM</td><td>Beggs &amp; Brill<br>Custom</td><td>Sour-PR<br>CPA</td><td>C-Inj2a</td></tr><tr><td>Three valves</td><td>UniSim®<br>GWellFM</td><td>Beggs &amp; Brill<br>Custom</td><td>Sour-PR<br>CPA</td><td>C-Inj2b</td></tr></tbody></table>

<table><caption>**Table 2**
Two-phase hydrodynamic models used in the numerical tools.</caption><thead><tr><th>Tool</th><th>Model</th><th>Type</th><th>Flow regime</th><th>Comments</th></tr></thead><tbody><tr><td>PIPESIM, OLGA</td><td>Hagedorn &amp;Brown</td><td>Empirical (α<sub>G</sub> &amp; ΔP)</td><td>Independent</td><td>Small diameter &amp; vertical pipes, water/oil &amp; air</td></tr><tr><td>UniSim®</td><td>Beggs &amp; Brill</td><td>Empirical (α<sub>G</sub> &amp; f<sub>TP</sub>)</td><td>Dependent</td><td>Horizontal &amp; upwards inclinations, small diameter pipes, water &amp; air</td></tr><tr><td>GWellFM</td><td>Custom</td><td>Drift flux (α<sub>G</sub>), semi-empirical (ΔP)</td><td>Independent</td><td>Large diameter pipes &amp; downwards flows</td></tr></tbody></table>

<table><caption>**Table 3**
Thermodynamic models used in the numerical tools.</caption><thead><tr><th>Tool</th><th>Model</th><th>EoS<sup>1</sup></th><th>Transport properties<sup>2</sup></th><th>Comments</th></tr></thead><tbody><tr><td>PIPESIM, OLGA</td><td>Black-oil</td><td>-</td><td>Equations</td><td>Not suitable for pure H<sub>2</sub>O &amp; multicomponent mixtures</td></tr><tr><td>PIPESIM</td><td>Compositional</td><td>PR</td><td>Equations for μ &amp; σ, constant λ</td><td>EoS not suitable for aqueous mixtures, equations for O&amp;G applications</td></tr><tr><td>UniSim®</td><td>Compositional</td><td>Sour-PR</td><td>Equations</td><td>EoS adapted for acid gases in contact with an aqueous phase, equations</td></tr><tr><td>Multiflash</td><td>Compositional</td><td>CPA</td><td>Equations</td><td>EoS suitable for aqueous mixtures, Specific equations</td></tr><tr><td>Carnot</td><td>Compositional</td><td>CPA</td><td>Equations per phase</td><td>EoS suitable for aqueous mixtures, Specific equations</td></tr></tbody></table>

switching from the second valve to the third one and so on for additional valves. The multiple injection solution was employed to ensure the entrainment of the gas phase to the liquid phase, and minimise the compression needs on the surface.

Two different deep mixing injection (Inj2) simulations were conducted for Castelnuovo; mixing of NCG and water through one deep point (single deep injection, C-Inj2a) and through three points (triple deep injection, C-Inj2b) located at different depths.

### 3. Methods

Three different commercial tools (PIPESIM steady-state flow simulator, OLGA dynamic flow simulator, and UniSim® Design modelling platform) and one in-house flow simulator (GWellFM) were used to model the injection techniques at the two geothermal sites. Calculations were performed under steady-state conditions with PIPESIM and GWellFM and transient conditions with OLGA and UniSim®. In the latter cases, only the final steady-state results were used for comparison purposes. All four codes discretised the domain and solved the mass, momentum and energy transport equations (see Appendix 1). Table 1 summarises performed simulations, and numerical tools applied for each study case.

For each numerical tool different types of two-phase hydrodynamic models (Table 2) and thermodynamic calculations (Table 3) were selected due to specific reasons and limitations. For thermodynamic calculations, additional equations were used to calculate transport properties. The black-oil (B-O) model was the only available option in OLGA, while in PIPESIM an EoS model with its default options was applied in addition to the B-O model. On the other hand, in UniSim® and GWellFM models for the specific type of mixtures were selected. In GWellFM, two options were available and the same EoS (with different equations for the transport properties) were used, for comparison purposes.

#### 3.1. Steady-state flow simulations with PIPESIM

PIPESIM is a steady-state flow simulator that can solve multiphase flow problems using physics-based governing equations (see Appendix 1). Two injection scenarios were investigated for a well injecting CO<sub>2</sub> or NCG and water in the Kizildere site (Table 1). Gas and water injection rates as well as reservoir injectivity sensitivity simulations were performed. For each scenario completion configurations were set up, one for injecting through the casing (Fig. 1a, K-Inj1) and one with a central tubing for mixing the two phases in a certain depth in the well (Fig. 1b, K-Inj2).

Although PIPESIM software does not have a direct interface and the tools to model separate flows in the annulus and the tubing (K-Inj2), the available features are enough to mimic such conditions. The fluid injected from the surface was defined as pure water. A pipe plug was inserted just two meters before the end of the tubing, as not to let the

<footnote n="1">¹ Phase composition, phase density.</footnote>

<footnote n="2">² Viscosity (μ), thermal conductivity (λ), surface tension (σ).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Tubular configurations and ambient temperature profiles in PIPESIM for (a) water or H₂O-CO₂ mix injection through the casing and (b) injection of water through the annulus and injection of the CO₂ through the inner pipe for the Kizildere case study.</figcaption></figure>

injected water flow back. A sliding sleeve was installed below the plug to ensure the water flow enters the tubing and a gas lift valve was placed at the end of the tubing to inject the CO₂/NCG and mix it with water. Moreover, once the simulation was performed and the pressure value at the valve location was obtained, the calculations for the gas-only (CO₂/NCG) flow through the tubing were done, and the pressure profile along the tubing was obtained. In the simulations, the diameter of the tubing was set at 2⁷/₈" and its initial length at 500 m. The domain was discretised with a maximum length segment of 10 m (similar to OLGA discretisation), and the tolerance for the pressure and the mass balances was set at 0.01 % of the relative error.

Several hydrodynamic models for two-phase flows are available. In this study, the Hagedorn & Brown model [35] was applied for the vertical downward multiphase flow. This model was developed from an experimental study of pressure gradients during continuous two-phase flows (water/air, oil/air) in vertical conduits with small diameters. It is used for calculating the pressure losses and the liquid holdup but does not predict the flow pattern. The Orkiszewski model [36] was additionally used in the study to define the two-phase flow regime.

PIPESIM model input scheme requires the reservoir depth, pressure, temperature, and injectivity. The reservoir depth was taken as the approximate midpoint of the open-hole section with a slotted liner, which occurs at 2200 m below ground. At this depth, the measured static pressure and temperature were used as the reservoir pressure and temperature, with values of 185 bar and 202 °C, respectively. For pure water injections an injectivity of 3.6 kg/s/bar was selected. A 27 bar pressure difference between the bottom hole and the reservoir was necessary for the water to penetrate the reservoir.

Another required parameter was the gas flow rate. The amount of gas to be injected with water was initially estimated based on the solubility of CO₂ in brine [22]. Solubility calculations suggested that the soluble CO₂ content may be as high as 66 g of CO₂ per kg of water at the reservoir pressure and at the injection temperature.

PIPESIM offers the options to perform fluid flow simulations either with a black-oil fluid model or with a compositional fluid model (see Appendix 3). Black-oil fluid modelling utilises correlations to simulate the key PVT fluid properties of oil/gas/water systems. These empirical correlations treat the oil/gas system as a simple two component system, unlike the more rigorous multi-component compositional methods which use an EoS to determine the behaviour and the thermodynamic properties. To observe differences between the two options, the simulations were performed first with the black-oil, and then with the compositional model.

The molar composition of the injected NCG is provided in Table 4. When injection of both water and NCGs were considered, the injectivity was set at 3.025 kg/s/bar to maintain the same pressure difference between the average pressure of the reservoir and the flowing bottom-hole pressure.

## 3.2. Transient flow simulations with OLGA

OLGA is a transient multiphase simulator for flows in pipes, flow lines and wells. Transient flow simulations were performed including a surface connection line. Water and NCGs were injected separately and mixed in a junction node representing the mixing of gas and water on the surface (K-Inj1). The rates and the injection pressures were specified at the injection nodes. The mixed fluid was then sent through the surface line to the wellhead and injected into the reservoir. The fluids were modelled with the black-oil model because it was the only active option in the available OLGA licence, and the same hydrodynamic models as in PIPESIM were applied. Despite this limitation, the transient simulations with OLGA were considered in the assessment of the numerical tools. In this paper, only the steady-state results inside the well are shown. The same cell size and convergence criteria as in PIPESIM were applied with OLGA. Such discretisation offered a good compromise between accuracy and computational cost.

<table><caption>**Table 4**
Molar and mass composition of NCG to be injected back into the reservoir for the Kizildere study case.</caption><thead><tr><th>Component</th><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>C<sub>3</sub></th><th>N<sub>2</sub></th><th>CO<sub>2</sub></th><th>H<sub>2</sub>S</th><th>Ar</th><th>O<sub>2</sub></th><th>H<sub>2</sub></th><th>He</th></tr></thead><tbody><tr><td>%mol</td><td>0.720</td><td>0.010</td><td>0.002</td><td>0.300</td><td>98.873</td><td>0.072</td><td>0.001</td><td>0.004</td><td>0.017</td><td>&lt; 10<sup>-3</sup></td></tr><tr><td>%wt</td><td>0.264</td><td>0.007</td><td>&lt; 10<sup>-4</sup></td><td>0.192</td><td>99.481</td><td>0.056</td><td>&lt; 10<sup>-4</sup></td><td>&lt; 10<sup>-4</sup></td><td>&lt; 10<sup>-4</sup></td><td>&lt; 10<sup>-4</sup></td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) Configuration of dynamic simulation and (b) simulation outline in UniSim® for the Castelnuovo case study.</figcaption></figure>

<table><caption>**Table 5**

Geometry configuration of UniSim® simulation for the Castelnuovo case.</caption><thead><tr><th>Vertical depth (m)</th><th>I.D. (mm)</th><th>O.D. (mm)</th><th>Roughness (µm)</th><th>U (W/m², K)</th></tr></thead><tbody><tr><td>0–2000</td><td>90</td><td>99</td><td>45.72</td><td>13.6</td></tr><tr><td>2000–2500</td><td>152.4</td><td>162</td><td>45.72</td><td>9.6</td></tr></tbody></table>

### 3.3. Transient flow simulations with UniSim®

UniSim® is a design suite software that allows the construction of dynamic and steady-state simulation models under variable environmental conditions. The dynamic simulations were conducted for the Castelnuovo case study, and the outline is shown in Fig. 2. The model consisted of four blocks: the well, the drain valve, the compressor, and the reinjection valve or valves (multi-point injections) located along the well. The well-block was the main part; the whole well depth of 2500 m was divided into 25 vertical segments, with 100 elements each, in order to calculate with a sufficient accuracy the dynamic momentum, the pressure gradient, the heat transfer, the solubility, and the void fraction. Refining further the mesh improved slightly the calculation accuracy, especially for the two-phase flows, but significantly increased the computational time. The upper part of the well was an annulus, whereas the lower part was a tube with a larger diameter, according to Table 5. Tolerance was set at a relative error of 0.1 % for energy, mass, and mole balances.

The Beggs & Brill [37] flow model was employed to determine the flow pattern, the liquid holdup (or void fraction), and the two-phase pressure losses. This model is applicable for any flow orientation in different inclination angles and includes correlations for the two-phase friction factor in various flow patterns: segregated (or annular), distributed (or dispersed), intermittent, and transition. However, its validity for vertical downward multiphase flows is limited.

Similar to PIPESIM, UniSim® cannot model flows in coaxial wells. Here, only the heat transfer in the upper dual compartment part was considered and the flow was not solved. The heat transfer of the upper part (0–2000 m) was considered to occur mainly between the NCGs (CO₂ and H₂S mixture) in the annulus section and the liquid-dominant fluid in the central tube, while the heat transfer of the lower part (2000–2500 m) occurred between the fluid and the rocks. The heat transfer coefficient and the heat transfer area were calculated based on the segment geometry and the temperature boundaries. The heat transfer coefficient was determined based on the local conditions, the inside film convection, and the outside conduction/convection with the gas/ground. The sour-type Peng-Robinson (PR) EoS was used for thermodynamic equilibrium and physical property calculations. The Sour-PR model combines the conventional PR EoS, and Wilson's API-Sour model to precisely determine the ionization of the H₂S and CO₂ in the aqueous phase (see Appendix 3).

The drain valve block adjusted the liquid flow and included a control valve, a controller with setpoint ramping, which enabled the valve to smoothly change the flow rate, and a vessel that left space when the liquid underwent a reduction in pressure, and a part of it flashed into vapor. The compressor block included a compressor, a bypass which was open only before the injection process, and a branch for distributing the NCGs between the injection valves. The simulation scenario consisted of the following stages: a) compressor start-up, b) shutting the compressor bypass, c) opening the injection valve, d) continuously checking the water level along the well, and e) opening the next injection valve in parallel of shutting the previous one (when multiple valves are present along the well).

The length of the initial hydrostatic column was around 747 m ($P_{res}/\rho g$) and the depth of the water surface in the well was about 1753 m. Thus, in the single injection configuration, the mixing point was placed at 1800 m depth (inside the liquid column), whereas in the three injection point configurations, the valves were located at depths of 1800 m.
<table><caption>**Table 6**
Process detail of UniSim® simulation for the Castelnuovo study case.</caption><thead><tr><th>Description</th><th>Drain side</th><th>NCG side</th></tr></thead><tbody><tr><td rowspan="3">Composition (%wt)</td><td>H<sub>2</sub>O</td><td>99.9</td><td>0</td></tr><tr><td>CO<sub>2</sub></td><td>0.1</td><td>97.5</td></tr><tr><td>H<sub>2</sub>S</td><td>0</td><td>2.5</td></tr><tr><td>Mass flow rate, m (kg/s)</td><td>16.56</td><td>1.44</td></tr><tr><td>Wellhead pressure, P<sub>wh</sub> (bar)</td><td>Variable</td><td>Variable</td></tr><tr><td>Wellhead temperature, T<sub>wh</sub> (°C)</td><td>89</td><td>89</td></tr></tbody></table>

<table><caption>**Table 7**
Reservoir conditions imposed in UniSim® for the Castelnuovo study case.</caption><thead><tr><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>Surface temperature, T<sub>surf</sub> (°C)</td><td>25</td></tr><tr><td>Geothermal gradient</td><td>Linear</td></tr><tr><td>Reservoir pressure, P<sub>res</sub> (bar)</td><td>70</td></tr><tr><td>Reservoir temperature, T<sub>res</sub> (°C)</td><td>280</td></tr><tr><td>Injectivity index, IP (kg/s/bar)</td><td>0.5</td></tr><tr><td>Bottom-hole pressure, P<sub>bh</sub> (bar)</td><td>106</td></tr></tbody></table>

m, 1400 m, and 1000 m (again the deeper valve was located inside the initial liquid column). The threshold for every valve to open was set to 50 m.

The model simulated the dynamic behaviour of the system from start-up to the final steady conditions. The final conditions were reached when the design NCGs flow rate was injected into the liquid through the uppermost injection valve, and both the water level and the compressor pressure were stabilised. It is important to note that the initial filling of pipes with gas was necessary for the simulations. Hence, a negligible amount of CO₂ was considered in the drain, which had a negligible impact on the injection. Further details and the boundary conditions are given in Tables 6 and 7.

### 3.4. Steady-state flow simulations with GWellFM

GWellFM (Geothermal Well Flow Model) is a steady-state 1D non-isothermal multicomponent and two-phase flow simulator. The model

<table><caption>**Table 8**
Summary of numerical tools.</caption><thead><tr><th>User</th><th>Tool</th><th>Limitations</th></tr></thead><tbody><tr><td>Middle East Technical University</td><td>PIPESIM</td><td><ul><li>Lack of interface to create coaxial well</li><li>Hydrodynamic model not adapted for water and NCG mixtures</li><li>Two-phase flow model not adapted for downward flows in large diameter pipes</li><li>Thermodynamic model not adapted to water and NCG mixtures</li></ul></td></tr><tr><td>Middle East Technical University</td><td>OLGA</td><td><ul><li>Hydrodynamic model not adapted for water and NCG mixtures</li><li>Two-phase flow model not adapted for downward flows in large diameter pipes</li><li>Thermodynamic model not adapted to water and NCG mixtures</li></ul></td></tr><tr><td>University of Florence</td><td>UniSim®</td><td><ul><li>Lack of interface to create coaxial well and simulate flow in the annulus</li><li>Hydrodynamic model not adapted for water and NCG mixtures</li><li>Two-phase flow model not adapted for downward flows in large diameter pipes</li></ul></td></tr><tr><td>IFP Energies Nouvelles</td><td>GWellFM</td><td><ul><li>Two-phase flow model not validated for water and NCG flows</li><li>Cannot simulate flow reversal or accumulation of gas phase for high gas velocity and/or low liquid velocity</li></ul></td></tr></tbody></table>

considers the single-phase flows of liquids and gases, the hydrodynamics of the two-phase downward and upward flows, constitutive laws for mixtures, the flow of fluids through orifices/valves, and finally the heat exchange between the well completion and the surrounding formation. The code was developed in Python and discretises and solves iteratively along a mesh the mass, the momentum, and the heat balances (see Appendix 1).

Hydrodynamic modelling consists of calculating the velocities of fluids and pressure drop along with the flow. For single-phase flows, the calculations are quite simple and well established. However, for two-phase flows the calculations are more complicated, and additional flow parameters, such as the void fraction (or liquid hold-up), the frictional pressure drop, and the phase distribution (flow regime), must be defined. Even though the two-phase gas-liquid flows attracted much attention over the last several decades, mainly for O&G applications and the nuclear industry, vertical downward flows were seldom considered and there is a lack of extensive experiments and data.

Several models have been proposed for the void fraction, but those based on the drift flux model have been proven the most accurate [38]. Yao et al. [39] developed a semi-empirical correlation for predicting the pressure drop due to friction of gas-liquid downward flows, adding a parameter to consider the effect of buoyancy. The above two models are independent of the flow configuration. The two phases can flow according to different flow patterns (regimes), which are determined by the interfacial structures between both phases and are known to vary with many factors, including the fluid properties, the flow channel size, geometry and orientation, the injection method, and the flow rates. Lokanathan & Hibiki [40] proposed criteria for all transitions; from dispersed (bubbly and cap-bubbly) to intermittent (slug and churn-turbulent) and annular patterns.

The thermal modelling (see Appendix 1) consists in calculating the heat transfer in the fluids and between the fluids in the wellbore and the surrounding bedrock system until the undisturbed formation temperature and it is based on the principle of the thermal resistances. Two mechanisms of energy transfer are considered: (i) heat exchange between the solid materials (i.e., rocks, steel) due to the conduction and (ii) convection transfer due to the fluid motion in the wellbore. When two phases are present, they are always in thermal equilibrium and their temperatures are equal.

The model is fully compositional and, in order to perform thermodynamic calculations, two advanced thermodynamic engines were integrated into the code. The commercial Multiflash (MF) software, which has a comprehensive PVT (pressure, volume, temperature) and physical properties package that quickly and reliably allows the complete modelling of the phase behaviour of complex mixtures and pure substances through equations of state; and the Carnot library, which is an in-house thermodynamic library, performing flash calculations using robust isothermal-isobaric algorithms coupled to various real-fluid EoS to compute the vapour-liquid equilibrium as well as the transport properties. Due to the presence of water and aqueous mixtures, the Cubic-Plus-Association (CPA) model was chosen to be used, in both models, since it explicitly considered the association of water molecules (see Appendix 3). This EoS can reproduce the properties of the vapour and water phases with the same model [41]. For viscosity and thermal conductivity, specific models of the aqueous and the vapour phases were used.

Simulations were performed for both sites and all injection scenarios, and the results were compared with the other numerical tools. The simulated GWellFM well geometries and conditions (inclination, diameter, depth, wall roughness, material thickness and thermal properties, injection pressure and temperature, fluid composition) were identical to those imposed in the other tools. The only difference was in the configuration of the coaxial well for the Castelnuovo site. In this case, the complete geometry consisting of the central tubing and the annulus side was constructed and the flows in both compartments were modelled. In UniSim®, on the other hand, a simple pipe (without
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Pressure and (b) temperature profiles calculated by PIPESIM, with the black-oil and the compositional options, and GWellFM, with the Carnot and the Multiflash options, for water-only injection (K-Inj0).</figcaption></figure>

annulus was considered and the fluid flow in the central tubing was only modelled. The annulus part was only considered in the heat transfer modelling through the overall heat transfer coefficient. In all simulations, the domain was uniformly discretised with a cell size of 10 m, and the convergence criteria of the pressure and the enthalpy were set at 10⁶. Simulations performed with finer meshes showed that the increase of accuracy was negligible.

### 3.5. Numerical solution

The simulation tools proceed similarly using iterative procedures. Starting from guess wellhead pressures, the pressure at the well outlet is calculated. The wellhead pressures correspond to the needs of pumping and/or compression on the surface for injecting the liquid and the gas phase or the two-phase mixture. The calculated pressure, $P_{calc}$, must always match the bottom-hole pressure, $P_{bh}$, considering the reservoir

pressure, $P_{res}$, the injectivity index, $IP$, and the total injection mass flow rate, $\dot{m}_{tot}$, according to Eq. (1). These parameters, in all simulations, are known and used as boundary conditions.

$$P_{calc} = P_{bh} = P_{res} + \frac{\dot{m}_{tot}}{IP} \qquad (1)$$

On the other hand, the temperatures of the injection fluids were considered constant, and the geothermal gradient (the undisturbed temperature of the formation) was also a boundary condition.

Table 8 summarises the limitations and the user of each code.

### 4. Results

During the injection of multiphase mixtures, the parameters of importance are the pressure, the temperature, and the two-phase flow regimes along the well. Moreover, some thermodynamic parameters are

<figure><img src="image_5.png" /><figcaption>Fig. 4. (a) Pressure and (b) temperature error analysis between the different models for water-only injection (K-Inj0).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of (a) pressure, (b) average velocity, and (c) mixture density between UniSim® and GWellFM, with the Carnot and the Multiflash options, for water-only injection (C-Inj0).</figcaption></figure>

also calculated. On the contrary, the mixture density, $ρ_m$, is calculated from the phase densities and the void fraction, and the average mixture velocity, $u_{av}$, based on the total mass flow and the total flowing cross-section of the well (see Appendix 2).

The error $δ$, Eq. (2), is introduced to compare the results between the different codes, where $Φ$ can be either the pressure or the temperature, and the subscripts 1 and 2 correspond to the case under investigation and to the reference case, respectively.

$$ \delta\Phi = \frac{\Phi_1 - \Phi_2}{\Phi_2} 100 \qquad (2) $$

## 4.1. Water-only injection

The first simulation set was carried out for the water-only injection; 350 t/h of water was injected into the Kizildere geothermal reservoir (K-Inj0) through the well having the configuration given in Fig. 1a. The pressure and temperature profiles from PIPESIM with the black-oil fluid representation and the compositional fluid models are shown in Fig. 3. As the injection of water into the reservoir requires a 27 bar pressure difference and the reservoir pressure is 185 bar, the pumping pressure at the wellhead should be around 10 bar at least. Both fluid models estimate very similar pressure profiles, with only a 3 bar discrepancy at the wellhead. The temperature profile shows that the temperature change from the surface injection temperature is not significant with 2 °C using the black-oil model and 3 °C using the compositional model.

GWellFM was compared with PIPESIM under the same single-phase conditions. In all cases (Fig. 3), the difference between the two simulators in pressure prediction is maximum at the wellhead considering the same bottom-hole pressure condition, whereas for the temperature the maximum is observed in the bottom-hole as injection temperature is imposed at wellhead. As expected, both thermodynamic options of GWellFM converge better closer to the compositional option of PIPESIM, even though the black-oil option also gives similar results.

The error analysis between the different models is given in Fig. 4. For the pressure, the graph shows only the first few hundred meters where the divergence is higher. Generally, GWellFM predicts slightly smaller pumping (wellhead) pressure. Even for the simple case of liquid water the different thermodynamic models do not predict the same conditions in the well. The black-oil model gives the larger divergence when

compared to the compositional models except for the comparison with the Multiflash option of GWellFM and the pressure evolution. The comparison of the two options of GWellFM with the compositional model of PIPESIM shows that the Carnot model agrees better in the pressure evolution, whereas the Multiflash option in the temperature profile. This is due to the different EoS, and the separated equations used for calculating the properties related to the pressure losses (density, viscosity) and the heat transfer calculations in the fluid (thermal conductivity). The same stands when comparing Carnot with Multiflash, since in this case all other parameters are identical and any discrepancy results exclusively from the thermodynamic calculations. More details will be given in comparison section.

Water-only (containing 0.1 %wt CO₂) injection tests were also conducted with UniSim® and GWellFM for the Castelnuovo geothermal case (C-Inj0); the results are compared in Fig. 5. As mentioned earlier, the water level in the well is located several meters below the surface due to the reservoir conditions. Even after the water injection, the surface rises but remains several hundred meters below the surface. This depth was found to be between 1347 and 1348 m by UniSim® and 1360–1370 m by GWellFM with the Multiflash option and between 1340 and 1350 m with the Carnot option.

In the part of the well above the water column, both models predict that the pressure is below 1 bar, and a two-phase flow of liquid and vapour mixture appears. A liquid flashing is occurring as soon as the fluid enters the well from the surface injection valve. This pressure is slightly different in the two tools, mainly because of the different thermodynamic models and the presence of the small quantity of CO₂ in the injected water; UniSim® gives an initial wellhead pressure of 0.90 bar and GWellFM 0.83 bar with the Multiflash option and 0.82 bar with the Carnot option. Further down, both models predict that the pressure remains practically constant (pressure losses due to friction are compensated by gravity) until the liquid water surface is reached. In this region, there is a divergence in the average mixture velocity, $u_{av}$, due to the differences in the hydrodynamic (void fraction, $α_G$) and thermodynamic (phase density, $ρ$) calculations. UniSim® predicts a void fraction of about 0.73 and a mixture density, $ρ_m$, of 263.17 kg/m³ (gas and liquid density 0.77 and 957.89 kg/m³), whereas with GWellFM these values are 0.66 and 326.17 kg/m³ ($ρ_G = 0.66 kg/m³$, $ρ_L = 971.46 kg/m³$) with Multiflash option and 0.66 and 322.99 kg/m³ ($ρ_G = 0.65 kg/m³$, $ρ_L = 958.30 kg/m³$) with Carnot option. Like the K-Inj0 case, the impact of
<figure><img src="image_7.png" /><figcaption>Fig. 6. (a) Pressure, (b) temperature, and (c) void fraction calculated by PIPESIM, with the black-oil and the compositional model, for injecting H₂O-CO₂ mixture from the surface (K-Inj1: solid lines, K-Inj0: dashed lines).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Impact of injectivity on (a) pressure, (b) temperature, and (c) void fraction calculated by PIPESIM, with the black-oil and the compositional options, for surface injection of H₂O-NCG mixture (K-Inj1).</figcaption></figure>

the thermodynamic models is obvious in the results.

The two-phase flow regime is always annular according to GWellFM, with the liquid phase flowing down attached to the walls and the gas phase remaining in the core of the central tubing. Then, and before reaching the water surface, the flow regime gradually passes to intermittent, and to dispersed pattern.

In the liquid dominant part of the well, both models give similar pressure, velocity, and density profiles since for single-phase flow the pressure losses are calculated with well-documented equations based on the friction factor and the estimation of the physical properties of the fluid (mainly water) is well-established. A modest impact of the CO₂ presence and its interaction with H₂O, however, can be still seen in the mixture density, similar to the upper part of the well.

## 4.2. Surface mixing

In the second simulation set with PIPESIM (K-Inj1), a mixture of 300 t/h water and 23.1 t/h CO₂ was injected from the surface. The pressure, temperature, and flowing gas volume fraction profiles from PIPESIM with the black-oil and the compositional fluid models are shown in Fig. 6. The water-only results are also included for comparison. The surface pressure requirement increases to approximately 50 bar when CO₂ is mixed with the water instead of 10 bar for the water-only injection, and the flowing temperature in the well changes slightly; the downhole temperature increases by 1 °C reaching 108 °C. The difference between the black-oil and the compositional fluid models for the temperature profiles is larger than in the water-only injection case. However, the difference in the pressure profiles is yet insignificant. The gas
<figure><img src="image_9.png" /><figcaption>Fig. 8. Impact of water flow rate on (a) pressure, (b) temperature, and (c) void fraction calculated by PIPESIM, with the compositional option, for surface injection of H₂O-NCG mixture (K-Inj1).</figcaption></figure>

fraction profiles from both fluid models are very similar; the gas volume fraction decreases with depth or increasing pressure. However, both fluid models estimate that the CO₂ amount does not dissolve completely in water, which is a contradicting result with the solubility calculations. Thus, the calculations of Duan & Sun [22] model are not in agreement neither with the black-oil model nor with the compositional model.

Similarly, simulations were performed by injecting a mixture of 300 t/h of water and 6 t/h of NCG, and the results of the pressure, temperature, and void fraction with the two flow models are presented in Fig. 7. Following the preliminary H₂O/CO₂ simulations, the complete composition of the NCG stream was used (Table 4) and a lower injection temperature was applied to increase the gas dissolution in water. Also, the quantity of NCG was reduced, compared with the above simulations, for achieving the complete dissolution of the gas.

Black-oil and compositional fluid models calculate almost the same pressure and temperature profiles. The void fraction in the upper section of the well is also very similar. However, the compositional model estimates that the complete dissolution of the gas phase in the water phase occurs below 800 m. In the black-oil model, the gas phase reaches a volumetric ratio of 3–4 % and then stays almost constant below 800 m. The complete dissolution of the gas phase impacts the fluid temperature. As soon as the compositional model predicts the complete dissolution of the gas phase, there is a difference of around 1 °C between the two calculated temperatures. This temperature jump corresponds to the heat release due to the exothermic nature of the gas dissolution reaction.

The effect of the injectivity on the pressure, the temperature, and the void fraction profiles is also shown in Fig. 7. The case was simulated with the compositional fluid model. Decreasing injectivity increases the

<figure><img src="image_10.png" /><figcaption>Fig. 9. Impact of gas flow rate on (a) pressure, (b) temperature, and (c) void fraction calculated by PIPESIM, with the compositional option, for surface injection of H₂O-NCG mixture (K-Inj1).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Comparison of (a, d) pressure, (b, e) temperature, and (c, f) void fraction calculated by PIPESIM, with the black-oil option, for deep mixing of (a-c) H₂O-CO₂ and (d-f) H₂O-NCG (K-Inj2: solid lines, K-Inj1: dashed lines).</figcaption></figure>

flowing bottom-hole pressure for the given reservoir pressure, therefore the pressure along the well and the required pumping pressure are increased to inject the same amount of water and NCG. On the other hand, the temperature remains practically constant. Only the location of the temperature jump moves upwards to follow the location of the complete gas dissolution which occurs closer to the surface, as can be seen from the void fraction profile. Due to the higher calculated pressure for the lower injectivity more gas is dissolved in the liquid phase and the void fraction is always lower.

Sensitivity studies were performed to investigate the effect of water and gas injection rates. As shown in Fig. 8, increasing the water rate (for $m_G = 5$ t/h) increases the pressure along the well because the higher injection rate increases the bottom-hole flowing pressure for the specified reservoir conditions ($P_{res} = 200$ bar, $IP = 2.135$ kg/s/bar). The faster the flow the less heat transfer is subject to, therefore the temperature along the well increases less (or the temperature profile shifts to lower temperature) with increasing water injection rate. Increasing water rate decreases the flowing gas volume (i.e., the void fraction). The gas phase disappears (i.e., the void fraction becomes zero) at a deeper level of the

well as the water injection rate increases because of the subsequent increase of the pressure along the well for the increasing rate. Besides, the expected flow regimes are shown in the figure; from the wellhead down to the level where the gas phase disappears, first the intermittent (slug) flow and then the dispersed (bubble) flow regimes are expected. Increasing the water rate shortens the length of the slug flow regime. The impact of higher liquid flow rates is the same as lower injectivity that tested previously.

The second sensitivity study examined the gas injection rate. The results are compared in Fig. 9. Increased gas injection rate (for $m_L = 200$ t/h) requires injecting a gas-water mixture with higher pressure at the wellhead. However, since the dominant phase for the pressure gradient is the water phase, the pressure profiles approach to each other as the water-gas mixture flows down the well, and the void fraction decreases until the gas phase disappears completely. The temperature profiles shift to higher values because higher gas content increases the heat transfer coefficient of the gas-water mixture due to the increase of the mixture velocity and the lower conductivity of the gas phases. Increasing water rate decreases the flowing gas volume. The gas phase disappears at
<figure><img src="image_12.png" /><figcaption>Fig. 11. (a) Well configuration, (b) pressure, (c) temperature, and (d) average velocity, and (e) flow configurations calculated by UniSim® for injection of H₂O-NCG and mixing occurring in different deep locations (C-Inj2a, C-Inj2b).</figcaption></figure>

almost the same depth of the well for all different gas injection rates. In addition, the expected flow regimes are also shown in the figure. From the wellhead down to the level where the gas phase disappears, first the slug flow and then the bubble flow regimes are expected. Decreasing the gas rate shortens the depth of the well where the slug flow regime can be observed. In both Figs. 8 and 9, the temperature jumps due to the complete dissolution of the gas are present in the profiles.

### 4.3. Deep mixing

The simulations with PIPESIM were repeated for the second well configuration of deep mixing (Fig. 1b). The gas phase, CO₂ or NCG, is injected through a separate tubing and the injection of water occurred through the tubing-casing annulus from the surface (K-Inj2). Gas and water are mixed downstream as CO₂/NCG leaves the tubing. The water and gas injection rates and the injection temperature were the same as before. The tubing length in the case of NCG injection was kept at 800 m, the depth at which the gas is completely dissolved in water according to the previous observations. The pressure, temperature, and flowing gas volume fraction profiles with the black-oil fluid model are shown in Fig. 10. The results are obtained solely with the black-oil fluid model, because the compositional fluid model simulations have failed.

Simulations with GWellFM revealed that part of the gas in the tubing transforms in its supercritical state. The profiles from surface mixing are also included for comparison purposes. As shown, the pumping pressure required for water injection decreases compared to the case of water-gas mixture injection from the surface. However, the compressor pressure for gas injection through the tubing increases. The temperature for this case stays slightly lower along the well. Soon after the mixing of the gas phase with the water (at 500 m for CO₂ and 800 m for NCG), the gas volume fraction follows the same profile as the water-gas mixture surface injection because the similar conditions in the well.

Simulations with UniSim® were performed for two different deep injection scenarios, mixing of the NCG and the water at 1800 m (C-Inj2a) and 1000 m (C-Inj2b) in the well (Fig. 11a). As described previously, mixing the NCG with water at 1000 m requires the water column to be previously rised above the level of the valve by injecting the NCG from valves located deeper in the well. Following the water-only injection, which results in an initial rise of the water column, initiation of the mixing of the two phases from the deeper point (1800 m) leads to a further rise towards the surface of the water column (Fig. 11b). Similarly, if the injection takes place from points that are located shallower in the well, further increase in the water level is observed. This corresponds to a transient simulation, but only the final steady-state results of

<figure><img src="image_13.png" /><figcaption>Fig. 12. Comparison of (a) pressure, (b) temperature, (c) void fraction, and (d) flow regimes (I: intermittent, D: dispersed, L: liquid) between GWellFM, with the Carnot and Multiflash options, PIPESIM, with the black-oil and the compositional options, and OLGA, with the black-oil, for H₂O-NCG surface injection (K-Inj1).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Comparison of (a, d) pressure, (b, e) void fraction, and (c, f) average velocity between UniSim® and GWellFM, with the Carnot and the Multiflash options, for H₂O-NCG deep mixing: (a-c): C-Inj2a, (d-f): C-Inj2b.</figcaption></figure>

injecting from the valve located at 1000 m are shown here so that the results can be compared with the corresponding results of the steady-state GWellFM simulator.

The impact of the well configuration and the mixing scenario can be seen on the profile of the mixture average velocity (Fig. 11d) and the flow configuration (Fig. 11e). Both in the upper part of the well and the liquid dominant area, the velocity is the same regardless of the injection scenario. In the lower part of the well, the reduction of the velocity due to the diameter change can be seen. When mixing of NCGs and H₂O takes place, part of the gas dissolves in the liquid phase. The dissolution increases downstream due to the more favourable dissolution local conditions and both the void fraction, and the mixture velocity reduce since less gas is present. After the deeper mixing point, the pressure is the same for both injection scenarios and therefore the behaviour of the two-phase mixture is also the same. In this part of the well, the velocity of the mixture is higher than the liquid velocity in the water-only injection case due to the presence of the extra injected amount of NCG fluid.

On the other hand, the impact of the flow regime modification is

visible on the temperature profiles (Fig. 11c). At the upper part of the well, when the flow passes from two-phase to liquid-only there is a sudden increase in fluid temperature which corresponds to the exothermic nature of the condensation. A second temperature jump is observed when the NCG is mixed with the liquid, which is caused by the enthalpy released due to the dissolution of the gas, similar to the observations with PIPESIM (K-Inj1).

### 4.4. Comparison of modelling tools

In addition to the comparisons between the results of the commercial tools and the in-house code for the water-only injection, simulations were performed for the mixing of water and NCGs. OLGA simulator with the black-oil option was also used to model one case from Fig. 9 ($m_L = 200$ t/h, $m_G = 5$ t/h) for the injection of a mixture of water and NCG from the surface (K-Inj1) and compared with PIPESIM and GWellFM (Fig. 12). Again, the behaviour of the mixture is different between the compositional and black-oil flow models.

The models with the compositional option predict that at one point
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** (a) Temperature profiles and (b) flow regimes calculated by GWellFM, with the Carnot and the Multiflash options, for deep mixing of H₂O-NCG (C-Inj2).</figcaption></figure>

in the well the gas is fully dissolved in the liquid, whereas the black-oil models estimate that the mixture is always under two-phase conditions. The complete dissolution of the gas can be also seen in the temperature profiles, where an increase in the local temperature is observed. The differences in the void fraction between the compositional GWellFM and PIPESIM can be due to either the thermodynamic calculations or the two-phase closure laws. On the other hand, the differences between Carnot and Multiphase are only due to the phase composition and properties calculations. As for the water-only injection (K-Inj0), GWellFM calculates lower pumping pressure when compared to PIPE-SIM, which can be attributed to the lower estimated void fraction which results in higher hydrostatic pressure and therefore leads to lower surface pressure. The pressure and the temperature profiles between the calculations of GWellFM with Carnot and Multiflash almost overlap, with an average error of 1 % and 3.26 %, respectively.

GWellFM was also used to simulate the two mixing scenarios of Castelnuovo (C-Inj2a and C-Inj2b). Pressure, void fraction, and average mixture velocity are compared in Fig. 13. At the upper part of the well, the results are the same as the water-only injection (C-Inj0). Generally, GWellFM always predicts a higher water level, and a lower void fraction and two-phase velocity when compared to UniSim®. On the other hand, the two thermodynamic options of GWellFM gave similar results, with the only difference occurring in the lower part of the well where the diameter of the well expands with the transition from casing zone to open hole zone.

In this section of the well, the Multiflash option predicts an increase in the void fraction and the dissolution of less NCGs. The calculations using Carnot model predict a flow regime that become dispersed when reaching the deepest part of the well (Fig. 14b). In the simulations, with Multiflash thermodynamic model, the temperature deep in the well is always lower when compared to simulations using Carnot thermodynamic model (Fig. 14a), which suggests that the dissolution of NCG reduces with the temperature. In addition to the temperature jumps observed with UniSim®, and reported previously, GWellFM also predicts a sharper increase of the temperature at the end of the coaxial completion (at 2000 m), where the fluid exchanges heat directly with the surrounding formation and a significant modification of the overall heat transfer coefficient occurs.

GWellFM also calculates the evolution of the pressure in the annulus space of the well completion (Fig. 13a, d), and the required wellhead gas pressure to inject the NCGs. The necessary gas compression depends on the mixing scenario; the shallower the mixing, the lower the NCG injection pressure. In fact, when the injection takes place from the deeper point in the well, the hydraulic head of the gas in the annulus must overcome the tubing pressure to pass from the annular space to the central tubing; around 600 m for the C-Inj2b instead of 800 m for C-Inj2a. This finding confirms the idea of the deep mixing of the NCGs and water through multiple injection points and suggests that the exact number and the locations of the valves can be further optimised to improve the operability and enhance the injection conditions.

To better understand the origin of the discrepancies in the pressure calculations between the two codes, the separated pressure losses (due to friction and gravity) are compared in Fig. 15 for the C-Inj2a deep mixing case. The results are similar for the C-Inj2b case. In the upper part of the well (H₂O + 0.1 % CO₂) both UniSim® and GWellFM give

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Comparison of (a) frictional, and (b) gravitational pressures losses between UniSim® and GWellFM, with the Carnot and the Multiflash options, for H₂O-NCG deep mixing (C-Inj2a).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Comparison of gas (G) and liquid (L) mixture compositions between UniSim® and GWellFM, with the Carnot and the Multiflash options, for H₂O-NCG deep mixing (C-Inj2a) at the upper (I) and lower (II) part of the well: (a) H₂O, (b) CO₂ and (c) H₂S (attention to the different scales).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Evolution of components concentration in the gas phase of H₂O-CO₂ mixture in the upper part of the well (GWellFM with Carnot, C-Inj2a).</figcaption></figure>

different values because of the differences in the velocity (or mixture density) and in the two-phase flow models, whereas in the lower part of the well (H₂O + NCG) the friction losses differ a lot which can be only attributed to the two-phase hydrodynamic models. In the liquid dominant regime, all models calculate similar losses.

Similarly, the impact of the thermodynamic models is evident in the mixture composition. Fig. 16 includes the gas and the liquid composition of each component for the parts of the well where a two-phase mixture is flowing. Due to the differences in the calculated pressure and temperature and the specificities of the thermodynamic models, the phase composition is not the same. All thermodynamic options show a peak of

the CO₂ concentration in the gas phase in the upper part of the well just before the liquid surface. It must be noted that the sum of all components' composition is always 100 %, thus this increase is followed by a decrease in the H₂O concentration, suggesting that the condensation of water is happening faster than CO₂ (Fig. 17). This peak was present in all injection scenarios (C-Inj0, C-Inj2a and C-Inj2b) and in the simulations with all tools.

GWellFM give access to all physical properties calculated with both thermodynamic options. Table 9 summarises the average error in the properties calculated by Multiflash and Carnot options for different mixtures and conditions. Important discrepancies are observed mainly for the liquid heat capacity, even for the water-only case. The Carnot option constantly underestimates the corresponding values. For the gas phase, even the presence of a small quantity of CO₂ mixed with water results in divergence in the estimation of the viscosity, and for high concentration of NCGs high differences are observed also for the thermal conductivity of the gas phase, the viscosity of the liquid phase and the surface tension.

A dimensional analysis is performed for the heat transfer phenomenon to correlate the overall heat transfer coefficient, *U*, with the heat capacity of the fluid (see Appendix 1). Indeed, the convective heat transfer coefficient, *h*, is proportional to *Cp* raised by a positive exponent, *n*₃. Then, *U* which is equal to the inverse sum of all thermal resistances, *R*₀, can be expressed by Eq. (3). The parameter *n*₂ includes all other physical properties. For constant heat conduction, *n*₁, and constant all other properties, *n*₂, the overall *U* value is lower for lower *Cp* values. Since Carnot underestimates the *Cp*, the estimated *U* is also lower when compared with the Multiflash value, which means that the fluid can absorb more heat from the hot surroundings and, thus, the fluid temperature is higher with the Carnot option, as presented previously. Similar impact is expected of the underestimation by Carnot of the thermal conductivity for the gas phase under two-phase conditions. Looking in Fig. 14, the fluid temperature with Carnot is always higher
<caption>**Table 9**</caption>

Average error, δ (%), in phase properties (density, viscosity, heat capacity, thermal conductivity, and surface tension) as calculated by GWellFM with the Multiflash and the Carnot options.

<table><thead><tr><th>Case label</th><th>Components</th><th>Phases</th><th>Phase</th><th>δρ</th><th>δμ</th><th>δCp</th><th>δλ</th><th>δσ</th></tr></thead><tbody><tr><td>K-Inj0</td><td>H<sub>2</sub>O</td><td>1</td><td>Liquid</td><td>-1.596</td><td>0.317</td><td><strong>-7.990</strong></td><td>-1.487</td><td>-</td></tr><tr><td rowspan="2">K-Inj1</td><td rowspan="2">H<sub>2</sub>O, NCG</td><td rowspan="2">2</td><td>Liquid</td><td>0.626</td><td>-0.253</td><td><strong>-11.100</strong></td><td>2.691</td><td rowspan="2">-0.580</td></tr><tr><td>Gas</td><td>0.462</td><td><strong>4.277</strong></td><td>-0.568</td><td><strong>-14.110</strong></td></tr><tr><td rowspan="2">C-Inj0</td><td rowspan="2">H<sub>2</sub>O, CO<sub>2</sub> (0.1 %)</td><td rowspan="2">2</td><td>Liquid</td><td>-1.359</td><td>0.108</td><td><strong>-7.809</strong></td><td>-1.299</td><td rowspan="2">2.440</td></tr><tr><td>Gas</td><td>-0.900</td><td><strong>-8.389</strong></td><td>-0.442</td><td>-0.962</td></tr><tr><td rowspan="2">C-Inj2a</td><td rowspan="2">H<sub>2</sub>O, CO<sub>2</sub>, H<sub>2</sub>S</td><td rowspan="2">1</td><td>Liquid</td><td>-1.458</td><td>-0.190</td><td><strong>-7.830</strong></td><td>-1.257</td><td>-</td></tr><tr><td>Gas</td><td>-1.465</td><td><strong>5.377</strong></td><td>-8.307</td><td>0.733</td><td><strong>7.860</strong></td></tr></tbody></table>

than the one calculated with Multiflash.

$$U \sim \frac{1}{R_{\theta,cond} + R_{\theta,conv}} \sim \frac{1}{n_1 + \frac{1}{h}} \sim \frac{1}{n_1 + \frac{1}{n_2 C_p^{\theta_3}}} \quad (3)$$

## 5. Conclusions

The present study sheds light on the feasibility of a novel approach for the simultaneous reinjection of water and NCGs in geothermal wells. While previous approaches examined the complete dissolution of CO₂ or NCGs in water and the single-phase reinjection, the solution proposed here benefits from the two-phase condition while the NCG mixture is sequentially introduced into the water column in multiple depths. The key drivers of complexity in this problem are the thermodynamics of the NCG mixture, which is near the critical condition, the uncertainty of gas mixtures solubility in water at high pressures, the hydrodynamics of downward two-phase flows in large diameter wells, which have not been analysed and characterised in detail, and the transient flow conditions. Due to a lack of experimental and field data, a set of numerical approaches were used to evaluate different reinjection scenarios. The tools set included UniSim®, OLGA, and PIPESIM, as well as an in-house code (GWellFM), which were all used to simulate different injection scenarios for two geothermal fields (see Table 1).

The main findings of the present study are:

a) Deep mixing helps fully entrain the gas mixture in the liquid with minimum work on the surface, driving the two-phase mixture down into the reservoir.

b) Although, for high-pressure reservoirs, deep mixing (K-Inj2) needs more NCG compression power on the surface than surface mixing (K-Inj1), it requires less pumping demand.

c) For specific geothermal conditions (as in the 2nd case study) with low reservoir pressure, limited water injection rate, or high NCG injection rate, the deep mixing is the only feasible solution.

d) Deep mixing through valves at different depths (C-Inj2b) is more technically promising compared to the single mixing point deep in the well (C-Inj2a). It was also found that the fluid temperature significantly impacts the pressure demand, and lowering the temperature leads to reducing the pressure required for both surface and deep mixing scenarios.

e) Flow patterns and regime transitioning along the well depend on liquid and gas flow rates. Higher flowrates require more pressure, which contributes to higher solubility; however, higher gas flow rates also result in a higher void fraction that offsets the benefits of the water column static head. There is, therefore, a need for a holistic model that optimises the reinjection process by changing the phase flow rates, pressure, solubility and mixing depths.

f) A wide range of models for solubility and density, as well as mixing rules, are utilised by the available tools. Despite the similarity and

the overall pattern of the dynamic behaviour, the deviations of the property profiles, and compositions along the well, are quite pronounced in the two-phase region. Particularly, the black-oil approach by PIPESIM and OLGA was found to not be flexible enough to consider both single-phase and two-phase properties in a single simulation. It was shown that a negligible amount of CO₂ (0.1 %wt) in water gave different results from UniSim® and GWellFM.

g) Existing pipe flow models are highly dependent on the phase fraction; therefore, the uncertainty of the phase fraction and physical properties will be projected on the flow model for downward flow on the reinjection well.

h) The effect of the thermodynamic and hydrodynamic models becomes more evident when the in-house code is compared to the commercial tools. A good agreement was found between UniSim® and the developed in-house code for property profiles, making them reliable to simulate the behaviour of the downward two-phase flows of water and NCG mixtures.

i) The discrepancies between the simulation tools are associated with: (i) the thermodynamic model responsible for phase interaction, particularly water and the NCG, (ii) the pressure loss model dedicated to the two-phase region, (iii) the void fraction model, (iv) the problem definition in terms of geometry and numerical approach implementation.

j) The difference and innovation of this work when compared to previous works, with the exception of the reinjection technique of deep mixing, is the use of two-phase modelling of NCGs and water instead of single-phase flows of dissolved CO₂ in water.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

The present research was funded by the European Union's Horizon 2020 Research and Innovation Program under grant agreement No 818169 (GECO Project).

The authors would like to acknowledge Zorlu Enerji and Graziella Green Power S.p.A, industrial partners of the GECO project and operators of the two geothermal sites, for providing the permission to use the field data for the simulations.
# Appendix 1

## Governing equations for transient flows in a geothermal wellbore

There are two general approaches for writing and solving the conservation equations: the two-fluid, or more generally the multi-fluid approach and the mixture formulation [42]. The mixture representation has been chosen here, and the conservation of mass, momentum, and energy in a 1D pipe (vertical direction) are governed by the below equations:

$$A \frac{\partial \rho_m}{\partial t} + \frac{\partial}{\partial z} (A \rho_m u_m) = 0 \qquad (4)$$

$$\frac{\partial}{\partial t}(\rho_m u_m) + \frac{1}{A} \frac{\partial}{\partial z} (A \rho_m u_m^2) = -\frac{\partial P}{\partial z} - \left(\frac{\partial P}{\partial z}\right)_f - \rho_m g \cos\theta \qquad (5)$$

$$\frac{\partial}{\partial t}(\rho_m \hat{h}) + \frac{1}{A} \frac{\partial}{\partial z} (A u_m \rho_m \hat{h}) = -\rho_m u_m g \cos\theta - \frac{q}{A} \qquad (6)$$

where $\rho_m$ is the mixture density, $u_m$ is the mixture velocity, $A$ is the cross-sectional area, $P$ is the pressure, $g$ the gravitational acceleration, $\theta$ is the inclination of the well, $\hat{h}$ is the specific enthalpy and $q$ are the heat exchange through the casing. Steady-state simulators do not consider the term of the temporal evolution ($\partial/\partial t$).

For single phase flows the pressure losses due to friction are [43]:

$$\left(\frac{\partial P}{\partial z}\right)_f = \frac{f \rho_k u_k |u_k|}{2A} \Delta z \qquad (7)$$

where $z$ is the depth, $k$ is the phase (gas or liquid) and $f$ is the Darcy friction factor calculated from the fluid Reynolds number, $Re$, the available diameter, $d$, and the wall roughness, $\epsilon$:

$$\left\{ \begin{array}{l} f = \frac{64}{Re}, \quad Re \le 2400 \\ \frac{1}{\sqrt{f}} = -2.0 \log \left( \frac{\epsilon}{3.7d} + \frac{2.51}{Re\sqrt{f}} \right), \quad Re > 2400 \end{array} \right. \qquad (8)$$

The concept of the thermal resistances, $R_\theta$, in series and the overall heat transfer coefficient, $U$, is applied for calculating the rate of heat thermal [44]:

$$q = \pi d U \Delta T \qquad (9)$$

where $\Delta T$ is the temperature difference between the fluid and a point of interest. In solid materials, the conductance transfer is expressed with the thermal conductivity, $\lambda$, of the material, whereas the convection, natural or forced, with the convective heat transfer coefficient, $h_b$, derived from Nusselt number correlations. Then, the overall heat transfer is calculated as:

$$\frac{1}{U} = \sum_{k} R_{\theta, \text{conv}} + \sum_{m} R_{\theta, \text{cond}} = \sum_{k} \frac{1}{h_{t,k}} + \sum_{m} \frac{d \ln \left( \frac{\tau_m}{\tau_{m-1}} \right)}{2\lambda} \qquad (10)$$

where $k$ is the phases (gas or liquid) and $r_m$ is the inner diameter of the m solid material (casing, cement, rocks, etc.).

# Appendix 2

## Two-phase flow equations

The two-phase mixture density, $\rho_m$, is calculated from the phase densities and the void fraction, $\alpha_G$, according to the mixture velocity, $u_m$, as the sum of the phase velocities. On the other hand, the average mixture velocity, $u_{av}$, is based on the total mass flow and the total flowing cross-section of the well [42].

$$\rho_m = \alpha_G \rho_G + (1 - \alpha_G) \rho_L \qquad (11)$$

$$u_m = u_{sG} + u_{sL} = \frac{\dot{m}_G}{\rho_G A} + \frac{\dot{m}_L}{\rho_L A} \qquad (12)$$

$$u_{av} = \frac{4 \dot{m}_{av}}{\pi d^2 \rho_m} \qquad (13)$$

where $u_{sG}$ and $u_{sL}$ are gas and liquid superficial velocities.

The two-phase models applied here propose equations for the gas void fraction and for calculating the pressure losses due to friction. For the latter, either a two-phase friction factor or an empirical two-phase pressure drop equation is proposed (see Table 2), which are also dependent on the fluid physical properties, the geometry and the wall roughness, like Eq. (8) for the single phase flow. However, these equations (or closure laws) are more complicated. For the complete two-phase equations used in each code the reader is referred to the corresponding sources.
# Appendix 3

## Property models

Many thermodynamics investigations adapted either empirical correlations or EoS compositional models for property modelling, and both approaches are employed in the present study, including the industry-standard black-oil correlation and cubic EoS compositional models. Selecting between correlations and EoS as the property model is a trade-off between accuracy and computational cost.

The black-oil model provides a wide range of viscosity correlations and fluid mixing rules developed for reservoir conditions. This model, which is commercially available in PIPESIM and OLGA, is applicable for oil and condensate, can also be used for simplified gases; however, this simplification may lead to significant errors, mainly in the gas phase composition.

The sour-PR, available in UniSim®, is based on the classical Peng-Robinson:

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

where the pressure of the fluid is a function of temperature, T, and molar volume, v, and parameters a and b are computed from van der Waals mixing rules [45]. However, it accounts for the ion balance for calculating phase composition, particularly the CO₂ and H₂S fraction in the aqueous phase, which increases the computational cost as well. The improvement is derived from the VLE data regression and the modified dissociation equilibrium constants assessed by Wilson [46]. The Wilson improvement is based on component partial pressure assessed by Henry's law, Eqs (15) & (16), where the Henry constant, H, is related to the temperature and undissociated CO₂ and H₂S in the liquid phase.

$$P_{\text{partial,CO}_2} = H_{\text{CO}_2} \cdot C_{\text{CO}_2} \qquad (15)$$

$$P_{\text{partial,H}_2S} = H_{\text{H}_2S} \cdot C_{\text{H}_2S} \qquad (16)$$

where $C_i$ is the liquid phase concentration of the components.

The other EoS is the CPA model which is developed based on the cubic equation of state, adapting the association terms. This model includes three terms [45]:

$$P(T, v) = \frac{RT}{v-b} - \frac{a(T)}{v(v+b)} - \frac{1}{2} \frac{RT}{v} \left(1 + \frac{1}{v} \frac{\partial \ln \beta}{\partial(1/v)}\right) \sum_{i} x_{i} \sum_{A_{i}} (1 - X_{A_{i}}) \qquad (17)$$

in which the first two terms are taken from the classic cubic models for the physical interaction, and the third accounts for the association interactions. Where $x_i$, $X_{A_i}$ and $\beta$ are the component mole fraction, the site monomer fraction and the radial distribution function.

## References

[1] A. Baldacci, M. Mannari, F. Sansone, Greening of geothermal power: An innovative technology for abatement of hydrogen sulphide and mercury emission, in: World Geothermal Congress 2005, Antalya, Turkey, 2005.

[2] A. Sbrana, P. Marianielli, M. Belgiorno, M. Sbrana, V. Ciani, Natural CO₂ degassing in the Mount Amiata volcanic-geothermal area, J. Volcanol. Geotherm. Res. 397 (2020), 106852, https://doi.org/10.1016/j.jvolgeores.2020.106852.

[3] S. Solomon, T. Flach, Carbon dioxide (CO₂) injection processes and technology, in: M.-V.M. Mercedes (Ed.), Developments and Innovation in Carbon Dioxide (CO₂) Capture and Storage Technology, Woodhead Publishing, 2010, pp. 435–466. http://doi.org/10.1533/9781845699574.4.435.

[4] A. Rivera Diaz, E. Kaya, S.J. Zarrouk, Reinjection in geothermal fields – a worldwide review update, Renewable Sustainable Energy Rev. 53 (2016) 105–162, https://doi.org/10.1016/j.rser.2015.07.151.

[5] E. Kaya, S.J. Zarrouk, M.J. O'Sullivan, Reinjection in geothermal fields: a review of worldwide experience, Renewable Sustainable Energy Rev. 15 (2011) 47–68, https://doi.org/10.1016/j.rser.2010.07.032.

[6] S.K. Sanyal, E.E. Granados, A.J. Menzies, Injection-related problems encountered in geothermal Projects and their mitigation: The United States experience, in: World Geothermal Congress 1995, Florence, Italy, 1995.

[7] E. Kaya, S.J. Zarrouk, Reinjection of greenhouse gases into geothermal reservoirs, Int. J. Greenhouse Gas Control 67 (2017) 111–129, https://doi.org/10.1016/j.ijggc.2017.10.015.

[8] M. Burton, S.L. Bryant, Surface dissolution: minimizing groundwater impact and leakage risk simultaneously, Energy Procedia 1 (2009) 3707–3714, https://doi.org/10.1016/j.egypro.2009.02.169.

[9] P.E. Eke, M. Naylor, S. Haszeldine, A. Curtis, CO₂/brine surface dissolution and injection: CO₂ storage enhancement, SPE Projects, Facilities & Construction 6 (2011) 41–53, https://doi.org/10.2118/124711-PA.

[10] H. Koide, Z. Xue, Carbon microbubbles sequestration: a novel technology for stable underground emplacement of greenhouse gases into wide variety of saline aquifers, fractured rocks and tight reservoirs, Energy Procedia 1 (2009) 3655–3662, https://doi.org/10.1016/j.egypro.2009.02.162.

[11] T.D. Rathnaweera, P.G. Ranjith, M.S.A. Perera, A. Haque, Influence of CO₂-brine co-injection on CO₂ storage capacity enhancement in deep saline aquifers: an experimental study on Hawkesbury Sandstone formation, Energy Fuels 30 (2016) 4229–4243, https://doi.org/10.1021/acs.energyfuels.6b00113.

[12] S.M. Shariatipour, E.J. Mackay, G.E. Pickup, An engineering solution for CO₂ injection in saline aquifers, Int. J. Greenhouse Gas Control 53 (2016) 98–105, https://doi.org/10.1016/j.ijggc.2016.06.006.

[13] K. Suzuki, H. Miida, H. Wada, S. Horikawa, T. Ebi, K. Inaba, Feasibility study on CO₂ micro-bubble storage (CMS), Energy Procedia 37 (2013) 6002–6009, https://doi.org/10.1016/j.egypro.2013.06.528.

[14] S. Zendehboudi, A. Khan, S. Carlisle, Y. Leonenko, Ex situ dissolution of CO₂: a new engineering methodology based on mass-transfer perspective for enhancement of CO₂ sequestration, Energy Fuels 25 (2011) 3323–3333, https://doi.org/10.1021/ef200199r.

[15] M. Zirrahi, H. Hassanzadeh, J. Abedi, The laboratory testing and scale-up of a downhole device for CO₂ dissolution acceleration, Int. J. Greenhouse Gas Control 16 (2013) 41–49, https://doi.org/10.1016/j.ijggc.2013.02.020.

[16] C. Kervévan, M.-H. Beddelem, K. O'Neil, CO₂-DISSOLVED: A novel concept coupling geological storage of dissolved CO₂ and geothermal heat recovery – Part 1: Assessment of the integration of an innovative low-cost, water- based CO₂ capture technology, Energy Procedia 63 (2014) 4508–4518, https://doi.org/10.1016/j.egypro.2014.11.485.

[17] B. Sigfusson, S.R. Gislason, J.M. Matter, M. Stute, E. Gunnlaugsson, I. Gunnarsson, E.S. Aradottir, H. Sigurdardóttir, K. Mesfin, H.A. Alfredsson, D. Wolff-Boenisch, M. T. Arnarsson, E.H. Oelkers, Solving the carbon-dioxide buoyancy challenge: the design and field testing of a dissolved CO₂ injection system, Int. J. Greenhouse Gas Control 37 (2015) 213–219, https://doi.org/10.1016/j.ijggc.2015.02.022.

[18] M.J. Shafaei, J. Abedi, H. Hassanzadeh, Z. Chen, Reverse gas-lift technology for CO₂ storage into deep saline aquifers, Energy 45 (2012) 840–849, https://doi.org/10.1016/j.energy.2012.07.007.

[19] R.W. Stacey, L. Norris, S. Lisi, OLGA modeling results for single well reinjection of non-condensable gases (NCGs) and water, GRC Transactions 40 (2016) 931–939.

[20] V. Leontidis, M. Gainville, L. Jeannin, M. Perreaux, S. Souque, Modelling of the non-condensable gases re-injection for geothermal emission control (GECO Project), in: World Geothermal Congress 2020+1, Reykjavik, Iceland, 2021.

[21] I. Yücetas, N. Ergiçay, S. Akin, Carbon dioxide injection field pilot in Umurlu geothermal field, Turkey, GRC Trans. 42 (2018).

[22] Z. Duan, R. Sun, An improved model calculating CO₂ solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar, Chem. Geol. 193 (2003) 257–271. https://doi.org/10.1016/S0009-2541(02)00263-2.

[23] E.S.P. Aradóttir, I. Gunnarsson, B. Sigfússon, G. Gunnarsson, B.M. Júliusson, E. Gunnlaugsson, H. Sigurdardóttir, M.T. Arnarson, E. Sonnenthal, Toward cleaner geothermal energy utilization: capturing and sequestering CO₂ and H₂S emissions
from geothermal power plants, Transp. Porous Media 108 (2015) 61–84, https://doi.org/10.1007/s11242-014-0316-5.

[24] I. Gunnarsson, E.S. Aradóttir, E.H. Oelkers, D.E. Clark, M.P. Arnarson, B. Sigfússon, S.Ö. Snæbjörnsdóttir, J.M. Matter, M. Stute, B.M. Júlfusson, S.R. Gíslason, The rapid and cost-effective capture and subsurface mineral storage of carbon and sulfur at the CarbFix2 site, Int. J. Greenhouse Gas Control 79 (2018) 117–126, https://doi.org/10.1016/j.ijggc.2018.08.014.

[25] A. Morin, M. Hjelstuen, B.T. Lovfall, E. Meese, Multiphase flow simulation of geothermal wells, in: 3rd Hydrocarbon Geothermal Cross Over Technology Workshop, Geneva, Switzerland, 2019.

[26] Z.P. Aunzo, G. Bjornsson, G.S. Bodvarsson, Wellbore Models GWELL, GWNACL, and HOLA User's Guide, California, USA.

[27] A. Abouie, A. Kazemi Nia Korrani, M. Shirdel, K. Sepehrnoori, Comprehensive modeling of scale deposition using a coupled geochemical and compositional wellbore simulator, in: Offshore Technology Conference, Houston, Texas, 2016. https://doi.org/10.4043/27072-MS.

[28] R.A. Tonkin, M.J. O'Sullivan, J.P. O'Sullivan, A review of mathematical models for geothermal wellbore simulation, Geothermics 97 (2021), 102255, https://doi.org/10.1016/j.geothermics.2021.102255.

[29] L. Cheng, Frontiers and Progress in Multiphase Flow I, in: L. Cheng (Ed.), Frontiers and Progress in Multiphase Flow I, Springer Cham, 2014, pp. 157–212. https://doi.org/10.1007/978-3-319-04358-6.

[30] P.H. Niknam, L. Talluri, D. Fiaschi, G. Manfrida, Improved solubility model for pure gas and binary mixture of CO₂-H₂S in water: a geothermal case study with total reinjection, Energies 13 (2020) 2883, https://doi.org/10.3390/en13112883.

[31] S. Erol, T. Akın, A. Başer, Ö. Saraçoğlu, S. Akın, Fluid-CO₂ injection impact in a geothermal reservoir: evaluation with 3-D reactive transport modeling, Geothermics 98 (2022), 102271, https://doi.org/10.1016/j.geothermics.2021.102271.

[32] P.H. Niknam, L. Talluri, D. Fiaschi, G. Manfrida, Sensitivity analysis and dynamic modelling of the reinjection process in a binary cycle geothermal power plant of Larderello area, Energy 214 (2021), 118869, https://doi.org/10.1016/j.energy.2020.118869.

[33] M. Vaccaro, F. Batini, M. Stolzuoli, S. Bianchi, R. Pizzoli, S. Lisi, Geothermal ORC plant case study in Italy: Castelnuovo Val di Cecina – Design and technical aspects, in: European Geothermal Congress 2016, Strasbourg, France, 2016.

[34] V. Leontidis, M. Gainville, Multiphase compositional well flow modelling: The case of injecting in low pressure reservoirs, in: 1st Geoscience and Engineering in Energy

Transition Conference, GET 2020, Strasbourg, France, 2020, pp. 1–5. https://doi.org/10.3997/2214-4609.202021008.

[35] A.R. Hagedorn, K.E. Brown, Experimental study of pressure gradients occurring during continuous two-phase flow in small-diameter vertical conduits, J. Petrol. Tech. 17 (1965) 475–484, https://doi.org/10.2118/940-PA.

[36] J. Orkiszewski, Predicting two-phase pressure drops in vertical pipe, J. Petrol. Tech. 19 (1967) 829–838, https://doi.org/10.2118/1546-PA.

[37] D.H. Beggs, J.P. Brill, A study of two-phase flow in inclined pipes, J. Petrol. Tech. 25 (1973) 607–617, https://doi.org/10.2118/4007-PA.

[38] S.M. Bhagwat, A.J. Ghajar, A flow pattern independent drift flux model based void fraction correlation for a wide range of gas-liquid two phase flow, Int. J. Multiphase Flow 59 (2014) 186–205, https://doi.org/10.1016/j.ijmultiphaseflow.2013.11.001.

[39] C. Yao, H. Li, Y. Xue, X. Liu, C. Hao, Investigation on the frictional pressure drop of gas liquid two-phase flows in vertical downward tubes, Int. Commun. Heat Mass Transfer 91 (2018) 138–149, https://doi.org/10.1016/j.icheatmastransfer.2017.11.015.

[40] M. Lokanathan, T. Hibiki, Flow regime transition criteria for co-current downward two-phase flow, Prog. Nucl. Energy 103 (2018) 165–175, https://doi.org/10.1016/j.pnucene.2017.11.021.

[41] X. Courtial, N. Ferrando, J.-C. de Hemptinne, P. Mougin, Electrolyte CPA equation of state for very high temperature and pressure reservoir and basin applications, Geochim. Cosmochim. Acta 142 (2014) 1–14, https://doi.org/10.1016/j.gca.2014.07.028.

[42] G. Yadigaroglu, G. Hewitt (Eds.), Introduction to Multiphase Flow: Basic Concepts, Applications and Modelling, first ed., Springer International Publishing, Cham, 2018.

[43] F.M. White, Fluid Mechanics, seventh ed., McGraw-Hill, New York, N.Y., 2010.

[44] M. Massoud, Engineering Thermofluids: Thermodynamics, Fluid Mechanics, and Heat Transfer, Springer-Verlag, Berlin Heidelberg, Berlin, Heidelberg, 2005.

[45] G.M. Kontogeorgis, M.L. Michelsen, G.K. Folas, S. Derawi, N. von Solms, E. H. Stenby, Ten years with the CPA (cubic-plus-association) equation of state. Part 1. Pure compounds and self-associating systems, Ind. Eng. Chem. Res. 45 (2006) 4855–4868, https://doi.org/10.1021/ie051305v.

[46] G.M. Wilson, A new correlation of NH₃, CO₂, and H₂S volatility data from aqueous sour water systems, American Petroleum Institute, Washington, USA, 1978.