Carbon dioxide emissions are a huge contributor to global warming and climate change. Human activities, such as burning fossil fuels, deforestation, and industrial activities, are significant sources of carbon dioxide emissions, which have several negative impacts on the environment [1–3]. Hence, reducing carbon dioxide emissions is a global concern as it is necessary for mitigating the impacts of climate change. Governments, businesses, and individuals all have a role to play in reducing carbon dioxide emissions through a variety of measures, including transitioning to clean energy sources, improving energy efficiency, and reducing reliance on fossil fuels.

Among the potential solutions, the capture, storage, utilization, and conversion of carbon dioxide into useful and highly value-added chemical products stand out, since carbon dioxide is a safe, cheap, and abundant carbon source; this represents a highly positive impact on the environment because it reintroduces carbon dioxide into the carbon cycle [4–8].

Some alternatives for carbon dioxide conversion into valuable chemical products are: electrochemical reduction, photochemical reduction, and hydrogenation of carbon dioxide [9]; here green hydrogen from water hydrolysis could play a key role. A common limitation of these approaches is the high thermodynamic stability of carbon dioxide, which means that for its chemical transformation, elevated external energy is needed.

One option to overcome carbon dioxide stability is through hydrothermal reduction, i.e., the conversion of CO₂ dissolved in high-temperature, high-pressure water, in a process that mimics the reactions that take place in the natural environment of submarine volcanoes [10]. This process offers the potential to turn a harmful greenhouse gas into valuable chemical products that can be used in a variety of industries, such as formic acid, methanol or methane. Furthermore, a useful characteristic of hydrothermal reduction is that it works with aqueous feeds, which means that the solutions obtained by absorption of CO₂ into alkaline water streams, which is the most technically developed carbon capture process available at the moment, can be directly used as feed for the hydrothermal reduction process. For example, if CO₂ is captured with aqueous solutions of sodium hydroxide, the resulting solution, rich in sodium bicarbonate, is a suitable inorganic carbon source for the hydrothermal reduction process [11–13].

The performance of the hydrothermal reduction process is based on the peculiar properties of water at high-temperature, high-pressure conditions, close to the critical point of water (T_C = 647.3 K and P_C = 221.2 bar). For example, properties such as dielectric constant reduce from 78 at ambient conditions (298 K and 1 bar) to 6 at the critical point, a value typical of a nonpolar solvent [14]. Besides, the solubility of some gases such as carbon dioxide increase close to the critical point. The ionic product of water (pKw) rises with temperature from 14 at ambient conditions until a maximum of around 11.2 close to 473 K, and then under supercritical conditions decreases until a value below 22. For that reason, pure water near 473 K behaves as if it were simultaneously more acid and more basic than pure water under normal conditions [15].

The hydrothermal reduction process can be carried out either using gaseous hydrogen as reductant, or using water itself as hydrogen source; in the second case, an additional reductant to promote the decomposition of water into hydrogen and oxygen, such as a metal [16, 17] or an organic reductant [18, 19], must be provided. Typically, formic acid is the main product yielded by this process [12, 18], but several products can be obtained from carbon dioxide under hydrothermal reduction using metal reductants and catalysts. For instance, methane was obtained from sodium bicarbonate using Raney Ni nanoparticles as a catalyst [20], and acetate was produced using a catalyst based on cobalt [21].

For the elucidation of the reaction mechanisms and the optimization of the process, a crucial aspect is to determine the ionic species present in the solution as a function of the reaction conditions. Indeed, upon dissolution in water, CO₂ undergoes transformations to carbonic acid, bicarbonate and carbonate according to its acid-base equilibrium. Of this species, it is considered that bicarbonate is the most reactive one, while carbonic acid and carbonate show a lower reactivity [11, 12]. Hence, it is of the highest significance to determine if the dissolution equilibrium of CO₂ is displaced towards bicarbonate.

This acid-base equilibrium is influenced by several factors, such as, of course, pressure and temperature, but also by the presence of other substances dissolved in the medium, and particularly ionic species such as those produced by the dissolution of salts. This equilibrium is also significant for geological studies, and as such has been studied by different authors [22], as discussed in detail in the following Section 2.

Based on these considerations, the objective of this work is to model the thermodynamic behavior of a mixture of water, carbon dioxide, and hydrogen under high temperature and pressure conditions, relevant for the development of the hydrothermal CO₂ reduction technology, in order to determine the molality of the species in reaction conditions for these process. In particular, conditions that favor the displacement of the equilibrium towards bicarbonate, which is the reactive species, are sought. The study aims to contribute to the field of carbon dioxide reduction transforming waste into economically valuable chemical products helping achieve the Sustainable Development Goal (SDG) due to the potential of hydrothermal reduction as a viable method for capturing and converting carbon dioxide into useful chemical products, with a focus on the production of formic acid.

The thermodynamic model considered in this work describes the speciation equilibrium of H⁺, Na⁺, OH⁻, Cl⁻, HCO₃ ⁻, CO₃ ^2-, CO₂, and H₂ in an aqueous solution of carbon dioxide, hydrogen, sodium chloride and sodium bicarbonate at hydrothermal conditions.

The equilibrium between the different species is a function of pressure, temperature, and composition. Two sorts of physicochemical equilibrium are essential to have a whole understanding of the systems; these are chemical equilibrium and phase equilibrium. Both are complex functions of pressure, temperature, and molality.

The model developed by [22] for the phase equilibrium of water, carbon dioxide, and sodium chloride was used to describe the behavior of the mixtures applied in the hydrothermal reduction of CO₂. This model considers a temperature range from 273.15 to 523.15 K and a pressure range from 0 to 1,000 bar; it is based on the previous model of [23] which considers a hydrate phase and an aqueous phase [24]. The parameters required to apply this model were taken from the literature, as detailed below. This model has been coupled with adequate Henry constants of the gas species, CO₂ and hydrogen, also collected from the literature [25].

The amount of carbon dioxide and hydrogen present in the aqueous dissolution of the system depends on its solubility in the vapor-liquid equilibrium according to Eqs 1, 2. This solubility was calculated as a function of temperature, pressure and gas composition with the corresponding Henry constants [25]. The equilibrium acid-base reactions involved in the aqueous dissolution are: C O 2   g ↔ C O 2   a q (1) H 2   g ↔ H 2   a q (2) H 2 O + C O 2   a q ↔ H C O 3 − + H + (3) H C O 3 − ↔ C O 3 2 − + H + (4) H 2 O ↔ H + + O H − (5)

Eqs 3, 4 represent the dissociation of the carbonic acid, and Eq 5 the autoionization of water. The respective equilibrium constants (K) are defined as follows: K 1 = m H + γ H + m H C O 3 − γ H C O 3 − m C O 2 γ C O 2 a H 2 O (6) K 2 = m H + γ H + m C O 3 2 − γ C O 3 2 − m H C O 3 − γ H C O 3 − (7) K w = m H + γ H + m O H − γ O H − a H 2 O (8) Here m represents molality, γ molal activity coefficient, and a H 2 O is molal activity of water.

The ionic product of water (8) was calculated following the model developed by [15], shown in Eq 9; this model was selected because of its wide T-P range of applicability. The parameters from A to G can be found in [22]. log K w = A + B T − 1 + C T − 2 + D T − 3 + E + F T − 1 + G T − 2 log ρ w (9)

The other equilibrium constants (6) and (7) were obtained with the following Eq (10): ln ⁡ K = a 1 + a 2 T + a 3 T − 1 + a 4 T − 2 + a 5 ⁡ ln ⁡ T + a 6 T − 1 + a 7 T − 2 + a 8 T − 1 ⁡ ln ⁡ T P − P s + a 9 T − 1 + a 10 T − 2 + a 11 T − 1 ⁡ ln ⁡ T P − P s 2 (10)

Where P_S is the saturation pressure of water calculated through the IAPWS-IF97 [26]. Parameters a₁ to a₁₁ were provided by [22]; these parameters were obtained by fitting experimental data.

To find the equilibrium of the different species: cations, anions, and neutral molecules of the system, it is necessary to calculate the activity coefficient ( γ ) in the electrolytic dissolutions. Thus, the activity coefficients of the different species were obtained with the model developed by Pitzer and co-workers [27–32]. The model is based on the Debye-Hückel model. For the activity coefficient of hydrogen, the Eq 20 listed below, presented by [33], was used, which estimates the molality of hydrogen in aqueous sodium chloride solutions (0–5 mol/kg) at temperatures and pressures of 273.15–373.15 K and 0–230 bars; this equations is also based on the Pitzer model.

All the different Virial Pitzer parameters ( β 0 , β 1 , C φ , λ , θ , and ψ ) for the mixture H₂O-CO₂-NaCl at different temperatures and pressures were obtained from the work of [22]; these parameters were used to solve the equilibrium compositions of the system H₂O-CO₂-NaHCO₃-NaCl-H₂. For the case of hydrogen, the parameters were obtained from [33]. Besides, experimental data on activity coefficients of hydrogen and other gases in different aqueous salt solutions can be found in [34].

The equations presented below describe the proposed model based in the model of Pitzer and co-workers [27–32] for the studied mixture. The Eqs 16–18 describe the activity coefficients, where M, X, and N in equations represent specific cations, anions, or neutral molecules, respectively. Eq. 19 calculates the water activity and with Eq. 20 the activity coefficient of hydrogen is obtained. Eqs 21–27 show how those parameters were obtained, the remaining parameters used to determine the activity coefficients were calculated with the polynomial functions of T and P developed by [22]. A φ = 1 3 2 π N A ρ w 1000 1 2 e 2 4 π ϵ 0 ϵ r k T 3 2 (11) I = 1 2 ∑ m i z i (12) Z = ∑ i m i z i (13) F = − A φ I 1 2 1 + 1.2 I 1 2 + 2 1.2 ln 1 + 1.2 I 1 2 + ∑ c = 1 N c ∑ a = 1 N a m c m a B c a ′ + ∑ c = 1 N c − 1 ∑ c ′ = c + 1 N c m c m c ′ Θ c c ′ ′ + ∑ a = 1 N a − 1 ∑ a ′ = a + 1 N a m a m a ′ Θ a a ′ ′ (14) ϕ − 1 = 2 ∑ i m i − A φ I 3 2 1 + 1.2 I 1 2 + ∑ c = 1 N c ∑ a = 1 N a m c m a B c a φ + Z C c a + ∑ c = 1 N c − 1 ∑ c ′ = c + 1 N c m c m c ′ Θ c c ′ ϕ + ∑ a = 1 N a m a ψ c c ′ a + ∑ a = 1 N a − 1 ∑ a ′ = a + 1 N a m a m a ′ Θ a a ′ ϕ + ∑ c = 1 N c m a ψ a a ′ c + ∑ n = 1 N n ∑ c = 1 N c m n m c λ n c + ∑ n = 1 N n ∑ a = 1 N a m n m a λ n a + ∑ n = 1 N n ∑ c = 1 N c ∑ a = 1 N a m n m c m a ζ n c a (15)

The equations for determining the activity coefficients of the various ions and neutral molecules are given below, in addition to the equation for calculating the activity of water. Ln ⁡ γ M = z M 2 F + ∑ a = 1 N a m a 2 B M a + Z C M a + ∑ c = 1 N c m c 2 Θ M c ϕ + ∑ a = 1 N a m a ψ M c a + ∑ a = 1 N a − 1 ∑ a ′ = a + 1 N a m a m a ′ ψ a a ′ M + z M ∑ c = 1 N c ∑ a = 1 N a m c m a C c a + 2 ∑ n = 1 N n m n λ n M + 6 ∑ n = 1 N n ∑ a = 1 N a m n m a ζ M n a (16) ln ⁡ γ X = z X 2 F + ∑ c = 1 N c m c 2 B c X + Z C c X + ∑ a = 1 N a m a 2 Θ X a + ∑ c = 1 N c m c ψ X a c + ∑ c = 1 N c − 1 ∑ c ′ = c + 1 N c m c m c ′ ψ c c ′ X + z X ∑ c = 1 N c ∑ a = 1 N a m c m a C c a + 2 ∑ n = 1 N n m n λ n X + 6 ∑ n = 1 N n ∑ c = 1 N c m n m c ζ n c X (17) ln ⁡ γ N = 2 ∑ n = 1 N n m n λ N n + 2 ∑ c = 1 N c m c λ N c + 2 ∑ a = 1 N a m a λ N a + ∑ c = 1 N c ∑ a = 1 N a m c m a ζ N c a (18) a H 2 O = exp − M H 2 O ϕ 1000 ∑ j m j (19) ln ⁡ γ H 2 = 2 ∑ c = 1 N c m c λ H 2 c + 2 ∑ a = 1 N a m a λ H 2 a + ∑ c = 1 N c ∑ a = 1 N a m c m a ζ H 2 c a (20)

The different Virial parameters of the Pitzer model are described in the following equations: B c a ′ = 2 β c a 1 α 2 I 2 − 1 + 1 + α I 1 2 + 1 2 α 2 I exp ⁡ ⁡ − α I 1 2 (21) B c a φ = β c a 0 + β c a 1 ⁡ exp ⁡ ⁡ − α I 1 2 (22) B M X = β c a 0 + 2 β c a 1 α 2 I 1 − 1 1 + α I 1 2 exp ⁡ ⁡ α I 1 2 (23) C c a = C c a φ 2 z M z X 1 2 (24) Θ i j ϕ = θ i j + θ i j . E I + I θ i j ′ . E I (25) Θ i j = θ i j + θ i j . E I (26) Θ i j ′ = θ i j ′ . E I (27)

Where θ i j is the second Virial parameter for each pair of cations or each pair of anions. The terms θ i j . E I and θ i j ′ . E I describes the electrostatic effects of unsymmetrical ion mixtures [35]. The Eqs 25–27 were obtained by [30]. The most important characteristic of θ i j . E I and θ i j ′ . E I is that they only depend on the charges of the ions and the total ionic strength ( I ). Furthermore, they do not represent an extra parameterization of the model. θ i j . E I and θ i j ′ . E I are set to zero when the ions have the same charge.

Some Virial parameters were set to zero; for instance, the second virial parameter ( β c a 0 , β c a 1 , and C c a φ ) of the following combinations: H⁺-OH^-, H⁺-HCO₃ ^-, and H⁺-CO₃ ^2-, because their influence on the activity coefficient is small. Furthermore, the third virial coefficients ( θ ′ s and ψ ′ s ) related to H⁺, OH⁻, HCO₃ ⁻, and CO₃ ^2- also can be fixed to zero due to their minor impact on free energy. Besides, the third-order interaction ( ζ ) of CO₂-Na⁺-HCO₃ ^- also can be set to zero because of the small contribution of the third virial coefficient. These assumptions were taken before by other researchers [22, 36].

It must be emphasized that the parametrization of the model used in this work, based on the work of [22], considered all species present in the hydrothermal conversion process, except hydrogen. Due to its nature and its relatively low concentration in the solution, it is not expected that hydrogen has a significant impact on the activity coefficients of other compounds in the solution. However, the inverse may not be true, and in particular the dissolved CO₂ species, that are not included in the parametrization of the activity coefficient model for hydrogen, may have a certain influence on this activity coefficient that is not considered in the model used.

Finally, the pH was calculated as the product between the concentration and the activity coefficient while the pOH as the difference between pK_w and pH. p H = − log m H + · γ H + (28) p O H = p K w − p H (29)

The speciation equilibrium was represented by a nonlinear problem formed by four equations (ecs. 30–33 listed below) and four variables: mH⁺, mOH⁻, mHCO₃ ^-, and mCO₃ ^2-. The thermodynamic system was calculated and solved using the commercial computing software MATLAB^® R2022b. The program was divided into three main parts: auxiliary functions, equilibrium functions, and a solution file. The flow diagram of the code can be seen in Figures 1, 2.

The auxiliary functions consist of five function codes that resolve different parts of the system: equilibrium constants, activity coefficients, Henry’s Law, saturation pressure of water and water density. The equilibrium constant function solves Eqs 6–8. The activity coefficient function computes the activity coefficient of every specie in the dissolution through the Pitzer model according to the Eqs 11–19. The density function returns the density of water according to the PC-SAFT equation of state [37], using an open-source program from [38].

The equilibrium functions file contains the four functions: obj_F (1), obj_F (2), obj_F (3), and obj_F (4), which are described in equations from (30) to (33); these are the function that must be solved by the computing software to calculate the thermodynamic equilibrium. The functions obj_F (1), obj_F (2), and obj_ F (3) represent the difference between the reference equilibrium constants calculated with Eqs 9, 10; and the named in Eqs 6–8 after dividing by its respective reference K to scale the function and facilitate the convergence of the numerical resolution method; and the fourth equation obj_F (4) describes the charge balance between the different ions, again scaled to facilitate convergence.

To solve the nonlinear equation system, formed by the four equilibrium functions named above, the fsolve function of MATLAB^® was used; this command takes an initial vector as initial guess and iterates until made zero the functions established in the equilibrium functions file by minimizing the sum of squares of the functions. Then the solution obtained is taken as the initial point for the next iterations. Some of the initial points were taken from the condition of pure water or water and carbon dioxide mixtures. o b j _ F 1 = 1 − K 1 c a l c K 1 (30) o b j _ F 2 = 1 − K 2   c a l c K 2 (31) o b j _ F 3 = 1 − K w   c a l c K w (32) o b j _ F 4 = m H + + m N a + − m O H − − m C l − − m H C O 3 − − 2 m C O 3 2 − / m H + + m N a + (33)

The algorithm used by the MATLAB^® function fsolve to find solutions to the nonlinear systems is the trust-region-dogleg algorithm: fsolve tries to solve a system of equations by minimizing the sum of squares of the components. If the sum of squares is zero, the system of equations is solved. The trust-region-dogleg algorithm is an iterative method for solving nonlinear systems of equations. The algorithm starts with an initial approximate solution and then uses a trust region to gradually adjust the solution until a satisfactory solution is found. The trust region is a set of possible solutions that are considered “reasonably close” to the true solution. At each iteration, the algorithm uses a root-finding method to find the closest solution within the trust region. If this solution is satisfactory, the process is stopped, and the solution is returned. If it is not satisfactory, the trust region is adjusted, and the process is repeated. The root-finding method used by the trust-region-dogleg algorithm is the dogleg method. This method combines two simpler methods: the Newton method and the gradient direction method.

The Newton method uses information about the curvature of the function to find a faster solution. However, this method can fail if the curvature of the function changes a lot from one iteration to the next. The gradient direction method, on the other hand, uses information about the slope of the function to slowly move in the direction of the solution. This method is more stable, but also slower. The dogleg method combines these two methods to achieve a balance between stability and speed. It uses the Newton method when possible, and the gradient direction method when necessary.

Several parameters of interest for the hydrothermal reduction of CO₂ were analyzed with the thermodynamic model: temperature, pressure, concentration of bases and concentration of sodium chloride (which would be applicable in case that CO₂ dissolved as bicarbonate in sea water is used as feedstock for the reduction process).

A thermodynamic model of the equilibrium species in the hydrothermal conversion of CO₂ dissolved in aqueous solutions has been presented. The influence of different process conditions on the equilibrium concentrations has been analyzed, including: the effect of pressure and temperature, sodium bicarbonate concentration, sodium hydroxide concentration and sodium chloride concentration, with a focus on the concentration of bicarbonate, which is the main reacting species, in the equilibrium achieved under each of these conditions. It has been observed that high concentrations of sodium bicarbonate, which in turn favor the performance of the reaction, are favoured by moderate temperatures, high initial concentrations of sodium bicarbonate and moderate initial concentrations of sodium hydroxide. The presence of sodium chloride in the range of typical concentrations in natural seawater has a negligible influence on the equilibrium concentrations of dissolved species. These results provide valuable guidelines for the development and optimization of hydrothermal CO₂ conversion processes.