Simultaneous Calculation of Chemical and Isotope Equilibria Using the GEOCHEQ_Isotope Software: Oxygen Isotopes

The GEOCHEQ_Isotope software package, elaborated previously for modeling chemical and carbon isotope equilibria in hydrothermal and hydrogeochemical systems by minimizing the Gibbs energy, is extended to the simultaneous calculation of carbon and oxygen isotopic effects. Similar to what was done for carbon, the β-factor formalism was used to develop algorithms and a database for calculating the isotopic effects of oxygen. According to the developed algorithm, the Gibbs energy of formation of a rare isotopologue, G*(P, T), is calculated through the Gibbs energy of formation of the main isotopologue, the value of the β18O factor of this substance, and the mass ratio of the rare (18O) and main (16O) isotopes. The isotope mixture is assumed to be ideal. The temperature dependence of the β-factor is unified as a polynomial in reciprocal absolute temperature. Necessary information on oxygen isotope equilibria involving important geochemical compounds was critically analyzed, and the available data were reconciled and modified. The temperature dependences of the β18O-factors were correspondingly optimized. The thermodynamic database was updated by adding information on the temperature dependence of β18O-factors specified by polynomial coefficients for each substance. The usage of the GEOCHEQ_Isotope software package and the corresponding database is demonstrated by modeling the dependence of oxygen and carbon isotope fractionation factors on the acidity of the solution (pH) in a carbonate hydrothermal system. The simulation results are in a good agreement with experimental data available from the literature. The enrichment of dissolved carbonates in the 18O heavy oxygen isotope relative to water decreases with increasing pH of the system. At the same time, a pH increase results in a decrease in the negative carbon isotope shift between calcite and dissolved carbonates. At high pH values (~11), the isotope shift inversion and the enrichment of the dissolved carbonate in the heavy carbon isotope relative to calcite are predicted.


INTRODUCTION
Studies of the isotopic effects of oxygen is of paramount importance in isotope geochemistry. The great majority of isotope geothermometers make use of oxygen isotopes. This paper is devoted to the simultaneous computation of chemical and oxygen isotope equilibria with the GEOCHEQ_Isotope software package. The principal approaches to and methods of calculating oxygen isotopic effects utilized in this study are generally analogous to those we have applied when studying carbon isotopic effects. The first ever simultaneous simulation of chemical and isotope equilibria by minimizing the Gibbs free energy (this method is referred to as the method of isotope-chemical system, and its idea was formulated by B.N. Ryzhenko) was carried out for sulfur isotopes (Bannikova et al., 1987;Grichuk, 1987Grichuk, , 1988Grichuk, , 2000Grichuk and Lein, 1991). The method enables calculating chemical and isotope equilibria simultaneously, in contrast to the usual approaches (Ohmoto, 1972), in which first chemical and then isotope equilibria are evaluated.
As was pointed out in (Mironenko et al., 2018), one of the main reasons why simultaneous simulations of chemical and isotope equilibria are still used inadequately little is the absence of internally consistent thermodynamic databases on rare isotopologues. One of the tasks of our study is to bridge this gap by developing a database for simulating geochemical processes with regard to oxygen isotopic effects, based on critical analysis of information on isotope equilibria.

APPROACHES AND METHODS
Similar to what has been done in studying carbon isotopic effects (Mironenko et al., 2018), the simulations used herein were conducted with a software package (Mironenko et al., 2020) that utilizes the SUPCRT92 database (Johnson et al., 1992), which has been later updated and corrected. The GEOCHEQ software package minimizes the Gibbs energy by the algorithm of the convex simplex (de Capitani and Brown, 1987). The isotopic effects of an element are calculated within the GEOCHEQ_Isotope extension of the GEOCHEQ software package by using isotopes of the element as independent components (instead of the chemical element as whole), and the list of substances is appended with the corresponding isotopologues (Mironenko et al., 2018). For oxygen, these are 16 O and 18 O and their isotopologues.
All isotopic effects are computed in GEO-CHEQ_Isotope in the approximation of an ideal mixture of isotopes, in which isotope substitutions in different sites do not affect one another, and the isotopic effects are independent of the multiplicity of the isotope substitution. For molecules containing several atoms of a given element, this approximation is commonly formulated as the rule of geometric mean (Galimov, 1973(Galimov, , 1974. The approximation by an ideal mixture is valid for all elements except hydrogen at temperatures above 100 K (Polyakov, 1993). This approximation is deemed to be valid by default in the overwhelming majority of geochemical simulations, if no clumped isotopic effects are computed (Grichuk, 2000;Hill et al., 2014Hill et al., , 2020Schauble and Young, 2020). Similar to carbon isotopic effects, the calculations in this study were carried out using the formalism of the β-factor. The concept of the β-factor was coined by Ya.M. Varshavskii and S.E. Vaisberg (1957) on the basis of the previously suggested (Bigeleisen and Mayer, 1947) reduced isotopic partition function ratio. However, it was E.M Galimov (1973Galimov ( , 1974Galimov ( , 1982 who introduced the concept of the β-factor into the practice of geochemical research. It is convenient to express the β-factor of a substance that contains atoms of the chemical element of interest in equivalent settings in a logarithmic form, in terms of the Gibbs energy (1) where β is the β-factor, G is the Gibbs energy, P is the pressure, T is the temperature, R is the absolute gas constant, z is the multiplicity of the isotope substitution, m is the mass of the isotope, and the heavy isotopes are marked with *. According to this, the equilibrium isotope fractionation factor (α A/B ) between two substances A and B is equal to the ratio of their β-factors or, in the logarithmic form for the isotope shift (Δ A/B ) (2b) In (2b), |β -1| 1. If atoms of an isotopically substituted element occur in nonequivalent settings, equations (1) and (2a) should be refined. Because natural isotope exchange processes involve all of the isotopologues, one should know how the β-factor of a compound at whole relates to the site-specific β-factors for isotope substitutions in certain nonequivalent positions. The theory of intermolecular isotopic effects developed by E.M. Galimov (1973E.M. Galimov ( , 1974E.M. Galimov ( , 1981E.M. Galimov ( , 1082E.M. Galimov ( , 2006 for carbon isotopes is universal and can be applied to the fractionation of oxygen isotopes, because geochemical studies often deal with substances containing several oxygen atoms and some of them are in nonequivalent positions 1 . In particular, E.M. Galimov has demonstrated that, in the case of rare isotope substitutions, the β-factor of the compound is generally related to the β-factors of monosubstituted isotopologues (β i ) by the relation (Galimov, 1973(Galimov, , 1981(Galimov, , 1982: Relation (3a) was referred to by E.M. Galimov as the first additivity law (Galimov, 1981;1982). 2 It has been later demonstrated that relation (3a) is valid not only at rare isotope substitutions but also when β i is insignificantly different from unity (Polyakov, 1987;Polyakov and Horita, 2021). Also, the logarithmic form is slightly more precise (Polyakov and Horita, 2021).
Because the intramolecular effects of oxygen are not computed in this GEOCHEQ_Isotope version, the use of the β-factor of the substance as a whole according to the Eq. (3b) makes it possible to minimize the increase in the database due to introduction of data on the isotopologues of various multiplicity and localization. Hence, in calculating an equilibrium fractionation of oxygen isotopes with the GEO-CHEQ_Isotope complex, it is sufficient to consider the 18 O isotopologue along with the most abundant 16 O one, and a change in the Gibbs energy per one atom of this 18 O isotopologue calculated according to (1), analogously to what has been earlier done for the fractionation of 13 C/ 12 C carbon isotopes (Mironenko et al., 2018). 1 The term intrastructural isotopic effect is sometimes used for intramolecular isotopic effects in solids (e.g., Ustinov et al., 1988). 2 The second additivity law is the method of the isotopic bond numbers (Galimov, 1973(Galimov, , 1981(Galimov, , 1982. The oxygen β-factors insignificantly depend on pressure (Polyakov and Kharlаshina, 1994), and hence, only their temperature dependences are incorporated into the GEOCHEQ_Isotope software. The temperature dependences of the β 18 O-factors are presented in the form of polynomials in the reciprocal absolute temperatures, as has been done for the carbon β-factors (Mironenko et al., 2018) (4) where , and A i is the coefficients of polynomial (4), which are given in Table 1.
Another additions to the computation of chemical and isotope equilibria in system involving the 18 O isotope is the following modification of the formulas for the molalities of the solute components in view of the occurrence of the H 2 18 O isotopologue (5) where m i is the molality of component i, and n i is the number of moles of component, and 55.51 is the number of moles corresponding to 1 kg of water.
β-FACTORS IN THE GEOCHEQ_ISOTOPE THERMODYNAMIC DATABASE Similar to what pertains to the carbon isotopic factors (Mironenko et al., 2018), the accuracy and justification of the calculations of the oxygen isotope equilibria depend, first of all, on the accuracy of data on the β-factors utilized in the GEOCHEQ_Isotope software package. Equilibrium oxygen isotope fractionation factors can be evaluated experimentally and/or theoretically. Equilibrium oxygen isotope fractionation factors are experimentally determined using the Northrop-Clayton partial isotope exchange method (Northrop and Clayton, 1966) and the three-isotope method (Matsuhisa et al., 1978). Oxygen isotope geothermometers are also calibrated using empirical techniques (Clayton et al., 1968;Ustinov et al., 1978;Macey and Harris, 2006). A review of these techniques and their advantages and disadvantages can be found in (Chacko et al., 2001). The theoretical calculations are conducted based on model vibrational spectra (Urey, 1947;Bigeleisen and Mayer, 1947;Kieffer, 1982), the interatomic potential (Patel et al., 1991), and ab initio principles, based on the density functional method (Schauble, 2004;Blanchard et al., 2009;Hill et al., 2014;Schauble and Young, 2021). Ab initio computations lately became the most widely used. The oxygen β-factors of many minerals were calculated with the application of the semiempirical method of increments (Schütze H., 1980;Zheng, 1991;1999), which, in essence, is E.M. Galimov's method of isotopic bond numbers (Galimov, 1973(Galimov, , 1974  have, however, to be treated with care, because they often come in conflict with both experimental data and theoretical calculations (Horita and Clayton, 2007).

Oxygen Isotope Fractionation Factors in the System H 2 O-CO 2 (gas)-dissolved Inorganic Carbon
The great majority of experimental determinations of the equilibrium oxygen isotope fractionation factors were conducted with the involvement of H 2 O and CO 2 (Chacko et al., 2001), and hence, an accurate choice between the β 18 O-factors of H 2 O and CO 2 is of principal importance for the calibration of the entire set of data on oxygen isotope equilibria.
β 18 O-factors of H 2 O. The oxygen β-factors of water vapor were calculated using various theoretical approaches. The computations on the basis of spectroscopic data and molecular constants were made with regard to the anharmonicity and interaction between oscillations and rotations (Bron et al., 1973;Richet et al., 1977). Recent ab initio computations were carried out based on the density functional (Hill et al., 2014(Hill et al., 2020Guo and Zhou, 2019;Schauble and Young, 2021). These calculations were conducted using diverse techniques, and their results are in reasonably good mutual agreement. The coefficients of approximating polynomial (4) used in the GEOCHEQ_Isotope software are listed in  (Horita and Wesolowski, 1994) β 18 O-factors of CO 2 and dissolved inorganic carbon. The β-factors of CO 2 were computed by analogy with what has been done for water: from spectroscopic data and molecular constants (Richet et al., 1977) and ab initio principles (Hill et al., 2014(Hill et al., , 2020Guo and Zhou, 2019;Schauble and Young, 2021). The data were then used to calculate the coefficients of polynomial (4) for O 2(gas) ( Table 1). Also, experimental data are available on the equilibrium oxygen isotope fractionation between CO 2(gas) and H 2 O (liq) at 25°C, because these values are used in determining water isotope composition. Numerous experiments (see references in Chacko et al., 2001) enabled the IAEA Consultants Group to recommend the value of 1.0412 at 25°C (Hoefs, 2018) as that of equilibrium fractionation coefficient in the system CO 2(gas) -H 2 O (liq) . Estimates are also now available for the temperature ranges of 5-100°C (Brenninkmeijer et al., 1983) and 130-350°C (Truesdell, 1974). Comparison of the equilibrium isotope fractionation factors in the system CO 2(gas) -H 2 O (liq) calculated using GEOCHEQ_Isotope data (Table 1) and corresponding experimental values are exhibited in Fig. 1; the experimental and calculation data are obviously in reasonably good agreement.
for the calculation of the temperature function of the carbon

Oxygen β-Factors of Carbonate Minerals
Fractionation of oxygen isotopes in the system calcite-CO 2 . From the standpoint of isotope geochemistry, calcite seems to be the most important carbonate mineral. The reason for this is that the fractionation − 3 HCO − 2 3 CO coefficients of oxygen isotopes experimentally determined in the anhydrous system (i.e., without H 2 O fluid) link to calcite (Chacko et al., 2001;Clayton and Kieffer, 1991). In view of this, it is very important to estimate equilibrium distribution in the system calcite-CO 2 because, as is demonstrated above, the β 18 O-factor of CO 2 has been reliably determined. For calcite, the ab initio calculations (Hill et al., 2014(Hill et al., , 2020Guo and Zhou, 2019;Schauble and Young, 2021) and those on the basis of the model oscillation spectra of the Kieffer type (Clayton and Kieffer, 1991;Chacko and Deines, 2008;Polyakov, 2008) agree well to each other. In GEOCHEQ_Isotope, the values of the β 18 O-factor are assumed equal to the ab initio computation results. At temperatures above 400 K, these data coincide with estimates in (Clayton and Kieffer, 1991). The difference between the values assumed for the β 18 O-factor of calcite in GEO-CHEQ_Isotope and calculated based on the Kieffer spectrum model in a quasiharmonic approximation (Polyakov, 2008) are also much smaller than the scatter of the experimental data (Fig. 2). A somewhat greater difference was found between the CaCO 3 /CO 2 fractionation curve (Chacko and Deines, 2008) and the fractionation curve calculated using the GEO-CHEQ_Isotope database. This difference is explained by the underestimation of the Einstein oscillator fre-    (Mironenko et al., 2018), which leads to underestimated values of the calcite β-factor. This effect is smaller for the β 18 O-factor of calcite than for the β 13 C-factor. Equilibrium fractionation in the system calcite-CO 2 that was calculated using the values of the β 18 O-factor for calcite from the GEO-CHEQ_Isotope database is generally in a good agreement with the experimentally measured equilibrium oxygen isotope fractionation factor for this system (Fig. 2). Effects of oxygen isotope fractionation between carbonate minerals. The method applied to calculating the coefficients of polynomial (4) for carbonate minerals was analogous to the method used for their β 13 Cfactors (Mironenko et al., 2018). Considering that the fractionation between carbonate minerals is well described by the model (Chacko and Deines, 2008) but the β 18 O-factors of calcite are underestimated compared to the values assumed in GEOCHEQ_Isotope (Fig. 2), the difference between the values of 10 3 lnβ of calcite accepted in the GEOCHEQ_Isotope database and those in the model (Chacko and Deines, 2008) was added to the values of 10 3 lnβ for carbonate minerals (6) where β carb and β calc are the β 18 O-factors for the carbonate mineral and calcite in the GEOCHEQ_Isotope database, and and are the same factors calculated by the model (Chacko and Deines, 2008). The coefficients of polynomial (4) for the β 18 Ofactors of carbonates calculated according to (6) are given in Table 1. Figure 3 makes it possible to compare the fractionation of oxygen isotopes in the system dolomite-water according to data in GEOCHEQ_Isotope and experimental calibrations of other authors. The experimental values of the oxygen equilibrium fractionation factors between dolomite and water and the calculation data found using the GEOCHEQ_Isotope database are obviously in a good agreement. GEOCHEQ_Isotope Polyakov, 2008Chacko, Deines, 2008Chacko et al., 1991Dove et al., 1992Coplen, 1959 1 0 10 6 /T -2 (K -2 ) 12 14 4 6 8 Jimenez-Lopez et al., 2001Kim, O'Neil, 1997Rosenbaum, 1994Zhou, Zheng 2003Tarutani et al., 1969 No. 11 2021 POLYAKOV et al.  25  50  100  175  300  450  700   12  3  6  9 Horita, 2014 Northrop, Clayton, 1966Mattew, Katz, 1977Vasconcelos et al., 2005Schmidt et al., 2005Clayton et al., 1968 Calcte-water GEOCHEQ_Isotope topes in silicates and oxides. The coefficients of polynomial (4) for the calculation of the β 18 O-factors of these minerals, which were incorporated into GEO-CHEQ_Isotope, are presented in Table 1. Most of these data were obtained by means of ab initio calculations. On the one hand, these results agree well with experimental and empirical data, and on the other hand, they enable extrapolations of the results to conditions inaccessible in the experiments. Some results obtained by various researchers with the use of different methods must be reconciled. An example of this is the calibrations of the β 18 O-factors of TiO 2 polymorphs: rutile, anatase, and brookite. The β 18 O-factors of these polymorphs computed ab initio with the application of the supercell method (Krylov and Kuznetsov, 2019) are in very good agreement with the empirical calibration of the geothermometers (Zack and Kooijman, 2017). However, the β 18 O-factors of rutile are in notably poorer agreement with experimental and theoretical estimates of the oxygen isotope fractionation factors between the rutile-quartz and rutile-calcite (Fig. 4) than the calibration for rutile calculated based on a Kieffer model spectrum (Clayton and Kieffer, 1991). We have reconciled the various calibrations using an approach analogous to that applied for carbonates. We estimated the difference between the β 18 O-factors obtained for rutile by these two methods and then corrected the β 18 O-factors of the TiO 2 polymorphs for this value and thus retained the calibration of the geothermometers and ensured reasonably good agreement with experimental and empirical estimates for the rutile-quartz and rutilecalcite isotope fractionation coefficients. When selecting the calibrations of the β 18 O-factors for the GEOCHEQ_Isotope database, we have paid much attention to the consistency of the data. It is the consistency of the estimates of the β 18 O-factors obtained by different techniques (or the possibility of making these data mutually consistent) that has been applied as the criterion of the reliability of the data.
The β 18 O-factors of oxides were determined, along with conventional techniques mentioned above, also by some new approaches. In this situation, it is particularly important to verify the results obtained by the new techniques by using conventional techniques, which have already proven to be reliable. For oxides whose cations possess a "Mössbauer" isotope, β 18 Ofactors can be estimated from experimental data on the gamma-resonant scattering and heat capacity. This method was applied to measurements of the β 18 O-factor for hematite (Polyakov et al., 2001) and cassiterite . For cassiterite, this test has been done experimentally, using the Norton-Clayton method  and by the method of empirical calibration with the use of other isotope geothermometers and temperature estimation from data on gas-liquid inclusions (Macey and Harris, 2006). Similar to what has been done for hematite, the test can be conducted by comparing the β 18 O-factors of hematite calculated ab initio and derived from experimental data on the temperature shift in the Mössbauer spectra and heat capacity. These estimates for the β 18 O-factor of hematite practically exactly coincide: even at room temperature, the difference in 10 3 lnβ is ~0.3 (0.4% of the 10 3 lnβ value), and this difference further diminishes with increasing temperature.
Analogously, data can be made consistent for grossular, for which two calculations are available (Krylov and Glebovitsky, 2017;Schauble and Young, 2020), both based on the density of states functional approach but obtained using different methods. The consistency of these data is slightly poorer (~3.3% of the 10 3 lnβ value) but is still acceptable (Fig. 5). Both computations are in a good agreement with the experimental fractionation coefficients for grossular-calcite (Rosenbaum andMattey, 1995, 1996) and grossular-aqueous supercritical fluid. The values of the β 18 O-factor of grossular accepted in GEOCHEQ_Isotope obtained by averaging the two calibrations (Table 1).
The usage of β-factors in databases requires the values of these factors computed for broad tempera-ture ranges. This is most commonly done in theoretical calculation studies, but sometimes no such calculation data are available. In such situations, data obtained by various authors can be made consistent. This situation occurs with the β 18 O-factors for magnetite. Reliable experimental data are available for the equilibrium fractionation of oxygen isotopes in the system magnetite-H 2 O at temperatures higher than 300 o C (Cole et al., 2004), and there are unsystematic scarce data on isotope fractionation in the system magnetite-H 2 O at low temperatures (Blattner et al., 1999;Zhang et al., 1995;Mandernack et al., 1999). We used these experimental data to plot an interpolation polynomial of form (4) for the magnetite-H 2 O fractionation curve (Fig. 6). To obtain a polynomial for the β 18 O-factor of magnetite, the coefficients of the polynomial for water (liquid or supercritical fluid, T > 647 K) from Table 1 were added to the coefficients of the polynomial for the fractionation curve. The polynomials thus obtained were used to plot the temperature dependence of the β 18 O-factor of magnetite. This curve was smoothed by approximating it again using a polynomial of form (4). The results are summarized in Table 4.
Above, we described the main principles, methods, and techniques used to reconcile the data and compile a base of natural isotopic data in GEOCHEQ_Isotope. Detailed discussions of particulars of some related issues is beyond the scope of this paper and Rutile-calcite, Chacko et al., 1996Rutile-calcite, Clayton and Kieffer, 1991Rutile-calcite, Krylov and Kuznetsov, 2019Rutile-quartz, Agrinier, 1991Rutile-calcite, Clayton and Kieffer, 1991Rutile-quartz, Krylov and Kuznetsov, 2019  shall be rather referred to a manual for the user of this software package.

SIMULTANEOUS COMPUTATIONS OF EQUILIBRIUM OXYGEN AND CARBON ISOTOPE EFFECTS AND CHEMICAL COMPOSITION IN A HYDROTHERMAL SYSTEM
The GEOCHEQ_Isotope software now enables one to simulate chemical processes with regard to oxygen and carbon isotope effects. As an illustrative example of this simulation, we chose a system, which is analogous to that described in (Mironenko et al., 2018). The calculation method of chemical equilibria with regard to both oxygen and carbon isotope equilibria is generally analogous to those applied to models that involved estimates only of equilibrium carbon isotopic effects (Mironenko et al., 2018(Mironenko et al., , 2021Sidkina et al., 2019). A natural change was the addition of the 18 O isotope as an independent component and 18 O isotopologues to the list of substances subject to Gibbs energy minimization. The molality of the components of aqueous solutions was calculated by formula (5).
As an example of the application of GEO-CHEQ_Isotope, we have computed chemical and isotope equilibria in CO 2 -bearing hydrothermal system of the bulk composition 1 mole CO 2 + 1 kg H 2 O + 1.4 moles CO 2 at 15, 25, and 40 o C and a pressure of 40 bar, depending on pH. The 18 O/ 16 O and 13 C/ 12 C isotope ratios were specified in the system according to the natural abundances of these isotopes. We considered 16 O-, 18 O-, 12 C-, and 13 C-isotopologues of H 2 O, CO 2 (aq), , ,CaCO 3 (aq), , CO 2 (gas), CH 4 (gas), and H 2 O(gas). For this system, the dependence of the composition of dissolved inorganic (carbonate) carbon on the equilibrium fractionation of oxygen isotopes between solute carbonate species and water depending on pH have been studied experimentally and theoretically (Beck et al., 2005). These experimental results only insignificantly different from the carbonate system discussed herein. In the GEO-CHEQ_Isotope simulations, we analyzed the possibility of formation of the solute CaCO 3 (aq) and species, whereas these species were ignored in (Beck et al., 2005). At the same time, theoretical modeling in (Beck et al., 2005) involved the solute H 2 CO 3 species. It should be stressed that the H 2 CO 3 concentration was, according to the modeling results, close to 0.2-0.3% of the total concentration of dissolved carbonate carbon at pH < 5 and drastically decreased at increasing pH. The results obtained by GEOCHEQ_Isotope simulation of the composition of the carbonate species are shown in Fig. 7. Our simulation results agree well with earlier results in (Beck et al., 2005). The dominance of one or another solute species of carbonate carbon depending on pH, with regard to the difference in the values of the β 18 O-and β 13 C-factors of the solute carbonate species (Table 1, Fig. 1), predetermines the dependence of the equilibrium fractionation coefficient of oxygen and carbon isotopes between total dissolved inorganic carbon and coexisting phases. Figure 8 shows examples of such dependences simulated with GEOCHEQ_Isotope. For oxygen isotopes, reliable experimental data (Beck et al., 2005) are currently available on the equilibrium fractionation of dissolved carbonate carbon and water at temperatures of 15-40°C. Our simulation results are in good agreement with these experimental data (Fig. 8). The decrease in dissolved carbonates in the heavy oxygen isotope relative to water at increasing pH (Fig. 8)  CO with that the equilibrium fractionation factors between these carbonate species and water (and hence, their β 18 O-factors) decrease in the same sequence (Fig. 1). The situation with carbon isotopes is opposite: the β 13 C-factors of the aforementioned carbonate species occur in the opposite order (Mironenko et al., 2021), and hence, the depletion of dissolved carbonate carbon in the 13 C heavy isotope is reduced as pH is increased. The isotope shift for carbon between dissolved carbonate carbon and calcite at high pH shows an inversion, because the β 13 C-factor for the dominant carbonate species is greater than that of calcite (Fig. 8).

CONCLUSIONS
We present an extension for the database of the GEOCHEQ_Isotope software package. This extension was designed to enable simultaneous simulations of chemical and isotope equilibria by means of minimizing the Gibbs free energy for the calculation of oxygen isotope equilibria. Similar to what has been done with carbon isotopes, this software version is underlain by the β-factor formalism in an approximation of an ideal mixture of isotopes. Retrograde curve at T > 300°C, Cole et al., 2005T > 300°C, Cole et al., 2005Zhang et al., 1997Blattner et al., 1983Mandernack et al., 1999 Interpolation polynomial 0 2 1 0 10 6 /T -2 K -2 12 14 4 6 8  4 5 6 7 8 9 10 11 12 pH 13 C/ 12 C Geocheq_Isotope, 25°C Geocheq_Isotope, 40°C Geocheq_Isotope, 200°C Beck et al., 2005, 15°C Beck et al., 2005, 25°C Beck et al., 2005 The presented version of GEOCHEQ_Isotope uses the GEOCHEQ database (Mironenko et al., 2000) to describes nonisotope thermodynamic properties of substances. The Gibbs energy of the isotopologues is calculated using known values of the β-factors according to Eq. (1), with regard to additivity relations between the lnβ of a compound as a whole and lnβ i of the monosubstituted isotopologues.
In this version, the computational algorithms are adapted to situations when oxygen isotopes are analyzed, and the accuracy of the computations is improved. In particular, the current equation for calculating the molality takes into account the presence of the H 2 18 O isotopologue. Available information on equilibrium fractionation factors of oxygen isotopes has been critically analyzed, and the data have been made mutually consistent. The temperature dependences of the oxygen β-factors are unified and written in the form of polynomials (4) in reciprocal absolute temperature.
The application of the GEOCHEQ_Isotope software package for simultaneous calculations of chemical, oxygen and carbon isotope equilibria was tested by simulating the dependence of the composition of solute carbonate species and the 18 O/ 16 O equilibrium fractionation coefficient between solute carbon and water, depending on the pH of the solution. The data are in a good agreement with experimental results (Beck et al., 2005). The enrichment of dissolved carbonate carbon in the 18 O heavy isotope decreases with increasing pH because of a systematic change of the dominant carbonate specie in the sequence CO 2 (aq) → → , which is characterized by a decrease in the values of the β 18 O-factor ( ). For carbon, a pH increase leads to an increase in the fractionation coefficient for dissolved carbonate carboncalcite because . At high pH (~11), the equilibrium isotope shift exhibits an inversion, and the dissolved inorganic carbon enriches in the 13 C isotope relative to calcite.

FUNDING
This study was carried out under government-financed research project 0137-2019-0015 for the Vernadsky Institute and was partly supported by Russian Foundation for Basic Research, project no. 19-05-00865a.
OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the