Introduction
Climate change is a global phenomenon strongly correlated with water, food, and energy changes. Thus, the energy-climate nexus is a key challenge for the sustainable development achievement4. Due to the increase of industrialisation, an increase in the released CO2 emissions in the atmosphere was observed, enhancing the global warming phenomenon as the main component of greenhouse gasses5. However, the continuous population development and economic growth require higher energy demands. Comparing with previous decades, this demand is slower increasing with an average of 0.7 % per year through 2050 than a percentage higher than 2 % average from 2000 to 20156. Thus, the European strategies aim for the energy policies’ rapid changes to achieve efficient energy demands combined with greenhouse gas reduction. The 2015 Paris Agreement proposed the maintenance of earth’s temperature at 1.5–2.0 °C and the mitigation of climate change consequences7. Moreover, the EU Energy Roadmap 2050 planed a reduction of greenhouse gas by 80–95 % by 2050 compared to the 1990 base. Hence, proper management of CO2 emissions aiming to the reduction of the quantities in the atmosphere and their use for energetical purposes is crucial for the sequence of the lines of the EU.
To achieve the CO2 emissions reduction, the CO2 capture and storage (CCS) technologies were developed. Two main categories were occurred; the storage in above ground facilities and in underground geological environments. The underground facilities are more favourable as they can store a higher amount of CO2, have a limited footprint as require lower land use for the aboveground installations, they are significantly less influenced to external factors and they are more safe due to the tightness of the caverns as well as the considerable distance from the biosphere and hydrosphere8. The CO2 underground storage was developed after understanding the different geological formations potential to contain tremendous amounts of hydrocarbons for million of years, which enhanced by the remarkable tightness and rock integrity of such formations. Potential geological environments for the carbon storage are the porous media and the hard rock caverns. The porous media are subdivided in depleted hydrocarbon reservoirs and in saline aquifers. The hard rock caverns separated in salt caverns and lined or unlined rock caverns9.
CO2 based electrothermal energy and geological storage system, known as CEEGS2, is one of novel technologies aiming to combine renewable energy storage systems, the transcritical CO2 cycle, the geothermal heat extraction and partial CO2 storage in geological formations. The concept is to integrate the transcritical CO2 cycles with underground energy storage through simultaneous CO2 geological storage and geothermal heat extraction. By this way, the advantages of Carbon Capture, Utilisation and Storage (CCUS) and renewable energy storage technologies with respect to efficiency and profit will be enhanced with a simultaneous minimal environmental impact.
Depleted hydrocarbon reservoirs are of high interest for CO2 storage for both scientific and industrial community compared to the other categories in terms of CO2 storage capacity, experience of reservoir characterisation, already established oil or gas well infrastructure, sealing performance as well as storage operability10–14. The CO2 storage in depleted hydrocarbon fields supposed to be one of the most realistic ways to reduce carbon emissions while it is the most economical-affordable method15. The known globally storage capacity of the available depleted hydrocarbon fields (DHF) is estimated to be close to 390–750 Gigatons which is approximately ten times the current annual CO2 emissions worldwide13. Despite the advantages of DHF, critical disadvantages influence the efficiency of CO2 storage and must be addressed prior its injection. More specifically, the problems can be summarised in; i. problems related to the evaluation method of storage potential and its applicability, ii. potential leakage of CO2 and description of its mechanism, and iii. the interactions between the injected CO2 and the remaining fluids10.
The purpose of this study is to investigate the potential interactions between the injected CO2 in a model DHF by a geochemical computer simulation method. As a model DHF, the depleted gas field of the Marismas fields in South Spain was selected, supposing that represents a range of depleted gas reservoirs with host rocks of carbonate-siliciclastic nature. Firstly, the existing geochemical interactions taken place in the Marismas field were specified, describing the host rock, remaining fluids (i.e. CH4), and remaining water interactions. Secondly, the interactions of the complex Marismas field system and the injected CO2 were examined. Moreover, the thermal-hydraulic-chemical simulation of the depleted gas reservoir using the CMG-GEM compositional simulator3, based on the features of the Marismas, was examined. The primary objective was to assess the technical feasibility of implementing a depleted hydrocarbons reservoir for the CEEGS technology. The models designed based on the needs of CEEGS concept, requiring a supercritical state of the injected CO2, while the results of the study may be helpful for other CO2 storage technologies too. The study represents a first attempt for Marismas field area conceptualization based on data from traditional geological mapping and primary field data, considering the selected DHF as a primary model of a carbonate-siliciclastic reservoir pending further research and data availability.
Methods
Study area
Among the suggested reservoir types for CEEGS technology implementation, the depleted reservoir (gas or oil) holds notable importance. These reservoirs are discarded remnants from the oil and gas industry. Rather than being left abandoned, they can be repurposed for specific applications, particularly within the realm of CEEGS technologies. The peculiarity of many of these reservoirs lies in the remnanent low pressure, which ensures optimal functioning of CEEGS technology. It is crucial that the depleted reservoir pressure does not fall below this critical minimum pressure. This ensures that the injected CO2 can be back-produced to the surface when needed (without requiring energy input that would decrease the efficiency of the system, and reach the wellhead under supercritical conditions.
An additional challenge is posed by the Joule-Thomson effect on the injection well. Low pressures at the reservoir can result in a strong temperature decrease of the CO2 as it expands at the bottomhole. This process can lead to damage to materials of the wells or to clogging of the system due to the formation of CO2 hydrates16,17.
The objective of the numerical model of the depleted hydrocarbon field is to assess the feasibility of this scenario to the CEEGS system from a technical point of view.
For the investigation of geochemical reactions during the CO2 injection into a depleted hydrocarbon reservoir, the Marismas fields were selected as a mock-up target area. Although it is not claimed to store CO2 for energy storage purposes at the Marismas field, some of its characteristics are fairly representative of DHFs where such could be considered. In this study we utilized only data in the public domain, and whenever there was a complete absence of information, expert-based assumptions were made regarded as realistic for the intended aims of the study.
The Marismas fields are operational gas fields belonging to Guadalquivir-Cadiz Miocene Foreland Basin in the South Spain, as given in the Figure 118. They are characterised by high porosity and permeability present in the onshore Guadalquivir and offshore Gadiz Gulf sandy turbiditic formations. The hydrocarbon traps of the area are mostly stratigraphic and the gas generation origin is mainly biogenic, produced by marine sediments and especially clays or shales which are located some meters under the reservoirs.

Figure 1. Study area of Marismas fields, showing the Marismas 3 field which is the selected model of depleted hydrocarbon field.
The Guadalquivir basin generated during Miocene due to the collision between the Iberian and African plates. The basin’s depth is shallow, ranging from 500 m depth in the east basement deepening to 1200 m depth to the west. The area divided in four main stratigraphic sequences19: 1. The basal transgressive calcarenite which is the base of the basin, 2. The Middle and 3. Late Miocene Guadalquivir formation consisted of shallow marine sand units continuing into deep-water turbidites, and 4. The Pliocene formation which overlies unconformably the Miocene unit. The lithology of the basin differentiates based on the geographic location. Marismas fields occur between Huelva and Sevilla provinces with thin sediments composed of calcarenite and siliciclastic deposits19.
The calcarenite composing the reservoir is a mixture of calcite (70 %), iron-rich dolomite (25 %), quartz (5 %) and clay minerals traces with a dominance of illite20. Neogene deposits of the lower-northwestern margin of Guadalquivir Basin are divided into four lithostratigraphic units21: 1. mixed carbonate–siliciclastic deposits of the Niebla Formation (bottom)22,23, 2. greenish–bluish clays of the Arcillas de Gibraleón Formation23, 3. fossiliferous sands and silts of the Arenas de Huelva Formation23, and 4. sands of the Arenas de Bonares Formation24.
As the specific mineralogical composition of the present area is not specified in the literature and the Marismas fields were used as a model depleted hydrocarbon field for the understanding of hydrogen interactions in similar environments, several assumptions were used for the specification of the main mineralogy of the reservoirs. The estimation of main mineralogical phases was performed taking into consideration the semiquantitative mineralogical analysis of the Gulf of Cadiz area, as it was specified by 25. Quartz, alkali-feldspars, illite-mica, kaolinite and Fe-oxides are reported in the area. However, as marine sand units form the Guadalquivir basin and the exact mineralogy of carbonate-siliciclastic deposits of the Niebla is not known, the percentage of carbonates supposed to be slightly higher than quartz. The Gibraleón Formation mainly consists by marls and clays with benthic remains and plankton22. At the base of both Gibraleón and Huelva Sands Formations, there are glauconite-rich horizons, indicating a paleoenvironment with shallow depths at the sediment-water interface. Amounts between 5 to 21 % of glauconite were examined in the areas. The glauconite-rich horizons, except of glauconite, contain also quartz, alkali feldspars and clay minerals as main mineralogical phases, while impurities of ilmenite, hematite, rutile, jarosite, and Fe-rich minerals were also present in the areas26. From the group of clay minerals, illite and kaolinite prevailed in the area, while smectite with its characteristic swelling properties, occurred only as an impurity. The Huelva Formation is rich in quartz, feldspars, phyllosilicates (clay minerals, and muscovite), illite, and chlorite. Mixed-layer clay minerals with the characteristic example of illite-muscovite are also present in the area27. The heavy fraction of sediments samples of the area verified the presence of magnetite, ilmenite, pyrite, hematite, limonite, goethite, andalusite, hornblende, epidote, sillimanite, kyanite, tourmaline, zircon, rutile, apatite, titanite and in smaller appearance granate, and galena. The heavy fraction covers a percentage of 5 to 15 % of the total samples and for this reason, the average of these percentages was assumed, which means the 10 %. Concerning the detection limits of the analytical techniques used for the identification of the minerals that are almost a 3 % of the sample, the minimum concentration of the above minerals specified to 0.3 %. The estimated mineralogical phases of the area are given in the Table 1. Estimated main mineralogical phases and impurities of the studied area with the empirical formula and the database used for the PHREEQC calculations., keeping into consideration the aforementioned knowledge of the percentages, mainly used for the mineralogical impurities. Some of them assumed to be present in higher concentrations related to their abundant presence or in lower due to the lack of their extend27.
Table 1. Estimated main mineralogical phases and impurities of the studied area with the empirical formula and the database used for the PHREEQC calculations.
| Minerals | Main | Impurities | Empirical formula | Estimated quantities % | Database |
|---|
Calcite or/and Aragonite | + | | CaCO3 | 26 | phreeqc.dat |
| Dolomite | + | | CaMg(CO3)2 | 9.2 | phreeqc.dat |
| Quartz | + | | SiO2 | 26 | phreeqc.dat |
| Alkali Feldspars | + | | KAlSi3O8 | 1.6 | phreeqc.dat |
| Muscovite (K-mica) | + | | Kal3Si3O10(OH)2 | 1.8 | phreeqc.dat |
| Illite | + | | K0.6Mg0.25Al2.3Si3.5o10(OH)2 | 3 | phreeqc.dat |
| Kaolinite | + | + | Al2Si2O5(OH)4 | 8.5 | phreeqc.dat |
| Chlorite | + | | Mg5Al2Si3O10(OH)8 | 2 | phreeqc.dat |
| Glauconite | + | | K0.75Mg0.25Fe1.5Al0.50Si3.75o10(OH)2 | 13 | sit.dat |
| Magnetite | | + | Fe3O4 | 0.6 | |
| Ilmenite | | + | Fe2+TiO3 | 1.0 | |
| Hematite | | + | Fe2O3 | 1.0 | |
| Goethite / Limonite | | + | FeOOH | 0.6 | |
| Pyrite | | + | FeS2 | 0.4 | |
| Hornblende | | + | Ca2(Fe2+4Al)(Si7Al)O22(OH)2 | 0.4 | |
| Tourmaline | | + | several | 0.4 | |
| Epidote | | + | (CaCa)(AlAlFe3+)O[Si2o7][SiO4](OH) | 0.4 | |
Sillimanite, Andalusite or Kyanite | | + | Al2(SiO4)O | 0.9 | |
| Zircon | | + | Zr(SiO4) | 0.3 | |
| Rutile | | + | TiO2 | 0.3 | |
| Jarosite | | + | Kfe3(SO4)2(OH)6 | 0.3 | |
Apatite (Hydroxyapatite) | | + | Ca5(PO4)3OH | 0.4 | |
| Titanite | | + | CaTi(SiO4)O | 0.4 | |
| Granate | | + | sevaral | 0.3 | |
| Galena | | + | PbS | 0.3 | |
| Smectite | | + | Na0.409K0.024Ca0.009(Si3.738Al0.262)(Al1.598Mg0.214Fe0.173Fe0.035)O10(OH)2 | 0.5 | |
| Ca-Vermiculite | | + | Ca0.43Mg3Si3.14Al0.86 bn10(OH)2 | 0.4 | |
For the aqueous-phase geochemical calculations, only the main mineralogical phases were taken into consideration, while in the Table 1. The impurities of the area are given for the sake of completeness. With this approach it is believed that the main geochemical interactions that could occur in the Marismas DHF are represented, however, as the heterogeneity of rock microstructure is possible even in microscale environments, the existence of additional geochemical reactions cannot be excluded. As calcite and aragonite are polymorphs, their presence was represented in the calculations by their common empirical formula given in the software as calcite mineral.
Reservoir features
Chevron discovered the Marismas fields in 1982, which produced gas since 1990. The average depth of the reservoirs ranges from 850 to 1000 m28. In this study, Marismas 3 field (Figure 2) was selected for further investigation with a depth of almost 1000 m29 and an estimated volume of 7 × 108 m3 30. Since this field is still under exploration, some data remain confidential. However, for the sake of the aqueous geochemical calculations presented in this report, some assumptions were made about the gas behavior.

Figure 2. Simplified schematic cross-section over the Guadalquivir Basin and the Gulf of Cádiz (based on 29,32).
The pressure response can change during the gas production from a field, separating the fields in three broad categories depending on the existence or absence of aquifer support as it is described below31:
1. Type 1 – No aquifer support: pressure depletes almost linearly with the increasing produced gas. The recovery can exceed 90 % while a very low abandonment pressure occurs.
2. Type 2 – Weak or limited aquifer support: pressure is reduced regarding the gas production and that is reduced by small amount of water influx. The recovery is almost 70 % with critically higher abandonment pressures.
3. Type 3 – Strong aquifer support: the pressure decrease during the production phase is quickly replenished from the aquifer. The recovery is around 40 %.
Type 1 and 3 are favourable for the CO2 storage, while in Type 2 problems may occur as some initial capacity will be present when compressing up the remaining hydrocarbon gas. Subsequently, the additional capacity will change at the rate at which the aquifer relaxes in terms of CO2 injection, which may be too low for practical application31.
In the Marismas area, the main Almonte-Marismas aquifer (Doñana aquifer) is at shallow depths ranging from 140 to 170 m; thus, strong aquifer support is lacking. As the depleted-gas reservoirs have, on average more than 15 % residual gas saturation as it was referred to in the literature, it was assumed that Marismas 3 belongs to Type 1 type with recovery close to 90 %. The abandonment pressure of such type depleted hydrocarbon reservoirs is assumed to be 2 MPa31.
For the investigation of boundary conditions, a depth of 1000 m was considered. For simplicity, the temperature inside the reservoir was calculated assuming a linear geothermal gradient and homogeneous thermal properties of the sedimentary cover based on the following equation:
where is the geothermal gradient and h is the reservoir’s depth. The mean land surface temperature over the year is almost 18 °C for the Marismas fields area33. As the geothermal gradient is known and equal to 30 °C/km, the expected average temperature according to Equation 1 was calculated as 48 °C.
The porosity of the field is 20 % on average and the residual water saturation is 25 %34. The residual hydrocarbon gas is composed by 98 % of CH4 and thus it was assumed to be the dominant hydrocarbon (100 %) for the run of the calculations. The initial water density was 1.04 and its pH was 7.02, while its temperature is assumed to be in equilibrium with the reservoir temperature. Based on these values, the activity of electrons (pe) in the investigated aqueous solutions was determined based on the following calculation:
where Eh is the redox potential and T the temperature (K). Thus, for a temperature of 48 °C = 321.15 K, the equation becomes:
Based on the graph in Figure 3, the Eh is close to -0.18 V for an environment of deeper ground-water with the specific features given above. Therefore, pe is equal to -2.68, a reasonable value between the typical range of -12 to 25 pe water values. Moreover, to represent the real water characteristics in Marismas field, the calculated salinity of the field was assumed to be the alkalinity of the solution, as given in the Table 2.

Figure 3. Characterisation of solutions by pH and Eh (figure has been reproduced with permission from Railsback37).
Table 2. Residual water characteristics of the studied depleted gas field.
| Properties | Water characteristics |
|---|
| Temperature (°C) | 48 |
| Pressure (MPa) | 2 |
| pH | 7.0 |
| pe | -2.68 |
| Density (kg/m3) | 1040 |
| Water mass (kg) | 364 × 108 |
| n (mol) | 2020507036249* |
| Alkalinity (ppm) | 35000 |
| Volume (m3) | 35 × 106 |
The concentration of the initial components in the reservoir was specified based on the Marismas field porosity (20 %), residual water saturation (25 %), and the reservoir volume which is 7 x 108 m3 30. Thus, the residual water was calculated as 35 x 106 m3 and the free volume of the reservoir was up to 105 × 106 m3. In Table 2, the water properties of the Marismas field used for the PHREEQC calculations are summarised.
Based on the available literature, the depleted-gas reservoirs have more than 15 % residual gas saturation on average5. For the calculation of the residual gas in the reservoir that percentage was adopted and as CH4 was the dominant residual gas (up to 98 %), 21 × 106 m3 CH4 assumed to occupy the reservoir. The presence of other residual components was assumed as negligible. The CH4 temperature should be equal to the reservoir temperature while the brine-methane mixture should be the abandonment pressure for depleted hydrocarbon gas reservoirs, equal to 2 MPa, as referred in the literature35. The CH4 density was found based on the Natural Gas Density Calculator and all the necessary CH4 properties for the run of the PHREEQC models are given in the Table 3.
Table 3. Residual gas (CH4) characteristics of the studied depleted gas field.
| Properties | CH4 characteristics |
|---|
| Temperature (°C) | 48 |
| Density (kg/m3) | 12.33 |
| Pressure (MPa) | 2 |
| Volume (m3) | 21 × 106 |
| m (kg) | 25.893 × 107 |
| n (mol) | 16142768080* |
During the CO2 injection in the DHF implemented for the CEEGS technology, the purity of CO2 was assumed to be 100 %. Based on CEEGS design, the injection of CO2 was planned in supercritical conditions, with temperature above 31°C, pressure above 7.38 MPa, and critical density 411.53 kg/m3 36. For the calculation of the CO2 injected moles in the reservoir, that critical density was used and the amount of moles kept constant in the reservoir despite the conditions’ changes. A deterministic modelling procedure was used to estimate CO2 storage capacity, taking into consideration the geological parameters of the selected area. Thus, the following equation referring to the Depleted Gas Hydrocarbon Fields was used4,38:
where
and Equation 5 becomes
where MCO2 is the mass of the CO2 storage capacity of a prospect field (kg), OGIP is the original gas in place at standard conditions, Bg is the formation volume factor of gas (%), DF is the depletion fraction as a percentage of OGIP, ρRES is the density of CO2 at reservoir storage conditions (kg/m3), EGAS is the storage efficiency factor (%), A is the total area of the reservoir (m2), H is the gross reservoir thickness (m), Φε is the effective porosity of the reservoir (%), and 1 – SW1 is the hydrocarbon saturation, i.e. 1 – water saturation (%).
In Marismas 3 case, the volume of the field is A × H = 7 × 108 m3 [7], Φε is 20 % and SW1 is 25 %. DF is 85 %, and the remaining hydrocarbon gas based on the literature is assumed to be 15 %. ρRES is equal to the CO2 density in the supercritical state, and the EGAS is typically in the range of 1.0 tο 5.0 %39. As the storage efficiency factor is unknown for the Marismas field, both extreme values were used for the estimation of the GCO2 creating to extreme scenarios. Finally:
, and
The above assumptions were used as representative and fed the PHREEQC software1 for subsequent calculations. The CO2 parameters are summarised in Table 4.
Table 4. Characteristics of the potential injected CO2 in the studied depleted gas field.
| Properties | CO2 characteristics |
|---|
| Temperature (°C) | 31.00 |
| Pressure (MPa) | 7.38 |
| Density (kg/m3) | 411.53 |
| Volume (Marismas 3) (m3) | 892500 – 4462500 |
| m (kg)* | 367290525 – 1836452625 |
| n (mol)* | 8345615201 – 41728076006 |
PHREEQC numerical modelling tool
PHREEQC1 is a computer simulation software written in C programming language that aims to perform a wide range of aqueous geochemical calculations based on PHREEQE Fortran program40 completely rewritten to give new capabilities. PHREEQC design is based on an ion-association aqueous model and has potential for:
i. speciation and saturation-index calculation, ii. Calculations for advective-transport and reaction-path, involving specified irreversible reactions, mixing of solutions, mineral and gas equilibria, ion-exchange reactions, and surface-complexation reactions, and iii. inverse modelling counting the transfers of a set of mineral and gas moles to investigate composition differences between aqueous solutions under specified compositional uncertainties41.
PHREEQC software output values include the elements’ concentrations, aqueous species molalities and activities, phase mole transfers to achieve the equilibrium, pH, pe, and saturation indices (see extended PHREEQC output raw data)41. For the relation of mass actions equations to actual solutions, known also as balanced chemical reactions, chemical thermodynamics, or specifically equilibrium thermodynamics are used. The basic principle to be achieved is that elements, molecules or compounds contain some internal energy and the whole system try to reach a state of minimum energy (equilibrium). Natural systems cannot always reach this state; however, they have this tendency42.
PHREEQC has been extensively applied for the investigation of underground gasses’ storage systems aiming to explain the geochemical interactions between reservoirs rocks, the injected gasses to be stored (e.g. H2 or CO2), and the aquatic phases present in the reservoir, e.g., brine and/or residual hydrocarbons43,44.
In the present study, the following assumptions were taken for the calculations:
1. Only the chemistry of the main minerals comprising the host rock (reservoir) was considered.
2. Water and methane found in the reservoir pore volume supposed to be in chemical equilibrium with the reservoir’s minerals.
3. The CO2 injection occurs in the supercritical state, reaching pressure equilibrium with the rock-water-methane system.
4. The injected CO2 is assumed to transmit to gas phase during i. its route to the reservoir and ii. in the reservoir when micing with methane, acting as an ideal gas.
The accuracy of the aqueous geochemical interactions was verified by the estimation of new mineralogical phases available in the phreeqc.dat that were in agreement with the known impurities of the area given in Table 1.
The dataset phreeqc.dat contains the default thermodynamic data of PHREEQC as was derived from PHREEQE40 with limitation regarding the smaller set of elements, aqueous species, and minerals45. The latest phreeqc.dat characteristics included;
1. the parameters needed for the Peng-Robinson equation of state to calculate gases’ fugacity coefficients in a gas mixture and their solubility in water,
2. Solids’ and aqueous species’ molar volumes to calculate the pressure dependence of log Ks (applicable up to ~100 MPa and 200 °C, where Ks is the solubility constant), and
3. the Redlich-type equation for temperature and pressure dependence of aqueous species’ volume46,47.
However, there is a limitation regarding the set of aqueous species, minerals, and gases that include the available parameters for corrections to ~200 °C. Some reactions do not include the parameters for calculating log Ks above 25 °C, while for others, only the enthalpy of reaction is available, and the estimation of logKs is performed by the use of van ‘t Hoff equation48 extrapolating from 25 °C to higher temperatures and assuming that enthalpy of reaction is invariant. This may result in relatively large uncertainties when higher temperatures are used49–51. In this study, only low temperatures were used and thus, the large uncertainties were avoided.
PHREEQC scenarios and input values
Four different scenarios were developed for the PHREEQC calculations, separated in two sets. The two scenarios of each set are based on the boundary conditions of the potential injected moles of CO2, as given in Table 4, reflecting the 1.0 tο 5.0 % range of the storage efficiency factor (SEFGAS)39. In the first set of scenarios, it is supposed that the CO2 is injected in the supercritical state and transmit to gas phase during its route to the reservoir and before interacting with the remaining hydrocarbons, while in the second set of scenarios it is supposed that the CO2 is in gas phase when gets into equilibrium with the remaining hydrocarbons in the reservoir. All the calculations were time-independently as the purpose was the systems to reach the chemical equilibrium between their initial components. For all developed scenarios, the phreeqc.dat database was used with addition of glauconite mineralalogical phase from sit.dat to allow for comprehensive main mineralogical phases’ consideration. The whole calculations were run with phreeqc.dat. Αs phreeq.dat take into consideration all the potential aqueous chemical reactions between the mineralogical phases of the database in the background of the software and thus they are not provided to the user. Table 5 provides details for the selected databases, fitting to the requirements of the present study. The sit.dat uses specific ion interaction theory (SIT) for calculations in concentrated solutions and is suitable for a range of conditions expected in both geological disposal facilities and near-surface conditions (i.e. pH 6–14, ionic strength up to SIT limit, redox potential (Eh) within the water stability fields, and 15 to 80 °C temperatures52.
Table 5. Thermodynamic data characteristics for each database used in this study, based on 45.
| Database | T-P range | Corrected P range | Aqueous activity coefficient model | Fugacity coefficients | Species’ number | Limitations | Reference |
|---|
| phreeqc.dat | <200 °C, <100 MPa | up to ~ 100 MPa | mixed WATEQ and Davies equation | Peng-Robinson | ~310 | • Relatively small species’ number • Ionic strength generally less than one molal | |
| sit.dat | 15–80 °C at 100 MPa | N.A. | Specific Ion Interaction Theory | Ideal gas law | ~2300 | • Small T-P range • C = logKs extrapolation using the van ‘t Hoff equation • Ideal gas law for gas fugacity coefficients • No pressure correction • Ionic strength generally less than one molal | Amphos 21, French Bureau of Research for Geology and Mining and HydrAsa for French Agency for the Management of Nuclear Waste based on ThermoChimie Version 9b045 |
All scenarios included a common first step of calculations, including the host rock/methane/water interactions, followed by the second step, which reflected the CO2 injection and thus the host rock/methane/water/CO2 interactions. The input values for each scenario are extensively described below.
Set 1-Scenario 1
As the main lithologies of the Marismas fields area are carbonate-silisiclastic rocks and sands, they have densities of 2.15-2.40 × 103 kg/m3 and 1.60-2.00 × 103 kg/m3, respectively53. Thus, an average density is almost 2.04 × 103 kg/m3. The reservoir volume is known by the literature and is estimated to be 7 x 108 m3 and the porosity of the reservoir is 20 %, a total volume of 80 % has a mass of 11424 × 108 kg (Table 6).
Table 6. Masses of components interacted in the studied depleted gas field.
| Components | Mass (kg × 108) |
|---|
| Rock | 11424 |
| Water | 364 |
| CH4 | 2.5893 |
| CO2 | 3.67290525 - 18.36452625 |
The mass of the main mineralogical phases was calculated based on the percentages’ estimation given in Table 1, presenting as raw data in the Table 7. The raw data were used for the calculation of moles of initial components provided as input values for PHREEQC calculations ensuring quality importance. The CH4 of the system was assumed to behave as an ideal gas for the purpose of simulations.
Table 7. Quantities of mineralogical phases in the host rock of depleted gas field.
Only the main mineralogical phases were used after normalisation of the amounts in 100 %.
| Mineral | Initial Estimation (in 100 %) | Mass in rock (× 106 kg) | Empirical formula based on the used databases | Molecular mass (g/mol) | Moles × 103
|
|---|
Calcite or/and Aragonite | 29 | 331296 | CaCO3 | 100 | 3312960000 |
| Dolomite | 10 | 114240 | CaMg(CO3)2 | 184 | 620869565 |
| Quartz | 29 | 331296 | SiO2 | 60 | 5521600000 |
| Alkali Feldspars | 2 | 22848 | KAlSi3O8 | 278 | 82187050 |
| Muscovite | 2 | 22848 | KAl3Si3O10(OH)2 | 398 | 57407035 |
| Illite | 3 | 34272 | K0.6Mg0.25Al2.3Si3.5O10(OH)2 | 384 | 89250000 |
| Kaolinite | 9 | 102816 | Al2Si2O5(OH)4 | 258 | 398511628 |
| Chlorite | 2 | 22848 | Mg5Al2Si3O10(OH)8 | 556 | 41093525 |
| Glauconite | 14 | 159936 | K0.75Mg0.25Fe1.5Al0.50Si3.75O10(OH)2 | 432 | 370222222 |
| Total | 100 | 1142400 | | | |
The rock mass and the stoichiometry of the main mineralogical phases given in Table 7 were used to specify the number of moles interacted in the system for the PHREEQC calculations.
To save computational time and avoid computational limitations posed by PHREEQC software, it was assumed that the whole system is composed of 10000 moles. The ratios between the system’s initial components are kept stable to represent reality. In this scenario, the minimum ratio of storage efficiency factor, which means 1.0 % was used (Table 8). For the gas phases (i.e., CH4 and CO2) used in the calculations, the log of gas partial pressure was used as a target saturation index. This may not be attained if the amount of the phase in the assemblage is insufficient. The Ideal Gas Law was used to calculate the partial pressure in the Partial Pressure Calculator54 as the gasses were assumed to behave as ideal. In this case, the CO2 supposed to transmit into gas phase during the route to the reservoir keeping its initial temperature (31°C). Moreover, it was supposed that the two gasses are not yet mixed and thus the partial pressure of each one as a separate gas was calculated.
Table 8. Quantities of initial components in the complex system after the injection of CO2 as estimated and used for PHREEQC calculations for scenario 1.
| Component | Moles × 103 | Moles in 10000 (used in PHREEQC) |
|---|
| Calcite or/and Aragonite | 3312960000 | 2642 |
| Dolomite | 620869565 | 495 |
| Quartz | 5521600000 | 4404 |
| Alkali Feldspars | 82187050 | 66 |
| Muscovite (K-mica) | 57407035 | 46 |
| Illite | 89250000 | 71 |
| Kaolinite | 398511628 | 318 |
| Chlorite | 41093525 | 33 |
| Glauconite | 370222222 | 295 |
| Water | 2020507036 | 1611 |
| CH4 | 16142768 | 13 |
| CO2 | 8345615 | 7 |
| Total | 12514608061 | 10000 |
Set 1 - Scenario 2
This scenario is similar with Set 1 – Scenario 1 with the only difference in the used maximum storage efficiency factor ratio, which was equal to 5.0 %. To achieve comparable results, the other initial components are kept in the same ratios; thus, the whole system is composed of 10026 moles (Table 9).
Table 9. Quantities of initial components in the complex system after the injection of CO2 as estimated and used for PHREEQC calculations for scenario 2.
| Component | Moles × 103 | Moles in 10026(used in PHREEQC) |
|---|
| Calcite or/and Aragonite | 3312960000 | 2642 |
| Dolomite | 620869565 | 495 |
| Quartz | 5521600000 | 4404 |
| Alkali Feldspars | 82187050 | 66 |
| Muscovite (K-mica) | 57407035 | 46 |
| Illite | 89250000 | 71 |
| Kaolinite | 398511628 | 318 |
| Chlorite | 41093525 | 33 |
| Glauconite | 370222222 | 295 |
| Water | 2020507036 | 1611 |
| CH4 | 16142768 | 13 |
| CO2 | 41728076006 | 33 |
| Total | 12572478906397 | 10026 |
Set 2 – Scenario 1
This scenario is the same as Set 1 – Scenario 2 with the only difference that the CO2 supposed to transmit into gas phase when interacting with the remaining hydrocarbons in the reservoir. Thus, the Ideal Gas Law was used to calculate the partial pressure of the gas mixture (CO2 and CH4) in the Partial Pressure Calculator55. In this case, the CO2 and the CH4 supposed to have a common temperature of 48 °C as the temperature of the reservoir.
Set 2 – Scenario 2
This scenario is the same as Set 2 – Scenario 1 with the only difference in the used maximum storage efficiency factor ratio, which was equal to 5.0 %.
CMG-GEM numerical modelling tool
CO2 injection and subsequent plume migration in the “Depleted Hydrocarbons Field (DHF)” scenario are modelled using the CMG-GEM compositional multi-phase flow subsurface simulator3. The modelling of CO2 injection and back production involves the component transport equations solution, the solution of equations for thermodynamic equilibrium between the gas (CO2 and methane) and aqueous (formation water) phase, and the geochemical reactions solutions. The solution adopted in GEM for CO2 transport-reaction simulations for porous media reservoirs is the flow simulation for adaptive-implicit multiphase multicomponent under phase and chemical equilibrium and rate-dependent mineral dissolution/ precipitation by the use of the fully-coupled approach56,57.
CMG-GEM model description and input values
A more generic, small-scale model of natural gas reservoir for a partially closed depleted gas reservoir with an underlying aquifer was developed. By introducing the element of an underlying aquifer in the system allowed to make a more complex approach possible to represent a wider range of areas of interest. According to Carneiro and Behnous guidelines58, the minimum pressure required for a reservoir at a depth of 2000 m to ensure the production of supercritical CO2 in the CEEGS technology, should exceed 17 MPa. In CMG-GEM simulation, the targeted reservoir is located at a depth of 2000 m in contrast to Marismas field depth which was 1000 m, starting with an initial pressure of 14 MPa and a temperature of 75°C (with a normal geothermal gradient of 30 °C/km). The reservoir exhibits a water saturation of 25%, and the remaining gas composition comprises 98% CH4, 1% C2H6, 0.5% C3H8, and 0.5% CO2. Water was assumed to have flow back to the reservoir, following production of the gas, and is virtually the single existing phase in the deepest parts of reservoir. The reservoir dimensions are specified in this more generic scenario as 1500m x 1500m x 50m (see Figure 4). Detailed information on other reservoir parameters are given in Table 10.

Figure 4. Water saturation in the depleted hydrocarbon reservoir (blue grids: water saturation 100%, green grids: water saturation is 25%).
Table 10. Parameters of the depleted reservoir simulation model.
| Parameter | Value |
|---|
| Reservoir size | 1500m x 1500m x 50m |
| Depth to reservoir top (m) | 2000 |
| Porosity (%) | 20 |
| Permeability (md) | 200 |
| Thickness (m) | 50 |
| Reservoir initial temperature (°C) | 75 |
| Reservoir initial pressure (MPa) | 14 |
| Residual water saturation (-) | 0.25 |
| Residual Gas saturation | 0.3 |
| Gas water contact (GWC) (m) | 2040 |
| Natural gas composition (mole fraction) | 98% CH4, 1% C2H6, 0.5% C3H8, and 0.5% CO2
|
| Injected gas rate | 33 kg/s of pure CO2
|
Top and bottom boundary conditions of the reservoir | No fluid flow and no heat flux |
The reservoir consists of two wells, designated as A and B. To inject CO2 into the depleted natural gas reservoir while minimizing water upconing at the production well, the process begins with a CO2 plume setup phase. During this phase, CO2 is continuously injected at a rate of 33 kg/s (equivalent to 1 Mtonne/year) for two years to establish a supercritical CO2 plume. Once the initial setup is complete, the second stage involves energy storage through charge-discharge cycles. In the charge phase, surplus renewable electricity is used to compress, heat, and inject CO2 into the reservoir via Well A. During the discharge phase, Well A produces supercritical CO2, which is utilized in the surface components of the CEEGS system to generate electricity. Meanwhile, Well B reinjects CO2 into the reservoir at the same flow rate as Well A but at a lower temperature. In brief, this scenario involves two years of continuous injection to establish the CO2 plume, followed by a one-month shut-in period. Additionally, six cycles of charge-discharges are included in the analysis.
The initial CO2 plume setup plays a critical role not only in the economic feasibility of the system but also in determining the relative contributions of CO2 permanent sequestration, energy storage, and geothermal heat extraction.
During the charge phase, Well A injects CO2 at a wellhead temperature of 60°C, while Well B remains inactive. In the discharge phase, Well A produces CO2, which is subsequently reinjected through Well B at the same flow rate but at a reduced temperature of 20°C. The injection process is subject to two key constraints: the bottomhole pressure must not exceed 20% of the initial hydrostatic pressure, and sufficient injection pressure is required to deliver CO2 to the bottom of the reservoir. The cycle of charge, shut-in, and discharge was repeated six times.
Results and discussion
PHREEQC calculations
Initial conditions
All the developed scenarios of the present study have a common first step including investigating interactions between the reservoir rock/residual water/CH4 of the Marismas 3 field. The pH of the aqueous solution in the reservoir before the injection changes from the initial water pH of 7.0 to 9.6, becoming more alkaline. This alkaline nature can be explained by the ions’ presence in the reservoir, i.e. Al, K, Ca, and Mg. The pH change is the result of a charge balance. The pe parameter was also changed to -7.26 adjusting to redox equilibrium. Among the main mineralogical phases of the reservoir host rock, i.e. calcite, dolomite, quartz, K-feldspars, K-mica, illite, kaolinite, chlorite, and glauconite, only the K-feldspar and illite were unstable. This can be explained by the gradual complete dissolution of K-feldspar, resulting in the precipitation of a variety of new mineralogical phases after long period of time59. In such conditions, gibbsite (Al(OH)3), K-mica, and kaolinite are some of the mineralogical phases that can occur after the dissolution of K-feldspar60. However, the precipitation rate of kaolinite is slower than the dissolution of the feldspar, leading to kaolinite-supersaturated solutions in the reservoir61. Illite dissolution is more selective and occurs only under alkaline pH conditions, at temperatures 100 °C and above, while when the pH decreases below 5, a rapid precipitation of the phase of aluminum oxy(hydroxide) is expected55. However, in this case study, the system’s newly formed minerals containing aluminum were undersaturated. This means that despite the lack of phreeqc.dat database in this phase and the computational limitations, there is no possibility of such phase precipitation. A change in the K-mica interacting moles (from 46 to 155) was observed due to the dissolution of K-feldspar and the precipitation of newly formed K-mica.
The other main mineralogical phases of the reservoir were in equilibrium with the system. The reservoir was supersaturated in hematite (Saturation Index-S.I. 11.48), goethite (S.I. 4.68), siderite (FeCO3) (S.I. 3.50), and talc (Mg3Si4O10(OH)2) (S.I. 1.48). The presence of hematite and goethite agreed with the impurities identified in the study area as given in Table 1, owing their origin in glauconite transformation62. Glauconite can also alter K-mica63 and precipitate siderite by diagenetic conversion64. A main difference among glauconite and siderite is the presence of mostly oxidised Fe3+ in the first one, and the presence of reduced Fe2+ to the second one63. Talc (Mg3Si4O10(OH)2) precipitation can be attributed to i. sedimentary magnesium carbonate rocks alteration at elevated pressure and temperature, and/or ii. magnesium-rich ultramafic rocks alteration under high pressures and temperatures up to 700 °C. The second origin is not possible in this study area, as the reservoir temperature is critically lower. Thus, talc can origin by chlorite + quartz = kyanite + talc + water65.
Set 1 - Scenario 1
The second step of PHREEQC calculations includes the investigation of interactions between rock/water/CH4 system after the injection of CO2. The pH critically decreases after the CO2 injection, changing from 9.6 to 7.6 due to the dissolution of CO2 in the water, which starts as soon as the CO2 is injected66. Due to buoyancy, the CO2 flows to the top of the reservoir. In this upward flow, an amount of CO2 partly reacts with the water becoming carbonic acid (H2CO3). This amount is generally small due to the small available storage capacity67 and the water availability in the system. The Total CO2 increased from a value of 0.46 mol/kg for the rock/water/CH4 system to 0.69 mol/kg due to the CO2 injection, as it was expected. This parameter represents the Dissolved Inorganic Carbon (DIC), without being restricted to CO2 but also including HCO3, , and other carbon-based complexes as these provided by each run of software’s simulation. The pe parameter was equal to -4.94, adjusting to redox equilibrium. The CO2 gas completely dissolved in the reservoir solution when reached the equilibrium, significantly decreasing its pressure from 7.38 MPa to 0.11 MPa, which is related to the change of the temperature and its equilibrium with the brine-methane system having an initial pressure of 2.00 MPa. As the system comes to equilibrium, the system pressure stabilises and the CH4 pressure will slightly decrease to 1.93 MPa.
In the case of rock/water/CH4 a complete dissolution of the illite was observed. In contrast, when CO2 is injected in the system, the illite not only is not dissolved but it is precipitated (S.I. 0.36). This behavior can be explained by the complete K-feldspar dissolution resulting in quartz and clay precipitation under abundant CO2 presence at underground P-T conditions, as it was verified by experimental data available in the literature68. Moreover, the CO2 injection in a reservoir containing K-feldspar increases significantly its dissolution rate as the initial equilibrium of the reservoir with the containing fluids is disrupted. However, the remaining fluids are then supersaturated to kaolinite approaching the previously established equilibrium, resulting in relatively low dissolution rates for the Al-Si minerals, as it was proved by experimental data Thus, the rhythm of secondary products precipitation is controlling the K-feldspar dissolution rate as it allows the fluid to approach saturation with K-feldspar than secondary products61. As the secondary minerals can control the dissolution of K-feldspar, its complete dissolution and the reach of the reservoir in chemical equilibrium can be a process that take a long period of time. The supersaturation of quartz (S.I. 0.49), chalcedony (S.I. 0.13), kaolinite (S.I. 5.04), K-mica (S.I. 5.58), and gibbsite (S.I. 1.57) can be also explained by this process. Another important differentiation in the system after the CO2 injection, which is in agreement with the literature68, is the carbonate minerals dissolution, with S.I. -0.68/-0.87 and -1.60 for calcite/aragonite and dolomite minerals, respectively. A series of reactions take place after the CO2 injection into the solution, starting from the CO2 dissolution in the water which forms H2CO3 acid. This acid is diluted in the aqueous solution of the reservoir causing a pH drop and an acid attack to the rock minerals. The following reactions with the carbonate minerals explain the less acidic nature of the brine:
The less acidic brine created by the CO2 injection, has the ability to dissolve minerals that contain divalent cations such as Mg, Ca, and Fe69. Common minerals containing divalent cations are chlorite, montmorillonite, talc and glauconite70. In this case, the dissolution of glauconite (S.I. -3.60), talc (-6.73), and chlorite (-11.42) were observed, while despite the partly dissolution of montmorillonite, other processes enhanced its precipitation. Thus, due to the dissolution of several mineralogical phases that extensively described above, as the rapidly dissolved calcite (S.I. -0.68) and to lesser extent dolomite (-1.60), and the prolonged dissolution of K-feldspar (-3.32) and chlorite (-11.42), liberated ions are present in the brine having the ability to react with each other resulting in clay minerals formulations, such as kaolinite (5.04), Ca-Montmorillonite (S.I. 3.85), and illite (S.I. 0.36). Hematite and goethite remain supersaturated also after the CO2 injection, with significantly lower values, i.e. S.I. 6.09 and 1.99, respectively. This decrease is owing to Fe2+ release from the glauconite and chlorite dissolve which is oxidised to Fe3+, resulting in such minerals’ precipitation. Amorphous phases of Fe(OH)3, Al(OH)3, and SiO2 are undersaturated in the system. Siderite (FeCO3) as a carbonate mineral is partly dissolutes, however, the dissolution of glauconite released Fe and the redox reaction of calcite and/or dolomite with H2CO3 acid, gives to the system the potential to form an amount of siderite (S.I. 3.08) which finally precipitates. Moreover, the dissolution rate of carbonates decreases in the following order: calcite > dolomite > siderite > magnesite71, which proves the slower siderite dissolution than calcite and dolomite.
Permeability reduction in porous reservoir media is widely investigated during the CO2 injection72,73. In Marismas field, a first increase of the permeability could be observed during the CO2 injection, owing to the carbonates dissolution as it was referred in the literature from studies held in similar environments74. However, potential decrease in the pressure result in the precipitation of the dissolute carbonate minerals or other secondary minerals, such as illite, kaolinite or Ca-montmorillonite that tend to migrate to pore throats decreasing the permeability75. Moreover, the existing clay minerals of the reservoir may cause permeability reduction due to the release of clay particles from pore walls of the reservoir and their subsequent redeposition in the pore throats downstream, which potentially have smaller diameters than the pores76. However, when chlorite is present in the reservoir, aluminum silicate minerals and chlorite itself are dissolved to some extent, minimising the precipitation of secondary carbonate minerals and improving the reservoir properties77. In this study, a significant amount of clay minerals in precipitating, including illite, K-mica, kaolinite, and Ca-montmorillonite which could transfer in the fluid flow path and potentially accumulate at pore throats, reducing the permeability. Significant changes in the permeability can be due to the swelling clays75. However, most clay minerals precipitated in the investigated system are non-swelling minerals, except Ca-montmorillonite. Ca-montmorillonite dissolution can be controlled by the control of temperature, pH, and time78. In this case, where low temperatures exist in the reservoir, a change of pH in more alkalic or more acidic conditions favors the dissolution of the mineral to avoid clogging problems, as in almost neutral environment the mineral is in equilibrium78. However, further investigation is needed.
The chemical precipitation of certain ions such as manganese (Mn) and iron (Fe) are extremely important in the process of chemical clogging. Simultaneously, it is strongly connected with the metal bacteria involvement that can enhance the chemical clogging. The partially oxidised and low crystalline intermediates of Fe (Fe2+, Fe3+) followed by stable crystalline Fe3+(hydro)oxides (e.g. goethite) are the end products of iron redox reaction79. The reduced precipitation of hematite minerals can ensure the goethite precipitation in a relatively stable form that is not correlated to the clogging of the pores80. Special attention is required for the siderite’s precipitation. Its small colloidal particles remain suspended in the fluid, and the ability to transport over long distances and eventually accumulate in certain areas may lead to the reservoir's clogging or corrosion81. The management of Fe-minerals precipitation can be easily controlled by changing the system's parameters and especially pH, pe, and oxygen concentrations to avoid failures in the reservoir79.
Despite the constant attention that permeability requires during injection, the porosity of the investigated system finally increases after the CO2 injection, as the solution occupies a volume of 30.92 L instead of 30.66 L of rock/water/CH4 system.
Set 1 - Scenario 2
Scenario 2 gave similar results to Scenario 1 as simulates the same system with only difference the higher concentration of CO2 reflecting the maximum ratio of storage efficiency factor. Partial differentiations occurred in the system relative to higher CO2 amount. Firstly, the pH of the aqueous solution is getting more acidic after the higher injected CO2, changing from 9.6 to 6.0, a value lower than the Scenario 1 (pH 7.6) related to the CO2 partially reaction with the water creating carbonic acid (H2CO3). Thus, a higher percentage of carbonate minerals is expected to be dissolved. This phenomenon is more intense when the storage capacity of CO2 or brine quantity are increasing. When the CO2 reaches the top of the reservoir, the diffusion is the dominant process determining the dissolution rate. As this process is slow, the resulting dissolution rate is also low. However, when the dissolved CO2 increases, an increase in the brine density is also observed, creating an unstable buoyant layer below the CO2. When this phenomenon occurs and the layer becomes sufficiently unstable, a downward convective flow with a fingers’ shape is developed, causing an upwelling transport of fresh brine to the interface of CO2/brine that has the ability to increase the CO2 dissolution rate66. Despite that in the Scenario 2 the storage capacity is significantly higher than in Scenario 1, as the CO2 was completely dissolved as a gas in the reservoir’s solution, it has not the ability to create an unstable buoyant layer that could significantly increase the carbonates dissolution causing potential stability problems of the reservoir. This complete CO2 dissolve in the reservoir’s solution as a gas caused a decrease in its pressure from 7.38 MPa to 4.53 MPa, critically higher than in Scenario 1 (0.11 MPa). The maintenance of pH values in values higher than 5, despite the increase of CO2 storage capacity, causing the dissolved carbon dioxide main transformation into bicarbonate ions, which will act as dissolution inhibitors later on, and based on the available literature82,83 may cause a greater reservoir stability than the one that was expected. The pe parameter was increased to -3.14. The volume of the aqueous solution occupied in this Scenario was 31.84 L, slightly higher than in Scenario 1 (30.92 L) corresponding to the higher dissolution of carbonate minerals due to the CO2 increase, leading also to the increase of free space. As a result, the Total CO2 was also increased from a value of 0.46 mol/kg for the rock/water/CH4 system to 0.69 mol/kg for the Scenario 1 and to an even higher value of 1.56 mol/kg for the Scenario 2 where a higher CO2 concentration was injected. As the pressure of CH4 after the system reached equilibrium remained constant as in the Scenario 1.
Concerning the mineralogical phases of the system, slight changes occurred when the concentration of CO2 increased. An amorphous Al(OH)3 supersaturation occurred (S.I. 0.44) due to chlorite and K-feldspar dissolution. This wasn’t observed in the Scenario 1 due to complete crystallisation of Al(OH)3 to gibbsite (Al(OH)3) and clay minerals. As the Si-based mineralogical phases remained almost stable, there was a lack of free silicon to potentially create other aluminosilicate minerals containing an Al(OH)3 amount. In the present Scenario, common clay minerals were precipitated as in the Scenario 1. A differentiation was occurred in the precipitation of Fe-minerals, i.e. goethite (S.I. -0.87), siderite (S.I. 1.51), and hematite (S.I. 0.37), related to the pH and pe changes that have the ability to dissolve the goethite that exist in the area as an impurity. This allows for the hematite and siderite precipitation to a lower extent. The less siderite percolates in the system, the better as it reduces the pore throat clogging effect.
Set 2 – Scenario 1
Scenario 1 of the second set gave almost the same results with Set 1 – Scenario 1, with only difference the CH4 pressure which was slightly lower than Set 1 -Scenario 1, equal to 1.80 MPa.
Set 2 – Scenario 2
Scenario 2 of the second set gave similar geochemical interactions with the Set 1 – Scenario 2. The general characteristics of the system kept constant. The saturation indices of the minerals were slightly differentiated, however, the minerals that were undersaturated kept undersaturated, e.g. anorthite (S.I. -2.65), aragonite (S.I. -2.28), calcite (S.I. -2.09), chlorite (S.I. -23.54), chrysotile (S.I. -19.84), dolomite (S.I. -4.46), glauconite (S.I. -8.81), goethite (S.I. -0.70), K-feldspar (S.I. -3.43), sepiolite (S.I. -14.02), and talc (S.I. -15.57), while the minerals that were precipitated had the same behavior, e.g. Ca-montmorillonite (S.I. 6.56), chalcedony (S.I. 0.15), gibbsite (S.I. 2.91), illite (S.I. 1.86), hematite (S.I. 0.72), K-mica (S.I. 8.15), kaolinite (S.I. 7.75), quartz (S.I. 0.51) and siderite (S.I. 1.59). In this case, the CO2 was not completely dissolved as a gas in the reservoir’s solution and its pressure was 3.59 MPa. The pressure of CH4 after the system reached equilibrium was 1.57 MPa.
CMG-GEM calculations
The development of the CO2 plume following two years of continuous injection is given in Figure 5. The plume takes on a distinct shape, taking to an upright cone. Notably, the maximum extent of the plume is noticeable in the lower section of the reservoir, near the point of separation between the aquifer and the gas. The plume's shape can be attributed to the higher density of CO2 in reservoir conditions (pressure and temperature) compared to impure methane (six times higher). A crucial finding is that for achieving optimal CO2 concentration, production is advised from beneath the horizontal centerline of the reservoir. In addition, as more CO2 is injected into the reservoir, a gas mixture of CO2 and CH4 is produced, where the CO2 fraction increases and the CH4 fraction decreases with time (Ezekiel et al., 2020). This may require additional fluid separation at the surface to obtain a high-purity supercritical CO2, making it suitable before entering surface components. However, this is disadvantageous as it could reduce the system’ efficiency.

Figure 5. CO2 Plume distribution (Gas saturation) after 2 Years of Continuous Injection.
Figure 6a illustrates the evolution of bottom-hole pressure in the injection well (Well A) over a two-year period. Initially, before starting injection, the bottom-hole pressure is assumed to be 14 MPa. This pressure, commonly occurring in depleted reservoirs due to long-term gas production, does not allow for the production of supercritical CO2. However, after approximately two years of injection, equivalent to roughly more than two million tons of CO2, the pressure increased to 21.3 MPa. Importantly, this pressure remains below the threshold assumed to avoid fracturing the reservoir and seal, that is 20% increase above the hydrostatic pressure, and enabling the production of supercritical CO2 for the effective operation of CEEGS technologies.

Figure 6.
a) Evolution of Bottom Hole Pressure in Well A injector over a 2-Year Period, b) Evolution of Bottom Hole Pressure during Two-Year Continuous CO2 Injection in Depleted Reservoir (8MPa, 2000m), c) CO2 Distribution and Storage Dynamics during Plume Setup Phase for the case of a highly depleted reservoir, the simulated conditions referred to a reservoir pressure of 8 MPa at 2000 meters.
The same scenario as the previous one was considered. After continuous CO2 injection for two years (2 million tonnes), as shown in Figure 6b, the bottom hole pressure increases as the injection proceeds and reaches a maximum of 15.5 MPa after two years.
As discussed in 58, this pressure is low and does not allow the production of supercritical CO2. The produced CO2 reaches the wellhead with a pressure of 5.11 MPa, which is below the supercritical state. Consequently, CO2 reaches the wellhead in the gaseous phase. To produce supercritical CO2 in this situation, it must either increase the CO2 injection rate or extend the duration of the plume stage to more than 2 years, with both possibilities resulting in an increase in the reservoir pressure.
Figure 6c illustrates the distribution of CO2 during the plume setup phase, highlighting the amounts trapped, dissolved, and stored in the supercritical phase. A substantial portion of CO2 is stored in the supercritical phase, with 12% becoming trapped, and 3.6% dissolving in brine (Table 11). Notably, the stored CO2 in the supercritical phase is essential for heat maintenance in CEEGS technologies.
Table 11. Amount of CO2 injected, supercritical CO2, residual trapped and dissolved in the reservoir during Plume setup phase.
| CO2 Amount | Moles | Fraction of injected CO2 |
|---|
| Total injected | 4.73873E+10 | - |
| CO2 Supercritical (mobile) phase | 4.60418E+10 | 96.4% |
| CO2 Trapped Sg < Sgc / Hysteresis | 5.98171E+09 | 12 % |
| CO2 dissolved in Water | 1.71166E+09 | 3.6% |
The variation in Bottom-hole Pressure (BHP) during charge-discharge cycles are illustrated in Figure 7. As evident during the discharge phase (indicated by the low-pressure segment in the curve), Well A is actively producing CO2. Notably, the BHP during this phase is 18.1 MPa. According to 58, this pressure level permits the production of supercritical CO2, a crucial aspect for the functioning of CEEGS technologies.

Figure 7.
a) Bottom Hole Pressure Variation during Charge-Discharge Cycles in Well A Producer, b) Water and Methane Dynamics in Production Well During Charge-Discharge Phase.
Figure 7b illustrates both the water rate and methane mass rate generated in the production well during the discharge phase. The water volume gradually increasing and reaches 750 litres per day in the last cycle, despite that, this amount remains very low compared to the produced CO2 (0.02 %wt). This emphasizes the significance of installing a water separator at the production well, to ensure good functioning of the supercritical CO2 surface installation. In contrast, the quantity of methane accompanying the CO2 is in the order of milligrams per day. This amount is sufficient for the effective operation of the CO2 surface installation.
According to the preliminary results from the hydrodynamic simulation of a generalized depleted gas reservoir, the establishment of a CO2 plume is crucial for restoring reservoir pressure to pressure levels compatible with the requirement of the CEEGS and effectively charging the reservoir with CO2. This process is essential for establishing a CO2 connection between the injection and production wells in the supercritical phase, as outlined by the CEEGS concept. Additionally, it is recommended to reduce the water content in the reservoir to minimize the potential interactions between acidic CO2-rich water and the host rock.
Conclusions
The CO2 storage in DHF is of high interest for the proper management of CO2 emissions. The storage capacity of DHF and the present infrastructure make them an attractive option for the underground CO2 storage. Geochemical interactions between the remaining fluids in the DHF (i.e. water and CH4), the mother rock of the reservoir and the injected CO2 are critical risks that could cause failure to implement this technology. The quasi steady-state simulations’ software PHREEQC was used to investigate the aquatic geochemical interactions in the present study. Marismas field selected as a model DHF for the CO2 injection under supercritical conditions required for the novel CO2 based electrothermal energy and geological storage system technology.
The reservoir rock/water/remaining gas (CH4) interactions proved that all the main mineralogical phases of the area were in equilibrium in the reservoir except K-feldspar and illite. Gibbsite (Al(OH)3), K-mica, and kaolinite were formed as new mineralogical phases after the dissolution of K-feldspar. The reservoir was supersaturated in hematite, goethite, siderite, and K-mica, which caused to glauconite transformation. Talc was also observed as a new phase by the geochemical calculations due to the carbonate rocks alteration. During the CO2 injection, it is partially dissolute in the water making the brine more acidic and causing an acid attack to the rock minerals. Carbonate minerals are strongly influenced by this attack while glauconite, talc, and chlorite dissolution were also observed. The prolonged K-feldspar dissolution resumed also after the CO2 injection. The liberated ions of these minerals were present in the brine, having the ability to react with each other, resulting in clay minerals formulations. Kaolinite, illite, and Ca-montmorillonite were the main newly formed clay minerals, while only Ca-montmorillonite threatens the permeability of the reservoir due to its swelling properties. In the case of the minimum potential CO2 injected concentration, hematite, goethite, and siderite were present in the reservoir. Goethite is a relatively stable phase of Fe-minerals that is not related to pore-clogging. However, during the increase of the CO2 concentration, the goethite dissolves and hematite as well as siderite precipitated in a lower extent. As siderite is composed by small colloidal particles, its presence could cause a pore throat clogging effect. Thus, the less siderites percolated in the system, the better.
Summarising, from a geochemical perspective, only Ca-montmorillonite and siderite could cause failures in the investigated reservoir, however, their presence can be easily controlled by anthropogenic changes in the reservoir parameters and especially pH.
This study highlights the technical feasibility of CEEGS (CO2-based Electrothermal Energy and Geological storage) technology, emphasizing the geochemical suitability of depleted natural gas reservoirs as underground storage candidates. Depleted gas reservoirs are considered ideal for CO2 storage due to their proven structural traps and caprock seals, which effectively prevent the lateral and vertical migration of gas, enabling long-term containment.
For successful CEEGS implementation, simulation results indicate that the establishment of a CO2 plume is a critical phase. This plume helps to stabilize or even increase reservoir pressure to levels conducive to CO2 storage, reducing the risk of water upconing. Establishing a stable CO2 plume between injection and production wells allows for continuous injection and withdrawal rates during the energy storage stage, optimizing the efficiency of the underground part.
Additionally, evaluating the geochemical integrity of depleted gas reservoirs is essential for safe, long-term storage. Key aspects of the storage complex—including the reservoir, caprock, and wells—must be carefully assessed. This assessment is best achieved through an integrated approach, combining field observations, laboratory experiments, and numerical modeling to gain a comprehensive understanding of geochemical stability.
The results of this study suggest that the proposed depleted gas reservoir for the CEEGS technology underground is promising for real-world applications. Future research should focus on exploring the combined effects of reservoir heterogeneity, CO2 dissolution, and rate-limited reactions to refine the understanding of this technology's long-term viability.