Immobilization of Ca²⁺ by OH⁻ ions on calcium silicate hydrates

Abstract

The interface of calcium silicate hydrate (C–S–H) plays a a crucial role in the mobility of ions through hardened cement paste and thus in the long term durability of concrete structures. Here, we combine experimental zeta potential measurements with non-equilibrium molecular dynamics simulations to elucidate the surface state of C–S–H under cement-relevant conditions. Our findings reveal that the co-adsorption of Ca²⁺ and OH⁻ ions overcompensates the negative surface charge of C–S–H, originating from deprotonated silanol groups, thereby forming a surface adsorption network. This mechanism is governed by the interplay between solution pH and calcium speciation, where higher pH values favor the formation of calcium–hydroxide solution complexes. Consequently, an increase in pH leads to a higher Ca/Si ratio in C–S–H and a greater adsorption of these complexes, further increasing the effective surface charge density. Conversely, a reduction in pH promotes ion desorption, leading to a progressive decrease in the Ca/Si ratio and decrease in the effective surface charge. These findings provide new insights into the long-term compositional evolution of C–S–H in cement paste, with implications for phase stability, mechanical performance, and concrete durability.

---

1. Introduction

Calcium silicate hydrate (C–S–H), the main hydration product of Portland cement, is the most abundant man-made material on the planet.¹ In hydrated systems, the interactions at the C–S–H-solution interface influence ion transport and consequently the lifetime of concrete structures.² To decrease the environmental impact of concrete and to increase its durability, a fundamental understanding of the interfacial phenomena at the C–S–H interface is needed.³

The C–S–H particles have a nanofoil morphology. The nanofoils are 100–200 nm long and below 5 nm thickness with a specific surface area of 200–300 m² g⁻¹.^4–7 At the atomic level, the C–S–H structure resembles the one of tobermorite 14 Å.^8,9 The nanofoils are constituted of calcium-silicate layers stacked in the c-axis direction. The mainlayer consists of a calcium oxide sheet flanked by defective silicate chains. The silicate chains consist mainly of dimers with some longer chain units. The mainlayers are separated by an interlayer containing water and ions. This arrangement results in C–S–H having dominant basal (001) surfaces in the ab-plane.⁴ While the bulk C–S–H structure is today well understood,^8,10,11 the surface structure is yet unexplored which hinders the interpretation and understanding of interfacial phenomena.^3,4,12

Casar et al. built an atomistic model of a high Ca/Si ratio C–S–H nanofoil with its basal (001) surfaces that matches experimentally measured properties (Ca/Si, proportion of silicate species, Si–OH/Si, and Ca–OH/Ca).^4,11 In order to fit the model to these experimentally measured properties, a substantial amount of Ca²⁺ ions needs to be adsorbed on the surface. Molecular dynamics (MD) simulations show that the amount of adsorbed Ca²⁺ is limited due to the strong hydration shell of Ca²⁺ and its preferential outer-sphere adsorption.^4,13 However, OH⁻ ions can act as a bridge between inner- and outer-sphere adsorbed Ca²⁺,⁴ which was previously observed on gypsum surfaces,¹⁴ effectively immobilizing Ca²⁺ close to the surface and consequently increasing the net surface charge density of the basal (001) surface. Enhanced calcium adsorption with increased pH has been previously reported for silicate materials^15–17 and hydroxyapatite (calcium apatite mineral).¹⁸ Shimanayashi and Nakagaki observed from adsorption isotherms a linear correlation between the adsorption of Ca²⁺ and OH⁻ on hydroxyapatite in high pH solutions.¹⁸

To investigate the surface structure, researchers have synthesized C–S–H through the reaction of CaO or Ca(OH)₂ with reactive SiO₂¹⁹ and measured the zeta potential to understand the structure at and near the C–S–H surface.^20–24 While C–S–H in early Portland cement paste reaches Ca/Si > 1.7, this synthesis method produces phase pure C–S–H with Ca/Si less than 1.4.

Zeta potential is an indirect technique to investigate the structure of the electric double layer (EDL) and consequently the surface composition and adsorbate coverage.^25–28 It is related to the charge state of the surface and infers the physical and chemical properties of the interfacial systems. The technique measures the potential-driven velocities close to the surface from which the electrophoretic or electroosmotic mobilities are calculated and converted to zeta or streaming potentials via various theoretical relations (such as the Helmholtz–Smoluchowski equation).^29,30

Since the exact solution conditions are rarely reported in experimental zeta potential studies, the interpretation of results is often difficult. However, some general trends for low (<1.4) Ca/Si C–S–H are seen.^20,23,24,31 The reported zeta potential values are consistently negative for low Ca/Si C–S–H, suggesting a negative surface charge. As calcium is added to the solution, the zeta potential gradually increases, which suggests immobilization of Ca²⁺ ions at the C–S–H surface.²⁵ When calcium is added as Ca(OH)₂ to the solution, the zeta potential becomes positive at approximately 2 mM of Ca(OH)₂ and increases further with the addition of Ca(OH)₂.^20,24,31 When CaCl₂ is added to the initial solution, the zeta potential again gradually increases, however, it remains negative.^24,32 This observation suggests a decisive role of OH⁻ ions in the mechanism of Ca²⁺ adsorption. Positive zeta potentials were exclusively reported for C–S–H in calcium rich solutions with pH above 11.³¹

The negative surface charge of the C–S–H surface originates from the deprotonation of surface silanol groups in the high pH solution.³³ Keen Koo et al. investigated the charge reversal at the C–S–H surface-electrolyte interface by a mean-field molecular theory approach.³⁴ By accounting for steric and electrostatic forces between ions, as well as for co-adsorption of Ca²⁺ and OH⁻ in the form of a calcium complex (CaOH⁺), they successfully predicted the experimentally determined interface electrostatic potential of Labbez et al. for low Ca/Si C–S–H.³³ Further, the results show a decrease in Ca²⁺ mobility at the interface with increasing pH values. It appears that an overcompensation of the negative surface charge by Ca²⁺ adsorption is possible when one accounts for co-adsorption of OH⁻ ions.⁴

Harris et al. standardized and published a method for the synthesis of phase pure high Ca/Si (>1.5) C–S–H.³⁵ This was the first methodological description of a synthesis protocol to produce synthetic phase pure C–S–H with Ca/Si similar to those in real cement pastes (Ca/Si > 1.6). Their C–S–H samples with nominal Ca/Si ratios from 1.8 to 2.0 exhibit higher Ca/Si ratios than expected from thermodynamic modeling. Further, they observed much lower concentrations of calcium in the supernatants than expected from thermodynamic modeling. Since an increase in Ca/Si ratio is achieved with an increase in pH,³⁶ again, a decisive role of OH⁻ ions in the formation of the C–S–H structure (bulk and surface) is suggested. The difference in low and high Ca/Si ratio C–S–H could result in completely different surface states as recently hypothesized by Casar et al.^4,11 According to this hypothesis at low Ca/Si a negatively charged, calcium-depleted (001) C–S–H surface is to be expected, while at high Ca/Si a calcium–hydroxide ion adsorption network is formed and creates a positively charged surface. If true, this implies that the interpretation of interfacial phenomena in real cementitious systems based on low Ca/Si synthetic C–S–H can be misleading.

In this study, phase pure high Ca/Si ratio (>1.5) C–S–H was synthesized following the protocol of Harris et al.³⁵ The zeta potential of samples was measured under various solution conditions. The solution conditions include supernatant solution, pure NaOH solution, and NaOH solution titrated with Ca(NO₃)₂·4H₂O. Since zeta potential and other bulk methods give little insight into local structure and ordering at interfaces, these studies were complemented by non-equilibrium molecular dynamics (NEMD) simulations. With combined experimental and simulation information this article clarifies the mechanism behind the positive zeta potential and gives atomic-level insights into the adsorption phenomena at the C–S–H interface in high pH solutions.

2. Results

2.1. C–S–H dispersed in supernatant

C–S–H samples with Ca/Si of 1.4, 1.75, and 2.0 were synthesized and vacuum filtered (unwashed), or vacuum filtered and washed with an ethanol–water mixture (washed).³⁵ Afterward, 3.1 wt% of the C–S–H samples were dispersed in their respective supernatants and the zeta potential was measured. Two clear trends (Fig. 1A) are observed for increasing Ca/Si ratios. Higher Ca/Si ratios correlate with higher pH of supernatants,^35,36 and with higher zeta potentials.³¹ While the zeta potential after washing decreases, it remains positive, even for Ca/Si = 1.4, suggesting a strong affinity for Ca²⁺ at the C–S–H surface. The difference between the zeta potential of unwashed and washed samples decreases with increasing Ca/Si. This suggests that C–S–H with higher Ca/Si has more ‘stable’ surface states, stronger adsorbed Ca²⁺, which are less affected by external factors, such as washing. However, it appears that there is always some small amount of loosely bound ions, possibly outer-sphere adsorbed, that can be easily removed by ethanol.

Fig. 1 Zeta potential of (A) 3.1 wt% C–S–H samples dispersed in reaction supernatant (pH 12.10, 12.69, and 12.80); (B) of 3.1 wt% C–S–H (unwashed) dispersed in reaction supernatant and titrated with HNO₃ (the annotations show the zeta potential of the respective samples at pH 12), and (C) measured during titration with Ca(NO₃)₂·4H₂O expressed as [Ca²⁺] addition (open symbols show the first measurement after filtration and dispersion in 10 mM NaOH, and full symbols show the measurement after the second filtration and dispersion). The reported Ca/Si is the nominal Ca/Si of precipitated C–S–H samples.

Next, the unwashed samples in the supernatant were titrated with nitric acid (HNO₃), and the pH and zeta potential were measured. Fig. 1B shows that a decrease in the pH results in a decrease of the zeta potential for all investigated C–S–H samples. This shows a correlation between the zeta potential and the solution pH. It seems that a decrease in pH results in Ca²⁺ desorption from the C–S–H surface. The similarity of Ca/Si = 1.75 and 2.0 samples in zeta potential change (slope) with pH suggests very similar initial surface structures that differ from the surface structure of Ca/Si = 1.4. The similar change in zeta potential suggests that C–S–H samples with Ca/Si = 1.75 and 2.0 desorb similar amounts of Ca²⁺ ions, which appear to be higher than the released amount of Ca²⁺ of C–S–H with Ca/Si = 1.4. While the zeta potentials continuously decreased with decreasing pH (Fig. 1B), they remained positive for all samples. The supernatants were titrated with nitric acid until pH 10, below this point C–S–H dissolves rapidly.

To gain a better understanding of the solution effect on the zeta potential, the supernatants at pH 12 (annotations in Fig. 1B) were analyzed by ICP-OES and its composition was used as input parameters for thermodynamic modeling (see SI).

As the solution pH (and, therefore, the Ca/Si ratio) increases, thermodynamic modeling predicts an increase in the concentration of Ca²⁺ and calcium solution complexes (Fig. S8). Ca²⁺ remains the dominant species over the entire range of Ca/Si ratios, which is consistent with the Porbaix diagram (log C-pH) of a 2.5 mM calcium solution (Fig. S9).³⁷ However, the relative amount of Ca(OH)⁺ in comparison to Ca²⁺ increases with an increase in Ca/Si (in pH). The relatively higher proportions of Ca(OH)⁺ at increased Ca/Si and the discrepancy between thermodynamically calculated (that does not account for surface complexation) and ICP-OES measured concentrations suggests increased co-adsorption rates of Ca²⁺ and OH⁻ at the surface.

2.2. Calcium titration in NaOH solution

Calcium titration experiments were carried out to investigate calcium adsorption onto C–S–H surfaces. Washed C–S–H samples with Ca/Si = 1.0 and 1.75 were investigated since they represent low and high Ca/Si C–S–H. The samples were dispersed in 10 mM NaOH solution and titrated with calcium nitrate tetrahydrate solution (Ca(NO₃)₂·4H₂O). Titration was carried out in 10-minute intervals to allow for equilibration before measurements (open symbols in Fig. 1C). After the final titration, the samples were allowed to sit for 15 hours before measuring the zeta potential. Afterward, the samples were again filtered and dispersed in a fresh 10 mM NaOH solution from which the zeta potential measurements were collected (full symbols in Fig. 1C).

Similar trends are observed for both samples. For both samples (Ca/Si = 1.0 and 1.75) after dispersion in 10 mM NaOH positive zeta potentials are measured (0.9 and 3.5 mV), and they gradually increase with the addition of calcium (Ca(NO₃)₂·4H₂O) until a concentration of 13 and 15 mM is reached. With an addition of 13 mM calcium the zeta potential of the Ca/Si = 1.0 sample increased by 1.2 mV, while the Ca/Si = 1.75 increased by 0.7 mV. The initial zeta potential of C–S–H with Ca/Si = 1.75 (after washing and dispersion in 10 mM NaOH pH 12.03) was ∼4 times lower than the measured zeta potential of the same unwashed C–S–H in the supernatant at pH 12 (3.5 mV versus 13.7 mV). Further, the zeta potential of the Ca/Si = 1.75 sample gradually increases until the calcium concentration of 15 mM, after which it shows a sharper, but less linear increase. At the final calcium concentration of 45 mM, the measured zeta potential was 9.5 mV, which is ∼30% less than the zeta potential of the unwashed sample in supernatant at pH 12. This difference could be due to calcium nitrate addition. When titrating with calcium nitrate the pH of the solution gradually decreases which lowers the zeta potential (Fig. 1B).^36,37

It appears that calcium adsorption is only partially reversible under these experimental conditions and that sodium promotes the removal of calcium from the C–S–H surface. These findings are further supported by the re-dispersion of the post-titration filtered C–S–H sample in 10 mM NaOH solution (full symbols in Fig. 1C), where the measured zeta potentials of the Ca/Si = 1.0 and 1.75 samples further decreased in comparison to the first dispersion in NaOH at the start of the titration experiment. After the second dispersion (full symbols in Fig. 1C) in NaOH, the measured zeta potential for Ca/Si = 1.0 sample was −1.1 mV, while the Ca/Si = 1.75 sample had 0.7 mV (after the first dispersion, open symbols in Fig. 1C, the measured zeta potentials were 0.9 and 3.5 mV).

To ensure the C–S–H samples were not altered during the titration experiments, part of the post-titration sample was freeze-dried, and X-ray diffraction (XRD) data was collected. XRD patterns (Fig. S5) of both samples are typical for C–S–H.^8,9,24,35

2.3. Molecular dynamics

To obtain an atomic-level insight into the phenomena leading to changes in the zeta potential of C–S–H, non-equilibrium molecular dynamics (NEMD) simulations were carried out.^25,26 The studied surface was a basal (001) C–S–H surface (2.8 nm × 4.5 nm).^4,38 93% of the surface silanol groups were deprotonated, which is in the expected range for C–S–H in high pH, calcium-rich solutions.³³ This level of deprotonation results in a surface charge density of −4.46 e nm⁻², i.e. −0.72 C m⁻². The negative surface charge was compensated by Ca²⁺ ions in the proximity of the surface (2.23 Ca²⁺ per nm² less than 1 nm away from the surface). Two such model surfaces were separated by a ∼6 nm slab which contained 2782 H₂O molecules and 7 CaCl₂ units (0.156 M).

In addition to the simulated systems described above, three more systems were investigated. In the additional systems, an integer number of one Ca²⁺ and two OH⁻ ions were initially randomly placed in the proximity of each (001) surface (Table S3). This addition is expressed as surface adsorbed OH⁻ ions per nm² in the following discussion. With this, a co-adsorption of Ca²⁺ and OH⁻ was approximated. In all simulations, the bulk solution concentration was kept constant (0.156 M CaCl₂, as explained in the SI) as an analogy to the experimental C–S–H samples where zeta potential were measured at pH 12 in the supernatant (vertical line in Fig. 1B). Simulation and model details are given in the SI.

The zeta potentials shown in Fig. 2A were calculated from the water mobilities (Fig. S11) as described by Předota et al.²⁵ It is seen that an increase in the OH⁻ ion surface concentration correlates with an increase in the zeta potential and the effective surface charge density. The effective surface charge density is calculated as the base surface charge due to deprotonation of silanol groups (−4.46 e nm⁻²) plus the charge due to immobile Ca²⁺ and surface adsorbed OH⁻ ions at the C–S–H surface. While both the zeta potential and the adsorbed amount of ions accounted in the effective surface charge density are physical and chemical properties of the interface related to the charge state of the surface, zeta potential is an electrokinetic property that depends on the mobility of charged species at the interface, not only their numbers.²⁵ Typically zeta potential reflects the sign of the surface charge, but exceptions can occur due to ion adsorption, solvent effects, and interfacial EDL structure. This can lead to counter-intuitive results as seen in Fig. 2A at OH⁻ ions per nm² = 0.32, where the zeta potential is positive (1.01 mV) while the effective surface charge density remains negative (−0.2 e nm⁻²). The here calculated effective surface charge was determined from the immobile ions defined by their MSD below a chosen threshold. Therefore, some Ca²⁺ ion that was treated as mobile (not contributing to the effective surface charge) under this assumption could have been immobile, possibly making the effective surface charge positive, in line with the positive zeta potential, which is determined directly from the relative motion of the surface and bulk solution, without the need of any assumptions on the (im)mobility of interfacial species. Regardless, the effective surface charge is a useful concept for the discussion of interfacial behavior and explaining the experimental trends.

Fig. 2 (A) Zeta potential and effective surface charge density as predicted by NEMD as a function of surface concentration of surface adsorbed hydroxide ions (OH⁻ ions per nm²). (B) Schematic representation of calcium sites at the CSH interface. Dark blue triangles are the surface silicate tetrahedra, the small red circle is an OH⁻ ion and water molecules are represented by red lines. The large circles are Ca²⁺ ions which are coordinated to silicate oxygens (I(O)), silicate oxygen and OH⁻ (I(O + OH)), OH⁻ (I(OH)), or are immobile outer-sphere adsorbed (I(⊘))), or mobile (M). (C) and (D) Show the Ca²⁺ concentration of different calcium sites. (E) Average distance (z_mean) of ions from the silicate chains (perpendicular to the surface).

To better understand the role of surface adsorbed OH⁻ ions at the C–S–H interface, an analysis of calcium sites and their relationship to the surface and OH⁻ ions was carried out. The traditional classification of ions in the electric double layer as inner-sphere, outer-sphere, and diffuse carries limited to no significance when studying the zeta potential since this classification is based on the structural and not dynamic (mobility) information.^25,26

Fig. 2B (Fig. S12 and S13) shows five unique calcium sites at the C–S–H surface which are clustered in two categories depending on their mobility. The mobility categories are established due to two distinct mean square displacement (MSD) values of the Ca²⁺ ions (Fig. S14). Ca²⁺ ions with high MSD values are fully mobile (M in Fig. 2B), similar to the traditional diffuse layer. Ions with low MSD are immobile (I) close to the surface and are further classified into 4 categories, where coordination of Ca²⁺ to a particular oxygen is considered if the distance between them is less than 3.15 Å:

• I(O): Ca²⁺ which coordinates to at least one silicate oxygen (traditional inner-sphere adsorption), but does not coordinate with an OH⁻ ion.

• I(OH): Ca²⁺ which coordinates with an OH⁻ ion, but not with a silicate oxygen.

• I(O + OH): Ca²⁺ which coordinates with a silicate oxygen and a OH⁻ ion.

• I(⊘): Ca²⁺: which coordinates neither with a silicate oxygen or OH⁻ ions (traditional outer-sphere adsorption).

The surface concentration of Ca²⁺ ions depending on their mobility is shown in Fig. 2C. At 0 OH⁻ ions per nm², approximately 20% of Ca²⁺ are mobile and 80% immobile. The charge contribution from immobile Ca²⁺ sites equals 4.1 e nm⁻², which is less than what is needed to charge balance the negative surface charge due to deprotonation of silanol groups (−4.46 e nm⁻²). This shows that the C–S–H surface at 93% deprotonation rate is not charge compensated by solely immobile Ca²⁺ ions. As OH⁻ per nm² increases, so does the proportion of immobile Ca²⁺ sites.

An interesting feature of the C–S–H surface emerges when the immobile calcium sites of the OH⁻ ions per nm² = 0 system are analyzed. Fig. 2D shows that roughly 50% of Ca²⁺ sites are I(⊘), traditionally considered outer-sphere adsorbed. While these sites are generally considered more mobile than inner-sphere adsorbed I(O) sites, they are immobile on the C–S–H surface. Such immobility of an outer-sphere adsorbed ion is presumably due to the strong affinity of Ca²⁺ to the negatively charged surface and the geometry of the uneven surface of C–S–H where silicate chains form open channels on the surface. In these channels, Ca²⁺ with a full H₂O hydration shell can be trapped without coordinating to a silicate oxygen (Fig. S12).⁴ This shows that an additional mechanism participates on immobilization of Ca²⁺ at the C–S–H interface to overcome the negative surface charge.

With an increase in OH⁻ ions per nm² the proportion of immobile sites (I) increases while the proportion of mobile (M) Ca²⁺ decreases (Fig. 2C). With an increase in the surface concentration of OH⁻ ions the I(O) site gradually transforms to I(O + OH). Simultaneously the mobile (M) and immobile I(⊘) sites are immobilized by the surface adsorbed OH⁻ ions of I(O + OH) sites and thus transforming to I(OH) sites. During this process the immobile ion region (I(O), I(OH), I(O + OH), I(⊘), and OH⁻) close to the surface grows, increasingly overcompensating the negative surface charge (Fig. 2A and Fig. S13), which repels the remaining mobile Ca²⁺ that are not immobilized further away from the surface (Fig. 2E). While doing so, the negatively charged anions (Cl⁻ in NEMD) occupy the space between the mobile and immobile calcium regions and become the ion species with dominant mobility close to the surface. Therefore, the mobile anions become the primary source of induced water mobilities which give rise to positive zeta potentials.²⁵

As seen in Fig. 2A, the effective surface charge density increases from −0.36 e nm⁻² to 1.01 e nm⁻². In all simulations, the base surface charge density due to deprotonation of surface silanol groups was kept at −4.46 e nm⁻². Immobile Ca²⁺ at the C–S–H surface (I(O) and I(O + OH) in Fig. 2B) contribute between 3.06 and 4.03 e nm⁻² to the effective surface charge density, which is less than the charge needed for charge compensation. Therefore, the remaining positive charge contribution comes from immobile ions that do not coordinate with the silicate chain of the C–S–H surface. Instead, they immobilize due to the bridging function of OH⁻ ions. This can happen through three different mechanisms:

• An OH⁻ adsorbs and coordinates to an immobile surface adsorbed Ca²⁺ and later immobilizes a mobile Ca²⁺;

• A calcium hydroxide ion complex (Ca(OH)⁺ or Ca(OH)_2,aq) adsorbs at the C–S–H surface, coordinating with its Ca²⁺ to an oxygen of the silicate chains. Its OH⁻ then immobilizes a mobile Ca²⁺;

• A calcium hydroxide ion complex is immobilized due to its OH⁻ coordinating to an immobile surface adsorbed Ca²⁺.

The majority of Ca²⁺ immobilizes (adsorbs) because of OH⁻ ions (Fig. S13 and S15) and not because of deprotonated silanol groups, as assumed in literature.^33,34 The inspection of the formed Ca²⁺–OH⁻ surface networks shows that one OH⁻ ion can simultaneously coordinate to up to three Ca²⁺ ions (Fig. S13). This underlines the strong influence of an OH⁻ ion in the immobilization of Ca²⁺, and consequently the formation of the Ca²⁺–OH⁻ surface network.^4,11 Through this co-adsorption the experimentally measured high zeta potential as well as the discrepancy between ICP-measured ion concentrations and those predicted by thermodynamic modeling, which does not include surface complexation, can be explained.

3. Discussion

The increase in OH⁻ ions per nm² (and the corresponding co-adsorption of Ca²⁺) of the investigated NEMD systems can be interpreted as an increase in the Ca/Si ratio of C–S–H.^4,11 Experimentally an increase in the Ca/Si ratio of precipitated C–S–H is achieved with an increase in solution pH due to the addition of NaOH which implies the role of OH⁻ on surface calcium complexation.^35,36 As seen from thermodynamical modeling (Fig. S8) and the Pourbaix diagram (Fig. S9),³⁷ the calcium speciation increasingly shifts from Ca²⁺ to calcium complexes (Ca(OH)⁺ and Ca(OH)_2,aq) as the pH increases.

NEMD predicted zeta potentials (Fig. 2A) show a correlation between the increase in the zeta potential with an increase in Ca²⁺ and OH⁻ co-adsorption, i.e. with an increase in Ca/Si ratio of C–S–H due to a modified surface structure.^4,11 One OH⁻ ion can coordinate and immobilize up to three Ca²⁺. This mechanism increases the effective surface charge density (Fig. 2A). As the effective surface charge density becomes more positive, by cation immobilization, the anions become the dominant species of the mobile, diffuse layer of the electric double layer (Fig. 2E). Due to the applied electric field in zeta potential measurements, the mobile anions drag the water molecules in their hydration shells with them in the opposite direction of the electric field. The movement of these water molecules induces electrophoretic water mobilities that give rise to positive zeta potential.^25,26 Therefore, NEMD simulations show that the coadsorption of Ca²⁺ and OH⁻ increases the effective surface charge density that causes an increase in the concentration of mobile anions close to the interface (in the electric double layer), consequently induces higher water mobilities and therefore higher measured zeta potentials.

The initially measured negative zeta potential of Yoshida et al. (Fig. S16) on C–S–H with Ca/Si = 1.34 in 20 mM NaOH solution (pH 12.3) can be explained as follows.²⁴ To synthesize a C–S–H sample with such low Ca/Si, a pH of 12 or less is sufficient.³⁵ At such pH, the calcium speciation strongly favors calcium as a Ca²⁺ ion over calcium solution complexes (Fig. S8 and S9). Therefore, the negative surface charge density due to the deprotonation of surface silanol groups is charge compensated by Ca²⁺. As seen from our NEMD simulation with OH⁻ ions per nm² = 0, the C–S–H surface cannot immobilize all of the charge compensating Ca²⁺. Therefore, the charge compensation extends into the mobile layer, which is cation dominated at these conditions. The cations with their hydration shells move in the direction of the external electric field, which results in positive water mobilities that translate to negative zeta potentials.²⁵ As the initial solution is titrated with a Ca(OH)₂ solution, the introduction of solution calcium leads to formation of CaOH⁺ and Ca(OH)_2,aq complexes. While mobile Ca²⁺ cannot approach the surface due to electrostatic repulsion, NEMD indicates that the hydroxide ion containing calcium complexes can overcome the energetic barrier and immobilize at the surface. While any calcium complex can immobilize through Ca²⁺ coordination to a deprotonated silanol group, it can also become immobilized by coordination of its OH⁻ to an already immobilized Ca²⁺, also observed on Portlandite surfaces.³⁹

Similarly, the increase in the zeta potential with the increase in Ca/Si of C–S–H in supernatant solution (Fig. 1B) is because of the increase in the effective surface charge density with increasing Ca/Si ratios. The increase in the Ca/Si ratio of precipitated C–S–H is achieved through the increase in the solution pH. Due to an increase in pH, proportionally more calcium complexes are present during nucleation and growth, which can become immobilized at the C–S–H surface and increasingly overcompensate the negative surface charge originating from deprotonated silanol groups.

The relevance of OH⁻ ions at the C–S–H interface is also observed from titration experiments. Reported titration of C–S–H–NaOH with a Ca(NO₃)₂ solution (Fig. 1C) is analogous to the CaCl₂ titration experiment of Yoshida et al.²⁴ In both cases, a slow increase in zeta potential with the addition of calcium solution is observed. This increase shows the sensitivity of the interfacial phenomena to the solution conditions. Titration of C–S–H–NaOH decreases the pH due to the mildly acidic nature of Ca(NO₃)₂,³⁶ while titration with Ca(OH)₂ increases it.⁴⁰ The rapid increase of the zeta potential from Yoshida et al. with the addition of Ca(OH)₂ can be explained by a higher abundance of initial and added OH⁻ ions,²⁴ coupled with a higher proportion of calcium solution complexes due to higher pH,³⁷ which preferentially adsorb at the C–S–H surface.

The dispersion of filtered C–S–H samples in a 10 mM NaOH (pH 12) solution revealed a decrease in the zeta potential in comparison to the values measured in the supernatant at pH 12 (3.5 mV versus 13.7 mV for Ca/Si = 1.75). Previous studies suggest competition of Ca²⁺ and Na⁺ ions at the C–S–H interface.⁴¹ It appears that calcium can be removed from the C–S–H interface in high sodium solution concentrations. Na⁺ has a hydration radius smaller than that of Ca²⁺.⁴² Therefore, the surface concentration of immobile (inner-sphere adsorbed) Na⁺ can theoretically be higher than that of Ca²⁺. However, the ionic charge of Na⁺ is half of that of Ca²⁺. Consequently, the C–S–H interface will always result in lower effective surface charge densities than can be reached solely by immobilization of Ca²⁺ ions. Sodium is also not known to form NaOH_aq solution complexes and therefore cannot increase the effective surface charge through co-adsorption.⁴³ Since a lower effective surface charge density results in more mobile cations at the interface, the water induced mobilities yield lower zeta potentials, as observed in Fig. 1C.

Some previously published models have considered co-adsorption of Ca²⁺ and OH⁻ ions and have successfully captured the experimentally measured surface properties (interface electrostatic potential and zeta potential) of low Ca/Si C–S–H.^24,34 However, the models went only as far as to consider the adsorption of Ca(OH)⁺ ions on deprotonated silanol groups. But to reach Ca/Si ratio of real cementitious systems (>1.5), higher degrees of adsorption are necessary.^4,11 Fig. 3 shows two nanofoils with an identical bulk structure (Ca/Si_bulk = 1.43) that is in agreement with experimentally measured properties.¹¹ When zero calcium adsorption at the (001) surface is considered this nanofoil has a Ca/Si = 1.33. Considering a calcium adsorption of 7.2 Ca²⁺ per nm² (and the corresponding co-adsorption of OH⁻), the Ca/Si ratio of the same nanofoil increases to 1.77. Such degree of calcium adsorption results in a nanofoil structure where 25% of all calcium atoms are found adsorbed (immobile) on the surfaces and the remaining 75% are located in the bulk structure (mainlayer and interlayer) of the C–S–H nanofoil. It was previously theorized that nanolayers of calcium hydroxide could be present in the C–S–H interlayer.^9,44 However, the layer spacings determined experimentally appear to be too small for the inclusion of a large number of Ca²⁺ and OH⁻ ions.^6,11 Therefore, it is more likely that these experimental observations, that lead to the hypothesis of calcium hydroxide nanolayers being intermixed with C–S–H, originate from the disordered Ca²⁺–OH⁻ surface adsorption network.

Fig. 3 The increase in the Ca/Si ratio of the C–S–H nanofoil model due to co-adsorption of Ca²⁺ and OH⁻ ions. The nanofoil model was constructed with a mean chain length of 3, full calcium occupancy of vacant interlayer bridging sites, and an additional packing of Ca_Int/dimer = 0.7. This configuration results in a bulk C–S–H structure with Ca/Si_bulk = 1.43.^4,11 Color legend: Si - dark blue, Ca - torques, O - red, H - white. Water molecules are hidden for clarity.

4. Conclusion

In summary, the here demonstrated results from NEMD in combination with experimentally measured zeta potentials on C–S–H with high Ca/Si ratios can be summarized as:

• Positive zeta potentials do not originate solely from the adsorption of Ca²⁺ to the negatively charged silicates of the C–SH surface.

• The effective surface charge density increases, and becomes positive, due to co-adsorption of Ca²⁺ and OH⁻ ions in the form of individual solution species or as calcium solution complexes (Ca(OH)⁺ and Ca(OH)_2,aq).

• The co-adsorption potential of these ionic species appears tremendous and can be a defying factor for the increase in Ca/Si ratio of C–S–H due to the formation of a surface complexation network (Fig. 3).⁴

Therefore, in real cement paste, where early hydration C–S–H reaches Ca/Si ratios above 1.7, high positive surface charge densities are expected. This could help explain why C–S–H initially grows as ‘low-density’ hydration product (outer C–S–H), and later as ‘high density’ (inner product).⁴⁵ While for both C–S–H phases a nanofoil morphology is observed, the higher porosity and lower packing density of the former could be partially caused by the high surface charge densities of the adsorbed surface network of Ca²⁺ and OH⁻ ions and complexes.^46,47 Such structural differences in C–S–H from the magnitude of the surface charge density could explain some of the observed difference between ordinary Portland cement (OPC) and blended cements. OPC's have high Ca/Si C–S–H whereas blended cements have lower Ca/Si in the C–A–S–H phase linked to the Al content.⁴⁸ C–A–S–H in blended cements reaches lower Ca/Si ratios and its alumina, incorporated at vacant Q sites, does not undergo deprotonation, as do silicate tetrahedra.⁴⁹ Therefore, the surface charge density of C–A–S–H is expected to be lower than that of C–S–H. This lower surface charge density could allow denser formations of the outer product and thus contribute to the differences of gel-pore sizes in OPC C–S–H and C–A–S–H in blended cements.⁵⁰ These differences in pore sizes and surface charge densities could also be expected to influence chloride ingress which is again different in OPC and blended cements.^2,51 With the surface adsorbed network or layer demonstrated in the current work further simulations of ion mobility (chloride, sulphate, and others) and interaction distances between surfaces can now confidently be further explored to enhance our knowledge of ion transport and formation mechanisms of C–S–H and eventually C–A–S–H.

5. Methods

5.1. Calcium silicate hydrate synthesis

C–S–H samples were synthesized in a segmented flow tubular reactor (SFTR),^52,53 with initial concentrations of reactants sodium metasilicate, calcium nitrate, and sodium hydroxide as specified in Table S1. Reactants were pumped into the reactor at a rate of 10 mL min⁻¹, reaction temperatures were maintained at 20 °C, and the solutions were allowed to mix for 3 h before collection. The total reactant volume was 800 mL. 30–40 g were produced per sample and precipitates were either filtered or washed with 60 mL of 1 : 1 pure de-ionized water-ethanol mixture and filtered. The filtration was carried out via a vacuum-filtration system on a 200 nm filter paper (Whatman ME 24/21 STL). To dry, 25 g of a filtered and washed sample were placed in a −80 °C freezer and afterwards freeze-dried for an additional 24 h in an Alpha 1–2 LDPlus freeze dryer operating at −50 °C and 0.01 mbar.

5.2. Acoustophoresis

Zeta potential was measured with electroacoustics by the AccoustoSizer II supplied by Colloidal Dynamics. The equipment was calibrated with potassium tungosilicate (KSiW). To prepare samples for acoustophoresis, a suspension was created by adding 5 g of filtered, wet C–S–H to 155 g of dispersion media, either filtered reaction supernatant or 10 mM NaOH solution. Suspensions were mixed for 15 min with a magnetic stirrer at 500 rpm, placed in an ultrasonic bath for 15 min, and allowed to mix at 500 rpm for another 15 min as the sample cooled to room temperature. Zeta potential measurements were taken immediately after the mixing steps. The measurements were taken in 10 min intervals. Each measurement in the calcium titration experiment (Fig. 1C) was taken three times, while all other measurements were taken only once. As seen in Fig. 1C, the standard deviation is negligible. This is consistent with other zeta potential measurements made using the same instrument on cementitious systems, where a standard deviation of around ±0.06 mV is expected.⁵⁴

5.3. Inductive coupled plasma-optical emission spectroscopy (ICP-OES)

To prepare C–S–H samples for ICP-OES analysis, 0.013 g of freshly prepared, wet C–S–H was dissolved in 1 mL of 65% nitric acid and placed in an ultrasonic bath for 30 min. 9 mL of 0.1 M diluted aqua regia was added to each sample. The solution was placed in the ultrasonic bath for an additional 15 min and then placed in a 50 °C oven for 24 hours. Desired dilutions were prepared with decarbonated water. Acidified precipitate samples exhibited maximum concentrations of Ca, Si, and Na at 200, 100, and 100 mg L⁻¹, respectively.³⁵

To prepare a reaction supernatant sample for ICP, 2 mL of supernatant was acidified in 3 mL of 65% nitric acid. The solution was placed in an ultrasonic bath for 30 min. 4-Fold, 10-fold, and 20-fold dilutions were made using decarbonated water. Acidified supernatant samples exhibited maximum concentrations of Ca, Si, and Na at 70, 5, and 2000 mg L⁻¹, respectively.

All diluted samples were stored in 15 mL polypropylene tubes and sent to Heidelberg Cement for analysis. ICP was carried out on an ICPE-9000 Shimadzu instrument in optical emission spectroscopy mode (ICP-OES).

5.4. Thermodynamic modeling

Thermodynamic modeling was done using the Gibbs free energy minimization software, GEM-Selektor, and PHREEQC Version 3.^55–57 An expanded Cemdata-18 database from Andalibi et al. was used including the CSHQ thermodynamic solid solution model, PSI-Nagra, redox-uncoupled gases from the PHREEQC database, and NO³⁻ & Na⁺ phases from the Lawrence Livermore National Laboratory (LLNL) database.^58–60

5.5. Molecular dynamics

The atomistic nanofoil model with the (001) basal surface was constructed with the brick code.^4,11,38 The lateral x × y size of the C–S–H nanofoil measured 2.79 × 4.46 nm. In z-axis direction, the nanofoil measured five mainlayers, four interlayers, approximately 6.5 nm. Half of the C–S–H nanofoil was placed at the bottom of the simulation box, while the other half was placed at the top of the simulation box in z-axis direction. The two (001) basal surfaces were separated by a 6 nm solution slab. The two surfaces had an identical structure and distribution of defects. The structure was constructed as a C–S–H with a mean chain length of 3, which equals to 75% vacant Q^2b sites. The Q^2b were randomly removed from the silicate chains and replaced by Ca²⁺ ions. 93% of surface silanol groups were deprotonated.³³ The negative surface charge was charge balanced by Ca²⁺ ions in the proximity of the surface. The bulk solution was modeled as a 0.156 M CaCl₂ solution. The solution concentration remained the same for all studied systems. To account for the effect of Ca²⁺ and OH⁻ co-adsorption at the C–S–H (001) surface, three variations of the above-described C–S–H nanofoil structure were made as described in the text and SI. The atomistic nanofoil model with the (001) basal surface was constructed with the brick code.^4,11,38 The lateral x × y size of the C–S–H nanofoil measured 2.79 × 4.46 nm. In z-axis direction, the nanofoil measured five mainlayers, four interlayers, approximately 6.5 nm. Half of the C–S–H nanofoil was placed at the bottom of the simulation box, while the other half was placed at the top of the simulation box in z-axis direction. The two (001) basal surfaces were separated by a 6 nm solution slab. The two surfaces had an identical structure and distribution of defects. The structure was constructed as a C–S–H with a mean chain length of 3, which equals to 75% vacant Q^2b sites. The Q^2b were randomly removed from the silicate chains and replaced by Ca²⁺ ions. 93% of surface silanol groups were deprotonated.³³ The negative surface charge was charge balanced by Ca²⁺ ions in the proximity of the surface. The bulk solution was modeled as a 0.156 M CaCl₂ solution. The solution concentration remained the same for all studied systems. To account for the effect of Ca²⁺ and OH⁻ co-adsorption at the C–S–H (001) surface, three variations of the above-described C–S–H nanofoil structure were made as described in the text and SI. The atomistic nanofoil model with the (001) basal surface was constructed with the brick code.^4,11,38 The lateral x × y size of the C–S–H nanofoil measured 2.79 × 4.46 nm. In z-axis direction, the nanofoil measured five mainlayers, four interlayers, approximately 6.5 nm. Half of the C–S–H nanofoil was placed at the bottom of the simulation box, while the other half was placed at the top of the simulation box in z-axis direction. The two (001) basal surfaces were separated by a 6 nm solution slab. The two surfaces had an identical structure and distribution of defects. The structure was constructed as a C–S–H with a mean chain length of 3, which equals to 75% vacant Q^2b sites. The Q^2b were randomly removed from the silicate chains and replaced by Ca²⁺ ions. 93% of surface silanol groups were deprotonated.³³ The negative surface charge was charge balanced by Ca²⁺ ions in the proximity of the surface. The bulk solution was modeled as a 0.156 M CaCl₂ solution. The solution concentration remained the same for all studied systems. To account for the effect of Ca²⁺ and OH⁻ co-adsorption at the C–S–H (001) surface, three variations of the above-described C–S–H nanofoil structure were made as described in the text and SI.

The molecular dynamics calculations were carried out with LAMMPS and the CementFF4 force field with a timestep of 0.28 fs.^61,62 First, the simulation box was relaxed in the isothermal–isobaric (NPT) ensemble for 10 ns at 297.15 K and 1 atm with periodic boundary conditions in all axes directions. Afterward, a fixed wall boundary was placed in z-axis direction to prevent self-induced polarizability effects,⁶³ and the atoms of the bottom and top two mainlayers and one interlayer were frozen. Only the upper half of the C–S–H nanofoil and solution were allowed to move. The simulation box was then further relaxed for 50 ns in the canonical (NVT) ensemble at 297.15 K.

The production run was carried out by non-equilibrium molecular dynamics for 11.2 ns. The simulation box was set up the same as during the equilibration NVT run. For the production run, an electric field of 2 × 10⁸ V m⁻¹ was prescribed in the y-axis direction affecting the solution H₂O molecules, Ca²⁺, Cl⁻ and OH⁻ ions, and the surface OH⁻ silanol groups. The simulation was carried out in the NVT ensemble at 297.15 K, however, the velocities were thermostated only in the x- and z-axis direction. For simulation details, we refer the reader to the works of Předota et al.^25,26

The electroosmotic mobilities of water in the y-axis direction, μ, were averaged over 5.6 ns on the fly during the calculations. The trajectories were saved every 14 ps and were analyzed with MDAnalysis and Travis.^64–67 MDAnalysis was used to calculate the mean square displacement of the individual Ca²⁺ and analyze their first shell environment. Travis was used to generate the density profiles along z-axis direction. All simulation and postprocessing files are given as part of the SI.

Author contributions

Z. C., L. M. and M. P. designed and carried out the simulations. M. H. carried out the experiments. Z. C. and L. M. analyzed the data. K. S., H. C. G., M. P. and P. B. supervised the work. All authors discussed the results and participated in the writing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

All LAMMPS input files and files for postprocessing are provided as part of supplementary information (SI). Supplementary information: a PDF with further details of experimental conditions and molecular dynamics simulations details is provided. A zip file containing the postprocessing scripts and simulation files is provided. See DOI: https://doi.org/10.1039/d5cp02821c.

Acknowledgements

The Swiss National Science Foundation is acknowledged for supporting the contributions of Z. C. (grant no. 179076 and 217756). The Horizon 2020 Marie Sklodowska-Curie Innovative Training Network Program (Grant Agreement No. 764691) is acknowledged for supporting the contribution of M. H. This work has made use of the Hamilton HPC Service of Durham University. Prof. Steve Parker, Prof. Franco Zunino and Prof. Aslam Kunhi Mohamed are acknowledged for the fruitful discussions. Dr Gabrielle Anne Sblendorio is acknowledged for her help with Python scripting.

References

1. K. L. Scrivener, T. Matschei, F. Georget, P. Juilland and A. Kunhi Mohamed, Cem. Concr. Res., 2023, 174, 107332.
2. K. De Weerdt, W. Wilson, A. Machner and F. Georget, Cem. Concr. Res., 2023, 173, 107287.
3. O. Heinz and H. Heinz, Langmuir, 2021, 37, 6347–6356.
4. Z. Casar, A. Kunhi Mohamed, P. Bowen and K. Scrivener, J. Phys. Chem. C, 2023, 127, 18652.
5. X. Zhang, W. Chang, T. Zhang and C. K. Ong, J. Am. Ceram. Soc., 2000, 83, 2600–2604.
6. I. Richardson, Cem. Concr. Res., 2004, 34, 1733–1777.
7. D. Viehland, J.-F. Li, L.-J. Yuan and Z. Xu, J. Am. Ceram. Soc., 1996, 79, 1731–1744.
8. A. Kumar, B. J. Walder, A. Kunhi Mohamed, A. Hofstetter, B. Srinivasan, A. J. Rossini, K. Scrivener, L. Emsley and P. Bowen, J. Phys. Chem. C, 2017, 121, 17188–17196.
9. I. G. Richardson, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2014, 70, 903–923.
10. E. Duque-Redondo, P. A. Bonnaud and H. Manzano, Cem. Concr. Res., 2022, 156, 106784.
11. Z. Casar, J. López-Zorrilla, H. Manzano, E. Duque-Redondo, A. Kunhi Mohamed, K. Scrivener and P. Bowen, Cem. Concr. Res., 2024, 183, 107593.
12. S. Simone, J. Phys. Chem., 2014, 8007–8013.
13. M. F. Döpke and R. Hartkamp, J. Chem. Phys., 2021, 154, 094701.
14. I. Siretanu, D. Ebeling, M. P. Andersson, S. L. Svane Stipp, A. Philipse, M. Cohen Stuart, D. van den Ende and F. Mugele, Sci. Rep., 2014, 4, 4956.
15. S. M. Ahmed and A. B. Van Cleave, Can. J. Chem. Eng., 1965, 43, 23–26.
16. K. Szymanek, R. Charmas and W. Piasecki, Adsorption, 2021, 27, 105–115.
17. C. Qin, R. Wang and W. Ma, Chem. Eng. J., 2010, 156, 540–545.
18. S. Shimanyashi and M. Nakagaki, Chem. Pharm. Bull., 1983, 31, 2976–2985.
19. R. Maddalena, K. Li, P. A. Chater, S. Michalik and A. Hamilton, Constr. Build. Mater., 2019, 223, 554–565.
20. H. Viallis-Terrisse, A. Nonat and J.-C. Petit, J. Colloid Interface Sci., 2001, 244, 58–65.
21. C. Labbez, A. Nonat, I. Pochard and B. Jönsson, J. Colloid Interface Sci., 2007, 309, 303–307.
22. G. Plusquellec and A. Nonat, Cem. Concr. Res., 2016, 90, 89–96.
23. E. L’Hôpital, B. Lothenbach, G. Le Saout, D. Kulik and K. Scrivener, Cem. Concr. Res., 2015, 75, 91–103.
24. S. Yoshida, Y. Elakneswaran and T. Nawa, Cem. Concr. Compos., 2021, 121, 104109.
25. M. Predota, M. L. Machesky and D. J. Wesolowski, Langmuir, 2016, 32, 10189–10198.
26. D. Biriukov, P. Fibich and M. Predota, J. Phys. Chem. C, 2020, 124, 3159–3170.
27. M. Kosmulski, Adv. Colloid Interface Sci., 2016, 238, 1–61.
28. M. Kosmulski, J. Colloid Interface Sci., 2002, 253, 77–87.
29. A. V. Delgado, F. González-Caballero, R. J. Hunter, L. K. Koopal and J. Lyklema, Pure Appl. Chem., 2005, 77, 1753–1805.
30. S. Qian and Y. Ai, Electrokinetic Particle Transport in Micro-/Nanofluidics, CRC Press, 1st edn, 2012.
31. J. Haas and A. Nonat, Cem. Concr. Res., 2015, 68, 124–138.
32. N. Krattiger, B. Lothenbach and S. V. Churakov, Cem. Concr. Res., 2021, 146, 106491.
33. C. Labbez, B. Jönsson, I. Pochard, A. Nonat and B. Cabane, J. Phys. Chem. B, 2006, 110, 9219–9230.
34. I. K. Koo, C. S. Yaw, Q. Liu, M. N. Chong and K. Goh, Cem. Concr. Res., 2024, 185, 107617.
35. M. Harris, G. Simpson, K. Scrivener and P. Bowen, Cem. Concr. Res., 2022, 151, 106623.
36. H. Matsuyama and J. F. Young, Adv. Cem. Res., 2000, 12, 29–33.
37. Y. Ruan, Z. Zhang, H. Luo, C. Xiao, F. Zhou and R. Chi, Minerals, 2018, 8(4), 141.
38. A. Kunhi Mohamed, S. C. Parker, P. Bowen and S. Galmarini, Cem. Concr. Res., 2018, 107, 221–235.
39. S. Galmarini, A. Kunhi Mohamed and P. Bowen, J. Phys. Chem. C, 2016, 120, 22407–22413.
40. R. G. Bates, V. E. Bower and E. R. Smith, J. Res. Natl Bur. Stand., 1956, 56, 305–312.
41. E. L'Hôpital, B. Lothenbach, K. Scrivener and D. Kulik, Cem. Concr. Res., 2016, 85, 122–136.
42. Y. Marcus, Chem. Rev., 1988, 88, 1475–1498.
43. T. Megyes, S. Bálint, T. Grósz, T. Radnai, I. Bakó and P. Sipos, J. Chem. Phys., 2008, 128, 044501.
44. S. Grangeon, A. Fernandez-Martinez, A. Baronnet, N. Marty, A. Poulain, E. Elkaïm, C. Roosz, S. Gaboreau, P. Henocq and F. Claret, J. Appl. Crystallogr., 2017, 50, 14–21.
45. C. Ruiz-Agudo and H. Cölfen, Chem. Rev., 2024, 124, 7538–7618.
46. S. Masoumi, S. Zare, H. Valipour and M. J. Abdolhosseini Qomi, J. Phys. Chem. C, 2019, 123, 4755–4766.
47. Y. Zhou, X. Sun, K. Zhong, D. G. Evans, Y. Lin and X. Duan, Ind. Eng. Chem. Res., 2012, 51, 4215–4221.
48. F. Zunino and K. Scrivener, Cem. Concr. Res., 2022, 153, 106693.
49. A. Kunhi Mohamed, P. Moutzouri, P. Berruyer, B. J. Walder, J. Siramanont, M. Harris, M. Negroni, S. C. Galmarini, S. C. Parker, K. L. Scrivener, L. Emsley and P. Bowen, J. Am. Chem. Soc., 2020, 142, 11060–11071.
50. E. Berodier and K. Scrivener, Cem. Concr. Res., 2015, 73, 25–35.
51. M. Wyrzykowski, A. M. Gajewicz-Jaromin, P. J. McDonald, D. J. Dunstan, K. L. Scrivener and P. Lura, J. Phys. Chem. C, 2019, 123, 16153–16163.
52. N. Jongen, M. Donnet, P. Bowen, J. Lemaître, H. Hofmann, R. Schenk, C. Hofmann, M. Aoun-Habbache, S. Guillemet-Fritsch, J. Sarrias, A. Rousset, M. Viviani, M. Buscaglia, V. Buscaglia, P. Nanni, A. Testino and J. Herguijuela, Chem. Eng. Technol., 2003, 26, 303–305.
53. A. Aimable, N. Jongen, A. Testino, M. Donnet, J. Lemaître, H. Hofmann and P. Bowen, Chem. Eng. Technol., 2011, 34, 344–352.
54. M. Palacios, M. Kappl, M. Stuer, U. Aschauer, P. Bowen, H. Butt, C. Pecharromán and F. Puertas, Mater. Constr., 2012, 62(308), 489–513.
55. D. A. Kulik, T. Wagner, S. V. Dmytrieva, G. Kosakowski, F. F. Hingerl, K. V. Chudnenko and U. R. Berner, Comput. Geosci., 2013, 17, 1–24.
56. T. Wagner, D. A. Kulik, F. F. Hingerl and S. V. Dmytrieva, Can. Mineral., 2012, 50, 1173–1195.
57. D. L. Parkhurst and C. A. J. Appelo, Description of input and examples for PHREEQC version 3-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, U.S. Geological Survey Techniques and Methods, 2013, book 6, ch. A43, p. 497, available at https://pubs.usgs.gov/tm/06/a43/.
58. M. R. Andalibi, A. Kumar, B. Srinivasan, P. Bowen, K. Scrivener, C. Ludwig and A. Testino, J. Mater. Chem. A, 2018, 6, 363–373.
59. B. Lothenbach, D. A. Kulik, T. Matschei, M. Balonis, L. Baquerizo, B. Dilnesa, G. D. Miron and R. J. Myers, Cem. Concr. Res., 2019, 115, 472–506.
60. W. Hummel, U. Berner, E. Curti, F. J. Pearson and T. Thoenen, Radiochim. Acta, 2002, 90, 805–813.
61. A. Thompson, H. Aktulga, R. Berger, D. Bolintineanu, W. Brown, P. Crozier, P. in't Veld, A. Kohlmayer, S. Moore, T. Nguyen, R. Shan, M. Stevens, J. Tranchida, C. Trott and S. Plimpton, Comput. Phys. Commun., 2022, 217, 108171.
62. Z. Casar, T. B. Montandon, M. Cordova, K. Scrivener, P. Bowen and A. K. Mohamed, Cem. Concr. Res., 2025, 187, 107708.
63. E. Galicia-Andrés, D. Petrov, M. H. Gerzabek, C. Oostenbrink and D. Tunega, Langmuir, 2019, 35, 15086–15099.
64. R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, J. Domanski, D. L. Dotson, S. Buchoux, I. M. Kenney and O. Beckstein, Proceedings of the 15th Python in Science Conference, 2016, 98–105.
65. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
66. M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, J. Chem. Phys., 2020, 152, 164105.
67. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.

---