Introduction

The long-term performance of the materials used in geological disposal of vitrified nuclear wasteforms raises numerous safety concerns.1 For example, over geological timescales, the exposure to groundwater of nuclear waste immobilization glasses may result in the precipitation on the glass surface of secondary crystalline phases, e.g., zeolites.2,3,4 This is significant since the formation of zeolitic secondary phases has been suggested to result in a resumption of glass corrosion, which can significantly reduce the long-term durability of nuclear wasteforms.5,6,7,8 Hence, accurate predictions of the propensity for formation and the stability of zeolitic phases are important to assess the long-term integrity of nuclear wasteforms in sub-surface conditions, e.g., < 250-m depth.9 But, assessing such propensities experimentally is challenging, especially since there is a myriad of factors that affect the formation of zeolites; which themselves feature a multitude of compositions. This is further complicated by the long equilibration times that are necessary to attain equilibrium at sub-boiling temperatures.10 The kinetics of zeolite formation may be accelerated with increasing temperature.11,12,13,14 However, some zeolites are metastable and may exhibit irreversible phase transformations within narrow temperature ranges or under specific synthesis conditions.15,16,17,18 Other zeolites can form new phases through ion-19,20,21,22,23 or water-exchange reactions,24 which increases the complexity of establishing thermodynamic data for these alterable phases. For these reasons, geochemical modeling may provide a comprehensive means to parametrically assess the diverse scenarios and conditions under which zeolites may form in geological repositories.

Geochemical modeling can be used to rapidly estimate the propensity for the formation of minerals in solutions, e.g., using Geochemist’s Workbench,25 PhreeqC,26,27 and Gibbs (free) energy minimization (GEM) solvers.28,29 However, the use of such platforms relies on the availability of comprehensive and accurate thermodynamic data.10 Although calorimetric data for some zeolites30,31,32,33,34,35 are available and empirical methods exist for estimating their thermodynamic properties,36,37,38,39,40,41 broadly speaking, available thermodynamic data remain sparse. This partially arises from the fact that the stability of some zeolites is limited to small domains of solution composition and temperature.42 As such, data obtained at high(er) temperature may not accurately describe low-temperature stability.10 Zeolitic minerals exhibit a wide range of compositions, structures, and hydration states.43 For example, while clinoptilolite and heulandite are isostructural compounds;44 other zeolites feature silicon-to-aluminum molar ratios ranging from 1 (e.g., thomsonite) to 5 (e.g., mordenite);45 while wairakite, leonhardite, and laumontite differ only in the number of zeolitic water molecules per formula unit (2, 3.5, and 4.5, respectively)46 and feature a Si/Ca and Si/Al ratios of 4 and 2, respectively, etc. Moreover, zeolites may coexist with other minerals, e.g., chabazite, mesolite/scolecite, stilbite, and laumontite,47 and other zeolites of differing composition including: ferrierite (a magnesium zeolite), mordenite, clinoptilolite,16 thomsonite, natrolite, gonnardite,15 erionite, and phillipsite.46 Zeolites can also form solid-solutions, e.g., wairakite–analcime,48 or stilbite–stellerite.49 As such, establishing and compiling the thermodynamic properties for the many different compositions of zeolite phases requires the ability to accurately quantify, and discern small differences in their phase stability, and equilibria which, for example, are evidenced by the formation, and co-existence of diverse zeolitic compositions in distinct geological environments as relevant to nuclear waste repositories in which glass corrosion may occur.

Currently, a few zeolite phases are included in thermodynamic databases, e.g., llnl.dat50 and sit.dat in PhreeqC,26,27 ThermoChimie,51 Thermoddem,52 and GEMS,53 which are integrated with solvers that compute stable phase equilibria based on the Law of Mass Action (LMA, e.g., PhreeqC26,27) or Gibbs Energy Minimization (GEM, e.g., GEMS-Selektor29). GEMS offers the advantage of automatically deriving a complete dataset of thermodynamic properties and extrapolating solubility constants at temperatures up to (and above) 100 °C.53 Therefore, by using the GEMS approach and integrating available calorimetry data,54,55,56,57,58,59,60,61 we compiled a unified, internally-consistent thermodynamic database to predict zeolite stability in neutral and alkaline solutions at low temperatures. The database when integrated with the GEMS solver is shown to correctly assess the propensity for zeolite formation, and predicted stability fields for a multiplicity of zeolite compositions viz-a-viz experimental observations.

Results and Discussion

The approach that we followed to develop the new database involved a comprehensive literature survey of measured thermodynamic and solubility data for zeolites (see Supplementary Discussion). In general, for any given zeolite, thermodynamic data; if unavailable experimentally was estimated using volume-based correlations or empirical equations. First, we undertake a detailed uncertainty analysis of the input (thermodynamic) data, and its implications on phase stability predictions. Furthermore, we evaluate the effect of the presumed, and widely indicated incongruent dissolution of zeolites, in terms of complexities that hinder use of measured solution data to assess either (precipitated) zeolite composition or the presence of secondary mineral phases. Second, the internal consistency of thermodynamic data, i.e., vis-à-vis the aqueous and gaseous species in the GEMS database, was examined. When integrated with the GEMS solver (see Methods section) that applies an interior points method (IPM)29 and Kuhn-Tucker conditions for energy minimization convergence,62,63 the data permits the construction of the stability diagrams for zeolites as relevant to nuclear waste disposal applications. Finally, we demonstrate the consistency of the database by constructing stability field diagrams for three systems of interest and comparing the predictions with experimental observations of zeolite precipitation.

Uncertainty in predicted solubility data

Geochemical modeling requires accurate thermodynamic data of solid, liquid, and vapor phases, including their associated dissolved species, complexes, and components. Such thermochemical data can be obtained via calorimetric measurements, solution-phase analysis, calculations, or by a combination of these methods. However, independent of how such data is obtained; experimental thermodynamic data present uncertainties, which impact the accuracy of the overall (phase equilibria) predictions. The sources of uncertainty may arise from the accuracy, repeatability, or reproducibility of a measurement, a calibration reference, equipment parameter, etc. In general, uncertainties in calorimetric data increase from 0.2 to 5% with decreasing temperature31,32,33,34,57,64 and their inter-sample variation can be as high as 7%.65 Also, inferring thermochemical data based on fitting equations can yield uncertainty ranging from 0.1 to 3%.66 Nevertheless, typically-reported uncertainties associated with calorimetric measurements of the entropy and enthalpy of formation range from 0.3 to 1.5% (see Supplementary Table 1). Hereafter, we explore the relationship between the uncertainty of each calorimetric data input and the impact on the predicted solubility (equilibrium) constant by conducting a sensitivity analysis for the example of analcime.

First, we focus on the effect of the Gibbs free energy (ΔfG°, kJ·mol−1) and enthalpy of formation (ΔfH°, kJ·mol−1). Based on the data presented in Supplementary Table 1, these quantities present an average and maximum uncertainty of ±0.3% and ±0.75%, respectively. As shown in Fig. 1a, b, an uncertainty of ± 0.3% in the Gibbs free energy or enthalpy of formation yields an uncertainty of about 1.2–1.7 log units in the predicted solubility constant for analcime. This uncertainty increases to 4.4 log units when the uncertainty of the Gibbs free energy or enthalpy of formation reaches ± 0.75%. Benning et al. noted that small variations in the Gibbs free energy of formation had a significant impact on the estimated solubility constant.67 They estimated that a variation in the Gibbs free energy of formation of just 5–10 kJ·mol−1 (i.e., which corresponds to a variation of 0.08–0.16% for sodium clinoptilolite) results in a variation of more than 1 log unit in the resulting solubility constant.

Fig. 1: A sensitivity analysis conducted for analcime to assess the effect of the uncertainty.
figure 1

Associated with the a Gibbs free energy of formation, b enthalpy of formation, c entropy, and d heat capacity on the estimated solubility constant. In each case, the solid line shows the predicted solubility constant as a function of temperature for an analcime phase with Si/Al = 2.0, while the dashed lines and dotted lines represent the confidence interval associated with the uncertainty of each thermodynamic input.

Full size image

In contrast, we find that uncertainties associated with the measured entropy (S°, J·K−1·mol−1) and heat capacity (Cp°, J·K−1·mol−1) only marginally affect the estimated solubility constant (see Fig. 1c, d). Indeed, reported uncertainties of ±1.5% and ± 3% in the entropy and heat capacity yield uncertainties of only 0.16 and 0.21 log units in the estimated solubility constant of analcime. Even increasing the uncertainty in the entropy to ± 10% yields an uncertainty of 1.25 log units in the solubility constant of analcime. Similarly, an uncertainty of ±10% in the heat capacity results in an uncertainty of about 0.2 log units in the solubility constant estimated for analcime. This a posteriori confirms that the temperature-dependence of the heat capacity can be neglected over small ranges of temperature. Overall, the sensitivity analysis shows that the uncertainties associated with the enthalpy and Gibbs free energy of formation dominantly affect the estimated solubility constants and stability fields of zeolite phases.

Uncertainty in experimental solubility data and the effects of incongruent dissolution

We now examine uncertainties associated with experimentally-obtained solubility data. The uncertainty associated with reaction data obtained from solution experiments (e.g., from under-saturation or over-saturation) is usually smaller than that associated with calorimetry measurements.10,52,67 For instance, the uncertainty associated with experimental estimates of the solubility constant of analcime at 25 °C was found to be around ± 0.3 log units; which decreases with increasing temperature.68 As such, if an uncertainty of ± 5% is estimated in the measured aqueous ion concentrations, this produces a consequent uncertainty of ± 0.1 units in log K and ± 0.5% in the Gibbs free energy of the dissolution reaction.

The determination of the solubility constants via solution-phase analysis requires congruent dissolution of the solid phase, whereby the molar ratios of the solid and liquid phase are similar. But, for some silicates, congruency may not be attained, whereas long equilibration times may be required before congruency is achieved for other compositions.69 Dissolution of zeolites tends to be nonstoichiometric (“incongruent”) with increasing Si/Al molar ratios.70,71 For example, nonstoichiometric molar ratios of aqueous [Si/Al] have been reported for analcime with Si/Al = 2.0 and 2.5, (Ca, K, or Na)-clinoptilolite, mordenite(Ca), natrolite, zeolites X, Y, and P(Ca).53,67,68

Indeed, careful analysis of the solution-phase compositions in Wilkin and Barnes showed that data obtained for analcime with Si/Al = 2.0 in the solid exhibited a Si/Al ≠ 2.0 in solution.68 Other solution data in Wilkin and Barnes68 showed stoichiometric aqueous Si/Alaq = 2.0 at 25 °C whose follow-on analysis yields log K = −16.55. This estimate of the solubility constant differs by 0.5 log K units as compared with solutions which showed nonstoichiometric ratios (Si/Alaq = 1.43, log K = −16.06 and Si/Alaq = 3.38, log K = −16.08); for the same analcime composition with Si/Al = 2.0. This difference in log K decreases further with increasing temperature. Surprisingly, comparison of the solubility constants for analcime that were estimated in spite of the use of nonstoichiometric dissolution data matched those derived from stoichiometric dissolution; within the uncertainty of the experimental measurements (i.e., of ion concentrations) on which such estimations are based.68

To better understand the origin of this unexpected outcome, the solution-phase data of Wilkin and Barnes were examined to assess how the different aqueous species impact log K, namely [Si], [Na], and [Al].68 To do so, three revised log K series were calculated at different temperatures for clinoptilolite and analcime under the assumption: [Al] = [Na], [Si] = x[Al], and [Si] = x[Na], where x is 2 for analcime (NaAlSi2O6·H2O) and 4.5 for clinoptilolite(Na) (Na1.1Al1.1Si4.9O12·3.5H2O), respectively, based on their [Si/Al] ratio. Using the solution data of Wilkin and Barnes,68 it was possible to recalculate log Ki for these compounds, assuming congruent dissolution of the zeolite of interest based on the original concentration of aqueous species, i. For instance, considering analcime, the log KNa series considers the dissolution of analcime to be congruent with respect to sodium, based on the assumption [Al] = [Na], [Si] = 2[Al], and [Si] = 2[Na] in accordance with the elemental ratios within solid analcime (NaAlSi2O6·H2O). Thus, Eq. (3) (see in Methods section) was used to calculate log Ki. The revised log Ki series were used in Eq. (1) as the KSi, KNa, and KAl terms. A fourth log K series was created to identify how the stoichiometric proportions of Si, Na, and Al in analcime and clinoptilolite(Na) impacted the solubility constant; i.e., as a means to moderate the Ki terms by the content of element “i” in the dissolving zeolite. Thus, based on the proportions of Si/Al and Si/Na in analcime and clinoptilolite(Na), the weighted log Kw was calculated as follows:

$$\log {K_{\mathrm{w}}} = \log \left( {{K_{\mathrm{Si}}}^{\frac{{\mathrm{Si}}}{{({\mathrm{Si}} + {\mathrm{Na}} + {\mathrm{Al}})}}}{K_{\mathrm{Na}}}^{\frac{{\mathrm{Na}}}{{({\mathrm{Si}} + {\mathrm{Na}} + {\mathrm{Al}})}}}{K_{\mathrm{Al}}}^{\frac{{\mathrm{Al}}}{{({\mathrm{Si}} + {\mathrm{Na}} + {\mathrm{Al}})}}}}\right)$$
(1)

where Si, Na, and Al reflect the elemental contents within analcime or clinoptilolite(Na). It should be noted that the log Ki values are expectedly calculated using solution-phase data of ion compositions as shown in Eq. (3) (see in Methods section).

The comparison of the three log Ki, and the weighted log Kw estimates is summarized in Fig. 2. Interestingly, the log Kw series corresponded well with the log KT curve (i.e., which is estimated from the measured solubility, and heat capacity data for a given zeolite; see Eq. 6), both for analcime and clinoptilolite(Na). This is because, via arithmetic compensation, the higher log KNa series at T < 90 °C was offset by the lower log KSi and log KAl for both analcime and clinoptilolite(Na). For analcime, at T > 100 °C, congruency in solutions was achieved leading to a convergence of log Ki values. In turn, estimation of the aqueous concentrations was, therefore, more accurate at T > 100 °C, as shown in Fig. 3a. A similar effect occurred for the clinoptilolite(Na) datasets, where the log Ki values became somewhat more consistent at T > 100 °C due to an approach towards congruent dissolution. However, since the log KSi values were higher than the log KAl and log KNa; it resulted in an under-prediction of [Si] abundance in solution, i.e., in terms of the amount of silicon that would be required to form clinoptilolite(Na) over the entire temperature range as shown in Fig. 3b.

Fig. 2: Recalculated solubility constants.
figure 2

Recalculated solubility constants of a analcime and b clinoptilolite(Na) at different temperatures using normalized congruent dissolution data.68 The solid line represents the extrapolated log KT of the respective phases based on Eq. (6). The circles represent the weighted averages of the recalculated log Kw series using Eq. (1).

Full size image
Fig. 3: Concentrations of aqueous species Al, Si, and Na as compared with those predicted with thermodynamic modeling.
figure 3

a analcime and b clinoptilolite(Na) dissolution for T < 300 °C. The dashed line represents the predicted concentrations of Si and the solid orange line represents the predicted concentrations of Al and Na.

Full size image

The analysis above offers a means to identify a propensity boundary for these zeolites to form. This approach however does not consider, e.g., nuances such as potential incongruency in dissolution, as a result of which, the (inferred) composition of the zeolite present, would be incorrect at lower temperatures. The simulations predict [Si] compositions, and relevant boundaries for analcime and clinoptilolite(Na) across the temperature range shown in Fig. 3; one of the key factors which determines the propensity for the formation of zeolites. While (alumino)silicate species are a critical framework component of zeolites, the extent of incongruency, which varies as a function of zeolite composition, temperature, etc., however, precludes accurate estimation of the zeolite composition from knowledge of aqueous silicon concentrations, alone. It should be noted that the lower levels of incongruency observed in the case of analcime—that features a lower [Si/Al] than that of clinoptilolite(Na)70,71—results in more accurate predictions of Na, Al, and Si (Fig. 3a) concentrations for this phase as compared with clinoptilolite(Na) (Fig. 3b).

Benning et al. observed that small variations in zeolite composition may have only a small impact on the Gibbs free energy of formation.67 For instance, data from Benning et al. suggests that a 2.5% variation in the [Si/Al] molar ratio results in a relative variation of 2% in the Gibbs free energy of formation,67 e.g., indicating a scaling that appears similar to Vegard’s Law.72 Moreover, Benning et al. noted that, at fixed Si/Al molar ratio, exchangeable cations have a limited effect on the Gibbs free energy of formation as they induced a relative variation of only around 0.5%.67 This is important as naturally formed zeolites can not only contain trace elements but their chemical composition (e.g., Si/Al molar ratio) can exhibit some variability. While the uncertainties noted above are negligible for calorimetric data; they can substantially affect the uncertainty that propagates into estimates of the solubility constant. Therefore, obtaining the thermodynamic properties of zeolite exhibiting varying compositions, e.g., as may form upon exposure to groundwater, is important for accurately ascertaining phase equilibria.

Chipera and Apps10 observed a relative difference between predicted and measured thermodynamic properties, namely, of 0.5% for the Gibbs free energy and enthalpy of formation (using data from Chermak and Rimstidt38) and 15% for the entropy (using data from Holland40). The method of Arthur et al.36 (i.e., based on a summation of hydroxides) typically offers Gibbs free energy and enthalpy of formation values that are within 0.2% and 0.3% of experimental data, respectively. On the other hand, the method of Vieillard,41 which is based on a summation of fictive oxides, offers entropy and heat capacity predictions that are within 3% of experimental data.36,41 Importantly, however, differences in molecular water content (i.e., water of crystallization) of zeolites can also affect their thermodynamic properties. For example, for analcime, we note that a 10% change in the number of zeolitic water molecules yields a variation of 0.75% in the Gibbs free energy of formation of the solid, which, in turn, yields a variation of 4 log K units in the predicted solubility constant. Moreover, differences in the thermodynamic properties of zeolitic water are also the origin of the systematic discrepancy that is observed between the predictive methods based on the summation of oxides or hydroxides. For example, differences in the reported entropy for water (namely, 52.0,41 59.0,39 and 69.95 J·K−1·mol−1 73) yield different estimates of the entropy of analcime (namely, 233.7, 240.7, and 251.65 J·K−1·mol−1, respectively); resulting in differences of up to 8% vis-à-vis experimental data.74

Data evaluation and completion of datasets

We review the zeolite phases included in the database and detail the method(s) that were used to estimate their thermodynamic properties. In general, thermodynamic data for the zeolite phases considered herein were sourced from published experimental data and, when such data was not available, were calculated via the predictive methods developed by: (a) Arthur et al.36 (for the enthalpy of formation), (b) Vieillard41 (for the heat capacity and entropy), and/or (c) volume-based correlations based on knowledge of the molar volume of the compound.75 Since predicted solubility data are very sensitive to small variations in thermodynamic inputs, the solubility constants at standard conditions (p = 1 bar, T = 25 °C) were calculated and assessed using reaction data from both the Gibbs free energy of formation (obtained from calorimetric or solubility measurements) and the enthalpy of formation and entropy (obtained by calorimetric measurements). In general, data offering good agreement between the solubility constants estimated via both approaches were prioritized for inclusion in the database.

New thermodynamic database for zeolite phases

The thermodynamic database developed herein, namely zeo19, comprises calcium- and sodium-dominant compositions and a minority of compositions containing other elements (see Fig. 4). Supplementary Table 2 presents the standard molar thermodynamic properties and the temperature-dependent coefficients used for the calculation of the solubility constants for both natural and synthetic zeolite phases considered herein (see Eq. (3) in Methods section) (for example, faujasite and gismondine are the natural variants, that contain impurity species, of the synthetic zeolites; zeolite X and Y, and zeolite P, respectively). An electronic copy of Supplementary Table 2 can be accessed via Supplementary Data 1. Supplementary Data 2 and 3 can be used to input the zeo19 database into the GEMS29 software directly.

Fig. 4
figure 4

Ternary phase diagram showing the chemical composition (in mole fraction) of the zeolite phases included in the thermodynamic database developed herein.

Full size image

Interestingly, previous studies have suggested the existence of empirical correlations among the various thermodynamic properties of solids (e.g., see the volume-based thermodynamics approach, and others76). In the following section, we examine such volume-based correlations within the thermodynamic database developed and such volume-based relationships are used to populate missing thermodynamic data (see below). Figure 5a shows the correlation between the molar volume and the molecular weight of the zeolitic phases considered in this database. We observe a positive correlation between molar volume and molecular weight, which arises from the fact that heavier atomic species often present larger radii. We obtain a coefficient of determination R2 = 0.94—which increases to R2 = 0.99 when zeolites of structurally-similar families are grouped together. These results are similar to those of Karzhavin,76 who also observed such a linear relationship for zeolites belonging to the natrolite group (with R2 = 0.9969). Figure 5b shows the correlation between the molar entropy and the molar volume of the zeolite phases. Again, following volume-based thermodynamics, we observe a strong linear correlation between these properties, in line with previous observations on fibrous zeolites.76 As shown in Fig. 5c, d, we also observe the existence of well-correlated empirical relations between the molar entropy, molar heat capacity, and molar enthalpy of formation of the zeolite phases.

Fig. 5: The relationships among the thermodynamic properties of the zeolite phases (all phases from Supplementary Table 2) considered in the new database.
figure 5

The new database includes a molar volume as a function of the molecular weight, b entropy as a function of molar volume, c entropy as a function of heat capacity, and d heat capacity as a function of enthalpy of formation. The solid lines show the best-fits. The associated equation and coefficient of determination (R2) are indicated in each panel. The red cross and purple star symbols indicate data for zeolite X and zeolite P(Ca), respectively, from Lothenbach et al.53

Full size image

The empirical relationships identified in Fig. 5 allow us to estimate the thermodynamic data of analcime, zeolite X, and zeolite P(Ca). Although the solubility of analcime for [Si/Al] = 2.5 was measured by Wilkin and Barnes,68 no calorimetry measurements for this composition are available. Neuhoff et al.77 derived a thermodynamic dataset using a solid-solution model for analcime with [Si/Al] = 1.9 and 3.0. We used this data for predicting the solubility of analcime with [Si/Al] = 2.5. In addition, we calculated the thermodynamic data for analcime with [Si/Al] = 2.5 using the methods developed by Arthur et al.36 and Vieillard.41 However, we find that these three datasets all tend to underestimate the solubility of analcime (viz-a-viz, experimental data)68 with [Si/Al] = 2.5 by ∼2 log units at temperatures up to 200 °C. To address this discrepancy, the relationships shown in Fig. 5 were used to estimate the thermodynamic properties for analcime with [Si/Al] = 2.5 (viz. Cp° = 158.12 J·K−1·mol−1, S° = 147.97 J·K−1·mol−1, Vm° = 88.43 cm3·mol−1). Notably, the predicted solubility values are greatly improved and agree with experimental data for temperatures up to 200 °C (see Fig. 5b). We note that the new entropy value that is predicted herein (viz. 147.97 J·K−1·mol−1) is significantly lower than those derived by Neuhoff et al. (viz. 226.98 J·K−1·mol−1 for [Si/Al] = 1.9 and 242.02 J·K−1·mol−1 for [Si/Al] = 3.0).77 La Iglesia and Aznar also estimated thermodynamic properties based on the summation of the constituent oxides for each zeolite—i.e., without considering any structural effect.78 The Gibbs free energy of formation (ΔfG° = −3043.5 kJ·mol−1) for analcime of [Si/Al] = 2.5 calculated using this method, agrees well with our estimate (viz. ΔfG° = −3042.65 kJ·mol−1) and the experimental measurement of Wilkin and Barnes.68

Careful inspection of Fig. 5b, d reveals that the molar entropy and heat capacity of zeolites X and P(Ca) sourced from Lothenbach et al. exhibit a significant deviation from the overall trend.53 While the reasons for this mismatch are not known, it may be on account of minor impurities or the formation of new phases through ion-exchange reactions or differences in the sources used for estimating the thermodynamic properties of the reference oxides. For example, zeolite P(Ca) was not observed to form at 20 °C due to slow reaction kinetics.53 However, a small quantity was detected at higher temperatures, i.e., 50 and 80 °C.53 The formation of zeolite P(Na) through ion-exchange reaction with the P(Ca) phase is unlikely, as the dissolved calcium concentration was found to be higher than that of sodium.53 Although the selectivity of zeolite P(Ca) for calcium is higher than for sodium, this process is reversible,79 which suggests that zeolite P(Na) and P(Ca) may coexist at 80 °C.80 Since calorimetric measurements are not available for these zeolites, Lothenbach et al.53 estimated the heat capacity of zeolite P(Ca) based on that of phillipsite(Ca),81 the oxide constituents, and water (from Helgeson).39 Similarly, the entropy and heat capacity of zeolites X were estimated based on data obtained for chabazite.53 As such, the large variability in reported values for the heat capacity of (zeolitic) water (viz. 47.7 J·K−1·mol−1 in ref.,39 75.4 J·K−1·mol−1 in the GEMS database,82 and 43.32 J·K−1·mol−1 based on the method developed by Vieillard41) may explain the discrepancy seen in Fig. 5b, d. Thus, using the methods of Arthur et al.36 and Vieillard,41 we calculate new values of the enthalpy of formation, entropy, and heat capacity for zeolite P(Ca) and zeolite X that coincide with the overall trends in Fig. 5b, d.

Predictions of zeolite solubility constants

We use the database developed herein and the GEMS approach to predict and validate the solubility of select zeolite phases as a function of temperature. To assess the effect of the new thermodynamic data calculated herein for zeolite P(Ca) and zeolite X, we first focus on these two systems and compare the solubility constants predicted by our new datasets with those estimated by Lothenbach et al.53 First, we observe that our estimations for zeolite P(Ca)’s solubility constant show reasonable agreement with experimental data from Lothenbach et al. at sub-boiling temperatures.53 In contrast with the predictions of Lothenbach et al.,53 we note that using our thermodynamic properties yields solubility constants that increase with increasing temperature, which agrees with the trend predicted by Blanc et al.52 Such an increasing trend is also observed for gismondine,83 which is the natural form of zeolite P(Ca). In addition, we observe that our data also yields a slightly improved prediction of the solubility constant of zeolite X.

Next, we focus on analcime. As shown in Fig. 6a, b, we observe that the estimated solubility for select analcime phases with [Si/Al] = 2.0 and 2.5 shows a good agreement with experimental data over the entire range of temperature.68 However, we observe a small mismatch between the predicted and experimental solubility values for T > 100 °C in the case of analcime with [Si/Al] = 2.5. Thereafter, for clinoptilolite, as shown in Fig. 6d–f, we once again find good agreement between predicted and experimental solubility data.67,68 We note that, in the case of clinoptilolite(K), the predicted solubility constant deviates slightly at low temperature from the experimental values reported in ref. 68 However, these experimental data were obtained by considering an equilibrium period of 213 days.68 In contrast, our predicted solubility data agree better with those of Benning et al.,67 which were obtained following a longer equilibrium time of 350 days. These results illustrate the critical need to sustain long equilibration times (i.e., especially at low temperature) to establish zeolite stability. Finally, as shown in Fig. 6c, we find that the solubility predicted for mordenite(Ca) agrees well with the experimental data of Benning et al.67 Moreover, and in general, the predictions of the solubilities of analcime, clinoptilolite, and mordenite derived herein agree within 1 log unit of the experimental solubility constants. This difference is consistent with the previously observed impact on solubility predictions resulting from an uncertainty of 0.3% in the Gibbs free energy of formation (or enthalpy of formation) or to an uncertainty of 10% in the entropy (or heat capacity). Overall, these results strongly support the ability of our approach to yield realistic predictions of the solubility constants of the zeolite phases considered herein.

Fig. 6: The temperature-dependence of the solubility constant (log KT) predicted using the zeo2019 database.
figure 6

Developed herein for: a analcime with Si/Al = 2.0, b analcime with Si/Al = 2.5, c calcium mordenite, d sodium clinoptilolite, e potassium clinoptilolite, and f calcium clinoptilolite. The predicted solubility data are compared with experimental data sourced from refs. 67,68

Full size image

Demonstrating the prediction capability of the database: estimating stability fields of zeolites

As a validation of our new database, we construct a series of equilibrium activity diagrams and compare the predicted domains of stability with experimental observations of zeolite precipitation in the context of glasses used for the vitrification of nuclear wastes. We simulate the formation of zeolitic compounds from a series of solutions at 25 °C (i.e., as relevant to ambient temperature glass dissolution) and 90 °C (i.e., as relevant to accelerated glass dissolution). To mimic conditions of the dissolution of Si-rich nuclear waste immobilization glasses,5 we focus on Si-saturated solutions (i.e., with the activities log a(SiO2) = −2.70 at 25 °C and log a(SiO2) = −2.23 at 90 °C82) and establish the associated stability fields for varying aluminum, calcium, potassium, and sodium activities. Figure 7 shows the predicted stability fields for solutions that present a silica activity that corresponds to a solution in equilibrium with amorphous silica. We observe that the boundaries of the stability domains exhibit a slope of −2 and −1 in log–log scale for alkaline earth- and alkali-ion activities, respectively.

Fig. 7: Predicted equilibrium activity diagrams for different Si-saturated solutions at 25 °C and 90 °C.
figure 7

Namely, a, d CaO–Al2O3–SiO2–H2O system, b, e K2O–Al2O3–SiO2–H2O system, and c, f Na2O–Al2O3–SiO2–H2O. The symbols indicates the solution (pH) conditions wherein the precipitation of select zeolite phases was observed, namely, (i) chabazite (cha, triangular symbols, log a(AlO2−) = −3.22, log a(Ca2+) = −4.59), (ii) merlinoite(K) (mer, star symbols, log a(AlO2−) = −4.26, log a(K+) = −0.95; log a(AlO2−) = −3.38, log a(K+) = −0.85), (iii) analcime (ana, red circular symbols, log a(AlO2−) = −4.26, log a(Na+) = −3.27), and (iv) zeolite X (x, green circular symbols log a(AlO2−) = −3.24, log a(Na+) = −2.19).86,87,89,90

Full size image

Nuclear waste containers may be disposed at depths ranging between 500 and 1000 m with storage temperatures of <100 °C.84 Under these disposal conditions, the hydrostatic pressure of water may reach up to 100 bar85 if the repository is fully saturated with groundwater. At these ranges of pressure, using Eqs. (6) and (7) to correct for temperature and pressure effects, respectively, at 100 °C the solubility constant of analcime at vapor pressure (1.01 bar) and 100 bar are −13.46 and −13.40, respectively, whereby the difference is far less than the experimental error of 0.5 log K units calculated by Wilkin and Barnes68 under such conditions. As the effect of pressure is far less than experimental error, the effect of pressure on solubility constants is considered insignificant under these circumstances. On the other hand, for example, the solubility constant of analcime increases from −16.10 to −13.46 from 25 to 100 °C, respectively, noting that temperature rather than pressure is the variable that affects phase equilibria.

First, in the CaO–Al2O3–SiO2–H2O system, at low aqueous activities of aluminum and calcium, several zeolite phases (e.g., heulandite, phillipsite, leonhardite, yugawaralite, stellerite, laumontite, clinoptilolite, mordenite, and scolecite) with small stability field boundaries are indicated to form and persist (see Fig. 7a). On the other hand, gismondine, wairakite, chabazite, and zeolite P(Ca) are found to be stable at higher aqueous cationic activities. However, we find that some of these phases do not persist at low temperature, namely, below 95 °C. For instance, wairakite, yugawaralite, leonhardite, and stellerite are encountered in high-temperature active geothermal areas,15,48 whereas laumontite and heulandite are found in the deepest zones of burial diagenesis.16,44 Moreover, clinoptilolite has been observed to transform partially to heulandite and, subsequently, heulandite to transform to laumontite with increasing depth since, in general, temperature increases with depth (25 °C for every 1 km of depth16). Note that, with increasing temperature, the stability field of zeolite P(Ca) predicted via the dataset of Lothenbach et al.53 actually falls into the wairakite stability field, whereas, with the new dataset, the stability field of zeolite P(Ca) is merged with that of chabazite.

Second, in the K2O–Al2O3–SiO2–H2O system, we find that merlinoite is the stable phase with respect to clinoptilolite and phillipsite in solutions having a high aqueous activity of potassium (see Fig. 7b). Finally, in the Na2O–Al2O3–SiO2–H2O system, phillipsite, heulandite, clinoptilolite, and natrolite phases are found to be stable at low aqueous activities of sodium compared with analcime and zeolites X and Y, which are found to be stable at higher aqueous sodium activity (see Fig. 7c). However, heulandite and natrolite have been observed in active geothermal systems,15,48 and so their formation is favored at high temperatures. On the other hand, at high aqueous activities of sodium, zeolites X and Y are found to be stable. Furthermore, with increasing temperature, the stability field of analcime is found to broaden and merge with that of zeolite Y. Thus, overall, we observe that the domains of stability of many zeolites are small and closely related to the activity of their constituent ions in solution.

We now compare the predicted stability fields with experimental observations of zeolite precipitation during the dissolution of select nuclear waste immobilization glasses.5,86 First, we model the solution chemistry of relevance upon the progressive dissolution of a complex glass (CG) at 90 °C after 564 days of corrosion.86 Herein, we find that, when the activity of aluminum increases to −3.22, the solution becomes super-saturated with respect to chabazite(Ca) (see Fig. 7d). In agreement with our modeling, the formation of chabazite(Ca) was indeed reported (e.g., see the triangular symbols in Fig. 7d) and was found to be associated with a resumption of glass corrosion (i.e., Stage III corrosion) after 508 days of dissolution.86 Second, we model the solution chemistry of relevance upon the progressive dissolution of an Advanced Fuel Cycle Initiative (AFCI) nuclear waste immobilization glass at 90 °C.87 We find that, after 600 days of corrosion, the solution becomes super-saturated with respect to chabazite(Ca) (see Fig. 7d). Under such conditions, the formation of chabazite(Na) was observed.87 Note that chabazite can accommodate both Na or Ca, but Ca is the most frequent extra-framework cation88—however, no thermodynamic dataset is presently available for chabazite(Na). Since chabazite can exchange Na with Ca cations19,23, our modeling results suggest that chabazite(Na) can be formed via an exchange reaction with Ca or via direct uptake of dissolved calcium cations present in the solution.

We now focus on the case of the progressive dissolution of the International Simple Glass (ISG) at 90 °C in K-rich solution.89 We find that the resulting aqueous activities of potassium and aluminum (log a(Na+) = −0.95 and log a(AlO2−) = −4.26)—in solution—lie in the stability field of merlinoite(K), which agrees with experimental observations (see the star symbols in Fig. 7e).89 Finally, we consider the case of the solutions forming upon the dissolution at 90 °C of an ISG glass following the data reported by Gin et al.89 We find that, after 1 day and 35 days after the adjustment of pH, the aqueous activities of sodium and aluminum fall in the stability fields of zeolite X and analcime (see Fig. 7f). In agreement with our predictions, the precipitation of zeolites was indeed observed by Gin et al.89 (see the circular symbols in Fig. 7f). In agreement with previous observations from Lothenbach et al.,53 we find that the formation of zeolite X is favored at high initial aqueous sodium and aluminum activities (log a(Na+) = −2.19 and log a(AlO2−) = −3.24). However, the precipitation of zeolite X tends to consume sodium and aluminum, which, in turn, reduces their activity (log a(Na+) = −3.27 and log a(AlO2−) = −4.26) and eventually favors the formation of analcime. Similarly, by modeling the dissolution of a WVUTh198 type IR glass at 90 °C at various test times,87 we find that the aqueous activities of sodium and aluminum fall in the stability field of zeolite X and progressively move toward that of analcime. These results suggest that analcime may become thermodynamically more favorable to form, on account of the consumption of solution species due to the precipitation of zeolite X. Overall, the striking harmony between the stability fields predicted in Fig. 7 and experimental observations, wherein zeolite precipitation occurs, offers a robust validation for the new database developed herein and strongly supports the ability of our approach to provide a realistic description of the stability of various zeolites, and the propensity for their formation across a wide range of solution chemistries, pH, and reaction temperatures.

In closing, we have developed, and validated a unified, comprehensive thermodynamic database for zeolite compounds for conditions representative of nuclear waste disposal, i.e., p = 1 bar and T < 95 °C. In addition to thermodynamic data sourced from the literature (viz. bikitaite, brewsterite, dachiardite, epistilbite, erionite, ferrierite, gonnardite, harmotome, leonhardite, paranatrolite, tetranatrolite, thomsonite, and yugawaralite), new thermodynamic inputs are derived for analcime, clinoptilolite, mordenite, zeolite X, and zeolite P(Ca) compounds, using volume-based or empirical correlations. When combined with the GEMS modeling approach, this new database yields realistic predictions of the propensity of the formation of zeolite phases as a function of temperature—which are in close agreement with experimental observations of zeolite formation.86,87,89,90 This analysis involved the construction of equilibrium activity diagrams and comparing them with experimental observations of zeolite formation (e.g., for analcime, merlinoite, and chabazite) as observed in glass dissolution studies. The good agreement between predictions and experiments offers a robust validation of our new database. Overall, the database provides a realistic description of the stability of various zeolite phases across a wide range of solution chemistry, pH, and temperature, which paves the way toward a robust prediction of propensity of zeolite formation (and the potential for the resumption of Stage III glass dissolution) in nuclear waste repository conditions.

Methods

Parameterization methods

GEM-Selektor has been extensively used to model geochemical systems including solid-solutions and aqueous complexes by using the principles of local and partial equilibrium.28,29 Herein, the thermochemical data and dissolution reaction for any new phase are developed using the reaction-defined components module (“ReacDC”) within the GEMS database mode.29 The reaction-induced change in thermodynamic properties for the dissolution reaction of zeolites is calculated using aqueous species and complexes from the PSI-Nagra 12/07 TDB91 and Cemdata18 databases.82 These databases are mutually consistent, namely, both contain the same primary aqueous, gaseous, and solid species.82,91 The PSI-Nagra 12/07 TDB comprises the same entries as the PSI-Nagra 01/01 TDB92 and includes improved thermodynamic properties for silicate and aluminosilicate species, e.g., Si4O8(OH)44− and AlSiO3(OH)43−.91 The Cemdata18 database encompasses hydrated phases commonly encountered in cementitious systems.82 The internal consistency of the thermodynamic database developed herein is maintained by using the same sources, and constituents for the aqueous species and complexes. Thus, the reaction data for a given phase is calculated with the primary aqueous species and complexes of these databases together with the data selected for the solid. Given a dissolution reaction of a zeolite with composition CcDdEe, species C, D, and E are generated in an aqueous solution as presented in (2):

$${\mathrm{C}}_{\mathrm{c}}{\mathrm{D}}_{\mathrm{d}}{\mathrm{E}}_{{\mathrm{e(s)}}} = {\mathrm{cC}}_{({\mathrm{aq}})} + {\mathrm{dD}}_{({\mathrm{aq}})} + {\mathrm{eE}}_{({\mathrm{aq}})}$$
(2)

The solubility constant, K, for this can be expressed as the ratio between the product of the activities {a} of the aqueous species raised to their stoichiometric coefficients and the activity of the reactant as written below:

$$K = \frac{{\left\{ {a_{{\mathrm{C}}({\mathrm{aq}})}} \right\}^{\mathrm{c}}\,\left\{ {a_{{\mathrm{D}}({\mathrm{aq}})}} \right\}^{\mathrm{d}}\,\left\{ {a_{{\mathrm{E}}({\mathrm{aq}})}} \right\}^{\mathrm{e}}}}{{\left\{ {a_{{\mathrm{C}}_{\mathrm{c}}{\mathrm{D}}_{\mathrm{d}}{\mathrm{E}}_{\mathrm{e}}({\mathrm{s}})}} \right\}}}$$
(3)

where the activity of the ith aqueous species is calculated as the product of molar concentration [X] and the activity coefficient (γ) as shown in \(\left\{ {a_i} \right\} = {\upgamma}_i\left[ {X_i} \right]\). The activity is assumed to be unity for pure minerals, therefore, \(\{ a_{{\mathrm{C}}_{\mathrm{c}}{\mathrm{D}}_{\mathrm{d}}{\mathrm{E}}_{\mathrm{e}}\left( {\mathrm{s}} \right)}\}\) is equal to 1. Standard state conditions are taken as 298.15 K and 1 bar. To consider the effects of solution non-ideality and water activity, the extended Debye–Hückel equation is used to calculate the activity coefficients of aqueous species as shown below:93

$${\mathrm{log}}\,\gamma _i = \frac{{ - A_\gamma z_i^2\sqrt I }}{{1 + \dot aB_\gamma \sqrt I }} + b_\gamma I + {\mathrm{log}}\frac{{X_{jw}}}{{X_w}}$$
(4)

where γi and zi are the activity coefficient and charge of the ith aqueous species, respectively, Aγ and Bγ are temperature- and pressure-dependent coefficients, I is the molar ionic strength, \(X_{jw}\) is the molar quantity of water, and Xw is the total molar amount of the aqueous phase. A common ion size parameter, \(\dot a\) (3.31 Å) and short-range interaction parameter, bγ (0.098 kg·mol−1), are used as constants for the background electrolyte. NaOH is selected herein to replicate alkaline pH required to cause the dissolution of glasses and precipitation of zeolites under typical geological repository conditions.

The solubility constant of the dissolution reaction at standard state conditions is related to the Gibbs free energy of reaction (ΔrG) as presented below:

$$\Delta _{\mathrm{r}}G^ \circ = - {\mathrm{RT}}\,\ln \left( K \right)$$
(5)

where R = 8.31451 J·mol−1·K−1 is the ideal gas constant and T is the temperature in Kelvin. The reaction data at standard conditions can be calculated with the expression \(\Delta _{\mathrm{r}}G^ \circ = \Delta _{\mathrm{r}}H^ \circ - {\mathrm{T}}\Delta _{\mathrm{r}}S^ \circ\), where \(\Delta _{\mathrm{r}}H\) and \(\Delta _{\mathrm{r}}S\) represent the reaction-induced changes in enthalpy and entropy, respectively. The reaction changes of these properties can be calculated using the temperature-dependence of the heat capacity of the solid and aqueous complexes. At low temperatures, lattice vibrations contribute mainly to the heat capacity of aluminosilicates.94 Furthermore, the thermal effect of hydrous silicates increases slightly with temperature as water species (hydroxyl and water molecules) impact the entropy but cause small changes in the heat capacity.95 Therefore, the heat capacity of a given phase is assumed to be constant, since, at room temperature, the influence of the temperature-dependence of the heat capacity on solubility constants is generally small.52

The Helgeson–Kirkham–Flowers (HKF) equation of state parameters, included in the PSI-Nagra 12/07 database,91 describe the effects of temperature and pressure.75,96 Predictions of the solubility constant along the liquid-vapor curve of water are influenced by the water vapor pressure, e.g., \(p({\mathrm{H}}_2{\mathrm{O}})\) = 30 bar at 300 °C.97 This is due to changes in the hydration state with temperature along the boiling curve of water. However, the solubility constants estimated from solubility data are assumed to be unaffected by hydration state since the system remains fully saturated.68,97 The total pressure of the system may be higher than that from the water vapor. However, extrapolation from standard state properties to pressure up to 5 kbar reproduces experimental data for aqueous species and complexes of aluminosilicate systems.96,98 Importantly, first, the phase equilibrium boundary of the NaAlSi3O8–NaAlSiO4–H2O system at 2 and 5 kbar from 150 to 900 °C was studied and the stability field of analcime was determined to be temperature-dependent but insignificantly affected by pressure.99 Second, the composition of an aluminum-rich analcime was noted to remain constant at T < 300 °C and p < 5.15 kbar.99 These conditions are rather dissimilar to those (100–300 °C) typically encountered in closed batch systems used for solubility experiments.67,68 Therefore, and in general, the effect of pressure at a constant temperature is described by volume and pressure induced changes of the reaction products; which being minimal do not affect the solubility constants at p < 5.15 kbar.100

The solubility constant at 298.15 K is fixed and its heat capacity coefficients are described by a 3-term approximation which is valid for temperatures below 100 °C (6):82

$${\mathrm{log}}\,K_{\mathrm{T}} = A_0 + A_2T^{ - 1} + A_3\,{\mathrm{In}}T$$
(6)

where log KT is the logarithm of the solubility constant for a given reaction at a given temperature, A0, A2, and A3 are temperature-dependent coefficients expressed as: \(A_0 = \left( {0.4343\,R^{ - 1}} \right)\left[ {{\mathrm{\Delta }}_{\mathrm{r}}S_{T_ \circ }^ \circ - {\mathrm{\Delta }}_{\mathrm{r}}Cp_{T_ \circ }^ \circ \left( {1 + {\mathrm{In}}{\it{T}}_ \circ } \right)} \right]\), \(A_2 = \left( {0.4343\,R^{ - 1}} \right)\left[ {{\mathrm{\Delta }}_{\mathrm{r}}H_{T_ \circ }^ \circ - {\mathrm{\Delta }}_{\mathrm{r}}Cp_{T_ \circ }^ \circ T_ \circ } \right]\), and \(A_3 = 0.4343\,R^{ - 1}\left[ {{\mathrm{\Delta }}_rCp_{T_ \circ }^ \circ } \right]\) S° is the entropy, and Cp° is the heat capacity at To = 298.15 K at constant pressure. The reaction-induced change in entropy and heat capacity at a temperature T are written as \({\mathrm{\Delta }}_{\mathrm{r}}S_T^ \circ = {\mathrm{\Delta }}_rS_{T_ \circ }^ \circ + {\mathrm{\Delta }}_{\mathrm{r}}Cp_{T_ \circ }^ \circ \,{\mathrm{In}}\left( {T/T_ \circ } \right)\) and \({\mathrm{\Delta }}_{\mathrm{r}}H_T^ \circ = {\mathrm{\Delta }}_{\mathrm{r}}H_{T_ \circ }^ \circ + {\mathrm{\Delta }}_{\mathrm{r}}Cp_{T_ \circ }^ \circ \,\left( {T - T_ \circ } \right)\).

At a fixed temperature, the following equation can be used to calculate the solubility constant as a function of temperature and pressure (KT,P)29:

$${\mathrm{log}}\,K_{{\mathrm{T,P}}} = {\mathrm{log}}\,K_{\mathrm{T}} - \frac{{0.4343}}{{{\mathrm{RT}}}}{\mathrm{\Delta }}_{Tr}V_{Tr}^ \circ \left( {P - P_ \circ } \right)$$
(7)

where \({\mathrm{\Delta }}_{Tr}V_{Tr}^ \circ\) is the volume change of reaction, KT is the solubility constant at temperature, T in kelvin, and P is pressure in bar.

Sensitivity analysis

The impact of uncertainties associated with the thermodynamic data of zeolites on estimations of the solubility constant at standard conditions is assessed via a sensitivity analysis. Herein, the uncertainty associated with the thermodynamic data of the aqueous species is assumed to be fixed. Thus, only the impact of uncertainty associated with the thermodynamic properties of the mineral, and its effect on the solubility constant are evaluated. We assume that the major source of uncertainty arises from calorimetric measurements of the solid’s properties. The sensitivity analysis is based on a local approach which involves varying one input parameter at a time, while all the other input parameter values remain fixed.101 The uncertainty range for each input parameter is estimated by averaging the experimental uncertainty of calorimetric measurements from Supplementary Table 1. For the entropy and heat capacity, average experimental uncertainty yields small impact, so the uncertainty ranges are increased to 10% and 20%, respectively. The magnitude of the impact of each input parameter is assessed by comparing the predicted solubility constant including its uncertainty with the corresponding value without considering any uncertainty.

Estimation of uncertainty for the thermodynamic properties of zeolites

The relative uncertainty of the thermodynamic properties (\(\Delta _fH^ \circ\), S°, Cp°, and \(V_m^ \circ\)) are based on experimental values from the literature. However, for some minerals, these uncertainties are not reported. In these cases, based on our previous sensitivity analysis, the uncertainties for these properties are assumed to be 0.3% for \(\Delta _fH^ \circ\), 1.5% for S°, 3% for Cp°, and 0.5% for \(V_m^ \circ\).

For a zeolite dataset that does not contain a Gibbs free energy value, this can be estimated using the relationship with entropy and enthalpy (\(\Delta _fG^ \circ = \Delta _fH^ \circ - T\Delta _fS^ \circ\)). The uncertainty (2μ) of the estimated Gibbs free energy value can be obtained from its standard deviation (μ), calculated using an algebraic sum formula based on a propagation of errors method:102

$$\mu _{\Delta _fG^ \circ } = \sqrt {\left( {\mu _{\Delta _fH^ \circ }} \right)^2 + \left( {\mu _{S^ \circ }} \right)^2}$$
(8)

Similarly, the uncertainty of the solubility constant is estimated from the Gibbs free energy of formation, whereby the Gibbs energy of reaction is subsequently calculated and used in Eq. (5). The uncertainty of the Gibbs energy of reaction and the molar gas constant (R = 8.3144598 J·mol−1·K−1, with a relative uncertainty of 2.7 × 10−7 103) are used to calculate the uncertainty of the solubility product.

Calculation of field boundary in equilibrium activity diagrams

Construction of equilibrium activity diagrams for four component systems is performed by balancing reactions at a given temperature.42 For instance, in the system Na2O–Al2O3–SiO2–H2O, a reaction between amorphous silica and analcime can be represented by:

$${\mathrm{SiO}}_{2({\mathrm{am}})} = {\mathrm{SiO}}_2$$
(9)
$${\mathrm{NaAlSi}}_2{\mathrm{O}}_6 \cdot {\mathrm{H}}_2{\mathrm{O}} = {\mathrm{Na}}^ + + {\mathrm{AlO}}_2^ - + 2\,{\mathrm{Si}}_2{\mathrm{O}} + \,{\mathrm{H}}_2{\mathrm{O}}$$
(10)

Eqs. (8) and (9) are combined by eliminating aqueous \({\mathrm{SiO}}_2\) and the resulting reaction is

$${\mathrm{NaAlSi}}_2{\mathrm{O}}_6 \cdot {\mathrm{H}}_2{\mathrm{O}} = {\mathrm{Na}}^ + + {\mathrm{AlO}}_2^ - + 2\,{\mathrm{SiO}}_{2({\mathrm{am}})} + {\mathrm{H}}_2{\mathrm{O}}$$
(11)

The solubility constant for the reaction (10) is expressed as

$$\log K = \log a_{{\mathrm{Na}}^ + } + \log a_{{\mathrm{AlO}}_2^ - } + \log a_{{\mathrm{H}}_2{\mathrm{O}}}$$
(12)

The boundary formed between amorphous silica and analcime in an activity diagram is determined by the log K value (which is the intercept in a diagram with x and y axes chosen as \(\log a_{{\mathrm{Na}}^ + }\) and \(\log a_{{\mathrm{AlO}}_2^ - }\)). The activity coefficient of the aqueous species at a given temperature is calculated with the GEMS code using the embedded PSI-Nagra 12/07 TDB database.91 The activity of water is assumed to be equal to 1. For the other Na-zeolites, similar calculations are performed to stablish the phase boundaries for the Na2O–Al2O3–SiO2–H2O system. An, a similar approach is once again used for determining the phase boundaries of the CaO–Al2O3–SiO2–H2O and K2O–Al2O3–SiO2–H2O systems.