All samples were prepared, stored and handled in an Ar-glove box (O₂ < 1 ppm) at T = (22 ± 2)°C. All solutions were prepared with purified water (Milli-Q academic, Millipore, 18.2 MΩ) purged with Ar for at least 1 h to remove CO₂(g). Beryllium sulfate tetrahydrate (BeSO₄∙H₂O, 99.99% purity) and calcium chloride (CaCl₂∙H₂O) EMSURE^® were purchased from Merck. Calcium hydroxide (Ca(OH)₂) was purchased from Sigma-Aldrich.

A combination pH electrode (ROSS Orion with 3.0 M KCl as filling solution) was used to determine the proton concentration ([H⁺] in mol∙kg^–1, with pH_m = –log [H⁺]). The pH electrode was calibrated using commercial pH buffers (pH = 1–12, Merck). In salt solutions of ionic strength I ≥ 0.1 mol kg^–1, the measured pH value (pH_exp) is an operational apparent value related to [H⁺] by pH_m = pH_exp + A_m. The values A_m entail both the activity coefficient of H⁺ and the liquid junction potential of the electrode at a given background electrolyte concentration and temperature. A_m values reported by Altmaier and co-workers for CaCl₂ systems were used for the determination of pH_m (Altmaier et al., 2008).

The solid phase used for undersaturation solubility experiments, α-Be(OH)₂(cr), was precipitated by slow addition of a 0.35 M BeSO₄ solution to a 2.0 M carbonate-free NaOH solution with a final pH ≈ 10.5. After an ageing time of ≈2 months, the solid phase was washed thoroughly with slightly alkaline water (pH ≈ 9.5) to remove the residues of the matrix solution (Na₂SO₄ and NaOH), and characterized by X-ray diffraction (XRD). Details on the preparation and characterization of the solid phase used in this work were provided in our previous study (Çevirim-Papaioannou et al., 2020).

A total of 30 independent batch samples were prepared by contacting 0.5–5 mg of α-Be(OH)₂(cr) to alkaline CaCl₂ solutions (0.05, 0.25, 1.0 and 3.5 M) with 8.9 ≤ pH_m ≤ 11.8. Be(II) concentrations and pH_m values were monitored regularly from 7 to 337 days. In each sampling, aliquots (100–500 µL) of the supernatant of each sample were pipetted into 10 kD filters (2–3 nm cut-off Nanosep^®, Pall Life Sciences), centrifuged for 10 min at 6000 g, and diluted (1:100 and 1:4000) with 2% ultrapure HNO₃ before the quantification by inductively coupled plasma mass spectrometry (ICP-MS, PerkinElmer ELAN 610). The detection limit of the technique ranged between 10^−6.5 and 10⁻⁸ M, depending upon the dilution factors applied at different CaCl₂ concentrations. Equilibrium conditions were assumed after repeated measurements over time with constant [Be] and pH_m values, typically after 372 days. After attaining equilibrium conditions, the solid phase of a selected sample (1.0 M CaCl₂, pH_m = 11.1) was characterized by XRD. For this purpose, an aliquot of the solid phase (≈1–2 mg) was washed 7 times with ethanol (0.5–1 mL) under Ar atmosphere in order to remove the background electrolyte, i.e., CaCl₂. After the last cleaning step, the washed solid phase was re-suspended in ethanol, deposited on the XRD sample holder and left drying for a few minutes inside the glovebox. The XRD diffractogram was collected using a Bruker D8 Advance X-Ray powder diffractometer (Cu anode) within the range 2θ = 2–80⁰, incremental steps of 0.02 and a measurement time of 0.4 s per step.

The simulation model was built in a box with dimensions of 43 Å × 43 Å × 49 Å, which contained 3,120 water molecules with 3 Be(OH)₄ ^2–, 50 Ca²⁺, and 84 Cl^–randomly distributed ions. This approximately corresponds to 1 M CaCl₂ and 0.05 M beryllium ion concentration. The high concentration of Be(II) is required in the simulations for better statistical sampling, and is in line with experimental observations in concentrated CaCl₂ solutions. The total electroneutrality of the simulation cell was ensured by balancing the number of counterions. Periodic boundary conditions were applied in three dimensions. Water was described using the extended simple point charge (SPC/E) model (Berendsen et al., 1987). The 12-6-4 Lennard-Jones-type nonbonded parameters for Ca²⁺ and Be²⁺ for SPC/E water, which include the contribution from the ion-induced dipole interaction, were taken from Li and Merz (Li and Merz, 2014; Li et al., 2020). Standard Lorenz-Berthelot mixing rules (Allen and Tildesley, 2017) were applied to calculate short-range Lennard-Jones interactions between the unlike atoms (with a cut-off distance of 1.4 nm). Long-range electrostatic forces were evaluated using the Ewald summation method. All molecular dynamics (MD) simulations were performed using the LAMMPS software package (3 March 2020 version) (Plimpton, 1995). Model equilibration was carefully monitored by assessing the temperature, pressure, kinetic and potential energy of the system, and dimensions of the simulation box, in order to confirm that these parameters reach their equilibrium steady state values on average (Braun et al., 2019). The Newtonian equations of the atomic motions were numerically integrated with a time step of 1 fs, and the model systems were initially equilibrated for 5 ns in the isobaric–isothermal statistical ensemble (NPT), then for 5 ns in the canonical ensemble (NVT). Temperature and pressure were constrained using the Nose-Hoover thermostat and barostat (Braun et al., 2019) at ambient conditions (T = 295 K, p = 0.1 MPa). The production run was performed in the NVT ensemble for 5 ns, and the trajectory was recorded every 500 fs. To describe local structural properties, radial distribution functions for pairs of atoms, and running coordination numbers were calculated from the recorded trajectory. The VMD software package (version 1.9.3) has been used for visualisation (Humphrey et al., 1996).

XRD patterns collected for the beryllium solid phase equilibrated in 1.0 M CaCl₂ at pH_m = 11.1 and 3.5 M CaCl₂ at pH_m = 10.6 are shown in Figure 1, together with the diffractograms previously reported for α-Be(OH)₂(cr) in Bear and Turnbull (1965) and in our previous study (Çevirim-Papaioannou et al., 2020). The main features of this sample are in excellent agreement with the patterns previously reported for α-Be(OH)₂(cr). This observation confirms that the starting material used in the experiments, α-Be(OH)₂(cr), is the solid phase controlling the solubility of beryllium also in CaCl₂ systems.

Figure 2 shows the experimental solubility data determined in alkaline CaCl₂ solutions with 8.9 ≤ pH_m ≤ 11.8, together with the solubility curve of α-Be(OH)₂(cr) at I = 0 calculated using the thermodynamic model reported in (Çevirim-Papaioannou et al., 2020) (see Supplementary Tables SI-S1, S2 in the Supplementary Material). The solubility of Be(II) in 0.05 M CaCl₂ shows a pH-independent behaviour within 9 ≤ pH_m ≤ 10.5, which is in good agreement with the calculated solubility line and can be assigned to the solubility reaction α-Be(OH)₂(cr) ⇔ Be(OH)₂(aq). Above pH ≈ 10.5, the solubility in 0.05 M CaCl₂ increases steadily with a slope of ≈ +2 (as log [Be(II] vs. pH_m), in disagreement with the slope of +1 observed for NaCl/KCl systems within the same pH-range and also predicted by the calculated solubility line. This observation reflects that in CaCl₂ solutions the solubility of the α-Be(OH)₂(cr) phase is dominated by an equilibrium reaction involving the release of 2 H⁺. A well-defined slope of ≈ +2 is also observed for all other CaCl₂ systems investigated in this work, with the highest solubility being achieved in 3.5 M CaCl₂ solutions. Our experimental observations in CaCl₂ systems and the direct comparison with solubility data in NaCl/KCl solutions support the earlier predominance of the [Be(OH)₄]^2– moiety in the former systems, expectedly as a result of its stabilization by Ca²⁺ ions. This is in line with previous observations reported in the literature for strongly hydrolyzing metal ions, e.g., M(III), M(IV) and M(V) (Altmaier et al., 2008; Neck et al., 2009; Fellhauer et al., 2010; Fellhauer, 2013; Yalcintas et al., 2016).

The local structure of a solution can be described by means of the radial distribution functions of the ion-ion and ion-water pairs. The RDF shows the probability of finding atoms at curtain distance during simulation time, and allows to calculate time-averaged coordination numbers. Figure 3 presents calculation results for Be²⁺ ions in CaCl₂ aqueous solution.

Figures 3A, B shows that the interaction of Be²⁺ with hydroxyl ions is strong and was stable during all the simulation time. There are 4 oxygen atoms of the aqueous hydroxyl ions in the first coordination sphere of Be²⁺ at r ≈ 1.7 Å. This result agrees well with the reported distances between Be²⁺ and hydroxyl ions calculated by DFT (Jin et al., 2015). The oxygen atoms of the first shell are coordinated with other solution molecules through hydrogen bonding (Hw peak at distances around 2–3 Å) to the second solvation shell (small Ow peak at r ≈ 3–4 Å). As it can be seen from Figures 3A, B, there is a low probability of binding between Be²⁺ and Cl⁻ions: no distinct peak of chloride ions was recorded, and Cl⁻presence was only found at r > 5 Å.

Three peaks for Ca²⁺ are found in the outer coordination sphere of Be²⁺: a high-intensity peak at r ≈ 3.1 Å, a small peak between 3 Å and 4 Å, and very broad peak of very low intensity around 5 Å. The first peak corresponds to binding of two calcium cations with oxygens of the hydroxyl ions in bidentate coordination, as can be seen in the snapshot from simulations (Figure 3C). The peak between 3 and 4 Å is the coordination of Ca²⁺ with a single hydroxyl group. The broad peak after 5 Å stands for Ca²⁺ ions in the solution that are not bound to the complex. We have previously shown that Ca²⁺ plays a crucial role in the sorption of Be(II) in cementitious systems (Çevirim-Papaioannou et al., 2021a; Çevirim-Papaioannou et al., 2021b). The surface of C-S-H provides a more structured environment, where both ions and water molecules are coordinated with the surface in the sorption layer increasing the probability for Be hydroxyl species to form complex with multiple Ca²⁺ ions (up to three). From the results of this study, it can be concluded that in the solution complexes of Be(OH)₄ ^2– with two Ca²⁺ are dominant, while binding of the third calcium cation remains possible, which is reflected in the slightly higher value of running coordination for the Be-Ca pair (n ≈ 2.2 at a distance of 4 Å).

The specific ion interaction theory (SIT) (Ciavatta, 1980) based on the Brønsted-Guggenheim-Scatchard model is used in this study to account for ion interactions in systems at I > 0. This is also the method adopted for ion strength corrections within the NEA-TDB project (Grenthe et al., 2020). In SIT, the activity coefficient of a given ion j (γ _j ) is calculated according to: l o g γ j = − z j 2 D + ∑ k ε j , k , I m m k (1) where z _j is the charge of the ion j, D is the Debye-Hückel term ( D = 0.509 I m 1 + 1.5 I m ), m _k is the molality of the oppositely charged ion k, and ε(j, k, I _m ) is the specific ion interaction parameter. The Debye Hückel term in Equation 1 accounts for the electrostatic, non-specific long-range interactions prevailing in dilute systems, whereas the term ∑ k ε j , k , I m m k accounts for short-range, non-electrostatic interactions that become relevant at higher ionic strength conditions. The SIT coefficient ε(j, k, I _m ) reflects also differences between ions of the same charge but different size. The applicability of the SIT approach is often considered limited to I _m ≤ 3.5 mol kg^–1, although previous studies have reported reliable results for 1:1 and 1:2 electrolytes (including CaCl₂ systems) up to I _m ≈ 13.5 mol kg^–1 (Neck et al., 2009; Fellhauer et al., 2010; Yalcintas et al., 2016; Altmaier et al., 2017).

Based on the slope analysis discussed in Section 3.1 and considering that α-Be(OH)₂(cr) is the solid phase controlling the solubility in all investigated systems, the equilibrium reaction (2) is proposed to describe the solubility of Be(II) in CaCl₂ solutions above pH_m ≈ 9–10 (depending upon salt concentration). As reflected in reaction (2), the number of Ca-atoms participating in the complexation reaction (n) is unknown and must be included in the optimization process. α − Be OH 2 cr + nCa 2 + + 2 H 2 O l ⇔ Ca n B e OH 4 2 n – 2 + 2 H + (2) with log * ⁡ K s , C a n B e O H 4 2 n − 2 ′ = log Ca n B e OH 4 2 n – 2 + 2 ⁡ log H + –  nlog  Ca 2 + (3) l o g *   K s , C a n B e O H 4 2 n − 2 ° = l o g *   K s , C a n B e O H 4 2 n − 2 ′ + log γ C a n B e O H 4 2 n − 2 + 2 log γ H + − n l o g γ C a 2 + − 2 log a w (4) l o g *   K s , C a n B e O H 4 2 n − 2 ′ − Δ z 2 D − 2 log a w = l o g *   K s , C a n B e O H 4 2 n − 2 ° − Δ ε m C l − (5)

Conditional solubility constants, log ^* K′_{s,CanBe(OH) 4} ⁽²ⁿ⁻²⁾, were determined from the experimental solubility data obtained in 0.05, 0.25, 1.0 and 3.5 M CaCl₂ using the minimization function ∑((log [Be]_exp–log [Be]_calc)²)^1/2. Three different chemical models were considered in the calculations, assuming the predominance of the ternary complexes Ca_nBe(OH)₄ ⁽²ⁿ⁻²⁾ with n = 1–3 in the aqueous phase. The SIT-plot was used to derive the solubility constant at the reference state, log ^* K°_{s,CanBe(OH) 4} ^(2n-2), and corresponding SIT ion interaction coefficient/s based on the linear regression (log ^* K′_{s,CanBe(OH) 4} ^(2n-2) + Δz²D—2 log a_w) vs. [Cl⁻] (in molal units). The SIT ion interaction coefficients of the ternary complexes Ca_nBe(OH)₄ ⁽²ⁿ⁻²⁾ were calculated from the slope of the linear regression (–Δε) in combination with ε(H⁺, Cl⁻) = (0.12 ± 0.01) kg∙mol^–1, ε(Ca²⁺, Cl⁻) = (0.14 ± 0.01) kg∙mol^–1 and a _w values as reported in the NEA-TDB (Grenthe et al., 2020). The criteria considered for the selection of the chemical model describing the solubility of Be(II) in CaCl₂ was based on:

(i) minimization of the quality parameter (Δ) calculated as ∑ l o g *   K s , C a n B e O H 4 2 n − 2 ′ , exp − log K s , C a n B e O H 4 2 n − 2 ′ , calc (Fellhauer et al., 2010; Fellhauer, 2013; Yalcintas et al., 2016),

Table 1 shows the conditional solubility constants, log ^* K′_{s,CanBe(OH) 4} ^(2n-2), determined for each chemical model and background electrolyte concentration, together with log ^* K°_{s,CanBe(OH) 4} ^(2n-2), ε(Ca_nBe(OH)₄ ⁽²ⁿ⁻²⁾, Ca²⁺/Cl⁻) and the quality parameter (Δ). The SIT-plot as well as the representation of log ^* K′_{s,CanBe(OH) 4} ^(2n-2) vs. [Cl⁻] for the three chemical models evaluated in this work are shown in Supplementary Figures SI-S1, S2 of the Supporting Information.

Table 1 shows very large values of the quality parameter (Δ) for the chemical models including the aqueous complexes CaBe(OH)₄(aq) and Ca₃Be(OH)₄ ⁴⁺, which reflect the significant deviations between experimental and calculated log *K′_{s,CanBe(OH) 4} ^(2n-2) for these systems (see also Supplementary Figures SI-S1 in the Supplementary Material). Note that although values of ε(i, j) = 0 for neutral species are considered per definition in SIT, a number of neutral species with ε(i, j) ≠ 0 are reported in the NEA-TDB books (see for instance Grenthe et al., 2020), and this cannot be considered as criteria to reject a specific chemical model. The chemical model including the complex Ca₂[Be(OH)₄]²⁺ shows the lowest value of Δ, a reasonable SIT coefficient, as well as a smooth shape of the plot log ^* K′_{s,CanBe(OH) 4} ^(2n-2) vs. [Cl⁻] (see Supplementary Figures SI-S2). This outcome is in agreement with the MD calculations summarized in Section 3.3, which predict the predominance of an aqueous complex with two Ca-atoms. The same model calculations described above but excluding solubility data in 3.5 M CaCl₂ were also performed in order to remain within the range of ionic strength normally considered for SIT. The best fit was obtained again for the chemical model including the complex Ca₂[Be(OH)₄]²⁺ (see Supplementary Figures SI-S3; Supplementary Tables SI-S4). Consistent values of log *K°_{s,Ca2[Be(OH) 4}]²⁺ and ε(Ca₂[Be(OH)₄]²⁺, Cl⁻) were obtained using both datasets, and thus the model derived including also solubility data in 3.5 M CaCl₂ was finally favored.

The chemical reaction (6) and corresponding thermodynamic functions are thus considered for the extension of the thermodynamic and activity models reported in Çevirim-Papaioannou et al. (2020) to the system Be²⁺–Ca²⁺–Na⁺–K⁺–H⁺–Cl^––OH^––H₂O(l)^. α − Be OH 2 cr + 2 Ca 2 + + 2 H 2 O l   ⇔   Ca 2 B e OH 4 2 + + 2 H + (6) log ⁡ * K s , C a n B e OH 4 2 n − 2 ° = – 25.13   ±   0.27 ε Ca 2 Be OH 4 2 + , Cl – = 0.00   ±   0.15  kg ∙ mol – 1

Figures 4A–D shows the excellent agreement between the experimental solubility data determined in 0.05–3.5 M CaCl₂ solutions with solubility calculations conducted using the thermodynamic and activity models derived in Çevirim-Papaioannou et al. (2020) and extended in this work. Future work will target the development of a Pitzer activity model for this system, which is clearly favored for the description of brine systems.

C-S-H phases of various compositions determine the most relevant properties of cement pastes in ordinary Portland and blended systems. The non-ideal multisite CASH + solid solution model (Miron et al., 2022a; Miron et al., 2022b; Kulik et al., 2022) can accurately model the C-S-H solubility, water content and elemental uptake accounting for the continuous change of solid and pore solution compositions as a function of pH, alkali concentration, and Ca/Si ratio in the system. In the present work, the CASH + model is used to compute the evolution of a pore solution composition assuming the degradation of a Portland cement in a repository setting defined by low ionic strength conditions. This is done by removing alkalis and of Ca from the system (Ca-Si-Na-K-H₂O) in equilibrium with C-S-H to produce the change of the pore solution composition (Figures 5A, B) as a function of [H⁺], present at different cement degradation stages. At the starting point the pore solution composition in equilibrium with a hydrated Portland cement is taken from (Vollpracht et al., 2016) and has a pH = 13.40 and [NaOH + KOH] = 0.555 M (with K/Na = 2). Alkalis are then removed and the pH decreases to a value of 12.4, when the system is in equilibrium with high Ca/Si C-S-H and portlandite. In natural systems, the pore solution composition is buffered at this condition until all portlandite is dissolved. By removing Ca from the system, after the full dissolution of portlandite, the composition of C-S-H goes towards lower Ca/Si ratios and the pH of the solution decreases to 9.8 at which point the system is in equilibrium with low Ca/Si C-S-H and amorphous silica. The evolution of the pore water composition obtained from calculations with the CASH + solid solution model in combination with the thermodynamic and activity models derived in this work and reported in Çevirim-Papaioannou et al. (2020) have been used to calculate the aqueous speciation of beryllium throughout the complete degradation process of cement (see Figure 5C).

Figure 5C shows that the hydrolysis species Be(OH)_3⁻ and (to less extent) Be(OH)₄ ²⁻ prevail within the degradation stage I of cement, characterized by very high pH and low Ca concentration. Due to the enhanced Ca concentration defined by the equilibrium with portlandite, the ternary complex Ca₂[Be(OH)₄]²⁺ dominates the aqueous speciation of Be(II) within the degradation stage II, as well as in the early steps of the degradation stage III at pH > 12.1. With the decrease in Ca concentration, the hydrolysis species Be(OH)_3⁻ becomes predominant again at 10.6 ≤ pH ≤ 12.1, in systems controlled by C-S-H phases with Ca:Si = 0.80–1.27. Below pH ≈ 10.6, the neutral species Be(OH)₂(aq) becomes predominant until the full degradation of cement.

Note that CASH + has not been yet verified for calculating systems with concentrated CaCl₂ solutions. In particular, the evolution of C-S-H and other cement phases at high Ca concentrations remains ill-defined. Although this contribution has focused on the SIT activity model, the development of the corresponding Pitzer activity model for the system Be²⁺–Ca²⁺–Na⁺–K⁺–H⁺–Cl^––OH^––H₂O(l) is foreseen in the next phase of this study. Pitzer formalism is clearly favored for the description of high-saline systems. Efforts dedicated to the thermodynamic modelling of cementitious systems in high saline conditions are also on-going in the context of the THEREDA Reference Database project (Moog et al., 2015).