<figure><img src="image_1.png" /></figure>

HAL
open science

# A review of coarse grained and mesoscale simulations of C–S–H

Katerina Ioannidou, Christophe Labbez, Enrico Masoero

► To cite this version:

Katerina Ioannidou, Christophe Labbez, Enrico Masoero. A review of coarse grained and mesoscale simulations of C–S–H. Cement and Concrete Research, 2022, 159, 10.1016/j.cemconres.2022.106857. hal-03698585

HAL Id: hal-03698585

https://hal.science/hal-03698585

Submitted on 18 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# A review of coarse grained and mesoscale simulations of C–S–H☆

Katerina Ioannidou<sup>a</sup>, Christophe Labbez<sup>b</sup>, Enrico Masoero<sup>c, *</sup>

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

### Keywords:

Mesoscale

Simulations

Modelling

Cement

## ABSTRACT

The nano-to-micro mesoscale is crucial for cementitious materials; here reactions and interactions between molecules produce complex mechanisms that determine the behavior of cement minerals, especially C-S-H. This manuscript reviews the current state of the art in coarse-grained and mesoscale simulations of C-S-H. These simulations leverage a rigorous statistical mechanical framework, linking atomistic description with coarse-grained modelling through several pivotal concepts: potential of mean force, ion-ion correlations between charged surfaces, and grand canonical reactive ensemble. The second part of the manuscript discusses the effective interaction potentials between C-S-H particles that are currently used, followed by methods to simulate C-S-H formation. Structural, physical and mechanical properties predicted by the existing simulations are then presented. Finally the manuscript highlights opportunities for future research, which are driving the multi-scale modelling of C-S-H but also of other mesostructured materials.

## 1. Introduction

Production of cement and concrete has a high environmental impact, with concerns linked to consumption of raw materials and global CO₂ emissions [1]. To mitigate these issues, new cements are being researched and tested, featuring large and sometimes complete substitution of traditional Portland cement with supplementary cementitious materials [2]. Predicting the properties of concretes with new chemistries is challenging, especially long-term performance for which there are scarce experimental data. Multi-scale modelling is being increasingly recognised for its potential to complement the experiments and accelerate innovation in cement and concrete science and technology [3].

Multiscale modelling [4] typically starts with ab-initio and atomistic simulations, which leverage first principles to predict thermodynamic, structural, mechanical, and kinetic data for various cement phases at the nanoscale. The effects of these fundamental properties are transferred across length and time scales to eventually inform macroscale engineering models, through constitutive laws that describe the properties of individual phases, interfaces, and mixtures. Atomistic and macroscale models of concrete and its phases are quite established and reviewed elsewhere [5–7]. However, there is clearly a need to connect the nanoscale with the macroscale. Such connection spans multiple scales in

length and time: from the nanometre and nanosecond of atomic processes, to at least millimetres and seconds at the engineering scales. This cover a wide range of mesoscales, and indeed different sub-communities within concrete science interpret the term differently. For example, in concrete modelling, the mesoscale of interest is typically in the $10^{-3}$ – $10^{-1}$ m range, where aggregates, cementitious matrix, and interfacial zones interact [6,7]. Instead, in the cement hydration community, the mesoscale is that of microstructure development: $10^{-7}$ – $10^{-3}$ m [8]. Our manuscript focuses on a mesoscale that is particularly important in the upscaling of properties from atomistic simulations: the nano-to-micro mesoscale, $10^{-9}$ – $10^{-6}$ m: see Fig. 1. This range of scales is implied whenever the word “mesoscale” appears in this manuscript without any further specification.

The nano-to-micro mesoscale is important because it is the first scale above atomistic simulations, so it is the first step in the upscaling effort. At this scale the (nano-) particles are still small enough for thermal fluctuations and thus for entropy to play an important role [9]. The use of statistical mechanics allows to rigorously describe the equilibrium state of mesoscale systems (resulting from the balance between internal energy and entropy) but also to coarse-grain and upscale information from the molecular scale below [10,11]. At larger scales, the motion of large particles (supermicronic) induced by thermal fluctuation is

<footnote><sup>a</sup> Laboratoire de Mécanique et Génie Civil, CNRS, Université de Montpellier, 34090 Montpellier, France</footnote>

<footnote><sup>b</sup> INTERFACES Dpt, ICB UMR 6303 CNRS, Univ. Bourgogne Franche-Comté, Dijon 21000, France</footnote>

<footnote><sup>c</sup> School of Engineering, Cardiff University, The Parade, Cardiff CF24 3AA, United Kingdom</footnote>

<footnote>* PACS: 0000, 11112000 MSC: 0000, 1111</footnote>

<footnote>* Corresponding author.
E-mail address: enrico.masoero@newcastle.ac.uk (E. Masoero).</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Multiscale modelling of C–S–H, from full atomistic simulations to coarse-grained and particle-based simulations, spanning length scales from nano to micro. Figures adapted with permission (a) Ref. [74] from the Royal Society of Chemistry (b) Ref. [72] from the American Chemical Society and (c) Ref. [90].</figcaption></figure>

negligible at room temperature; moreover intergrain dissipation and friction quickly drain the kinetic energy from the system. The nano-to-micro mesoscale is challenging too; at larger scales, one often assumes soft coupling between chemical transformations and mechanical deformation [12–14]. However, at the nano-to-micro mesoscale, these processes contribute rather commensurately to the same free energy landscape, which governs the evolution and properties of the material. In other words, chemistry and mechanics are fully coupled at the mesoscale [15].

This manuscript focuses on mesoscale simulations of calcium-silicate-hydrate (C–S–H) only. The reason is that almost all the existing mesoscale simulations of cementitious systems focus indeed on C–S–H, which is the main binding phase in traditional Portland cement. However, the techniques and results that will be presented here for C–S–H can be extended and applied to other phases too. Indeed, one of the objectives of this review is to stimulate a broader uptake of mesoscale simulations, to encompass other cement phases and processes as well as possibly other materials beyond concrete. The manuscript is structured as follows:

* Section 2 focuses on solid-solution interfaces, starting with the charging process of silica and C–S–H. This is followed by a presentation and discussion of the adsorption of inorganic ions and large organic molecules (e.g. superplasticizers) at the interface between C–S–H and cement pore solution. The section then illustrates the application of the grand canonical reactive ensemble to the liquid and solid speciation of the $SiO_2-CaO-H_2O$ system. This provides the starting physico-chemical properties to model and rationalize the effective interaction potentials, or Potentials of Mean Force (PMF), between C–S–H nanoparticles used in particle-based mesoscale simulations.

* Section 3 builds on the knowledge on the interface between C–S–H and pore solution, showing how these translate into PMF between neighboring C–S–H surfaces and particles. The PMF determined from atomic force microscopy experiments and coarsened grained simulations are presented first, and their limitations are discussed. Results from full atomistic simulations are presented next. The main effective potentials currently in use in particle-based mesoscale simulations are then presented, discussed, and compared.

* Section 4 presents the techniques to simulate the mesoscale formation of C–S–H using interacting nano-particles. These techniques are based on variations of the Monte Carlo approach, including Grand Canonical and Kinetic Monte Carlo. Results from the literature are presented, focusing in particular on the experimental hydration rate curve of cement from calorimetry tests. The section also features an initial introduction on previous models and simulations to describe

the same curve based on geometric rules, without mechanical interactions.

* Section 5 collects literature results from mesoscale simulations of C–S–H and cement paste. The results encompass structural, physical, and mechanical properties. In detail, these include: (i) volume fraction, cluster, connectivity, and local order analyses, (ii) spatial correlations as measured in small angle scattering; (iii) characterization of the pore network and, related to it, vapor sorption isotherms and capillary stress; (iv) formation eigenstresses; (v) elastic moduli; (vi) mechanical strength; (vii) plasticity and fracture; (viii) creep behavior. The presentation is accompanied by brief explanations of the main methods employed and discussion of the relevance and limitations of the current results.

The manuscript ends with a summary of main achievements and conclusions to date, as well as an outlook on some open challenges in the field of coarse grained and mesoscale simulations of C–S–H and other cement hydrates. The manuscript is accompanied by a brief introduction to the statistical mechanics background and tools used in the various mesoscale simulations presently reviewed, see the supplementary material S1. In particular, it is shown how the essential features of the free energy landscape of a complex molecular system can be coarse-grained using Potentials of Mean Force (PMF) between interacting species or nanoparticles. The statistical mechanical description of systems whose evolution is dictated by chemical reactions used in Sections 2 and 4 is introduced too.

## 2. Solid-solution interfaces

Solid-solution interfaces play a central role in many physico-chemical processes. In particular, adsorption and transport of solutes, interparticle interactions, paste rheology and microstructure, and precipitation and dissolution of solids, all fundamentally depend on the physico-chemical properties of solid-solution interfaces. In this section, we review the key features and principles of solid-liquid interfaces and illustrate their relevance to adsorption. The aim is to introduce fundamental concepts that are pivotal to interparticle interactions and precipitation of hydrates described in later sections.

Most minerals present an electric charge on their surface when in contact with an aqueous solution. The surface charge is either permanent/structural, originating from the isomorphic ionic substitution in the mineral structure, like in the case of ettringite [16], or amphoteric (pH dependent), that is due to the ionisation of surface groups, like e.g. C–S–H particles [17]. On top of this, due to electrostatic forces and thermal agitation, the ions in the solution form a cloud at the solid-solution interface. At equilibrium, this ionic cloud neutralizes the mineral surface charge. The electrostatic interactions in the ionic cloud, best known as the electrical double layer (EDL), control the properties of individual mineral particles, such as their precipitation mechanisms and rates, as well as the interactions between different particles. We shall illustrate this in Sections 3 and 4.1. It is therefore essential to understand and accurately model charge formation at mineral surfaces, the focus of this section.

The Gouy-Chapman-Stern model (GCSM) is the most used approach to describe the EDL; it dates back to the beginning of the last century and it is based on experimental observations with a mercury drop electrode in aqueous solutions (see Refs. [18,19] and also the work of Grahame [20]). In this model, the EDL is described as two adjacent layers, a compact (Stern) layer close to the surface and a diffuse layer that extends into the bulk solution. The GCSM is based on the Poisson-Boltzmann equation, a mean field approximation of the primitive model (PM) where ions, considered as point charges within an implicit solvent, interact with each other and with charged surfaces through a mean electrostatic potential, used as an approximation for the PMF. In other words, calling $\psi(x)$ the mean electrostatic potential at a normal distance $x$ from a charged surface, the PMF of an ion $i$ of charge $z_i$ with
the surface is approximated within the GCSM as $w_{is}(x) = z_i e \psi(x)$, and its concentration profile as $c_i(x) = c_i^{bulk} \exp[-\beta w_{is}(x)]$. The ion and potential profile in the diffuse layer are obtained by solving the Poisson-Boltzmann equation, see e.g. [21] The charging mechanism at mineral surfaces is self consistently described by the pK Stern model. Let us illustrate this taking silica as an example, starting with some generalities.

The surface charge of silica, like C–S–H, arises from the deprotonation of the silanol groups according to the reaction:

$$Si-OH \rightleftharpoons Si-O^{-} + H^{+} \quad (1)$$

governed by the intrinsic mass balance constant $K_a$, defined by the activity product:

$$K_a = \frac{a_{Si-O^{-}} - a_{H^{+}}}{a_{Si-OH}} \quad (2)$$

Introducing the chemical potential of the surface sites $\mu$ as well as the definition for pH ($pH = -\log a_H$) and pK, Eq. (2) can be rewritten as:

$$\ln(10)pK_a = \ln(10)pH - \ln\frac{\alpha}{1-\alpha} - \beta(\mu_{Si-O^{-}}^{ex} - \mu_{Si-OH}^{ex}) \quad (3)$$

where $\alpha$ is the fraction of deprotonated silanols. Recognizing that the excess chemical potential term $\beta(\mu_{Si-O^{-}}^{ex} - \mu_{Si-OH}^{ex})$ in Eq. (2), is nothing but the excess free energy change (PMF) to ionize a silanol group, $w_{ss}^{ex}(0)$, one can then write the change in Gibbs free energy for charging up a silanol group as

$$\Delta G = kT \ln \frac{\alpha}{1-\alpha} = kT \ln(10)(pH - pK_a) - w_{ss}^{ex}(0) \quad (4)$$

from which the surface charge density is obtained from the relation

$$\sigma_0 = \alpha z_s e \Gamma_s \quad (5)$$

where $\Gamma_s$ is the total number per unit area of silanol groups and $z_s e$ their charge in the deprotonated state, $z_s e = -1e$ here. So far, no approximations have been used: Eqs. (3), (4) are general and exact. Within the GCSM, the PMF for the surface groups is approximated by the mean electrostatic potential, here the surface potential, i.e. $w_{ss}(0) = -e\psi(0)$. The drawback of this approximation is that the ion-ion correlations between charged species (ions as well as surface groups) are neglected. Nevertheless, the GCSM has been successful in predicting the surface charging process and the colloidal stability of many experimental systems, but limited to weakly charged surfaces in contact with electrolyte solutions of monovalent ions at low concentrations [22]. For highly coupled systems (high concentration, multivalent ions, and high surface charge), where ion-ion correlations are important, the GCSM leads to quantitatively and qualitatively wrong results [23,24]. This is typically the case for cementitious systems and C–S–H in particular, as we shall see in the next section.

## 2.1. Surface charge

Starting with the charging process of minerals, Labbez et al. in their work on silica [25] found that both ion–ion and ion–surface-site correlations are important. This is illustrated in Fig. 2, which compares predictions from Monte Carlo (MC) simulations and GCSM with experimental measurements of the charging process of silica in sodium and calcium salt solutions, varying the pH and salt concentrations. The MC simulations were performed in the grand canonical reactive ensemble (see Ref. [17] and Supplement S1.4) and provide an exact solution of Eqs. (3), (4). The silica-solution interface was modelled in the framework of the PM as an infinite surface featuring explicit sites (silanol groups) next to an electrolyte solution, made of explicit ions embedded in a dielectric continuum, see Ref. [25] for more details. Both the GCSM and the MC simulations considered only electrostatic interactions. Identical parameters and calculation approaches were used

<figure><img src="image_6.png" /><figcaption>Fig. 2. Surface charge density of silica as a function of pH, obtained from potentiometric titration experiments [26], mean field GCSM, and Monte Carlo simulations in the reactive ensemble [25], see also Supplement S1.4.2. For the calculations, the intrinsic dissociation constant and the surface density of the silanol groups were set to $pK_a = 7.7$ and $\Gamma_s = 4.8$ sites/nm², respectively. In the MC simulations, where ions and surface sites are explicit, the minimum separation between an ion and a surface site was set to $d_{is} = 3.5Å$. In the GCSM, where ions are point charges and the surface sites are smeared out, a gap between the surface and the solution (Stern layer) of thickness comparable to $d_{is}$ was introduced to facilitate the comparison with the MC simulations. Figure adapted from Ref. [17]. Copyright 2009 American Chemical Society.</figcaption></figure>

for both the Na-silica and Ca-silica systems. Fig. 2 shows that both the GCSM and the MC simulations reproduce the qualitative trend of the experimental data. The surface charge density $\sigma_0$ is found to increase with pH, as expected from Eq. (1), but it is always smaller than in the ideal case due to strong electrostatic repulsion between neighboring sites. A change in the salt concentration modulates the screening of the electrostatic interactions which, in turn, changes $\sigma_0$. Fig. 2 further shows that the predictions from the GCSM are only accurate at low pH ($< < pK_a$) and salt concentration, where $\sigma_0$ is small. This gets more pronounced when considering divalent $Ca^{2+}$ instead of monovalent $Na^+$ ions. On the other hand, the MC simulations describe very well the charging process of silica also at high pH and for both monovalent and divalent counterions. This shows that ion-ion correlations can be pivotal in the generation of electrostatic charge at mineral surfaces.

The combination of high pH and calcium concentration in the pore solution of a typical Portland cement paste is thus expected to create a strong negative charge on the C–S–H surface. We recall here that the surface charge of C–S–H, like silica, arises from the ionisation of silanol groups. Using the same simulation strategy as above Labbez et al. showed, indeed, that the surface silanol groups of C–S–H are fully ionized in saturated and supersaturated lime solutions. These simulations used the structure of Tobermorite as a model for C–S–H, thus featuring a surface site density of 4.82 / nm², and $pK_a = 9.8$ for the Si-OH ⇌ Si-O⁻ + H⁺ reaction: see Ref. [17,27] for more details. In these conditions, C–S–H particles possess negative charge density as strong as –772 mC/m², see Fig. 3. To date, this is the strongest charge (in absolute value) ever reported for a mineral surface in aqueous solutions. The MC-simulated charge densities, arising from the deprotonation of the silanol groups, are in line with those calculated by Nonat et al. [28] using a classical surface complexation model. Unlike Monte Carlo simulations, however, the classical complexation model necessitates the use of four complexation reactions of the calcium ions with the C–S–H silanol groups in addition to the ionisation reaction in Eq. (1). This is due to the neglect of the ion–ion correlations [27].
<figure><img src="image_7.png" /></figure>

**Fig. 3.** Simulated surface charge density of C–S–H and solution pH as a function of the equilibrium calcium concentration in the bulk solution for the pure CaO-SiO₂-H₂O system. The simulations were performed with Reactive Monte Carlo (RMC) in the grand canonical ensemble in the framework of the primitive model, see Supplement S1.4.2, taking the Tobermorite as a model for C–S–H. The surface density of the silanol groups was set accordingly to 4.82 sites/nm². Their intrinsic dissociation constant was adjusted to pKa = 9.8 based on a combination of independent experimental data including potentiometric titration, Na⁺ adsorption and electrokinetic potential. The minimum separation between an ion and a surface site was set to dᵢs = 2Å, see main text and Ref. [17] for more details. Figure adapted from Ref. [27]. Copyright 2011 with permission from Elsevier.

All these modelling results illustrate the prominent role played by the divalent calcium ions in the understanding of the charging process of C–S–H and, more generally, that played by the ion–ion and site–ion correlations. In the following, we will further illustrate the importance of these correlations in the understanding of electrokinetic behavior and ion adsorption.

## 2.2. Electrokinetic potential

Electrokinetic measurements, like electrophoresis, are simple and rapid methods to gather qualitative information on the adsorbed state and amount of chemical species (e.g. ions and polyelectrolytes) at particle surfaces. They provide the electrostatic potential (electrokinetic potential or zeta potential) a few angströms away from the particle surface, resulting from the bare surface charge modulated by the adsorbed species (in its close vicinity) but also by the ionic strength of the solution. A clear example of the impact of ions is the continuous raise and change in sign, from negative to positive, of the electrokinetic potential of C–S–H particles, when equilibrated with a solution containing an increasing amount of Ca(OH)₂. This experimental behavior is perfectly reproduced by the primitive model of electrolyte solution and titrating surfaces, see Fig. 4. It can thus be explained on the basis of electrostatic interactions and correlations. Upon increasing both pH and calcium concentration, the negative bare charge of C–S–H becomes stronger and is compensated by an increasing amount of Ca²⁺ at the C–S–H/aqueous-solution interface. Simultaneously, the surface/solution calcium ion concentration profile deviates progressively from a classical diffuse layer (Poisson-Boltzmann distribution); see Fig. 5. The calcium ions accumulate in a thin layer close to the C–S–H surface and as much so as the ion-ion (and ion-site) correlations strengthen and, thus, as the charge of C–S–H becomes important. Above a threshold value of ≈ 2 mM of Ca(OH)₂ (pH 11.7) the build-up of the divalent calcium ions in this layer is so important that it more than compensates the bare charge and reverses the sign of the electrokinetic potential. In these conditions, a charged species more than a few angströms away from the surface,

<figure><img src="image_8.png" /></figure>

**Fig. 4.** Comparison between experimental (points) and simulated (lines) electrokinetic potentials (ζ) of C–S–H dispersions in equilibrium with various electrolyte solutions. The experimental and simulation data are from [17]. Ca(OH)₂: pure CaO-SiO₂-H₂O system where pH is increased with Ca(OH)₂; Ca 1 mM: pH raised with NaOH starting from a C–S–H dispersion at C/S 0.7; Ca 20 mM: pH increased with Ca(OH)₂ and Ca²⁺ concentration adjusted to 20 mM with CaCl₂. The simulations were performed with Reactive Monte Carlo in the grand canonical ensemble, see Supplement S1.4.2. Figure adapted from Ref. [17]. Copyright 2006 American Chemical Society.

<figure><img src="image_9.png" /></figure>

**Fig. 5.** Simulated ion profiles at the C–S–H/solution interface for increasing equilibrium bulk concentration of calcium (and pH) in the pure CaO-SiO₂-H₂O system. Note, that above 2.5 mM of Ca²⁺ (C/S 0.8) the surface charge of C–S–H is overcompensated by the calcium ions and the co-ions, here mainly OH⁻, start to adsorb as a second layer.

perceives an *apparent* charge whose sign is opposite to the bare charge of C–S–H. A process that is strongly impacted by this, and that we shall develop below, is the adsorption of anionic species by C–S–H, see e.g. the *OH⁻* profile in Fig. 5.

Multivalent cations such as cationic oligomers and polymers as well as metal complexes, were also found to overcompensate the surface charge of C–S–H. Conversely, addition of monovalent cations, e.g. alkalies (Na₂SO₄, NaOH, KCl, etc…), decrease the zeta potential of C–S–H [17,27,29–31]. The electrokinetic behavior of C–S–H in the presence of cationic oligomers and alkalies were shown to be well predicted by MC simulations at the level of the primitive model, see [17,27,29]. All these results confirm the generalities of the mechanisms described above and, in particular, the prevalence of ion-ion
correlations in explaining the charge reversal of C–S–H, which has significant impact on interactions at the mesoscale.

## 2.3. Ion and polyelectrolyte adsorption

Adsorption of ions on cement hydrates is one of the key factors which controls the ion composition of the pore solution and, thus, the saturation index of the different phases which dissolve and precipitate during cement hydration. It is also a determinant factor for the understanding of ion transport, pollution leakage, durability of the material and of inter-particle interactions. Due to its large specific surface area and ubiquity in hydrated cement paste, a large number of ion adsorption studies has been reported on C–S–H, with particular emphasis on alkali, sulfate and chloride ions.

Classical modelling of ion sorption based on the GCSM and/or geochemical modelling have been unsatisfactory so far, despite using numerous complexation constants. By contrast, the PM for electrolyte solutions and titrating surfaces was shown to provide accurate predictions for the adsorption isotherm of a number of charged solutes, including Na⁺, K⁺, Ca²⁺, SO₄²⁻ and small cationic oligomers, see e.g. [17,27,29].

As an example, Fig. 6 compares experimental and simulated adsorption isotherms of sodium and sulfate ions on C–S–H. The measurements were performed on C–S–H dispersions, with different initial calcium to silicon stoichiometric ratios, and an increasing amount of Na₂SO₄ salt added to the equilibrium solution [32,33]. The simulations used a short-range attractive potential between calcium and sulfate ions, w<sub>Ca, SO₄</sub>, in addition to the Coulomb pair potential, to accurately describe the Ca²⁺–SO₄²⁻ ion pair formation.<fnref n="1" /> w<sub>Ca, SO₄</sub> was adjusted to reproduce experimental conductivities and chemical activities of pure CaSO₄ aqueous solutions. The following functional was used:

$$ \beta w_{Ca,SO_4}(r) = A_{CaSO_4} \exp(-r/\tau) u_{disp}(r) \quad (6) $$

where A<sub>CaSO₄</sub> = 2.15kT and τ = 2Å are parameters that regulate the magnitude and range of the ion-pair potential and:

$$ \beta u_{disp}(r) = 2\pi \left[ \frac{2}{45(\sigma/r)^9} - \frac{1}{3(\sigma/r)^3} \right] \quad (7) $$

where σ = 4Å is the ion diameter. w(r) may be regarded as an effective potential of mean force, including dispersion forces and solvent-induced effects.

As shown in Fig. 6, the simulations perfectly capture the simultaneous increase in the sulfate adsorption and drop of the sodium adsorption upon increasing the Ca/Si stoichiometric ratio, or equivalently, upon raising the equilibrium calcium concentration and pH of the bulk solution. The adsorption of sulfate ions on the negatively charged surface of C–S–H is found to be mediated by the Ca²⁺ ions. It is a consequence of the build up of the overcompensating Ca²⁺ layer covering the C–S–H surface, see e.g. Figs. 4 and 5, promoted by the strong (negative) surface charge of C–S–H. The sulfate adsorption is further found to be strengthened by the Ca–SO₄ ion pair formation.

On the other hand, the behavior of sodium adsorption is the result of the Ca²⁺/Na⁺ competition with the negatively charged surface of C–S–H, which, for electrostatic reasons, largely favors the divalent calcium ion. In a similar way, Viallis et al. [34] observed an increased sodium adsorption when decreasing the C/S ratio. More interestingly, using NaCl instead of NaSO₄ to raise the bulk sodium concentration up to 1 M, the measured adsorption isotherm was not observed to saturate, but instead to increase linearly in logΓ<sub>Na</sub> - logc<sub>bulk</sub> plot (Freundlich type). Simulations at the level of PM were shown to perfectly predict this behavior which again is explained by the electrostatic competition

between Ca²⁺ and Na⁺ ions with the C–S–H surface [27].

The identified mechanism (ion pairing) for the adsorption of SO₄²⁻ on the negatively charged C–S–H surface provides the basis for the understanding of the anion affinity with C–S–H but also the link between the ionic interactions in bulk solutions and at solid/liquid interfaces. In practice, this means that the affinity of anions with C–S–H can be directly inferred from measurements or databases of thermodynamic properties of CaₓAnion_y salt solutions, e.g. solubility, complexation constant or ion activity. As an example, based on the solubility of calcium salts (the magnitude of the Ca–anion pair potential being inversely proportional to the calcium salt solubility) the following order of anion affinity with C–S–H can be deduced: CO₃²⁻ >> SO₄²⁻ and OH⁻ > IO₃⁻ >> I⁻ > Cl⁻ > NO₃⁻ > Br⁻ > MnO₄⁻, in good agreement with experimental observations. The implication of this finding is further illustrated below for anionic polymer/superplasticizer.

Organic admixtures are commonly used in the formulation of concrete and other cement-based construction materials for different purposes, from adjusting reactivity (e.g. retarders, accelerators) to controlling workability and mechanical performance (e.g. superplasticizers). Despite their common usage and the important empirical knowledge acquired, the understanding of the mechanisms involved is still in its infancy, in particular when it deals with reactivity (see Ref. [36] for a recent review on molecular modelling of admixtures). During the last decade, however, Turesson and co-workers [35,37] brought an important contribution to the rationalization of the cement-hydrate-superplasticizer interactions, starting with the adsorption mechanism. The primitive model for polyelectrolyte solutions was used to simulate bulk solution mixtures of polyelectrolyte and ions in contact with a planar charged surface, as depicted in Fig. 7a. Building on the knowledge acquired on simple ions, Turesson et al. generalized the concept of ion-pairing to bulky superplasticizers and polyanions. They showed that the strength of ion-pairing between calcium ions and polyanion functions is a key parameter to characterize the affinity of superplasticizers with C–S–H (or hydroxylated allite) and by extension their sulfate tolerance. This is illustrated in Fig. 7b which shows the increased adsorption of a comb copolymer model on a negatively charged surface (C–S–H or hydroxylated allite) upon deepening the attraction of the Ca²⁺–organic-function pair potential.

In line with experimental observations, the model predicts a lower adsorption for comb polyanions having longer neutral side chains (PEO), higher PEO grafting density and smaller polymerization degree for the charged backbone. This behavior is best explained by the balance between attractive surface energy (ion pairing, overcharging) and the repulsive contribution of the internal entropy of side chains. Similarly to simple coions, polyanion adsorption is also found to increase with the magnitude of negative charge density of C–S–H.

From a practical point of view, Turesson et al. proposed a simple experimental method of evaluating ion pairing and thus to explain and predict the marked difference observed in adsorption isotherms of superplasticizers with different geometries and functionalities. The ion pairing was quantified from measurements of potentiometric titration and calcium binding isotherms of superplasticizers. Although a proper quantification of ion-pairing also requires the use of simulations, it should be noted that a qualitative evaluation of the ion pairing strength can easily be obtained from the slope of the calcium binding isotherm at low Ca²⁺ concentrations. Such strategy was successfully used by Nalet et al. [38] to rationalize the adsorption of various negatively charged organic molecules onto C–S–H surface.

The adsorption of polyanions on C–S–H was also predicted to deviate from the common Langmuir isotherm characterized by a plateau when saturation is reached. The adsorption isotherm was, instead, found to first rise, reach a maximum, and then decrease upon further increasing the polymer content, as exemplified in Fig. 8. This non-Langmuir type of adsorption is in fact quite general for ion-mediated adsorption of polyelectrolytes and simple molecules/ions. This has been observed for various systems [39–45]. Turesson et al. further showed that a typical

<footnote n="1">¹ The latter is, in part, responsible for the low solubility of CaSO₄ (18 mM in pure water and at room temperature)</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** (a) Simulation snapshot of the calcium mediated adsorption of comb copolyanions onto a liked-charged mineral surface (C–S–H or hydroxylated allite). The polyelectrolyte salt solution was described at the level of the primitive model while the C–S–H (or hydroxylated allite) was modelled by a uniformly charged planar surface, see [35] for more details. (b) Simulated monomer density profile (full line: charged backbone monomers, dashed line: neutral side chain monomers) at the C–S–H/solution interface for a model superplasticizer (comb copolyanion) varying the specific calcium-charged-monomer pair potential, $\beta w = C/r^6$. $\beta w_{min}$ indicates the maximum strength of this non-electrostatic potential at the minimum separation (4 Å) between charged backbone monomers and calcium ions. The bulk concentration of charged monomers and calcium was ≈30 mM and 20 mM, respectively. The charge of C–S–H was set to –4e/nm². The polymerization degree of the comb polyanion backbone was set to 13 and that of the PEO side chains to 8 with a grafting density of 23%. Reprinted from Ref. [35]. Copyright 2011 American Chemical Society.</figcaption></figure>

Langmuir adsorption isotherm is, however, recovered when the solution is buffered with an excess of lime. These conditions, which are typical of cement pastes, are equivalent to having a very large reservoir of calcium with constant chemical potential, see Fig. 8. The experimental demonstration of this behavior was recently obtained by Bouzouaid et al. [46], studying the adsorption of gluconate on C–S–H.

## 2.4. Speciation of C–S–H

Modelling the speciation of solid and liquid phases allows predicting phase assemblage, which is key for the long-term durability of cementitious materials and for the development of new cement blends with lower CO₂ impact. Despite lot of efforts, there are still remaining issues due to the complexity of the cementitious systems, but also due to the limitations of thermodynamic models (geochemical models) based on

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Simulated adsorption isotherm of comb-polyanion on a negatively charged surface (C–S–H or allite) under conditions of constant calcium content and constant calcium chemical potential (solution buffered with portlandite). $\sigma = -4 \text{ e.m}^2$ and $\beta w_{min} = 0$. The same comb polyanion model and simulations as in Fig. 7 was used. Reprinted from Ref. [35]. Copyright 2011 American Chemical Society.</figcaption></figure>

solid-solution or on surface complexation approaches.

The system complexity can be illustrated by the main hydrate phase, namely C–S–H. Its composition, expressed in terms of Ca/Si and H₂O/Si stoechiometric ratio, varies with the composition (ion concentrations) of the pore solution with which it is in equilibrium. Due to its charge, C–S–H particles can physically adsorb both anions and cations, see Section 2.3. Furthermore, C–S–H can incorporate aluminate ions in its structure through a Si/Al substitution reaction. Those variations are accompanied by changes in the silicate chain length of C–S–H and in the ion concentrations of the pore solution (silicate from few mM to a few µM, 1 mM < Ca < 30 mM, 10 < pH < 13).

The limitations of the classical geochemical models have in part been discussed in the previous Section 2. They are principally due to the need of introducing a large number of surface complexation constants (for the surface complexation models) or end-members (for the solid solution models) to compensate for the many microscopic physical interactions which are not explicitly accounted for at this level of approximation. For example, the overcompensation of the C–S–H surface charge by calcium ions is modelled by a surface complexation constant or, equivalently, by a new end-member, whereas it appears naturally as the result of ionic correlations when the latter are taken into account, c.f. Section 2.2. Furthermore, the activity calculation of electrolyte solutions is generally
performed with extended Debye-Huckel theories whose range of validity in ionic strength is relatively small, typically up to 100 mM.

Alternatively, the speciation calculation can be performed at the level of the PM with RMC simulations, see Supplement S1.4.2. As an illustration, Fig. 9 compares experimental and simulated variation of the Ca/Si ratio and silicate mean chain length of C-S-H in the pure CaO-SiO₂-H₂O system. In the simulations the same model as in Section 2 was used to which was added the reaction of depolymerization of the silicate chains at the C-S-H surface and the two main speciation reactions of the silicates in solution (namely the first and second ionisation of the Si(OH)₄ species). In other words, beside the solubility constant of C-S-H, the speciation of C-S-H is completely described by two surface reactions: the ionisation of the surface silanols and the depolymerization of the silicate chains. For comparison, the classical surface complexation model for C-S-H of Nonat [28] requires the use of five surface complexation reactions. The RMC was also employed by Churakov et al. [51] to simulate with an equal success the aluminum incorporation in C-S-H.

### 3. Interactions between C-S-H particles

The previous sections have provided the fundamental concepts to understand the physical interactions between C-S-H surfaces and cement pore solution. This section applies those concepts to the quantification and modelling of effective interaction between coarse-grained units of C-S-H, viz. nanoparticles. The main concept to develop such interactions is the Potential of Mean Force (PMF) between C-S-H surfaces and particles. Such PMF is pivotal in mesoscale simulations as it allows to reduce a complex molecular system to a one-component system featuring only coarse-grained particles in a rigorous manner; see the Supplemental Information S1.3 for more details on the derivation and description of PMF. In the field of cement modelling, the development of PMFs from molecular simulations has progressed in parallel with the development of first coarse-grained nanoparticle models. As a result, most of the existing particle models of C-S-H use simplified interaction potentials that have been built using the emerging information from atomistic simulations, but without being constructed around a complete PMF. Ref. [56] was indeed the first work using a PMF derived from atomistic simulations as an effective interaction potential in mesoscale, particle-based simulations. Applications and results of particle-based simulations of C-S-H, using both simplified interaction potentials and those directly derived from a full PMF, are reviewed in Sections 4 and 5.

The section is structured as follows. Experimental data from atomic force microscopy (AFM) on the interaction between C-S-H surfaces are presented first, and compared to simulation results that consider slit

pores between infinite flat surfaces of C-S-H. Meaningful quantities, such as adhesive force and interfacial tension, are computed and discussed. The limitations of the AFM experiments in quantifying pairwise particle-particle interactions are discussed next. Results from molecular simulations are then presented and various effective interaction potentials currently in use are reviewed and discussed.

#### 3.1. Interactions from AFM experiments and slit-pore simulations

Interactions between C-S-H particles in various solutions representative of cement pore solutions were measured experimentally by atomic force microscopy (AFM). The force measurement was performed between an AFM silica tip covered with C-S-H particles and a large and atomically flat surface of C-S-H, grown by Ostwald ripening on the surface of a calcite monocrystal. A medium-range attractive force was observed, see negative values in Fig. 10-a. The adhesion force (pull-off force) was further observed to increase in magnitude with the pH and calcium content of the solution [52], see Fig. 10-b, that is with the charge (in absolute value) of the C-S-H particles, c.f. Section 2. These results show that C-S-H surface forces are essential components of cement cohesion, highlighting also the fundamental role played by the physico-chemical properties of the C-S-H/solution interface in controlling those forces.

Before discussing the physical origin of the cohesive interaction between C-S-H particles in solutions, let us address the relationship between the measured force, $F(h)$, and the interaction free energy per unit area, $W(h)$, between two infinite flat surfaces, as well as their link to the surface tension, $\gamma_s$. $W(h)$ is a PMF, in that it measures the change in free energy when two surfaces approach each other from infinity down to a distance $h$, see Supplement S1.3. $W$ can be estimated from the measured force between the C-S-H single crystal and the AFM tip covered with C-S-H, whose effective curvature is $1/R_{eff}$, applying the Dejarguin approximation:

$$W(h) = \frac{F(h)}{2\pi R_{eff}} \qquad (8)$$

The relation between the surface tension (C-S-H/solution) and the adhesion free energy $W_{adh}$ then follows from the contact theorem, $\gamma_{sL} = -1/2W_{adh}$ which when applied to Eq. (8) provides an estimation of the surface tension of C-S-H in aqueous solution,

$$\gamma_{sL} \approx \frac{F_{adh}}{4\pi R_{eff}} \qquad (9)$$

where $F_{adh}$ is the measured adhesion force. Fig. 10-b gives the C-S-H/solution surface tension estimated from AFM force measurements with

<figure><img src="image_12.png" /></figure>

grand canonical ensemble, see Supplement S1.4.

**Fig. 9.** (a) Comparison of the simulated (lines) and experimental (points) stoichiometric calcium to silicon ratio (C/S) of C-S-H as a function of the equilibrium calcium concentration in the pure CaO-SiO₂-H₂O system. Experimental points from Haas et al. [28] as obtained from the difference method which provides a measurement encompassing both structural and adsorbed calcium. The simulations distinguish the structural calcium in the octahedral sites of the C-S-H layers from the physically adsorbed calcium ions at the surfaces of the layers. (b) Simulated (lines) and experimental (points) mean length of the silicate chains of C-S-H. The experimental data are from [47–50]. The simulated mean chain length is calculated for both an infinite C-S-H crystal (black line) and a finite C-S-H crystallite with a maximum mean chain length of 40 silica units. In the simulations, the speciation of the silicate species in solution is accounted for and an equilibrium constant of pK = 6.9 is assumed for the depolymerization of the silicate chains. The simulations were performed with Reactive Monte Carlo in the
<figure><img src="image_13.png" /><figcaption>Fig. 10. (a) Measured attractive force between an AFM tip covered with C–S–H particles and an atomically flat C–S–H single crystal in equilibrium with aqueous solutions at various amount of Ca(OH)₂ (b) Adhesion free energy (W<sub>adh</sub>) between C–S–H flat surfaces in aqueous solutions of Ca(OH)₂ at various concentrations. W<sub>adh</sub> was determined from MC simulations in the framework of primitive model, green triangles, as well as from the adhesion force (F<sub>adh</sub>) measured by AFM, normalized by the effective radius of the AFM tip using Eq. (8), black squares. Each point is the average of 100 measurements, see [52] for more details. The right hand scale gives the C–S–H/solution surface tension γ<sub>el</sub> calculated using Eq. (9). The red circle provides the surface tension as obtained from homogeneous nucleation of C–S–H, see [53]. Figure (a) reprinted from Refs. [52,54]. Copyright 2004 American Chemical Society. Figure (b) was produced with Monte Carlo simulations in the grand canonical simulations as in Ref. [55] including the charge formation of C–S–H following Refs. [17,27], that is using the reactive method, see Supplement S1.4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$R_{eff} = 20 \pm 10$ nm. $R_{eff}$ was evaluated independently from AFM force measurements in dry air using various substrates of reference, see [52] for more details. The C–S–H/solution surface tension is observed to increase from 1 to 13 mJ/m² with the Ca(OH)₂ concentration of the equilibrium solution; higher Ca(OH)₂ concentrations also imply higher absolute values of the surface charge density of C–S–H, c.f. Fig. 3. γ<sub>el</sub> of C–S–H was also estimated from measurements of nucleation rates and applying the classical nucleation theory [53]; this gave a similar value of 12 mJ/m², see Fig. 10-b. It should be noted that the obtained surface tension in water γ<sub>el</sub> for C–S–H is significantly lower than the surface tension measured, also by AFM, in dry air, γ₀= 185 mJ/m² [52], and than that determined from atomistic simulations, see e.g. [56] and Section 3.4 thereafter. Similar values of surface tension were obtained for mica in water (γ<sub>el</sub> = 3–11 mJ.m⁻²) and in dry air (γ₀ = 130–160 mJ.m⁻²) from force measurements, see e.g. [57]. The great difference in values of the C–S–H surface tension as obtained from AFM measurements in solutions and atomistic simulations may indicate that the free energy of interaction between C–S–H particles (or mica surfaces) in water is composed of two free energy minima separated by a sharp free energy barrier induced by hydration forces, i.e. the work needed to partially dehydrate the adsorbed calcium ions. In other words, this difference may indicate the presence of two regimes, i.e. the colloidal and the crystalline regimes [58]. We will further develop this point in the next sections.

$W(h)$ (and $\gamma_{el}$) was also calculated with MC simulations in the framework of the primitive model and compared with experiments [54,55,59] using a slit geometry formed by two infinite and parallel C–S–H surfaces separated by an aqueous solution in equilibrium with a bulk. With this geometry, $W(h)$ is obtained by integrating the osmotic pressure of the interstitial solution, $p_{osm}(h)$, with respect to $h$

$$W(h) = - \int_{\infty}^{h} (p_{osm}(h') - p_{osm}^{bulk}) dh' \quad (10)$$

where $p_{osm}^{bulk}$ is the osmotic pressure of the bulk solution. In the framework of the primitive model, we can write $p_{osm}(h)$ as the sum of three contributions,

$$p_{osm}(h) = k_B T \sum_i \{ \rho_i(h/2) + p_i^{hc} + p_i^{corr} \} \quad (11)$$

where $\rho_i$ is the ideal contribution to the osmotic pressure (repulsive), $p_i^{hc}$

is the contribution from the hard core (repulsive), and $p_{osm}^{corr}$ accounts for the electrostatic correlations. The latter, generically called ion-ion correlations, arises from the ion interactions induced by the fluctuation of ion densities over the mid-plane which, when strongly correlated, generates a net attraction, $p_{osm}(h) - p_{osm}^{bulk} < 0$ [60–62]. The ion correlations grow with the magnitude of the surface charge, with the counterion valence, and with decreasing dielectric permittivity of the solvent. In some cases, low salt concentration and monovalent counterions or small surface charge densities, $p_{osm}^{hc} (>0)$ and $p_{osm}^{corr} (<0)$ compensate each other and the mean field Poisson-Boltzmann equation is recovered. That is, the osmotic pressure of the slit solution reduces to $p_{osm}(h) = k_B T \sum_i \rho_i(h/2)$. In this case, the osmotic pressure is always found to be repulsive [63]. Conversely, ion-ion correlations prevail in the presence of multivalent ions and relatively high surface charge densities; this generates a medium range attraction which increases with the magnitude of the surface charge density, as observed in the experiments. The predicted adhesion forces of C–S–H modelled in the framework of the primitive model are indeed found to compare well with those measured by AFM as illustrated in Fig. 10b.

Goyal et al. developed a primitive model for C–S–H based on the slit geometry with explicit water molecules interacting with SPC/E potential [64] in salt free conditions. This semi-atomistic model highlighted the importance of water-ions interlocking in reducing the dielectric screening and enhancing the ionic correlations in the confinement of the interfoliar space of charged C–S–H sheets. This provided net cohesive forces in very good agreement with molecular simulations, but much stronger than measured by AFM in model experimental systems.

Within the framework of the primitive model, Turesson et al. were able to study the interaction free energy between C–S–H particles in various anionic co-polymer solutions, including superplasticizers [37]. In the presence of superplasticizers, a long range repulsion was found as due to overlapping brushes of neutral side chains supplemented with a short range attractive ion-ion correlation interaction. Unlike copolymers, block copolymers were found to intercalate in between charged C–S–H nanoplatelets, forming hybrid tactoids, as due to a bridging mechanism. In these conditions, the interaction free energy shows two free energy minima separated by a strong free energy barrier. These simulations successfully helped in the development of hybrid C–S–H mesocrystals [65], with highly aligned C–S–H nanoplatelets interspaced with polymers. This hybrid C–S–H material has been shown to have a bending strength similar to nacre, outperforming all C–S–H-
based materials known to date.

## 3.2. Limitations of AFM in measuring particle-particle interactions

AFM experiments helped clarifying the nature and the evolution with the solution composition of the interaction potential between C–S–H particles. However, they may not provide a quantitative measurement of the PMF between pairs of interacting, anisotropic C–S–H particles within a concentrated cement paste.

A first issue is that the pair PMF between charged particles in a dispersion (or paste) depends on the particle density. In the AFM setup, the interactions are measured between two single C–S–H surfaces immersed in a large solution bath. In a concentrated paste, instead, each particle is surrounded by neighbours, a situation which can greatly affect the pair particle PMF due to multi-body interactions [66,67]. In this case, the solution surrounding each of the particles (solvent and ions) is not solely under the influence of electric fields generated by themselves but also by their near neighbours. Therefore, the mean ionic concentration of a concentrated aqueous particle paste is very different from that of the reservoir solution with which it is in equilibrium, and strongly depends on the particle density. In the strong coupling regime (high surface charge density and multivalent counterions), literature results [68,69] showed, indeed, that the cohesion forces from ion-ion correlations increase with the density of particles. This strengthening of the cohesion force with the particle density was however found to be altered with increasing the equilibrium calcium salt concentration (or calcium activity) [70,71]. Indeed, in the limit of very large salt concentrations in the bulk solution, the mean ion concentrations of the paste pore solution becomes equal to that of the bulk.

A second issue is that the pair PMF between particles depends on the size and shape of the particles [73]. As an example, Fig. 11 shows the simulated PMF between two charged C–S–H platelets of various sizes, immersed in a 10 mM calcium salt solution. The simulations were run in the framework of the primitive model [72]. Attraction from ion-ion correlations is only observed for sufficiently large platelets, and it grows with the size of the platelets. The minimum in the PMF scales approximately with the surface area. This becomes more precise for large platelets. A long range repulsion is also observed, which increases with the platelet size. The repulsion is purely entropic in origin and is due to the restricted rotation of the platelets at close contact. Two freely rotating platelets with their centers of mass at some fixed separation will possess a rotational entropy $S_{rot}$. As the platelets approach, this entropy decreases and leads a repulsive interaction free energy ($A = U - TS$). For platelets with very large radius $R_p$, this entropic contribution becomes negligible. Indeed, it can be shown that

$$S_{rot}(r/R_p) \approx S_{rot}(r) - S_{rot}(2R_p) \quad \text{for } r < 2R_p \quad (12)$$

Therefore, increasing the radius of the platelet decreases the rotational entropy by a constant amount for all center of mass separations $r$, whereas the attraction grows as the surface area of the platelet, i.e. $\sim R_p^2$. Finally, as already addressed in the previous section, the measured AFM force may not be representative of the cohesion energy of a C–S–H crystal [58], that is the interaction free energy between inner layers of C–S–H at very small separations, but rather of the colloidal behavior of C–S–H particles at medium-range separations.

## 3.3. Interaction potentials from molecular simulations of C–S–H particles

The previous sections presented experimental and theoretical results on interactions between C–S–H surfaces from AFM and MC simulations in the framework of the primitive model. The contribution and limitations of these measurements have been discussed, also in relation to simulation results for slit-pore geometries and for interacting platelets. The next step in model accuracy is to perform atomistic simulations of C–S–H particles with full atomistic details interacting within an aqueous solution also described at the molecular level.

Bonnaud et al. [74] applied this method first, considering prismatic C–S–H nano-blocks of size ~5 nm at 10% relative humidity, causing a few layers of water to be adsorbed on their surfaces and filling the space between interacting blocks. The resulting potentials, computed in various relative orientations of the blocks, returned a set of curves of grand potential energy $\Omega$ vs. interparticle distance. Curve-fits using orientation-dependent functions were proposed too: see green curves in Fig. 12. However these results and fitted potentials have not been used yet in particle-based simulations at the larger mesoscale, and their computed values of interaction forces and energy may be affected by capillary pressure due to the low relative humidity imposed.

The first adoption of a PMF from molecular simulations into mesoscale particle-based simulations is more recent (Masoumi et al. [56]) and will be presented in the next section. This PMF was obtained using a free energy perturbation approach from a previous article [75]. The same method was also used to compute tangential components of the interactions [76], viz. shear forces due to sliding between C–S–H particles and nanolayers. The importance of shear forces in determining the strength of C–S–H at the molecular scale has been appreciated in various works, e.g. [77–79], however to date a rigorous inclusion of tangential interaction components in mesoscale models is still missing. Atomistic simulations have also investigated how the Ca/Si ratio of C–S–H [80] and how temperature and pressure [81], all impact the PMF between C–S–H nanoparticles; also these effects, however, are still to be considered in the effective interaction potentials for mesoscale simulations.

## 3.4. Interaction potentials used in mesoscale simulations

The first mesoscale simulations of C–S–H using interacting particles

<figure><img src="image_14.png" /><figcaption>Fig. 11. (a) Schematic picture of two C–S–H platelets in a salt solution and confined to a cylinder (not shown). The center of mass of the platelets are constrained to the cylinder axis. The sites of the platelets are shown as red spheres while the divalent counterions as blue spheres. (b) The potential of mean force between two freely rotating C–S–H platelets with $\sigma = -1.2 \text{ e/nm}^2$ in a 10 mM calcium salt solution at a volume fraction of 0.023. The number of sites as indicated in the legend from 50 to 91, corresponds to 5 to 11 nm in diameter. Figures reproduced from Ref. [72], with permission from the American Chemical Society. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Particle-particle interaction potentials from full atomistic simulations. (a) Typical configuration of interacting C–S–H particles, covered by adsorbed water, at a surface-surface distance $d_{ij}$; (b,c) PMF, viz. the grand potential of the system Ω, as a function of $d_{ij}$, for different relative orientations between the particles. The green curves are analytical fits proposed in Bonnaud et al. [74], from where all the images in this figure have been adapted with permission from the Royal Society of Chemistry. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

date back to the early 2010s, before PMFs were computed from full atomistic simulations (Bonnaud et al.'s work from 2016 [74], discussed in the previous section). The first effective interaction potentials between C–S–H particles were thus developed using both logical and empirical considerations, inspired and supported by results from modelling and experiments.

The first effective potential has been proposed by Masoero et al. [82] in 2012. The functional form of the potential was chosen a priori to feature a strong repulsion when two particles are pushed together, an energy minimum at short distances to ensure cohesion, and a fading interaction at longer distances. A pairwise, spherical, Lennard-Jones potential was chosen, but with two modifications: (1) the exponents were not fixed to 12–6, thus leading in fact to a Mie potential [83]; (2) the energy scale was set to depend on the volume of the interacting particles, $\epsilon \sim \sigma_{ij}^3$, to account for particle size polydispersity. $\sigma_{ij} = \frac{1}{2}(\sigma_i + \sigma_j)$ is a measure of the average size of two interacting particles i, j. This led to:

$$U_{ij} = 4\epsilon(\sigma_{ij}^3) \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{2\alpha} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{\alpha} \right] \quad (13)$$

This potential features three parameters: $\sigma$, $\alpha$ and $\epsilon$. The first one, $\sigma$, can be set for each particle i based on its diameter, which is identified as the equilibrium distance $r_0$, i where U is minimum for two particles with same size; specifically, $\sigma_i = r_{0,i}/\sqrt{2}$. The other two parameters, $\alpha$ and $\epsilon$ (with its dependence on $\sigma_{ij}^3$), were calibrated to capture the elastic modulus $E_{eq}$ and the strain at failure $\epsilon_u$ of C–S–H from molecular simulations [82]. Typical potential curves are shown in Fig. 13.a. A detailed explanation of the calibration is provided in Ref. [84].

In 2014, Ioannidou et al. proposed attracto-repulsive spherical potentials [85]. The potentials feature a modified Lennard-Jones term akin to the $\alpha$ exponent of Eq. (13) but with energy well depth $\epsilon_0$ constant and independent of the particle size $\sigma$. The Lennard-Jones term ensured

short-range attraction, but the potential also included a novel feature: a medium-range repulsion described by a Yukawa term. All this led to:

$$U_{ij} = 4\epsilon_0 \left[ \left( \frac{\sigma}{r_{ij}} \right)^{2\alpha} - \left( \frac{\sigma}{r_{ij}} \right)^{\alpha} \right] + A \frac{e^{-\kappa r_{ij}}}{r_{ij}} \quad (14)$$

where $\kappa$ is the inverse screening length and A the Yukawa energy constant. The screening length was set to $\kappa^{-1} = 0.5\sigma$ and different ratios of the energy constants $\epsilon_0$ and A were considered in order to qualitative explore the forces measured from AFM in different Ca(OH)$_2$ solutions, and here previously shown in Fig. 10.a [87]. Fig. 13.b. shows two curves for a high and low Ca(OH)$_2$ concentrations. The energy constants were chosen so that the attractive energy wells has the same depth. In ref. [88], the phase diagrams of such interactions has been computed. Goyal et al. used such potentials for C–S–H gelation at early (low Ca(OH)$_2$ concentration) and later (low Ca(OH)$_2$) stages of hydration [89]. The Yukawa term reduces during hydration and disappears ($A = 0$ at Eq. (14)) at the setting of C–S–H [90] as the cohesion dominates due to elevated surface charges [55].

The potentials in Eqs. (13) and (14) became the bases for most interaction potentials in subsequent particle-based simulations of C–S–H. In some cases the potentials and parametrizations have been used directly and individually, e.g. in Refs. [90–92]. In other cases they were combined and slightly modified, such as in Refs. [93–95].

The potentials in Eqs. (13) and (14) do not capture the orientation-dependence of the interactions, as shown by the full atomistic simulations in Fig. 12. An orientation-dependent potential has been proposed and used in 2015 by Yu and Lau [96], who considered interactions between disk-shaped C–S–H particles. The chosen functional form for the interaction was Gay-Berne [97], calibrated using experimental data as well as adhesion energy between C–S–H interfaces from atomistic simulations. Nevertheless, potentials of Eq. (14) at low Ca(OH)$_2$ concentration have shown the growth of anisotropic elongated clusters akin to

<figure><img src="image_16.png" /><figcaption>Fig. 13. Particle-particle interaction potentials first proposed in (a) Masoero et al. 2012 [82] (Eq. (13)) and (b) Ioannidou et al. 2014 [85]. In (a), the two curves refer to different combinations of particle sizes. In (b), the two curves refer to same particle size but different relative strength of attraction and repulsion. This was inspired by AFM experiments of low and high Ca(OH)$_2$ contents (see Fig. 10). Figures reproduced from Refs. [85,86], with permission from the Royal Society of Chemistry.</figcaption></figure>
Bernal spiral [85,89].

In 2017, Shvab et al. [98] adopted the Lennard-Jones potential in Eq. (13) to simulate heterogeneous precipitation of C–S–H. For this study, they found it necessary to calibrate the potential in such a way that the cohesion energy $\epsilon_0$ be consistent with the interfacial energy of C–S–H in water, taken as 87.6 mJ/m² from Ref. [99]. This calibration of the potential, however, changed the values of elastic modulus compared to the previously discussed calibration, where the Young modulus was determining the value of $\epsilon_0$, not the interfacial energy. To enable a calibration based on both interfacial energy and mechanical properties, Coopamootoo et al. [100] proposed a 4-parameter (instead of 3 parameters in Eq. (13)) shifted harmonic potential. This potential was proposed for tricalcium silicate and later applied also to C–S–H particles [101]:

$$U_{ij} = \frac{1}{2} k_{ij} (r_{ij} - r_{0,ij})^2 - \epsilon_{ij} \quad \text{when } r_{ij} < r_{u,ij} \quad (15)$$

$$U_{ij} = 0 \quad \text{when } r_{ij} \ge r_{u,ij}$$

The four parameters in this potential are the stiffness $k_{ij}$, the equilibrium distance between interacting particles $r_{0,ij}$, the cohesive energy scale $\epsilon_{ij}$, and the interaction cutoff distance $r_{u,ij}$. In all the parameters, the subscript $ij$ indicates that they can account for interacting particles of different sizes (see Ref. [101] for more details). For simplicity, Eq. (15) considers a sharp interaction cutoff, but extension to a smooth cutoff is straightforward (e.g. by multiplying $U_{ij}$ times a logistic function). Section 3.5 below will discuss how the parameters in Eq. (15) can be linked to interfacial energy and mechanical properties.

A last interaction potential has been proposed and used by Masoumi et al., in 2020 [56]. The potential is an orientation dependent Gay-Berne fitted to PMFs from full atomistic simulations. These latter were performed using an approach similar to Bonnaud et al. [74] and considering various relative orientations of the C–S–H particles. However, unlike Bonnaud et al., Masoumi et al. considered thin particles, akin to platelets, and focused their subsequent mesoscale simulations on the agglomeration and stacking of such platelets. This focus is somewhat similar to the aforementioned work of Yu and Lau, who used disk-shaped particles [96]).

### 3.5. Comparison between effective potentials

Comparing the potentials presented in the previous section is not straightforward, as they were developed for different purposes and include unique parameters. There are however some general features that are common to all or most the potentials above. This section will present and discuss these features and use them to attempt a qualitative and somewhat quantitative comparison between the potentials. The purpose is not to extract precise quantities, but rather to provide the reader with an intuition of how different potentials and parametrization are likely to impact the resulting properties from mesoscale simulations featuring multiple particles. The detailed results of actual mesoscale simulations will be then presented and discussed in the next Sections 4 and 5.

Fig. 14 depicts a curve whose features are, in full or in part, displayed by most of the interaction potentials in the previous section. The qualitative resemblance is complete with the well-shoulder Lennard-Jones-Yukawa potential in Fig. 13.b, as well as with the Gay-Berne potential fitting PMFs from atomistic simulations for parallel particles in Fig. 12.b. The curve in Fig. 14 can also be made to resemble the Lennard-Jones potential in Fig. 13.a and the Gay-Berne potential for perpendicular particles in Fig. 13.c, as long as the medium-range shoulder is removed. Consequently, one can also consider the harmonic potential in Eq. (15) as a local approximation of the curve in Fig. 14 around its energy minimum.

The suggested approach to compare the potentials in the previous section is to extract from each of them the descriptor of the curve in Fig. 14 that are indicated on the curve itself. Namely:

<figure><img src="image_17.png" /><figcaption>Fig. 14. Generic pair interaction potential between nano-units. The highlighted quantities are used in the text to compare different potentials from the literature.</figcaption></figure>

* The equilibrium distance $r_0$, which is linked to the diameter of the interacting particles;

* The curvature of the potential $\rho^{-1}$ near its minimum, which controls the linear elastic mechanical response of the C–S–H if most particles in a mesoscale arrangement sit relatively close to their equilibrium positions. Following Ref. [84], the elastic modulus $E_{eq}$ of the material at zero porosity and $\rho$ may be linked as: $E_{eq} = \frac{4}{\pi r_0 \rho}$;

* The change in interaction energy $\epsilon$ when particles are separated from $r_{ij} = r_0$ to $r_{ij} = \infty$. This can be linked to the interfacial energy $\gamma$ between particles and surrounding medium, as previously discussed in Section 3.1. To this end, one must assume an effective contact surface area between interacting particles; one possibility is to link such area to the cross sectional area of two neighboring particles, which for particles of same size would give $\gamma = 2\epsilon/(\pi r_0^2)$;

* The activation energy to separate two particles, $\epsilon + \Delta\epsilon$, whose magnitude relative to the thermal energy $k_B T$ controls the entity of the energy fluctuations that may induce structural rearrangements. This is important during structural formation but also for plastic and viscous rearrangements in the mature solid;

* The force to separate two interacting particles, $F_{max}$, which is the maximum positive $dU_{ij}/dr_{ij}$. It can be linked to the failure stress $\sigma_u$ of the material at zero porosity, e.g. by $\sigma_u = 4F_{max}/(\pi r_0^2)$;

* The interparticle distance $r_u$ at which $F_{max}$ is reached, which is linked to the strain at failure of the material, $\epsilon_u = (r_u - r_0)/r_0$

Table 1 compares these quantities for some of the interaction potentials from the previous section. Depending on the definition and parametrization of each potential, the terms in Table 1 may end up depending (physically or not) on particle diameter. Therefore, in the table we indicated the diameter $D$ for which the terms were computed. For example, the equivalent modulus $E_{eq}$ is independent of $D$ for the Lennard-Jones potential in Ref. [82], because the potential features $\epsilon \sim$

Table 1

Physical quantities linked to the parameters highlighted in Fig. 14, for reference: particle diameter $D$, equivalent elastic modulus $E_{eq}$, solid-solution interfacial energy $\gamma$, energy well depth vs. thermal fluctuations $(\Delta\epsilon + \epsilon)/k_B T$ (for room temperature $T = 298$ K), stress $\sigma_u$ and strain $\epsilon_u$ at failure.

<table><thead><tr><th>Potential description</th><th>D (nm)</th><th>E<sub>eq</sub> (GPa)</th><th>γ (mJ/m<sup>2</sup>)</th><th>Δε + ε / k<sub>B</sub>T</th><th>σ<sub>u</sub> (MPa)</th><th>ε<sub>u</sub> -</th></tr></thead><tbody><tr><td>Lennard-Jones in Eq. (13) [82,90]</td><td>10</td><td>63.6</td><td>811</td><td>31,000</td><td>1080</td><td>0.05</td></tr><tr><td>Harmonic in Eq. (15) [101]</td><td>10</td><td>63.6</td><td>87.6</td><td>3090</td><td>1289</td><td>0.05</td></tr><tr><td>Bonnaud et al.'s PMF [74]</td><td>5</td><td>141</td><td>349</td><td>12,900</td><td>5900</td><td>0.066</td></tr><tr><td>Masoumi et al.'s PMF [56]</td><td>5</td><td>425</td><td>417</td><td>4000</td><td>5800</td><td>0.061</td></tr></tbody></table>
$\sigma^3$. However, the same potential does not ensure the size-independence of the interfacial energy $\gamma$. Furthermore, for all the potentials in the table, the activation energy $\Delta\epsilon + \epsilon$ changes with particle size.

The Lennard-Jones potential in Ref. [82] was parametrized to capture $E_{eq}$ and $\epsilon_u$ from atomistic simulations. Therefore the potential has been used successfully to describe elastic properties of the C–S–H gel at the nano-to-micro mesoscale, such as elastic moduli, hardness [86,90], and creep modulus [102] from nanoindention experiments. This potential does not include a repulsive shoulder, $\Delta\epsilon = 0$, and entails a relationship between energy scale $\epsilon$ and interaction stiffness $\rho^{-1}$. Consequently, when calibrated to capture modulus and failure strain, as in Ref. [82], the potential yields $\epsilon$ that scales as $D^3$ and thus an interfacial energy $\gamma$ that is reasonable for particles of ca. 5 nm diameter (101 mJ/m$^2$, compared to 87.6 mJ/m$^2$ in Ref. [99]), but overestimated for larger particles (e.g. 811 mJ/m$^2$ for particles with $D = 10$, as shown in Table 1). In reality, the interfacial energy should scale as the area of the interface between particles, ~$D^2$, whereas the curvature of the potential at equilibrium, which is linked to the elastic modulus, should scale as $D^3$. In the Lennard-Jones potential in Eq. (13), everything scales with $\epsilon$ which in turn scales as $D^3$, hence the potential is appropriate for elastic properties but not to describe interfacial energies. This limitation was overcome in the harmonic potential in Ref. [101], which separately fits $E_{eq}$, $\gamma$, and $\epsilon_u$.

Both the Lennard-Jones and the harmonic potential feature very high activation energies $\Delta\epsilon + \epsilon$. This makes structural rearrangements from energy fluctuations very unlikely. The LJ-Yukawa potential in Ref. [85], instead, was parametrized using AFM data that featured a weaker cohesion. As a result, its $\Delta\epsilon + \epsilon$ is now in the order of magnitude of $k_B T$, leaving room for interesting structural evolution that may be particularly important for the C–S–H gel during early hydration. On the other hand, the mechanical properties and interfacial energies associated to this potential are orders of magnitude lower than for C–S–H during setting.

Table 1 also includes quantities computed from the atomistic simulations in Refs. [56,74]. They predict high values for all the quantities. Nevertheless, one must consider that the AFM experiments in Fig. 10.a indicated a much lower energy scale, as fitted by most of the LJ-Yukawa potentials in Eq. (14). Lower energy scales, although greater than in the AFM experiments, were also predicted by the simulations in Fig. 10.a, where the solvent was treated implicitly. All this suggests that the interactions between C–S–H particles may feature two different regimes: a strong cohesion at very short range (crystalline regime) and softer colloidal interactions in the medium range. These regimes would be separated by a high energy barrier but, thus far, experimental observations of these regimes are still lacking as well as attempts to consider them both in mesoscale simulations.

# 4. Formation of C–S–H mesostructures

This section presents nanoparticle-based simulations of C–S–H at the nano-to-micro mesoscale. The focus is kept on simulations of C–S–H formation, as creating model structures is the first step to then compute structural features and mechanical properties. Structure and mechanics will be covered in the next section.

## 4.1. Nucleation and growth models

Nucleation of minerals is a complex topic, as many non-classical mechanisms are possible [103]. The experimental literature on nucleation mechanisms of C–S–H is limited. First data came from Garrault-Gauffinet et al. [53]; by testing homogeneous nucleation in supersaturated solutions, they found that Classical Nucleation Theory could be used to interpret the results, returning critical nuclei of C–S–H as small as 2 molecules only. This was attributed to a low (liquid-solid) surface tension, i.e. 12 mJ/m$^2$ see Section 4.1 and Fig. 10. Recent experiments found that the homogeneous nucleation of C–S–H may actually follow a

non-classical pathway involving cations other than calcium [104]. Recent in situ experiments of C₃S hydration in a wet cell by soft X-ray spectromicroscopy reported Ca-containing nanoparticles in the solution around the C₃S particles and agglomeration of them by the C₃S surface [105].

While the experiments have tested fundamental aspects in homogeneous nucleation regimes, traditional semi-analytical approaches to modelling C–S–H formation [8] have typically assumed heterogeneous nucleation for two reasons: first, heterogeneous nucleation is usually thermodynamically favourable, hence it is believed to take over in a cement paste where the dissolving cement grains provide surfaces on which C–S–H can form; second, experimental results could be better fitted using semi-analytical models that assume heterogeneous nucleation and growth rather than homogeneous ones [106,107]. Coarse-grained simulations have shown that growth of C–S–H should be strongly limited by electrostatics [108], hence what is described as “growth” in nucleation and growth models, is more likely to reflect a mechanism of secondary nucleation of nanoparticles, or of locally homogeneous nucleation followed by aggregation. Mesoscale nanoparticle-based simulations of C–S–H formation have thus emerged as a way of bridging the molecular with the continuum scales.

A usual target of nucleation and growth models is the early hydration rate curve of cement paste, shown in Fig. 15. The curve describes the temporal evolution of the heat release rate, originating from the exothermic reaction of cement dissolution. Based on an underlying assumption that C–S–H precipitation is rate-limiting, the rate curve is commonly used as a signature of C–S–H formation. In reality, multiple processes may contribute to the various regimes in Fig. 15 [99,109,110]; a review of these processes is beyond the scope of this manuscript, but the interested reader can start from Ref. [111]. What is important here is that the same curve has become a target for nanoparticle-based simulations too.

The first nano-particle based simulations of C–S–H nucleation and growth were proposed by González-Teresa et al. [112], but these simulations did not include interaction potentials and did not aim at a statistical mechanical description of the system. Instead, the authors implemented a series of rules, listed in Fig. 16, which led to C–S–H nuclei forming and growing into clusters that eventually filled the space. The simulations produced rate curves for C–S–H precipitation similar to the experimental hydration rate curve, and they were also used to discuss cement setting [113].

Two other mesoscale models have been proposed more recently, which did not feature mechanical interactions between particles: (i) a model based on the attachment of triangular nanoparticles mimicking the growth of C–S–H sheets [115], and (ii) a model based on the Phase Field Method, for now only used in 2D due to its high computational cost [116]. The focus of these models is increasingly toward an accurate description of the morphology-dependent chemical kinetics of the system. On the other hand, since 2012 a thread of mesoscale simulations of

<figure><img src="image_18.png" /><figcaption>Fig. 15. Schematic of a typical cement hydration curve [111].</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Schematic of the nucleation and growth simulations in Ref. [114], including predicted heat release rate curve. Some of the images in the figure are reprinted from Ref. [114] with the permission of AIP Publishing.</figcaption></figure>

C–S–H based on mechanically interacting nanoparticles has developed: these will be discussed in the rest of this section.

## 4.2. Mesoscale simulations of C–S–H formation using interaction potentials

To date, most nanoparticle-based models of C–S–H at the mesoscale employ coarse-grained units that interact with effective potentials. As discussed in Supplement S1, a key point of such approaches is that if the coarse-graining unit is well represented, access to the physics depends only on having appropriately modelled the interactions (here the PMF) within the statistical mechanics framework. The simulation techniques consist of Molecular Dynamics (MD) and Monte Carlo (MC) where fluctuations and heterogeneity emerge as results, rather than being imposed. A strength of these simulations is that both structural evolution and mechanical stress are captured into the same, coupled, free energy landscape. In other models, such as those based on continuum micro-mechanics, stresses and chemical evolution of the various phases are rarely coupled at the constitutive level. In these mesoscale simulations, instead, stress and chemical transformations contribute at the same time to the same free energy landscape, hence both are accounted for when minimizing the free energy and exploring the landscape.

Mesoscale models of C–S–H aim to capture the reactive solidification process that takes place during cement hydration and setting. The appropriate statistical ensemble in order to account for the creation of C–S–H particles due to the dissolution-precipitation reactions, is the Grand Canonical (µVT) statistical ensemble, where µ is the chemical potential, V the volume and T the temperature (see Supplement S1.2). This defines an open system that can exchange particles with a reservoir of fixed chemical potential µ or activity z (as see Supplement S1.4.1). In the rest of this section, results from different particle-based simulations of C–S–H formation are compared with reference, where possible, to the hydration rate curve of cement paste in Fig. 15.

### 4.2.1. Monte Carlo simulations

Monte Carlo simulations of C–S–H formation, using spherical particles polydisperse in size and cohesive interactions, were first published in 2012 [82]. The interaction potential was the generalized Lennard-Jones in Eq. (13), with $\epsilon$ proportional to particle volumes and with exponent $\alpha = 14$ [82]. The simulations started with an empty box and particles were inserted or deleted with probability $P_{ins/del} \approx \min\{\exp\left[\frac{\pm\mu-\Delta U}{k_B T}\right], 1\}$, where $\Delta U$ is the change in total interaction energy in the system between the two states (see Supplement S1.4.1). The choice of $\mu = 0$ led to indefinitely increase the number of particles within the simulation box, but the Lennard-Jones interactions prevented the insertion of new particles overlapping with already existing ones (i.e. for which $\Delta U \gg 0$). In between insertion moves, the algorithm allowed for a number $N$ of Monte Carlo particle displacements, stimulating aggregation. Monte Carlo simulations are not meant to simulate time evolution of the system. Because of this limitation, the

simulation results were not compared with the hydration rate curve of cement. Different mesoscale morphologies were thus obtained by varying $N$ and the allowed polydispersity of the particles. The simulations produced meso-structures with different packing fractions and particle size polydispersity that were used to compute mechanical properties.

### 4.2.2. Monte Carlo coupled to Molecular Dynamics simulations

Ioannidou et al. [85] proposed a hybrid scheme of GCMC insertions/deletions followed by Molecular Dynamic in the Canonical NVT ensemble in order to introduce a time scale. The energy dissipation in the canonical ensemble was ensured by Nosé-Hoover thermostat [117,118]. A family of interparticle potentials motivated by AFM measurements (Eq. (14)) have been investigated using this simulation approach. This simulation approach with the aforementioned interparticle potentials successfully reproduced the hydration kinetics of C–S–H. Fig. 17(d) shows the comparison of the simulation volume fraction with the degree of hydration from experiments. Fig. 17 (a-c) shows snapshots of the GCMC/MD simulations while the system progressively densified. In the beginning of the acceleration regime the system forms clusters (Fig. 17 (a)) that eventually grew in elongated structures that formed a percolating gel network (Fig. 17 (b)). The gel formation indicates the inflection point between acceleration and deceleration. Fig. 17 (c) shows the formation of a porous solid C–S–H from the gel precursor.

The number of GCMC attempts over the number of steps of MD defined a kinetic rate $R = N_{GCMC}/N_{MD}$. The kinetic rate had an effect on the hydration curve and on the resulting meso-structures. The effect on the hydration curves is a shift in time whereas, on the meso-structures, higher kinetic rates induce more structural disorder [85]. The produced meso-structures were validated through a large set of available experimental results ranging from electron-microscopy imaging [90], nanoindentation tests [90], small angle X-rays and neutron scattering (SAXS and SANS) [90] and adsorption/desorption of nitrogen and water [119] that are discussed in Section 5.

Goyal et al. recently extended this model to account for heterogeneous precipitation of C–S–H particles using a spatial gradient of chemical potentials corresponding to higher nucleation of C–S–H near the surface of dissolving cement grains [89]. This enhances the effect of different interparticle potentials on the hydration curves and the structural and mechanical heterogeneities.

### 4.2.3. Kinetic Monte Carlo simulations

Another way to introduce temporal evolution is offered by the Kinetic Monte Carlo method that is based on transition rates among states. In 2017, Shvab et al. [98] published mesoscale simulations of C–S–H formation based on an off-lattice Kinetic Monte Carlo (KMC) scheme that explicitly accounted for the rates of C–S–H precipitation and dissolution. The KMC algorithm lets C–S–H particles appear or disappear, deriving the characteristic time of the events from coarse-grained rate equations. To compute the rates, Shvab et al. considered a
<figure><img src="image_20.png" /><figcaption>Fig. 17. Snapshot from the GCMG and MD precipitation simulation with increasing volume fraction of solid (a)η = 0.08, (b)η = 0.16 and (c)η = 0.26. Only the bonds between the particles are displayed, using lines. The colour code corresponds to the number of nearest neighbours. Red is 10–12, blue 6–9, and gray 0–5 neighbours; (b) Comparison between solid volume fractions from the simulations and degree of hydration from experiments. The simulations capture an initial accelerating phase and the subsequent deceleration. Figures adapted from Ref. [85] with permission of the Royal Society of Chemistry [85]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

simulation box containing two phases: a solid represented explicitly by mechanically interacting spherical particles with diameter D = 10 nm, and a surrounding implicit solution represented by the average concentrations of ions in it.

The per-particle rate equations were obtained starting from the rate of forming or dissolving a single molecule of solid C–S–H. Employing the law of mass action and assuming one-step reactions, the unimolecular rates were expressed as:

$$r_{\text{diss}} = ka^2 \exp \left[ (1-\chi) \frac{-\Delta U_{ij,\text{diss}}(r) + \gamma \Delta \Omega}{k_B T} \right] \quad (16)$$

$$r_{\text{prec}} = ka^2 \exp \left[ \chi \frac{-\Delta U_{ij,\text{prec}}(r) - \gamma \Delta \Omega}{k_B T} \right] \beta \quad (17)$$

k is the rate constant of C–S–H dissolution. a = 0.65 nm is the linear size of a C–S–H molecule. ΔΩ is the change in surface area of a particle when the molecule under consideration is dissolved or precipitated. γ is the solid-solution interfacial energy of C–S–H. ΔUij is the change in interaction energy in the system upon dissolving or precipitating the molecule (this changes the size of the pertaining particle, hence the interactions between particles if their size-dependence is accounted for). Together, the ΔUij + γΔΩ term represents the excess enthalpy of the solid compared to its standard state. Transition State Theory assumes that all this enthalpy is dissipated upon dissolution, hence does not feature in the precipitation rate. Instead, Shvab et al. allowed for such enthalpy to also contribute to the precipitation rate, assuming that a fraction 0 < χ < 1 is still present in the activated complex (see [98] for more details). Transition State Theory is recovered for χ = 0. Finally, in Eq. (17), β is the saturation index of the solution with respect to C–S–H precipitation, viz. the ratio of the ion activity product and the equilibrium constant for the reaction.

The formation (or dissolution) of a D = 10 nm C–S–H particle, containing thousands of molecules, was then modelled as the serial precipitation (or dissolution) of n = D/(2a) C–S–H molecules in radial direction:

$$R^{\text{del}} = \left( \sum_{i=1}^{n} r_{\text{diss},i}^{-1} \right)^{-1} \quad (18)$$

$$R^{\text{in}} = \kappa \left( \sum_{i=1}^{n} r_{\text{prec},i}^{-1} \right)^{-1} \quad (19)$$

κ is a parameter related to the concentration of calcium and silicon ions in solution, which are all assumed to provide initial nucleation sites (from the aforementioned result that critical C–S–H nuclei are very small). In implementing Eqs. (18) and (19), the ΔUij term contained in r was expressed as a function of the entire interaction energy between a particle and its neighbour, for which the authors employed the generalized Lennard-Jones potential in Eq. (13).

Fig. 18 shows two simulation results for heterogeneous C–S–H precipitation from solutions with assigned β(t). Predicted rate curves are compared to a typical one from Chan’s model (BNG). Irrespective of β(t), the KMC simulations predict initial acceleration reflecting the increasing surface of the hemispherical domains formed by precipitated particles. A peak rate is also always predicted, when neighboring domains of precipitates impinge laterally on each other: this was originally discussed in the seminal work and model of Garrault and Nonat [120]. The time-dependence of β becomes important after the peak rate, as deceleration is only captured when ion depletion in solution is considered, leading to the decreasing β(t) curve in Fig. 18.a (the curve was taken from Ref. [99], where this mechanism for deceleration was originally modelled at the microscale).

Recently, these KMC simulations have been extended to use net reaction rates [100], although thus far only for the dissolution of tricalcium silicate and considering that each particle represents one molecule of solid only.

# 5. Properties of C–S–H gel and cement paste in mesoscale simulations

## 5.1. Structural properties

Particle-based simulations can be used to explore the growth of C–S–H gels and cement paste. The simulated particle configurations can be further analysed in terms of structural characteristics and subjected to numerical mechanical testing. This allows for comparison with experimental data and analysis of structural evolution of the growth
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Kinetic Monte Carlo simulation results from Shvab et al. [98]. (a) Temporal evolution of saturation index from Ref. [99]; (b) Rate curves predicted by KMC simulations with different assumptions on $\beta(t)$, and comparison with a typical curve from Chan's model (BNG); (c) Snapshots of the evolving mesoscale structure of C–S–H during a typical KMC simulation (the numbers correspond to different points on the rate curves in (b)). Reprinted with permission from Ref. [98]. Copyright 2017 American Chemical Society.</figcaption></figure>

models.

### 5.1.1. Volume fraction

The configurations of particle mesoscale simulations are characterized by their volume or packing fraction, which depends on the shapes, sizes and interactions of/between the particles. A crystalline system of monodisperse spheres can pack up to 74%. C–S–H exhibits areas of crystalline packing but overall is amorphous. Hence, models with spherical particles, monodisperse in size, aim for volume fractions up to the close random packing of 64%. Size polydispersity of spheres allows to increase the random close packing [121]. Masoero et al. investigated size polydispersity of C–S–H spherical particles and reached volume fractions up to 0.72 [82]. Masoumi et al. reached similar volume fractions using ellipsoid particles of polydisperse sizes [56] and Yu et al. by varying the aspect ratio [96]. Masoumi et al. have summarised results on volume fraction and particle asphericity showing that aspherical particles percolate at a lower volume fraction [56]. The final volume fraction of the configurations can be obtained by different preparation methods, such as MD in isothermal ensemble (NVT) [88], MD in isobaric-isothermal ensemble (NPT) [56] or with GCMC precipitation models (see Section 4).

The volume of C–S–H increases during hydration following a sigmoidal evolution with time [122]. This experimental result has been used to connect to the time evolution of the volume fraction in simulated configurations. An example of the time evolution of the volume fraction in precipitation models is shown in Fig. 17(d) [85]. Goyal et al. investigated the evolution of volume fraction in heterogeneous surface growth models [89]. In heterogeneous growth models, the probability of C–S–H precipitation is higher near the surface of other solids, therefore computing local volume fractions (lvf) allows to rationalize the C–S–H

aggregation process [89]. The lvf distributions show that the local volume fraction near a surface can be 1.5 to 2.5 times higher than away from it, depending on the particle interactions.

The distribution of local volume fractions is also useful in homogeneous precipitation models as it allows to explore heterogeneities within the configuration. Such distributions have been computed in different models [56,89,90,92,123]. The definition of the local volume for the computation of the lvf differs between models. It is either defined in a spherical volume of few particles diameter around each particle [90,92] or based on a 3D grid [56]. Distributions of local volume fractions for spherical and ellipsoidal particle simulations are shown in Fig. 19. a and b respectively [56,92]. Overall, simulation data for both spherical or ellipsoidal particle models agree that, as the average volume fraction of the configuration increases, the distribution of local volume fractions gets skewed with the mode moving to higher values and a long tail developing toward the lower volume fractions [56,90,92,123]. This indicates that the areas of highly densely packed C–S–H grow locally within the configuration.

The solid volume fraction is interconnected to porosity, which is one minus the volume fraction. The Powers-Brownyard (PB) model of C–S–H estimates the capillary pore volume in the hardened cement paste from the water-to-cement ratio [124]. Ioannidou et al. used the PB model to connect the volume fraction of the mesoscale simulation configurations to the water-to-cement (w/c) ratio in order to compare structural and mechanical properties with experimental data [90,92]. Fig. 19.c shows a comparison of local volume fraction distributions of C–S–H estimated from nano-indentation experiments in cement paste samples and from simulations [90]. The nano-indentation measurements have been coupled with chemistry information of elemental maps to select only the C–S–H volume fractions.

### 5.1.2. Clusters, connectivity and local order

Particles simulations give the opportunity to analyse further the positions of particles and compute spatial or temporal correlations. Common quantities to characterize the evolution of the microstructure during hydration and setting is the cluster size distribution, coordination number, gyration radius, bond orientational order parameter, pair distribution functions and structure factor. The cluster size distribution has been used to characterize the stacking of ellipsoids or the gelation process of C–S–H. The cutoff distance for two particles to belong in the same cluster is determined by the position of the first minimum after the main peak in the pair distribution function. This is also the first coordination shell. The gyration radius indicates the morphology of the clusters and the bond orientational order parameter the type of local packing. The pair distribution function describes how density varies as a function of distance from a reference particle.

Masoumi et al. computed the cluster size distribution of ellipsoidal particles for configurations of different volume fractions [56]. Interestingly, the estimated average thickness is between 3 and 6 nm for all volume fractions and ellipsoid sizes, in agreement with SANS observations [125].

For the precipitation model described in Section 4.2.2 the cluster size distributions, coordination number, pair distribution function, gyration radius and bond orientational parameters were computed for low and high Ca(OH)₂ concentration potentials and different kinetic rates [85,126]. The combined analysis of these structural quantities suggested that the low Ca(OH)₂ concentration potential promotes the growth of elongated clusters with local packing akin to Bernal spirals that branch and impinge during precipitation. This potential has thermodynamically stable clusters at low volume fraction that promote quick filling of the simulation box and geometric percolation. On the contrary, the high Ca (OH)₂ concentration potential grows large, compact clusters with crystalline local packing that favor densification. Comparison of the same structural properties with MD equilibrium simulations of the same interactions potentials showed that the local packings are thermodynamically stable, hence persist through a broad range of kinetic rates
<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Distributions of local volume fractions η<sub>local</sub> for simulations with spherical (a) and ellipsoidal (b) particles for different overall volume fractions η. (c) Comparison of packing fractions estimated from nano-indentation experiments and simulations. Reproduced from (a) ref. [92] with permission from ASCE, (b) ref. [56] with permission from Wiley and (c) ref. [90].</figcaption></figure>

[85,88]. Fig. 20 a and b show heatmaps of bond orientational parame-
ters distribution for low and high Ca(OH)₂ concentration potentials
respectively [88]. Low Ca(OH)₂ potential shows a high intensity at
Bernal spiral packing whereas high Ca(OH)₂ potential shows high in-
tensity at face-centered cubic and hexagonal close packed crystals.

The geometric percolation for a precipitation model as a function of
the kinetic rate (as defined in Section 4.2.2) was investigated in
Ref. [127]. The authors computed size, number, shape and percolation
volume fraction for configurations of about 1700 spherical particles of 5
nm in a cubic box of 60 nm, with particles interacting via a Mie po-
tential. Despite the small system size, they reported that high kinetic rate
simulations produced more clusters of smaller size and elongated shape
and pushed geometric percolation to higher volume fractions.

Cluster analysis in surface growth simulations of the low Ca(OH)₂ concentration potential resulted in percolating clusters with low density whereas high Ca(OH)₂ concentration potential resulted in locally dense and crystalline crust on the cement grain surfaces and an anisotropic network in between as shown in Fig. 20 c [89].

5.1.3. *Spatial correlations*

Non-invasive experimental techniques to investigate correlations at
different length scales are Small Angle X-ray and Neutron scattering
(SAXS and SANS). Scattering originates from density fluctuations. In
porous materials, such as cement paste, the fluctuations occur at the
inter interfaces between void (pore) and matter (solid) and depend on the
orientation and roughness of interfaces. The scattering intensity I(q) is
computed as the Fourier transform of the auto-correlation of fluctua-
tions in local density [128]. When the scattering intensity is computed in
simulation configurations, the system size must be large enough to not
induce cutoff effects in real space, which would cause artifacts in Fourier

space.

Small angle scattering is also connected to chord length distribution
functions, providing a link between imaging techniques and scattering
[129]. Chords are line segments that have both ends at a pore-solid
interface. Such distributions can also provide average pore size, pore
size distribution, surface roughness and spatial correlation between
points belonging to the solid, the pore network, or the interface between
them. Chord length distributions also provide a way to estimate the
specific surface area of porous materials. For sufficient sampling of the
materials, the specific surface area is $S_{sp} = 4\phi/[\langle l_p > \rho(1-\phi)\rangle]$, where
$\langle l_p \rangle$ is the average of the pore chord length distribution, $\rho$ is the density
of a nanosized C-S-H crystal, and $\phi$ is the porosity of the porous sample.
Specific surface areas were computed for different cutoffs of the smallest
chord length in Ref. [90] and used to rationalize experimental mea-
surements of specific surface area that vary broadly upon drying con-
ditions and choice of adsorbate. Goyal et al. computed the evolution of
specific surface area during hydration [89].

Cement pastes have been extensively investigated with SANS/SAXS. An interesting result is that the scattering intensity invariably behaves as a power of the *q* vector over a wide range of *q* and this result persist regardless the composition or age of hardened paste. Fig. 21 shows experimental measurements and lines of powers of *q*. The characteristic *q*⁻³ dependence is shown to extend over about two orders of magnitude. Such extended spatial correlations of cement pastes could not be obtained with dense packings of spheres polydisperse in size [128]. Fig. 21. b shows that the *I*(*q*) of a very dense, highly polydisperse packing, simulated with the process in Section 4.2.1 [82], has a *q*⁻² dependence. Scattering intensities computed from ptychographic imaging show that this is relevant for inner product C-S-H over a limited range of *q* [130].

Fig. 21.b shows that the I(q) of a low density packing simulated with

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Bond orientational parameters heatmaps for (a) low Ca(OH)₂ concentration potential showing high intensity at Bernal spiral (BS) packing and (b) high Ca (OH)₂ concentration potential showing high intensity in hcp or fcc crystal both at volume fraction 0.052. (c) Snapshot of the surface growth model for high Ca(OH)₂ concentration potential. The dark colour indicates higher density crust close to cement grain surfaces. Reproduced (a) and (b) from ref. [88] and (c) reprinted with permission from ref. [89]. Copyright 2020 American Chemical Society.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Scattering intensity computed from simulations and compared with experiments. The configurations show (a) an extended $q^{-3}$ power law and (b) limited $q^{-3}$ dependence for $\phi = 0.36$. The dense configuration of $\phi = 0.76$ has $q^{-2}$ dependence. $\eta$ and $\phi$ both indicate volume fraction. Reproduced from (a) ref. [90] and (b) ref. [131] with permission from ASCE.</figcaption></figure>

the process in Section 4.2.2 [85] has a very limited $q^{-3}$ dependence [131]. The first simulations data that could attain an extended spatial correlations of $q^{-3}$ are shown in Fig. 21.a. The $q^{-3}$ dependence persist for about two orders of magnitude for both a very porous and a denser system (volume fraction 0.33 and 0.52 respectively). These simulated configurations were also produced with the process in Section 4.2.2 [85], however to obtain long range correlations over a wide range of $q$, system sizes of about 0.6µm had to be simulated, using about half a million particles [90]. In the computation of $I(q)$, finite size effects and boundary conditions must be handled with care. In Fig. 21.a the plateau at low $q$, where the simulation data are restricted by the system size, is associated to a correlation length of about $\pi/q = 40$ nm, ten times smaller than the box size 0.6µm.

Combining the scattering intensity with chord length distribution analysis, Ioannidou et al. found that the extended $q^{-3}$ behavior of the small angle scattering originates from the coexistence of surface roughness of the capillary pores at the scale of 3 - 40nm with a fractal pore network [90,132]. The authors suggested that the out-of-equilibrium process that was used to produce the growth of the structure (see Section 4.2.2 produced such a complex mesopore network [85].

Other growth models, not based on interaction but rather on geometrical triangulation in 2D with a prescribed set of rules, have been developed to estimate the scattering intensity [115]. However, the system size of such structures was small and the computed scattering intensity was significantly affected by finite box effects.

Goyal et al. computed the scattering intensity of heterogeneously grown structures [89]. In these configurations, the packing and morphology is anisotropic along the z direction; high density packed areas exist in boundaries and more porous in the middle of the simulation box. The scattering intensity calculation of anisotropic structures is more complex as it requires the computation of autocorrelation of density fluctuations in each orientation before performing the 3D Fourier transform. To overcome this, the authors used a method based on 2D projections [133]. They found that at high $q$ the signal in all directions is oscillating due to monodisperse size of spherical particles. At low $q$, the comparison of $I(q)$ along the z anisotropic direction and the isotropic directions suggested that the differences are due to high differences in density along the z direction. However, for high Ca(OH)$_2$ concentration potential the anisotropy persists also in the middle zone, suggesting that is inherent to the potential.

### 5.1.4. Pore network

At the sub-micrometer scale, hardened cement paste is composed of gel pores ($< 5$ nm), small capillary pores (10–100 nm) and large

capillary pores (100–1000 nm). Pores of sizes larger than 50 nm cannot be captured easily in mesoscale simulations. In particle-based simulations, the pore size distribution (PSD) can be directly calculated from the mesoscale structure of the aggregated particles. The PSD can be obtained using Monte Carlo (MC) sampling of the geometric pore space in simulation configurations [134]. The pore volume is probabilistically probed with either 3D shapes, such as spheres, or with 1D chords. In the first case, the PSD is obtained from the diameter of the largest sphere that fits at a given point [134]. Sufficient sampling provides a representative PSD. This method was used in Ref. [89,90,92] to estimate the PSDs from the simulations.

The PSDs of the simulation samples in Refs. [90,92] showed a first peak of gel nanopores and wide range distribution of small capillary mesopores in the 4–30 nm range (see Fig. 22c). The majority of these pores were connected with few of them being isolated. These PSDs were in agreement with pore size distribution estimated from nuclear magnetic resonance spectroscopy (NMR) [135]. In Ref. [92] the results were compared with the nitrogen adsorption/desorption experimental results of Mikhail et al. [136] for different ordinary Portland cement pastes, prepared with different initial water-to-cement ratios.

In the case of chords, a line segment starting from a random point grows in random directions until it reaches solid boundaries [129]. The same methods can be applied to digitised tomography images [129]. This provides a direct comparison of the pore structure between experiments and simulations. In Ref. [90] the distribution of chord lengths in the pore obtained from particle simulation configurations were compared to ones obtained from experimental SEM images.

Yu et al. computed PSD for ellipsoidal packings based on a 3D grid algorithm that grades different levels of voids and that clusters them based on a breadth-first search (BFS) algorithm, assigning to connected clusters a spherical volume [137]. The PSD computed for models with ellipsoid of different aspect ratios follow a normal distribution, with high aspect ratios leading to larger pores. Masoumi et al. provided a 2D voxel map of the pore volume in ellipsoidal particle configurations in box sizes 75 nm and various volume fractions [56]. Goyal et al. computed the PSD using MC geometrical sampling [134] for early and late stages of hydration of C–S–H and showed the evolution of the average pore size for C–S–H close to clinker surface and away from it [89].

Another common way to assess pore size distribution (PSD) are physisorption experiments. In physical adsorption the molecules of an adsorbate (vapor) attach to the surface of the adsorbent (solid), similar to condensation. In capillary pores of 2–50 nm, multilayer adsorption from the vapor into a porous medium proceeds to the point at which pore spaces become filled with condensed liquid from the vapor. The
<figure><img src="image_25.png" /><figcaption>**Fig. 22.** Isotherms and pore size distributions. Adsorption/desorption isotherms for (a) nitrogen and (b) water in cement paste. Symbols show the mass ratio of wet over dry sample, $M_{adsorbent}/M_{dry}$ (blue squares and red dots for adsorption and desorption respectively) from experiments of Ref. [136] for (a) and Ref. [142] for (b). Black lines show simulated isotherms on cement model configurations with realistic 3D pore network of water to cement ratio (w/c) 0.52. For water sorption (b) the degree of saturation (DOS) estimated from the simulations is shown at the right y axis. (c) and (d) show pore size distributions calculated by Monte Carlo method on the three cement model configuration, and from applying the Kelvin's equation to the isotherms of nitrogen and water respectively. Y and X-axis values are normalized probability density and pore radius respectively. The terms hardened and aged referred to samples with and without eigenstresses. Reproduced from ref. [119]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

measurement of the quantity of adsorbent as a function of its vapor pressure at constant temperature provides the most common experimental observation: the adsorption/desorption isotherm. The amount of adsorbed vapor at a given temperature and pressure depends on the interaction between the vapor and the solid. Note that the isotherm of a solid-gas system should be reversible if there is no irreversible change in the nature of the solid and/or gas, since there is one state of equilibrium at a given vapor pressure. At lower pressures this is what is usually observed, but at higher pressures a hysteresis arises. This is commonly related to the pore structure of the solid and to the formation of a meniscus; the difference in the process of formation and recession of this meniscus causes the hysteresis. The distribution of capillary bridges and menisci is key to understand the capillary stresses induced on a porous solid matrix.

PSD experimentally derived from adsorption/desorption isotherms is based on the Kelvin's equation and surface adsorption models that assume that pores are cylinders or spheres [138,139]. A recent theoretical model for adsorption and PSD estimation in cement paste was based on such hypotheses for pore sizes and connectivity [140]. Zhou et al. developed a simulation scheme based on mesoscale particle configurations of C–S–H and lattice gas Density Functional Theory [141]. Fig. 22a and b show adsorption/desorption isotherms obtained from the fluid distribution around C–S–H particles, compared with experimental data for nitrogen and water respectively [119]. The sorption curve for nitrogen displays only minimal hysteresis, whereas for water displays a significant hysteresis, both in agreement to experimental observations in Refs. [136,142]. The water hysteresis arises from ink-bottleneck states upon desorption due to constrictions in the pore network of C–S–H

[141]. Comparison of PSDs obtained from the Kelvin equation and the geometrical Monte Carlo sampling is shown in Fig. 22 c and d for nitrogen and water respectively. The Kelvin equation PSD for nitrogen show similar features to the geometrical Monte Carlo PSD, whereas for water, the Kelvin equation features a bimodal distribution for gel pores and lack of large capillary pores. This may be because water molecules interact strong with C–S–H particles and the interlayer spaces within the C–S–H [143], which nitrogen does not. All this suggests that nitrogen sorption may be more appropriate for characterizing the mesopore network of C–S–H [119].

Capillary stresses due to changes in relative humidity cause material deformation, as observed in drying shrinkage experiments [144–146]. The state of the art in atomistic modelling of fundamental mechanisms driving drying shrinkage has been recently reviewed [147]. At the larger mesoscale, in the model proposed by Zhou et al., capillary stresses due to imbalanced fluid distributions around the particles can be computed and used as input to MD simulations [141]. Zhou et al. computed the drying shrinkage for various relative humidities and compared those to experimental values [141].

## 5.2. Mechanical properties

Various mechanical properties have been computed using mesoscale simulations of C–S–H. This section starts with the eigenstresses generated during the precipitation of the gel, then will move to short-term properties (elastic moduli and strength), and will finally address failure and long-term viscous behaviours.
### 5.2.1. Eigenstresses during C–S–H formation

During its formation, the C–S–H phase develops autogeneous strain generating self-equilibrated eigenstresses. If a sample is unconfined and free to undergo volume changes, the average eigenstresses will be null but still the sample will feature a heterogeneous field of local eigenstresses. In 1997, Bažant et al. [148] indicated these eigenstresses as the promoters of long-term basic creep in concrete, arguing how they could originate from processes such as competitive growth of C–S–H nano-particles and other crystals (e.g. calcium hydroxide) as well as heterogeneous drop of pore water chemical potential during hydration. Based on these mechanisms, Bažant et al. predicted an order of magnitude of 150–200 MPa for the average intensity of the local eigenstresses at sub-micrometer length scales. Samples whose deformations are constrained will still develop a field of local eigenstresses and, on top of that, they will also feature non-zero average eigenstresses, experimentally computed in the 0–40 MPa range with non-monotonic evolution during hydration [149].

In a particle simulation, the components of the eigenstress tensor $\sigma$ on particle $i$ can be quantified using the virial expression [150,151]. For pairwise interactions only and negligible kinetic energy, this is:

$$ \sigma_{ab,i} = \frac{1}{2V_l} \sum_{i \in V_l} \left[ \sum_{j \neq i} (r_{a,i} F_{ij,b,i} + r_{a,j} F_{ij,b,j}) \right] \quad (20) $$

where $V_l$ is the volume over which the contributions to the eigenstresses are averaged. $a, b$ are Cartesian coordinates, $j$ are the neighbours interacting with particle $i$, $r_a$ is the $a$-component of the position vector, and $F_{ij, b, i}$ is the $b$-component of the force on particle $i$ coming from its interaction with particle $j$. If the focus of the simulation is on the average eigenstresses in the whole system, then $V_l$ should equal the volume of the entire simulations box $V_{box}$. Instead, if one is after meaningful local stresses, one should use $V_l \ll V_{box}$ but, at the same time, $V_l$ should be sufficiently large to capture a representative elementary volume, RVE.

Sometimes a related but different quantity is computed: the per-particle contribution to the average stress, whereby only one particle $i$ is considered at a time. In this case taking $V_l = V_{box}/N$, where $V_{box}$ is the volume of the simulation box and $N$ is the number of particles in it, leads to a quantity whose average over the $N$ particles coincides with the average eigenstress over the whole simulation box. A variation of the per-particle stress contribution uses the per-particle Voronoi volume instead of just $V_{box}/N$.

In a simulation, only very ordered structures of particles can be free from local eigenstresses, viz. have all the particles sitting in local minima for all the interaction potentials, such that $F_{ij} = 0$ everywhere: see the “crystalline” data-set in Fig. 23.b. All the C–S–H simulations presented in Section 4 produce disordered structures, hence all would feature local eigenstresses. However, such eigenstresses would not result only from the physical process of C–S–H formation; in part they would also be intrinsic to arbitrary steric constraints imposed by model-specific assumptions, e.g. on the size, shape, and polydispersity of the particles. Therefore, when simulating formation eigenstresses one should take care of controlling and possibly minimizing the impact of such steric constraints on the predicted stresses.

Ioannidou et al. [92] computed local eigenstresses using a version of the potential in Eq. (14), without Yukawa shoulder and with $\epsilon = 0.15$ k_BT. This low energy scale and the employed Grand Canonical Monte Carlo formation algorithm (see Section 4.2.2) ensure that the system can explore the potential energy landscape quite efficiently, leading to better equilibrated structures and thus minimizing the effect of steric constraints on the resulting eigenstresses. The eigenstresses were then computed fixing the structure and applying the pressure unit $\epsilon/\sigma^3 = 0.55$ GPa that the authors estimated for hardened C–S–H gel. This led to eigenstresses in the $\pm 0.1$ GPa range, with intensity increasing with the local packing density of the structure: see the increasing standard deviation of the curves in Fig. 23.a.

Masoero and Di Luzio [101] used a different approach to compute

<figure><img src="image_26.png" /><figcaption>Fig. 23. Simulated fields of eigenstress and per-particle stresses from C–S–H formation on structures, all relaxed to zero average axial stresses. (a) Snapshot of local packing densities $\eta$, which correlate with local eigenstress intensity, and distributions of hydrostatic eigenstress (pressure) intensity $a$ for different average $\eta$ of the structures (reproduced from Ref. [92], with permission from ASCE); (b) Snapshot of per-particle shear stresses straight from a packing algorithm akin to Section 4.2.1, and distribution compared to a crystalline structure featuring no eigenstress (reprinted from Ref. [91], Copyright (2019), with permission from Elsevier); (c) Snapshot of tensile (red) and compressive (blue) hydrostatic component of per-particle stresses, and distribution of intensity before and after particle inflation as explained in the main text (reprinted from Ref. [101], Copyright (2020), with permission from Elsevier). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
per-particle contributions to the eigenstresses. Initial C–S–H structures were created using a packing algorithm similar to that described in Section 4.2.1, with monodisperse particles interacting via harmonic potentials of the type in Eq. (15). Strong per-particle stresses of ±1 GPa emerged as a result of the steric constraint: see the dataset named “after monodisperse packing” in Fig. 23.c. Similar intensities were previously obtained by Liu et al. [91], see Fig. 23.b, using a similar packing algorithm and a slightly different, but similarly intense, interaction potential. Masoero and Di Luzio [101] found that these per-particle stresses depended indeed on arbitrary assumptions in the packing algorithm, which did not represent a physical process of C–S–H formation. To correct for this, a fraction $\delta$ of the particles in the structure were inflated from a diameter of 5 nm to 7 nm, and the average stress on the whole simulation volume was relaxed back to zero by energy minimisation, while letting the simulation box expand to remove average eigenstresses and preserve local stresses only. This mimicked physical processes taking place during C–S–H formation and generating local pressure, e.g. crystallization pressure or disjoining pressure in the pore solution. The result was an even stronger field of per-particle stresses, but consistent across different packing schemes to obtain the initial monodisperse structure: see Fig. 23.b. The resulting per-particle stresses were capped by local plastic deformations to ±2 GPa, with average positive and negative intensity in the 150–200 MPa range, consistent with the aforementioned estimation by Bațant.

The simulated stresses in Fig. 23 differ by one order of magnitude, but this is because Fig. 23.a shows eigenstresses averaged over a sufficiently large RVE (see details in Ref. [92]) whereas Fig. 23.b,c show per-particle stresses, with $V_l$ set to the tributary volume of a single-particle. Another contribution to the difference in stress intensity may be that the structures in Fig. 23.a had maximum packing density $\eta = 0.52$, whereas the structures in Fig. 23.b,c had $\eta \approx 0.64$; the intensity of the eigenstresses increases rapidly when approaching the close packing limit, which induces strong steric constraints [92].

The presented simulations have mitigated but not completely removed the issue of intrinsic steric constraints causing artificial contributions to the eigenstresses. A step to further reduce such artifacts would be to model partial dissolution/reprecipitation of C–S–H by letting particles change their shape anisotropically, hence simulating physical interactions between growing and shrinking nano-domains. Another improvement would be to couple particle simulations with pore saturation models (e.g. those shown later in this section), which would help account for the now-overlooked contribution of capillary stresses during C–S–H formation.

### 5.2.2. Elastic moduli

Most nanoparticle-based simulations of C–S–H in the literature report computed elastic moduli. González-Teresa et al. [112] were first to estimate elastic moduli from mesoscale models of C–S–H. In particular, they used a purely geometric (no interaction forces) particle-based model to compute the density $\rho$ of so-called “low-density” and “high-density” C–S–H structures. Subsequently they used linear functions $E(\rho)$, separately obtained from previous molecular simulations [152], to estimate various elastic moduli (here collectively called $E$). The results agreed well with experimental data on C–S–H.

With the advent of interaction-based nanoparticle simulations of C–S–H, the most common method to compute elastic moduli became the finite deformation method. This consists in starting from a C–S–H structure at zero average stress, applying an increasing finite deformation to the simulation box in one or more directions, and computing the average stress components during the deformation process (using Eq. (20) with $V_l = V_{box}$, for all $a, b = x, y, z$ components). The slopes of the resulting stress-strain curves in their initial linear regime quantify the elastic moduli. Depending on the applied deformations, different moduli can be extracted, and the full elastic tensor can be constructed from multiple simulations. An alternative way to compute moduli is by monitoring the thermal fluctuations of the average stress components at

constant volume [153], as done for example in Refs. [82,86]. The advantage of the fluctuations method is that the whole elastic tensor is obtained from a single simulation. The drawback is that the reference configuration must be in a rather stable energy minimum, because if fluctuations induce even small irreversible particle rearrangements, the error on the computed moduli can be large. Instead, small irreversible rearrangements have only a limited impact when the finite deformation approach is used, causing small drops in the stress-strain curves which can be easily neglected when computing the initial slope.

When the finite deformation method is used, an important variable is the applied strain rate. Most nanoparticle simulations of C–S–H perform energy minimisation between subsequent deformation steps, because the energy barriers in the simulations are usually so large that the impact of thermal fluctuations can be neglected: see Table 1, where in most cases $\epsilon + \Delta\epsilon \gg k_B T$.<fnref n="2" />

The indentation modulus $M$ is typically computed, because it allows direct comparison with experimental results from nanoindentation. These results also relate local $M$ with local packing density $\eta$ of the C–S–H. The numerical calculation of $M$ is usually carried out in one of two ways: (i) via finite deformations in uniaxial strain, while keeping the simulation box fixed in the perpendicular directions, which returns $M$ as the slope of the axial stress-strain curve, or (ii) by computing two other moduli, e.g. the Young modulus $E$, the shear modulus $G$, the bulk modulus $K$ or the Poisson coefficient $\nu$ (or the full elastic tensor with the stress fluctuations method) and combining them via linear elastic relationships, such as $M = \frac{G(4G-E)}{3G-E}$. In fact, the ability to predict $M$ and its relationship to $\eta$ has become the basic validation for nanoparticle models of C–S–H that aim at reproducing mechanical properties, also going beyond the linear elastic regime.

Fig. 24 collects various simulated relationships between elastic moduli and packing density. The first such relationship was presented in 2012 by Masoero et al. [82] and is reproduced here in Fig. 24.a. The structures were obtained with the packing algorithm described in Section 4.2.1, with modified Lennard-Jones potential in Eq. (13). Particle size polydispersity was used to reach $\eta > 0.64$. Similar results are shown in Fig. 24.b, from a slightly improved way of computing the contact surface and thus the bond energy between adjacent particles [84,86]. The results in Fig. 24.c [90] extended the range down to very low $\eta$, for which nanoindentation experiment cannot provide reliable data. The simulation method here featured an initial packing stage using the Grand Canonical Monte Carlo method in Section 4.2.2, with the softer Lennard-Jones-Yukawa potential in Eq. (14). Subsequently, when the desired $\eta$ had been achieved, the interaction potential was replaced with the stronger one in Eq. (13) before computing the elastic moduli. The results in Fig. 24.d were obtained with the same approach in Fig. 24a,b, but a thorough investigation on the impact of polydipsersity was conducted in the manuscript [91], as well as comparing the moduli of disordered C–S–H structures to those of crystalline structures featuring the same interaction potentials. This showed that $M$ does not depend exclusively on $\eta$, or on the polydispersity level, or on the structural order or disorder, but on all of the above. Fig. 24.e shows recent results on Young modulus, which is quite similar to $M$ but does not coincide with it. The underlying simulations [93] used the same packing algorithm as in Fig. 24.c and then switched to peridynamics simulations parametrised on strong mechanical properties similar to those leading to the entries in Table 1 for the Lennard-Jones and the Harmonic potentials. The results in Fig. 24.f instead were obtained from simulations [94] employing the Lennard-Jones-Yukawa potential in Eq. (14) but parametrises similarly

<footnote n="2">² Sometimes the simulations involving energy minimisation are called “quasi static”; this may be misleading, because energy minimisation takes the system to the nearest metastable energy minimum, whereas simulations at finite temperature and finite strain rate might better explore the energy landscape and find deeper local minima. Using energy minimisation, therefore, might be more representative of loading regimes with large strain rates.</footnote>
<figure><img src="image_27.png" /><figcaption>**Fig. 24.** Simulated indentation moduli *M* [82,86,90,91] and Young moduli *E* [93,94] as functions of C–S–H packing density *η* and, in some cases, comparison with experimental results from [90,154] Permissions to reproduce images: (a) with permission from the American Physical Society [82], (b) with permission from the Royal Society of Chemistry [86], (d,e,f) from Ref. [91], Copyright (2019), from Ref. [93], Copyright (2021), and from Ref. [94], Copyright (2021), with permission from Elsevier.</figcaption></figure>

to the Lennard-Jones potential in Eq. (13). The computed Young moduli in Fig. 24.f are in line with those from the other simulations and with the experimental values. However the physical interpretation for the resulting Yukawa shoulder, with intensity of hundreds of $k_B T$, is not clear.

### 5.2.3. Strength

Computing strength from particle-based simulations is not straightforward. A first issue is that various types of strength can be computed, e.g. tensile, compressive, shear, or indentation hardness, and the relationships between them are complex, when they exist at all. A second issue is that irreversible deformations accumulate progressively while transitioning from linear elasticity to full plasticity, and eventually to fracture (and not all these regimes necessarily emerge from a simulation). Therefore conventional decisions must be made when defining strength, e.g. it could be the stress at which the very first irreversible deformation occurs, or the absolute peak stress on a stress-strain curve, or other more convoluted definitions. A notable example is the indentation hardness, whose calculation requires constructing several Mohr

circles at failure (which itself requires a definition of failure) and finding the strength envelope that they describe, so that shear cohesion at zero confinement can be obtained [155]. In presenting the literature results below we will only indicate the type of strength being computed, hence providing just an intuitive appreciation of the data; the interested reader can find all the methodological details in the referenced works.

Among the various types of strength being computed, indentation hardness *H* has attracted most attention. This is because, like for the indentation modulus *M* in the previous section, also *H* can be directly compared with experimental results from nanoindentation, which also relate it to the local packing density *η* of the C–S–H. Ioannidou et al. [88] were first to report simulated *H* values: see Fig. 25.a. Liu et al. [91] obtained similar values, in Fig. 25.b, which was expected as they used the same interaction potential as in Ref. [88], although not the same packing algorithm. They also compared the results to simulated *H* values for crystalline structures made with similar nanoparticles, as shown in Fig. 25.b. Recently, Yaphary et al. [94] also simulated *H* values, using a potential that featured both a cohesive well and a repulsive shoulder, with similar bond energy as in the previous works. Their *H* values,

<figure><img src="image_28.png" /><figcaption>**Fig. 25.** Simulated indentation hardness *H* [90,91,94] as functions of C–S–H packing density *η* and, in some cases, comparison with experimental results from [90,154]. Images (b,c) reprinted from Ref. [91], Copyright (2019) and from Ref. [94], Copyright (2021), with permission from Elsevier.</figcaption></figure>
shown in Fig. 25.c, are in a similar range as previous ones, although their $H-\eta$ relationship is less representative of the experimental trend.

A few other simulated values of strength have been reported in the literature. Masoero et al. [86] computed strength under pure shear, using energy minimisation after increments of shear strain and keeping the box dimensions fixed (hence positive or negative pressures built up during the deformation). The obtained shear strength was between 80 and 320 MPa, increasing with packing densities from $\eta = 0.64$ to 0.71 (different $\eta$'s were obtained via random close packing with different particle size polydispersity). Recently Liu et al. [91] reported a single, smaller, simulated value of shear strength: ca. 50 MPa. Their simulations were similar in packing algorithm and potentials as those in Ref. [86], but the tests were simple shear, not pure shear.

Uniaxial tensile strength was first simulated by Davie and Masoero [156], using energy minimisation between increments of applied axial strain and keeping zero the axial stress perpendicular to the loading axis. The computed tensile strength was between 550 and 750 MPa, also in this case increasing from $\eta = 0.64$ to 0.76. Recently Hou et al. [93] computed tensile strength using a Grand Canonical Monte Carlo to create the C–S–H structures, followed by a peridynamic algorithm with interaction energies in the same order as in Ref. [156]. They computed tensile strengths for structures with $\eta \le 0.53$, but extrapolating their results to a $\eta = 0.64$ would lead to a strength of ca. 600 MPa, which is similar to the previous result in Ref. [156].

Compressive strength has been computed only recently by Yaphary et al. [94]. In the simulations, the axial stresses perpendicular to the loading directions were kept null by adapting the dimensions of the simulation box; without this feature, a simulation under compression would not be able to predict a transition to plasticity hence a compressive strength (differently from tensile simulations, where lack of transversal relaxation would still allow breaking the material and computing a corresponding strength). Their computed values are 280–550 MPa for C–S–H configurations that they associate to a low-density phase, whose $\eta$ is between 0.58 and 0.66. For high-density C–S–H structure, with $\eta = 0.72 – 0.75$, they computes strengths in the 500 – 900 MPa range. These values suggest a bigger impact of C–S–H packing density $\eta$ on compressive strength that on tensile strength (cf. previous paragraph). However, more results on compressive strength are needed before drawing strong conclusions.

### 5.2.4. Plasticity and fracture

The current literature on this property is very thin, with very few works having reported on simulated plastic deformations in some depth [86], and none having quantified fracture energy or toughness yet.

Masoero et al. [86] computed the accumulation of irreversible strain under pure shear, using the same simulations that led to the shear strengths in the previous section. Their study considered different levels of particle size polydispersity, inducing different packing densities $\eta$. An interesting result is that, upon large shear deformations, less dense structures ($\eta \le 0.7$) initially developed a negative pressure, i.e. they tended to densify. This tendency reversed with the accumulation of plastic strain, leading to a buildup of positive pressure indicating a tendency to dilate. By contrast, denser structures ($\eta \ge 0.7$) always tended to dilate, i.e. they developed a positive pressure that kept increasing as the plastic strain grew. Irreversible deformations were linked to non-affine displacements showing that, in polydisperse systems at large strain levels, local deformations would only develop in the matrix of smaller particles, whereas larger particles would behave as inclusions drifted by the average field of affine strain. Masoero et al. [86] argued that the energy to realise this matrix-inclusion behavior would control the activation of a fully developed plastic regime. Later simulations monitored the strain energy $\delta U$ to activate plastic deformations during quasi-static simulations of simple shear strain $\epsilon$; the computed $\Delta U$ were increasing logarithmically with $\epsilon$ until full plasticity was attained, at which point $\Delta U$ starbilled around 3 – 4 eV. [157,158].

Regarding fracture energy or toughness, one may wonder why those

have not been reported in the literature, since these could be easily obtained from published stress-strain curves under uniaxial tension (with or without applying a pre-existing notch). The difficulty here is that fracture propagation is inherently a multi-scale process, where governing defects are likely to exist at length scales greater than the hundreds of nanometers of particle-based simulations of C–S–H. As a result, correspondence between simulated and experimental results is unlikely to be retrieved (even for micro-scale experiments, such as scratch tests [159,160]). If correspondence be retrieved, one should then make sure that toughening and softening mechanisms at larger scales are considered in the validation [161]. A recent attempt to compute fracture properties from a PMF-based is in Refs. [162,163], where the lattice-element methods has been employed to address fracture of heterogeneous solids and has captured indeed some toughening mechanisms in two-phase layered composite materials. These challenges are still to be addressed in the field of mesoscale simulations of C–S–H.

### 5.2.5. Creep

Particle-based simulations of C–S–H have been used to predict long-term basic creep, viz. slow accumulation of irreversible strain $\epsilon$ under constant stress $\sigma$ and without concurring chemical reactions. The basic creep strain of concrete evolves logarithmically, $\epsilon = \sigma C \log(\frac{t}{t_0} + 1)$, where C and $t_0$ are two material parameters: creep modulus and characteristic time scale. C–S–H is believed to control the basic creep of concrete, and nanoindentation experiments have indeed shown that C–S–H enters a logarithmic creep regime already a few seconds after loading [164] (as opposed to macroscale creep tests on concrete, which take longer to establish a similar regime [165]). Simulating creep is challenging because the high energy barriers in hardened C–S–H (see Table 1 and Fig. 26) imply that irreversible strain accumulates over time scales that are too large for molecular dynamics.

A first approach to simulate creep was proposed by Masoero et al. [157]. They used the activation energies $\Delta U$ from the plastic strain simulations in the previous section, where each $\Delta U$ accompanied a finite increment of plastic strain $\Delta \epsilon$. The Arrhenius law was used to convert $\Delta U$ into a corresponding time increment $\Delta t = A \exp(\frac{\Delta U}{k_B T})$. The resulting evolution of $\epsilon(t)$ displayed indeed an initial logarithmic regime until full plasticity developed: see Fig. 27.a. An analytical explanation as to why $\Delta U \sim \log(\epsilon)$ entails logarithmic creep was provided later [158]. The result in Fig. 27.a, however, predicted an excessively large timescale originating from two features of the model: (i) the steric constraint imposed by not resolving molecular displacements below individual particles, and (ii) the association of creep-driving deformations, which are local with only few particles rearranging, with plastic deformations under stress, which instead involve an increasing number of particles eventually spanning the whole system. To remove the second issue, new simulations were carried out where irreversible deformations under constant shear stress were activated via small oscillations of shear stress or strain. Each oscillation entailed an approximately constant change in strain energy, which defined the magnitude of the energy barriers sampled by the simulation. Based on this, each cycle of shear oscillation was considered to be proportional to a constant time increment, so that the computed $\epsilon \sim \log(n)$ relationship (e.g. in Fig. 27.b), where n is the number of shear oscillations, effectively corresponded to logarithmic creep in time. The oscillatory approach was first proposed [158,166] and then fully deployed [102,167,168] for C–S–H simulations both at the molecular and at the mesoscale. A highlight of these simulations is their ability to predict the creep modulus C as a function of the packing density $\eta$ of the C–S–H, providing direct comparison with results from nanoindentation creep experiments [164]: see Fig. 27.c.

Recently, Masoero and Di Luzio [101] addressed another aspect of the basic creep of C–S–H: microprestress relaxation. In 1997, Bažant et al. [148] proposed a theoretical model that predicted logarithmic creep provided that the microprestresses in C–S–H, which are closely
<figure><img src="image_29.png" /><figcaption>**Fig. 26.** Simulated plastic behavior of C–S–H under shear. (a) Monodisperse C–S–H model structure; (b) stress-strain curve under simple shear; (c) evolution of energy barriers $\Delta U$ to activate irreversible deformation (plastic) events. Images reproduced from Ref. [158], with permission from ASCE).</figcaption></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>**Fig. 27.** Mesoscale simulation of C–S–H creep: (a) Creep curve constructed from activation energies during monotonic shear tests [157], such as in Fig. 26 above; (b) creep curves from the oscillatory approach, for C–S–H structures with different packing densities $\eta$ [102]; (c) creep modulus as a function of $\eta$ [102] from the creep curves in (b), and comparison with experimental results from nanoindentation [164]; (d) simulated power-law relaxation of eigenstresses during oscillatory creep simulations [101]. Permissions to reproduce images: (a) from Ref. [157] with permission from ASCE, (b,c) from Ref. [102], Copyright (2019), with permission from Elsevier, (d) from [101], Copyright (2020), with permission from Elsevier.</figcaption></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

related to the formation eigenstresses discussed previously, relax as a power law of time. Using oscillatory shear simulations and temperature-accelerated molecular dynamics on a C–S–H model with harmonic interactions as in Eq. (15), Masoero and Di Luzio found that indeed microprestresses relaxed as a power law of the number $n$ of activation attempts, and thus approximately of time: see Fig. 27.d. This confirmed Bažant et al.'s assumption, which was originally based purely on self-similarity arguments and on it being necessary for capturing logarithmic creep.

## 6. Summary and outlook

The development of mesoscale simulations and methods for materials is an active field of research that has significant challenges and

practical applications. In this manuscript we have provided an overview of computational models and methods used in mesoscale simulations of cementitious materials, focusing on C–S–H, which is the main cement hydrate. We have illustrated with several examples how mesoscale simulations can bridge the gap of length and time scales between molecular simulations and continuous modelling, to predict macroscopic (thermodynamic, kinetic and mechanical) properties and behaviours from fundamental molecular processes using a rigorous statistical mechanics treatment. The chosen examples ranged from adsorption isotherms of simple ions and complex polymers up to the growth, structure and mechanical properties of C–S–H gels through liquid and solid phase composition (speciation) and interparticle interactions, always comparing the results with experimental observations.

Within the statistical mechanics framework we showed how to
develop a hierarchy of coarse grained models where models at different levels of coarse graining are related via the calculation of the interaction free energy, that is the pair potential of mean force (PMF), between the chemical species or nano-particles of interests. Using such a method we showed, for example, that the many details of a chemical reaction (electronic, kinetic pathway etc...) can be reduced to an equilibrium constant, linked to the Gibbs free energy. When embedded into an appropriate Monte Carlo move and thermodynamic ensemble, the speciation of liquid and solid phases can be studied while keeping a microscopic (albeit coarse grained) description of the main physical and chemical processes. This was illustrated in the simple CaO-SiO₂-H₂O system, but with some efforts can be extended to the prediction of the phase assemblage and pore solution composition of a cement paste. This approach may support the development of low clinker cementitious materials, which increasingly rely on organic admixtures (e.g. super-plasticizers) and, for the same reason, are difficult to treat with classical speciation models.

Mesoscopic modelling can also help select and develop organic additives. It can be used to predict the affinity of the organic molecules with the different phases of the cement (anhydrous and hydrates) according to their structure and functionality, but also to predict the changes in inter-particle interactions induced by the organics adsorption. This research field is still in its infancy, especially for what concerns the rationalization of the dissolution processes of anhydrous phases and nucleation/growth of hydrates, while they interact with organic additives. This is indeed a field that can significantly benefit from mesoscopic simulations such as those reviewed in this manuscript.

The hierarchical coarse graining approach reviewed in this manuscript, which relies on PMF calculations, can be powerful but some aspects require caution. An example was illustrated here with regard to the interaction potential between C-S-H particles. Experimental measurements from atomic force microscopy and simulations in the framework of the primitive model describe similar interaction potentials, whereas atomistic simulations employing empirical force fields predict interaction free energies that are one to two orders of magnitude greater. On one hand, this shows a limitation of the primitive model, which does not capture the many-body effects that strengthen the interactions in the extreme confinement of the inter-lamellar spacing of C-S-H crystals. On the other hand, atomistic simulations do not capture the colloidal behavior of C-S-H particles at larger distances, nor there is an available method for them to predict the changes in interaction potentials in response to a change in pH or ion concentration in the cement pore solution. Moreover, the PMFs depend not only on the geometry and charge of the particles but also on the density of the cement paste. Therefore, the tools must then be chosen and used with care according to their limitations and to the properties that one wants to model.

At the coarser mesoscale of hundreds of nanometers, we have shown how a few simulations have started to incorporate PMFs coming from smaller scales, but still most of the existing particle-based simulations use empirical potential which are built rationally to target certain properties. Such empirical potentials miss some aspects of atomistic processes; this is evidenced, for example, by the practice of using certain effective potentials to represent the early stages of C-S-H formation, and others to study its mature behaviours. Empirical potentials capturing both regimes of colloidal and crystalline interactions are still to be developed and could stimulate interesting future studies. Despite such limitations, this manuscript has shown how particle-based mesoscale simulations to date have already successfully captured a remarkable range of C-S-H gel properties, both structural and mechanical ones. However, there is still room for significant new developments, such as:

* Modelling particles that can partially grow and dissolve. The existing simulations model growth as secondary nucleation and aggregation, but C-S-H often forms elongated crystals, such as foils and needles [169]. Aggregation of anisotropic nanoparticles could lead to such shapes, but a more efficient way to simulate the formation of such

morphologies would be to let particles grow and dissolve partially. This would also enable simulations of processes where small changes in particle sizes might enable significant structural rearrangements or stress relaxation, such as in creep or crystallization pressure.

* Detailed account of particle nucleation and growth mechanics. In the current simulations, particles usually appear or disappear without specifying the details of the underlying mechanisms. Such details are not important in simulations where the target is the thermodynamic limit, but they are important for the kinetics. Shvab et al. [98] have assumed a simple classical growth mechanism in their Kinetic Monte Carlo simulations of C-S-H formation. Nonclassical mechanisms may better describe the formation of C-S-H [104,105] but are still to be included in mesoscale simulations.

* Coupling particle precipitation and dissolution with ion diffusion in solution. Most of the current particle-based simulations of C-S-H formation assume uniform composition of the solution surrounding the C-S-H particles. Recently, Goyal et al. considered a gradient of chemical potential to model a solution that is increasingly saturated as a surface of unhydrated cement is approached [89]. Modelling also diffusive processes in solution might lead to capturing concentration gradients near the surfaces of individual C-S-H particles, where local ionic concentration can significantly affect the strength of the interactions. Furthermore, limitation of reactions by diffusion might become important during C-S-H formation, especially as the packing density of the solid increases. Incorporation of diffusion might be through fine grid calculations, e.g. using Lattice Boltzmann [170], but for the sake of computational efficiency one might instead use a coarser grid and consider the mesostructure of C-S-H in each grid cell as an effective porous medium (e.g. as done in Ref. [171], where the porous solid is a film of agglomerated bacteria).

* Addressing long timescales. To date, particle-based simulations have used heuristic methods to explore the long-term behavior of the materials: see e.g. the simulations of creep reviewed in this manuscript. The challenge of simulating long timescales at low length scales is common across computational material science. However, state-of-the-art activation-relaxation techniques are not yet much applied to C-S-H, for example metadynamics [172], ABC [173], or ART [174]. ART has been recently used to access long term diffusion of cesium ions in the C-S-H matrix and combined with mesoscale simulations provided useful information for nuclear waste repositories [175]. In another multiscale assessment of alkali-silica reaction, hours long diffusion times of alkali in the C-S-H interlayer were reported using ART [176]. These techniques can be very computationally expensive and, to date, there is no definite way of ensuring that all the important events underlying the slow evolution of the system will be captured. However, incorporating such techniques is worth exploring, as they could improve our understanding of slow processes such as those underlying creep or the rheology of fresh cement pastes.

* Simulating cement hydrates other than C-S-H. Undoubtedly, C-S-H is a key phase in cement paste and concrete. However, various degradation mechanisms involve other phases too, if not primarily; for example, calcium hydroxide is key to carbonation, whereas sulfoaluminates underlie delayed ettringite formation. Furthermore, new phases may largely supplant C-S-H in future, more sustainable concrete, e.g. geopolymers. Including such phases in mesoscale simulations would require detailed understanding of relevant fundamental processes at the atomic scale, as well as coarse-graining through the mesoscales, as shown here for C-S-H. Appropriate interaction potentials will be needed too, which can either be derived rigorously as PMFs from smaller-scale simulations, or developed empirically but rationally directly at the particle scale: such an attempt for interfaces between C-S-H and fly ash can be found in Ref. [177]. The approaches and results reviewed in this manuscript will hopefully provide good guidance to the researchers who want to contribute to this exciting field.
CRediT authorship contribution statement

**Katerina Ioannidou:** Conceptualization, Data curation, Formal analysis, Methodology, Project Administration, Visualization, Writing – original draft, Writing – review & editing.
**Christophe Labbez:** Conceptualization, Data curation, Formal analysis, Methodology, Project administration, Visualization, Writing – original draft, Writing – review & editing.
**Enrico Masoero:** Conceptualization, Data curation, Formal analysis, Methodology, Project Administration, Visualization, Writing – original draft, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cemconres.2022.106857.

References

[1] G. Habert, S. Miller, V. John, J. Provis, A. Favier, A. Horvath, K. Scrivener, Environmental impacts and decarbonization strategies in the cement and concrete industries, Nat. Rev. Earth Environ. 1 (11) (2020) 559–573.

[2] U. Environment, K.L. Scrivener, V.M. John, E.M. Gartner, Eco-efficient cements: potential economically viable solutions for a low-co2 cement-based materials industry, Cem. Concr. Res. 114 (2018) 2–26.

[3] J.J. Biernacki, J.W. Bullard, G. Sant, K. Brown, F.P. Glasser, S. Jones, T. Ley, R. Livingston, L. Nicoleau, J. Olek, et al., Cements in the 21st century: challenges, perspectives, and opportunities, J. Am. Ceram. Soc. 100 (7) (2017) 2746–2773.

[4] W. Andreoni, S. Yip, Handbook of Materials Modeling: Applications: Current and Emerging Materials, Springer, 2020.

[5] R.K. Mishra, A.K. Mohamed, D. Geissbühler, H. Manzano, T. Jamil, R. Shahsavari, A.G. Kalinichev, S. Galmarini, L. Tao, H. Heinz, et al., Cemff: a force field database for cementitious materials including validations, applications and opportunities, Cem. Concr. Res. 102 (2017) 68–89.

[6] N. Kringos, B. Birgisson, D. Frost, L. Wang, Multi-scale modeling and characterization of infrastructure materials, in: Proceedings of the International RILEM Symposium: Estocolmo, Springer, Springer, 2013.

[7] P. Thilakarathna, K.K. Baduge, P. Mendis, V. Vimosatit, H. Lee, Mesoscale modelling of concrete—a review of geometry generation, placing algorithms, constitutive relations and applications, Eng. Fract. Mech. 231 (2020), 106974.

[8] J.J. Thomas, J.J. Biernacki, J.W. Bullard, S. Bishnoi, J.S. Dolado, G.W. Scherer, A. Luttge, Modeling and simulation of cement hydration kinetics and microstructure development, Cem. Concr. Res. 41 (12) (2011) 1257–1278.

[9] E. Del Gado, K. Ioannidou, E. Masoero, A. Baronnet, R.J.M. Pellenq, F.-J. Ulm, S. Yip, A soft matter in construction – statistical physics approach to formation and mechanics of C-S-H gels in cement, Eur. Phys. J. Spec. Top. 223 (2014) 2285–2295.

[10] E.Del Gado, K. Martens, R.J.M. Pellenq, in: From Microscopic Insight to Constitutive Models: Bridging Length Scales in Soft and Hard Materials, Springer International Publishing, Cham, 2020, pp. 1–19.

[11] K. Ioannidou, B. Carrier, M. Vandamme, R. Pellenq, The potential of mean force concept for bridging (length and time) scales in the modeling of complex porous materials, in: EPJ Web of Conferences Vol. 140, EDP Sciences, 2017, p. 01009.

[12] D. Gawin, F. Pesavento, B. Schrefler, Simulation of damage–permeability coupling in hygro-thermo-mechanical analysis of concrete at high temperature, Commun. Numer. Methods Eng. 18 (2) (2002) 113–119.

[13] C.T. Davie, C.J. Pearce, N. Bičanić, A fully generalised, coupled, multi-phase, hygro-thermo-mechanical model for concrete, Mater. Struct. 43 (1) (2010) 13–33.

[14] M. Pathirage, D. Bentz, G. Di Luzio, E. Masoero, G. Cusatis, The onix model: a parameter-free multiscale framework for the prediction of self-desiccation in concrete, Cem. Concr. Compos. 103 (2019) 36–48.

[15] D. Rodney, A. Tanguy, D. Vandembroucq, Modeling the mechanics of amorphous solids at different length scale and time scale, Model. Simul. Mater. Sci. Eng. 19 (8) (2011), 083001.

[16] M. Medala, C. Labbez, I. Pochard, A. Nonat, Ettringite surface chemistry: interplay of electrostatic and ion specificity, J. Colloid Interface Sci. 354 (2) (2011) 765–770, https://doi.org/10.1016/j.cis.2010.11.031.

[17] C. Labbez, B. Jönsson, I. Pochard, A. Nonat, B. Cabane, Surface charge density and electrokinetic potential of highly charged minerals: experiments and Monte Carlo simulations on calcium silicate hydrate, J. Phys. Chem. B 110 (18) (2006) 9219–9230, https://doi.org/10.1021/jp057096+.

[18] G. Gouy, Sur la fonction électrocapillaire, Ann. Chim. Phys. 29 (7) (1903) 145.

[19] D.L. Chapman, LI. A contribution to the theory of electrocapillarity, Lond. Edinb. Dublin Philos. Mag. J. Sci. 25 (148) (1913) 475–481, https://doi.org/10.1080/14786440408634187.

[20] D.C. Grahame, Properties of the electrical double layer at a mercury surface. I. Methods of measurement and interpretation of results, J. Am. Chem. Soc. 63 (5) (1941) 1207–1215, https://doi.org/10.1021/ja01850a014.

[21] J.N. Israelachvili, Intermolecular and Surface Forces, 2nd Ed., Academic Press, London, 1991.

[22] B. Jönsson, H. Wennerström, When ion-ion correlations are important in charged colloidal systems, in: C. Holm, P. Kekicheff, R. Podgornik (Eds.), Electrostatic Effects in Soft Matter and Biophysics, Kluwer Academic Publishers, 2001.

[23] G.M. Torrie, J.P. Valleau, Electrical double layers. 4. Limitations of the Gouy-Chapman theory, J. Phys. Chem. 86 (16) (1982) 3251–3257, https://doi.org/10.1021/j100213a035.

[24] L. Guldbrand, B. Jönsson, H. Wennerström, P. Linse, Electrical double layer forces. A Monte Carlo study, J. Chem. Phys. 80 (5) (1984) 2221–2228, https://doi.org/10.1063/1.446912.

[25] C. Labbez, B. Jonsson, M. Skarba, M. Borkovec, Ion- ion correlation and charge reversal at titrating solid interfaces, Langmuir 25 (13) (2009) 7209–7213.

[26] P.M. Dove, C.M. Craven, Surface charge density on silica in alkali and alkaline earth chloride electrolyte solutions, Geochim. Cosmochim. Acta 69 (2005) 4963–4970.

[27] C. Labbez, I. Pochard, B. Jönsson, A. Nonat, C-S-H/solution interface: experimental and Monte Carlo studies, Cem. Concr. Res. 41 (2) (2011) 161–168, https://doi.org/10.1016/j.cemconres.2010.10.002.

[28] J. Haas, A. Nonat, From C-S-H to C-A-S-H: experimental study and thermodynamic modelling, Cem. Concr. Res. 68 (2015) 124–138, https://doi.org/10.1016/j.cemconres.2014.10.020.

[29] I. Pochard, C. Labbez, A. Nonat, H. Vija, B. Jönsson, The effect of polycations on early cement paste, Cem. Concr. Res. 40 (10) (2010) 1488–1494, https://doi.org/10.1016/j.cemconres.2010.06.002.

[30] F. Brunel, I. Pochard, M. Turesson, S. Gauffinet, C. Labbez, Elastic response of cementitious gels to polycation addition, ACS Omega 2 (5) (2017) 2148–2158, https://doi.org/10.1021/acsomega.6b00445.

[31] E. Bernard, Y. Yan, B. Lothenbach, Effective cation exchange capacity of calcium silicate hydrates (CSH), Cem. Concr. Res. 143 (2021) 106393.

[32] C. Labbez, M. Medala, I. Pochard, A. Nonat, Adsorption of multivalent ions in cementitious materials : importance of electrostatics, URL, in: 2nd International Workshop on Mechanisms and Modelling of Waste/cement Interactions, Le Croisic, France, 2008, https://hal.archives-ouvertes.fr/hal-00415128.

[33] T. Matschei, R. Skapa, B. Lothenbach, F. Glasser, 12th International Congress on the Chemistry of Cement, Vol. 3008, Montréal, Canada, 2007. https://www.dora.libri.ch/empa/islandora/object/empa.

[34] H. Viallis, P. Faucon, J.-C. Petit, A. Nonat, Interaction between salts (NaCl, CsCl) and calcium silicate hydrates (C-S-H), J. Phys. Chem. B 103 (25) (1999) 5212–5219, june 24, 1999XC 16.

[35] M. Turesson, C. Labbez, A. Nonat, Calcium mediated polyelectrolyte adsorption on like-charged surfaces, Langmuir 27 (22) (2011) 13572–13581, https://doi.org/10.1021/la2030846.

[36] A.K. Mohamed, S.A. Weckwerth, R.K. Mishra, H. Heinz, R.J. Flatt, Molecular modeling of chemical admixtures; opportunities and challenges, Cem. Concr. Res. 156 (2022), 106783.

[37] M. Turesson, A. Nonat, C. Labbez, Stability of negatively charged platelets in calcium-rich anionic copolymer solutions, Langmuir 30 (23) (2014) 6713–6720, https://doi.org/10.1021/la501228w.

[38] C. Nalet, A. Nonat, Ionic complexation and adsorption of small organic molecules on calcium silicate hydrate: relation with their retarding effect on the hydration of C3S, Cem. Concr. Res. 89 (2016) 97–108, https://doi.org/10.1016/j.cemconres.2016.08.012.

[39] L. Dupont, A. Foissy, R. Mercier, B. Mottet, Effect of calcium ions on the adsorption of polyacrylic acid onto alumina, J. Colloid Interface Sci. 161 (2) (1993) 455–464, https://doi.org/10.1006/jcis.1993.1489.

[40] A.K. Bajpai, Interface behaviour of ionic polymers, Prog. Polym. Sci. 22 (3) (1997) 523–564, https://doi.org/10.1016/S0079-6700(96)00003-2.

[41] D.N. Misra, Adsorption of low molecular weight poly(acrylic acid) on hydroxyapatite: role of molecular association and apatite dissolution, Langmuir 7 (11) (1991) 2422–2424, https://doi.org/10.1021/la00059a003.

[42] J. Sun, L. Bergstroöm, L. Gao, Effect of magnesium ions on the adsorption of poly (acrylic acid) onto alumina, J. Am. Ceram. Soc. 84 (11) (2001) 2710–2712, https://doi.org/10.1111/j.1151-2916.2001.tb01078.x.

[43] J.M. Lamarche, J. Persello, A. Foissy, Influence of molecular weight of sodium polyacrylate in calcium carbonate aqueous dispersions, Ind. Eng. Chem. Prod. Res. Dev. 22 (1) (1983) 123–126, https://doi.org/10.1021/i300009a028.

[44] A. Foissy, A. El Attar, J.M. Lamarche, Adsorption of polyacrylic acid on titanium dioxide, J. Colloid Interface Sci. 96 (1) (1983) 275–287, https://doi.org/10.1016/0021-9797(83)90029-2.

[45] L. Jaärnstööm, P. Stenius, Adsorption of polyacrylate and carboxy methyl cellulose on kaolinite: salt effects and competitive adsorption, Colloids Surf. 50 (1990) 47–73, https://doi.org/10.1016/0166-6622(90)80253-7.

[46] L. Bouzouaid, B. Lothenbach, A. Fernandez-Martinez, C. Labbez, Gluconate and Hexitols Effects on C-S-H Solubility, arXiv:2111.05152. URL, arXiv, Nov 2021, https://arxiv.org/abs/2111.05152v1.

[47] M. Grutzeck, A. Benesi, B. Fanning, Silicon29 magic angle spinning nuclear magnetic resonance study of calcium silicate hydrates, J. Am. Ceram. Soc. 72 (4) (1989) 665–668.
[48] X. Cong, R.J. Kirkpatrick, 29Si MAS NMR study of the structure of calcium silicate hydrate, Adv. Cem. Mater. 3 (3-4) (1996) 144-156, 1996/0.

[49] I. Klur, B. Pollet, J. Virlet, A. Nonat, C-S-H structure evolution with calcium content by multinuclear NMR, in: P. Colombet, H. Zanni, P. Sozzani, A.-R. Grimmer (Eds.), Nuclear Magnetic Resonance Spectroscopy of Cement-Based Materials, Springer, Berlin, 1998, pp. 119-141, fG 27.

[50] J.J. Chen, J.J. Thomas, H.F.W. Taylor, H.M. Jennings, Solubility and structure of calcium silicate hydrate, Cem. Concr. Res. 34 (9) (2004) 1499-1519.

[51] S.V. Churakov, C. Labbez, Thermodynamics and molecular mechanism of Al incorporation in calcium silicate hydrates, J. Phys. Chem. C 121 (8) (2017) 4412-4419, https://doi.org/10.1021/acs.jpc6.b12850.

[52] Ceérd Plassard, Eric Lesniewska, Isabelle Pochard, Andrée Nonat, Nanoscale experimental investigation of particle interactions at the origin of the cohesion of cement, Langmuir 2005 (21) (2005) 16, https://doi.org/10.1021/la050440+.

[53] S. Garrault-Gauffinet, A. Nonat, Experimental investigation of calcium silicate hydrate (csh) nucleation, J. Cryst. Growth 200 (3-4) (1999) 565-574.

[54] B. Joönnsson, H. Wennerstroöm, A. Nonat, B. Cabane, Onset of cohesion in cement paste, Langmuir 20 (16) (2004) 6702-6709, https://doi.org/10.1021/la0498760.

[55] B. Joönnsson, A. Nonat, C. Labbez, B. Cabane, H. Wennerstroöm, Controlling the cohesion of cement paste, Langmuir 21 (20) (2005) 9211-9221, https://doi.org/10.1021/la051048z.

[56] S. Masoumi, D. Ebrahimi, H. Valipour, M.J.Abdolhosseini Qomi, Nanolayered attributes of calcium-silicate-hydrate gels, J. Am. Ceram. Soc. 103 (1) (2020) 541-557.

[57] H.K. Christenson, Adhesion and surface energy of mica in air and water, J. Phys. Chem. 97 (46) (1993) 12034-12041, https://doi.org/10.1021/j100148a032.

[58] R.J.-M. Pellenq, H. Van Damme, Why does concrete Set?: the nature of cohesion forces in hardened cement-based materials, MRS Bull. 29 (5) (2004) 319-323, https://doi.org/10.1557/mrs2004.97.

[59] R.J.-M. Pellenq, J.M. Caillol, A. Delville, Electrostatic attraction between two charged surfaces: a (N, V, T) Monte Carlo simulation, J. Phys. Chem. B 101 (42) (1997) 8584-8594, https://doi.org/10.1021/jp971273s.

[60] L. Guldbrand, B. Joönnsson, H. Wennerstroöm, P. Linse, Electrical double layer forces. A Monte Carlo study, J. Chem. Phys. 80 (5) (1984) 2221-2228, https://doi.org/10.1063/1.446912.

[61] R. Kjellander, S. Marcełja, Correlation and image charge effects in electric double layers, Chem. Phys. Lett. 112 (1) (1984) 49-53, https://doi.org/10.1016/0009-2614(84)87039-6.

[62] B. Joönnsson, H. Wennerstroöm, ION-ION correlations in liquid dispersions, J. Adhes. 80 (5) (2004) 339-364, https://doi.org/10.1080/00218460490465551.

[63] B. Joönnsson, H. Wennerstroöm, B. Halle, Ion distributions in lamellar liquid crystals. A comparison between results from Monte Carlo simulations and solutions of the Poisson-Boltzmann equation, J. Phys. Chem. 84 (17) (1980) 2179-2185, https://doi.org/10.1021/j100454a014.

[64] A. Goyal, I. Palaia, K. Ioannidou, F.-J. Ulm, H. van Damme, R.J.-M. Pellenq, E. Trizac, E.Del Gado, The physics of cement cohesion, Sci. Adv. 7 (32) (2021), https://doi.org/10.1126/sciadv.abg5882.

[65] A. Picker, L. Nicoleau, Z. Burghard, J. Bill, I. Zlotnikov, C. Labbez, A. Nonat, H. Coölfen, Mesocrystalline calcium silicate hydrate: a bioinspired route toward elastic concrete materials, Sci. Adv. 3 (11) (Nov 2017) e1701216. https://www.science.org/doi/10.1126/sciadv.1701216.

[66] S. Alexander, P.M. Chaikin, P. Grant, G.J. Morales, P. Pincus, D. Hone, Charge renormalization, osmotic pressure, and bulk modulus of colloidal crystals: theory, J. Chem. Phys. 80 (11) (1984) 5776-5781, https://doi.org/10.1063/1.446600.

[67] A.A. Louis, Beware of density dependent pair potentials, J. Phys.: Condens. Matter 14 (40) (2002) 9187-9206, https://doi.org/10.1088/0953-8984/14/40/311.

[68] M. Turesson, B. Joönnsson, C. Labbez, Coarse-graining intermolecular interactions in dispersions of highly charged colloids, Langmuir 28 (11) (2012) 4926-4930, https://doi.org/10.1021/la3005008.

[69] G. Bareigts, C. Labbez, Effective pair potential between charged nanoparticles at high volume fractions, Phys. Chem. Chem. Phys. 19 (6) (2017) 4787-4792, https://doi.org/10.1039/C6CP08056A.

[70] A. Thresson, M. Ullner, T. Åkesson, C. Labbez, B. Joönnsson, Monte Carlo simulations of parallel charged platelets as an approach to tactoid formation in clay, Langmuir 29 (29) (2013) 9216-9223, https://doi.org/10.1021/la401272u.

[71] A. Thresson, M. Ullner, M. Turesson, Interaction and aggregation of charged platelets in electrolyte solutions: a coarse-graining approach, J. Phys. Chem. B 118 (26) (2014) 7405-7413, https://doi.org/10.1021/jp502015g.

[72] M. Delhorme, C. Labbez, M. Turesson, E. Lesniewska, C.E. Woodward, B. Joönnsson, Aggregation of calcium silicate hydrate nanoplatelets, Langmuir 32 (8) (2016) 2058-2066, https://doi.org/10.1021/acs.langmuir.5b03846.

[73] T.L. Hill, An Introduction to Statistical Thermodynamics, Dover Publications Inc., New York, 1986.

[74] P.A. Bonnaud, C. Labbez, R. Miura, A. Suzuki, N. Miyamoto, N. Hatakeyama, A. Miyamoto, K.J. Van Vliet, Interaction grand potential between calcium-silicate-hydrate nanoparticles at the molecular level, Nanoscale 8 (7) (2016) 4160-4172.

[75] S. Masoumi, H. Valipour, M.J.Abdolhosseini Qomi, Intermolecular forces between nanolayers of crystalline calcium-silicate-hydrates in aqueous medium, J. Phys. Chem. C 121 (10) (2017) 5565-5572.

[76] S. Masoumi, H. Valipour, M.J.Abdolhosseini Qomi, Interparticle interactions in colloidal systems: toward a comprehensive mesoscale model, ACS Appl. Mater. Interfaces 9 (32) (2017) 27338-27349.

[77] H. Manzano, E. Masoero, I. Lopez-Arbela, H.M. Jennings, Shear deformations in calcium silicate hydrates, Soft Matter 9 (30) (2013) 7333-7341.

[78] M. Kai, L. Zhang, K. Liew, New insights into creep characteristics of calcium silicate hydrates at molecular level, Cem. Concr. Res. 142 (2012) 106366.

[79] E. Duque-Redondo, E. Masoero, H. Manzano, Nanoscale shear cohesion between cement hydrates: the role of water diffusivity under structural and electrostatic confinement, Cem. Concr. Res. 154 (2022), 106716.

[80] S. Masoumi, S. Zare, H. Valipour, M.J.Abdolhosseini Qomi, Effective interactions between calcium-silicate-hydrate nanolayers, J. Phys. Chem. C 123 (8) (2019) 4755-4766.

[81] T. Honorio, Monte Carlo molecular modeling of temperature and pressure effects on the interactions between crystalline calcium silicate hydrate layers, Langmuir 35 (11) (2019) 3907-3916.

[82] E. Masoero, E. Del Gado, R.-M. Pellenq, F.-J. Ulm, S. Yip, Nanostructure and nanomechanics of cement: polydisperse colloidal packing, Phys. Rev. Lett. 109 (15) (2012), 155503.

[83] G. Mie, Zur kinetischen theorie der einatomigen körper, Ann. Phys. 316 (8) (1903) 657-697.

[84] E. Masoero, H. Jennings, F. Ulm, E. Del Gado, H. Manzano, R. Pellenq, S. Yip, Modelling cement at fundamental scales: from atoms to engineering strength and durability, Comput. Model. Concr. Struct. 1 (2014) 139-148.

[85] K. Ioannidou, R.J.-M. Pellenq, E. Del Gado, Controlling local packing and growth in calcium-silicate-hydrate gels, Soft Matter 10 (8) (2014) 1121-1133.

[86] E. Masoero, E. Del Gado, R.J.-M. Pellenq, S. Yip, F.-J. Ulm, Nano-scale mechanics of colloidal c-s-h gels, Soft Matter 10 (3) (2014) 491-499.

[87] C. Plassard, E. Lesniewska, I. Pochard, A. Nonat, Nanoscale experimental investigation of particle interactions at the origin of the cohesion of cement, Langmuir 21 (16) (2005) 7263-7270.

[88] K. Ioannidou, M. Kandúc, L. Li, D. Frenkel, J. Dobnikar, E. Del Gado, The crucial effect of early-stage gelation on the mechanical properties of cement hydrates, Nat. Commun. 7 (1) (2016) 1-9.

[89] A. Goyal, K. Ioannidou, C. Tiede, P. Levitz, R.J.-M. Pellenq, E.Del Gado, Heterogeneous surface growth and gelation of cement hydrates, J. Phys. Chem. C 124 (28) (2020) 15500-15510, https://doi.org/10.1021/acs.jpc6.0c02944, arXiv:doi:10.1021/acs.jpc6.0c02944, doi:10.1021/acs.jpc6.0c02944, URL.

[90] K. Ioannidou, K.J. Krakowiak, M. Bauchy, C.G. Hoover, E. Masoero, S. Yip, F.-J. Ulm, P. Levitz, R.J.-M. Pellenq, E.Del Gado, Mesoscale texture of cement hydrates, Proc. Natl. Acad. Sci. 113 (8) (2016) 2029-2034.

[91] H. Liu, S. Dong, L. Tang, N.A. Krishnan, G. Sant, M. Bauchy, Effects of polydispersity and disorder on the mechanical properties of hydrated silicate gels, J. Mech. Phys. Solids 122 (2019) 555-565.

[92] K. Ioannidou, E. Del Gado, F.-J. Ulm, R.J.-M. Pellenq, Inhomogeneity in cement hydrates: linking local packing to local pressure, J. Nanomech. Micromech. 7 (2) (2017), 04017003.

[93] D. Hou, W. Zhang, P. Wang, M. Wang, H. Zhang, Mesoscale insights on the structure, mechanical performances and the damage process of calcium-silicate-hydrate, Constr. Build. Mater. 287 (2021), 123031.

[94] Y.L. Yaphary, F. Sanchez, D. Lau, C.S. Poon, Mechanical properties of colloidal calcium-silicate-hydrate gel with different gel-pore ionic solutions: a mesoscale study, Microporous Mesoporous Mater. 110944 (2021).

[95] A. Yadav, N.M.A. Krishnan, Role of steric repulsions on the precipitation kinetics and the structure of calcium-silicate-hydrate gels, Soft Matter 17 (2021) 8902-8914.

[96] Z. Yu, D. Lau, Nano-and mesoscale modeling of cement matrix, Nanoscale Res. Lett. 10 (1) (2015) 1-6.

[97] J. Gay, B. Berne, Modification of the overlap potential to mimic a linear site-site potential, J. Chem. Phys. 74 (6) (1981) 3316-3319.

[98] I. Shvab, L. Brochard, H. Manzano, E. Masoero, Precipitation mechanisms of mesoporous nanoparticle aggregates: off-lattice, coarse-grained, kinetic simulations, Cryst. Growth Des. 17 (3) (2017) 1316-1327.

[99] J.W. Bullard, G.W. Scherer, J.J. Thomas, Time dependent driving forces and the kinetics of tricalcium silicate hydration, Cem. Concr. Res. 74 (2015) 26-34.

[100] K. Coopamootoo, E. Masoero, Simulations of crystal dissolution using interacting particles: prediction of stress evolution and rates at defects and application to tricalcium silicate, J. Phys. Chem. C 124 (36) (2020) 19603-19615.

[101] E. Masoero, G. Di Luzio, Nanoparticle simulations of logarithmic creep and micropreedom relaxation in concrete and other disordered solids, Cem. Concr. Res. 137 (2020), 106181.

[102] H. Liu, S. Dong, L. Tang, N.A. Krishnan, E. Masoero, G. Sant, M. Bauchy, Long-term creep deformations in colloidal calcium-silicate-hydrate gels by accelerated aging simulations, J. Colloid Interface Sci. 542 (2019) 339-346.

[103] P.G. Vekilov, Nonclassical nucleation, in: Crystallization via Nonclassical Pathways Volume 1: Nucleation, Assembly, Observation & Application, ACS Publications, 2020, pp. 19-46.

[104] N. Krautwurst, L. Nicoleau, M. Dietzsch, I. Lieberwirth, C. Labbez, A. Fernandez-Martinez, A.E. Van Driessche, B. Barton, S. Leukel, W. Tremel, Two-step nucleation process of calcium silicate hydrate, the nanobrick of cement, Chem. Mater. 30 (9) (2018) 2895-2904.

[105] S. Bae, M. Kanematsu, D. Hernández-Cruz, J. Moon, D. Kilcoyne, P.J.M. Monteiro, In situ soft x-ray spectromicroscopy of early tricalcium silicate hydration, Materials 9 (12) (2016), https://doi.org/10.3390/ma9120976. https://www.mdpi.com/1996-1944/9/12/976.

[106] J.J. Thomas, A new approach to modeling the growth and growth kinetics of tricalcium silicate hydration, J. Am. Ceram. Soc. 90 (10) (2007) 3282-3288.

[107] G.W. Scherer, J. Zhang, J.J. Thomas, Nucleation and growth models for hydration of cement, Cem. Concr. Res. 42 (7) (2012) 982-993.

[108] C. Labbez, B. Jönsson, C. Woodward, A. Nonat, M. Delhorme, The growth of charged platelets, Phys. Chem. Chem. Phys. 16 (43) (2014) 23800-23808.
[109] P. Juilland, L. Nicoleau, R.S. Arvidson, E. Gallucci, Advances in dissolution understanding and their implications for cement hydration, RILEM Tech. Lett. 2 (2017) 90-98.

[110] K. Scrivener, A. Ouzia, P. Juilland, A.K. Mohamed, Advances in understanding cement hydration mechanisms, Cem. Concr. Res. 124 (2019), 105823.

[111] J.W. Bullard, H.M. Jennings, R.A. Livingston, A. Nonat, G.W. Scherer, J. S. Schweitzer, K.L. Scrivener, J.J. Thomas, Mechanisms of cement hydration, Cem. Concr. Res. 41 (12) (2011) 1208-1223.

[112] R. Gonzalez-Teresa, V. Morales-Florez, H. Manzano, J.S. Dolado, Structural models of randomly packed tobermorite-like spherical particles: a simple computational approach, Mater. Constr. 60 (298) (2010) 7-15.

[113] A. Prabhu, J.-C. Gimel, A. Ayuela, S. Arrese-Igor, J.J. Gaitero, J.S. Dolado, A multi-scale approach for percolation transition and its application to cement setting, Sci. Rep. 8 (1) (2018) 1-11.

[114] R. Gonzalez-Teresa, J.S. Dolado, A. Ayuela, J.-C. Gimel, Nanoscale texture development of csh gel: a computational model for nucleation and growth, Appl. Phys. Lett. 103 (23) (2013), 234105.

[115] M.A. Etzold, P.J. McDonald, A.R. Fouth, Growth of sheets in 3d confinements—a model for the c-s-h meso structure, Cem. Concr. Res. 63 (2014) 137-142.

[116] T. Petersen, P.-L. Valdenaire, R. Pellenq, F.-J. Ulm, A reaction model for cement solidification: evolving the c-s-h packing density at the micrometer-scale, J. Mech. Phys. Solids 118 (2018) 58-73.

[117] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1) (1984) 511-519, https://doi.org/10.1063/1.447334, arXiv:10.1063/1.447334, doi:10.1063/1.447334. URL.

[118] W.G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A 31 (1985) 1695-1697, https://doi.org/10.1103/PhysRevA.31.1695. https://link.aps.org/doi/10.1103/PhysRevA.31.1695.

[119] T. Zhou, K. Ioannidou, F.-J. Ulm, M.Z. Bazant, R.J.-M. Pellenq, Multiscale poromechanics of wet cement paste, URL, Proc. Natl. Acad. Sci. 116 (22) (2019) 10652-10657, https://doi.org/10.1073/pnas.1901160116, https://www.pnas. org/content/116/22/10652.

[120] S. Garrault, A. Nonat, Hydrated layer formation on tricalcium and dicalcium silicate surfaces: experimental study and numerical simulations, Langmuir 17 (26) (2001) 8131-8138.

[121] P. Mutabaruka, M. Taiebat, R.J.-M. Pellenq, F. Radjai, Effects of size polydispersity on random close-packed configurations of spherical particles, Phys. Rev. E 100 (2019), 042906, https://doi.org/10.1103/PhysRevE.100.042906. https://link.aps.org/doi/10.1103/PhysRevE.100.042906.

[122] S. Garrault, E. Finot, E. Lesniewska, A. Nonat, Study of C-S-H growth on C3S surface during its early hydration, Mater. Struct. 38 (4) (2005) 435-442, https:// doi.org/10.1007/BF02482139. URL doi.org/10.1007/BK-2018-1296. ch018.

[123] K. Ioannidou, Gels and Other Soft Amorphous Solids, in: Ch. Heterogeneity in Cement Hydrates, 2018, pp. 357-371, https://doi.org/10.1021/bk-2018-1296. ch018. https://pubs.acs.org/doi/abs/10.1021/bk-2018-1296.ch018.

[124] T. Powers, T. Brownyard, Studies of the physical properties of hardened cement paste, PCA Bull. 22 (1948).

[125] A.J. Allen, Time-resolved phenomena in cements, clays and porous rocks, URL, J. Appl. Crystallogr. 24 (5) (1991) 624-634, https://doi.org/10.1107/ S0021889890012237, https://onlinelibrary.wiley.com/doi/abs/10.1107/S00 21889890012237.

[126] E.Del Gado, K. Ioannidou, Nanoscale numerical study of C-S-H precipitation and gelation, in: ASCE, 9th International Conference on Mechanics and Physics of Creep, Shrinkage, and Durability of Concrete and Concrete Structures, 2013, pp. 30-39, https://doi.org/10.1061/9780784413111.003.

[127] H. Liu, L. Tang, N.M.A. Krishnan, G. Sant, M. Bauchy, Structural percolation controls the precipitation kinetics of colloidal calcium-silicate-hydrate gels, J. Phys. D. Appl. Phys. 52 (31) (2019), 315301, https://doi.org/10.1088/1361- 6463/ab217b. URL doi.org/10.1088/1361-6463/ab217b.

[128] S. Brisard, P. Levitz, Small-angle scattering of dense, polydisperse granular porous media: computation free of size effects, Phys. Rev. E 87 (1) (2013), 013305.

[129] P. Levitz, D. Tchoubar, Disordered porous solids: from chord distributions to small angle scattering, J. Phys. I France 2 (6) (1992) 771-790, https://doi.org/ 10.1051/jpl:1992174. URL doi.org/10.1051/jpl:1992174.

[130] S. Bae, R. Taylor, D. Shapiro, P. Denes, J. Joseph, R. Celestre, S. Marchesini, H. Padmore, T. Tyliszczak, T. Warwick, D. Kilcoyne, P. Levitz, P.J.M. Monteiro, Soft x-ray ptychographic imaging and morphological quantification of calcium silicate hydrates (c-s-h), J. Am. Ceram. Soc. 98 (12) (2015) 4090-4095, https:// doi.org/10.1111/jace.13808.

[131] K. Ioannidou, E. Masoero, P. Levitz, R.J.-M. Pellenq, E.Del Gado, Hydration kinetics and gel morphology of C-S-H, ASCE, URL, in: 10th International Conference on Mechanics and Physics of Creep, Shrinkage, and Durability of Concrete and Concrete Structures, 2015, pp. 565-573, https://doi.org/10.1061/ 9780784479346.069.

[132] K. Ioannidou, F.-J. Ulm, P. Levitz, E. Del Gado, R.J.-M. Pellenq, Nano-granular texture of cement hydrates, EPJ Web Conf. 140 (2017) 15027, https://doi.org/ 10.1051/epjconf/201714015027. URL doi.org/10.1051/epjconf/201714015027.

[133] S. Brisard, R.S. Chae, I. Bihannic, L. Michot, P. Guttmann, G. Thieme, G. Schneider, P.J. Monteiro, P. Levitz, Morphological quantification of hierarchical geomaterials by X-ray nano-CT bridges the gap from nano to micro length scales, Am. Mineral. 97 (2-3) (2012) 480-483.

[134] S. Bhattacharya, K.E. Gubbins, Fast method for computing pore size distributions of model materials, Langmuir 22 (18) (2006) 7726-7731.

[135] R. Valckenborg, L. Pel, K. Kopinga, Combined nmr cryoporometry and relaxometry, J. Phys. D. Appl. Phys. 35 (3) (2002) 249.

[136] R.S. Mikhail, L.E. Copeland, S. Brunauer, Pore structures and surface areas of hardened Portland cement pastes by nitrogen adsorption, Can. J. Chem. 42 (2) (1964) 426-438.

[137] Z. Yu, A. Zhou, D. Lau, Mesoscopic packing of disk-like building blocks in calcium silicate hydrate, Sci. Rep. 6 (1) (2016) 36967, https://doi.org/10.1038/ srep36967. URL doi.org/10.1038/srep36967.

[138] E.P. Barrett, L.G. Joyner, P.P. Halenda, The determination of pore volume and area distributions in porous substances. i. Computations from nitrogen isotherms, J. Am. Chem. Soc. 73 (1) (1951) 373-380.

[139] J. Hagmasy Jr., S. Brunauer, R.S. Mikhail, Pore structure analysis by water vapor adsorption: I. T-curves for water vapor, J. Colloid Interface Sci. 29 (3) (1969) 485-491.

[140] M.B. Pinson, E. Masoero, P.A. Bonnaud, H. Manzano, Q. Ji, S. Yip, J.J. Thomas, M. Z. Bazant, K.J. Van Vliet, H.M. Jennings, Hysteresis from multiscale porosity: modeling water sorption and shrinkage in cement paste, Phys. Rev. Appl. 3 (6) (2015), 064009.

[141] T. Zhou, K. Ioannidou, E. Masoero, M. Mirzadeh, R.J.-M. Pellenq, M.Z. Bazant, Capillary stress and structural relaxation in moist granular materials, Langmuir 35 (12) (2019) 4397-4402.

[142] V. Baroghel-Bouny, Water vapour sorption experiments on hardened cementitious materials: part i: essential tool for analysis of hygral behaviour and its relation to pore structure, Cem. Concr. Res. 37 (3) (2007) 414-437.

[143] E. Masoero, G. Cusatis, G. Di Luzio, C-s-h gel densification: the impact of the nanoscale on self-desiccation and sorption isotherms, Cem. Concr. Res. 109 (2018) 103-119.

[144] R. Feldman, Sorption and length change, scanning isotherms of methanol and water on hydrated portland cement, in: Proc. Fifth Int. Symp. Chem. Cement, Tokyo, 1968, 1968.

[145] J.J. Thomas, A.J. Allen, H.M. Jennings, Structural changes to the calcium-silicate-hydrate gel phase of hydrated cement with age, drying, and resaturation, J. Am. Ceram. Soc. 91 (10) (2008) 3362-3369.

[146] I. Maruyama, G. Igarashi, Y. Nishioka, Bimodal behavior of cSH interpreted from short-term length change and water vapor sorption isotherms of hardened cement paste, Cem. Concr. Res. 73 (2015) 158-168.

[147] M.J.A. Qomi, L. Brochard, T. Honorio, I. Maruyama, M. Vandamme, Advances in atomistic modeling and understanding of drying shrinkage in cementitious materials, Cem. Concr. Res. 148 (2021), 106536.

[148] Z.P. Bazant, A.B. Hauggaard, S. Baweja, F.-J. Ulm, Microprestress-solidification theory for concrete creep. I: aging and drying effects, J. Eng. Mech. 123 (11) (1997) 1188-1194.

[149] M. Abuhaikal, K. Ioannidou, T. Petersen, R.J.-M. Pellenq, F.-J. Ulm, Le châtelier's conjecture: measurement of colloidal eigenstresses in chemically reactive materials, J. Mech. Phys. Solids 112 (2018) 334-344.

[150] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford university press, 2017.

[151] A.P. Thompson, S.J. Plimpton, W. Mattson, General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions, J. Chem. Phys. 131 (15) (2009), 154107.

[152] H.Manzano Moro, Atomistic Simulation Studies of the Cement Paste Components, Servicio Editorial de la Universidad del Pas Vasco/Euskal Herriko, 2009.

[153] J. Lutsko, Generalized expressions for the calculation of elastic constants by computer simulation, J. Appl. Phys. 65 (8) (1989) 2991-2997.

[154] G. Constantinides, F.-J. Ulm, The nanogranular nature of c-s-h, J. Mech. Phys. Solids 55 (1) (2007) 64-90.

[155] M.A. Qomi, K. Krakowiak, M. Bauchy, K. Stewart, R. Shahsavari, D. Jagannathan, D.B. Brommer, A. Baronnet, M.J. Buehler, S. Yip, et al., Combinatorial molecular optimization of cement hydrates, Nat. Commun. 5 (1) (2014) 1-10.

[156] C. Davie, E. Masoero, Modelling damage from the nano-scale up, in: CONCREEP 10, 2015, pp. 613-620.

[157] E. Masoero, H. Manzano, E.Del Gado, R.-M. Pellenq, F.-J. Ulm, S. Yip, Kinetic simulation of the logarithmic creep of cement, in: Mechanics and Physics of Creep, Shrinkage, and Durability of Concrete: A Tribute to Zdeněk P. Bazant, 2013, pp. 166-173, FG 27.

[158] E. Masoero, M. Bauchy, E.Del Gado, H. Manzano, R. Pellenq, F.-J. Ulm, S. Yip, Kinetic simulations of cement creep: mechanisms from shear deformations of glasses, in: CONCREEP 10, 2015, pp. 555-564.

[159] A.-T. Akono, F.-J. Ulm, Scratch test model for the determination of fracture toughness, Eng. Fract. Mech. 78 (2) (2011) 334-342.

[160] K.J. Krakowiak, J.J. Thomas, S. Musso, S. James, A.-T. Akono, F.-J. Ulm, Nano- chemo-mechanical signature of conventional oil-well cement systems: effects of elevated temperature and curing time, Cem. Concr. Res. 67 (2015) 103-121.

[161] L. Brochard, G. Hantal, H. Laubie, F.J. Ulm, R.J.-M. Pellenq, Fracture mechanisms in organic-rich shales: role of kerogen, in: Poromechanics V: Proceedings of the Fifth Biot Conference on Poromechanics, 2013, pp. 2471-2480.

[162] S. Monfrared, H. Laubie, F. Radjai, R.J.-M. Pellenq, F.J. Ulm, Mesoscale poroelasticity of heterogeneous media, J. Nanomech. Micromech. 7 (4) (2017) 04017016.

[163] H. Laubie, F. Radjai, R.J.-M. Pellenq, F.J. Ulm, A potential-of-mean-force approach for fracture morphologies of heterogeneous materials using the lattice element method, J. Mech. Phys. Solids 105 (2017) 116-130.

[164] M. Vandamme, F.-J. Ulm, Nanogranular origin of concrete creep, Proc. Natl. Acad. Sci. 106 (26) (2009) 10552-10557.

[165] M. Königsberger, M. Ifan-ul Hassan, B. Pichler, C. Hellmich, Downscaling based identification of nonaging power-law creep of cement hydrates, J. Eng. Mech. 142 (12) (2016), 04016106.
[166] M. Bauchy, E. Masoero, F.-J. Ulm, R. Pellenq, Creep of bulk csh: insights from molecular dynamics simulations, in: CONCREEP 10, 2015, pp. 511–516.

[167] M. Bauchy, M. Wang, Y. Yu, B. Wang, N.A. Krishnan, E. Masoero, F.-J. Ulm, R. Pellenq, Topological control on the structural relaxation of atomic networks under stress, Phys. Rev. Lett. 119 (3) (2017), 035502.

[168] A. Morshedifard, S. Masoumi, M.A. Qomi, Nanoscale origins of creep in calcium silicate hydrates, Nat. Commun. 9 (1) (2018) 1–10.

[169] I. Richardson, The nature of the hydration products in hardened cement pastes, Cem. Concr. Compos. 22 (2) (2000) 97–113.

[170] M. Zalzale, P. McDonald, K. Scrivener, A 3d lattice boltzmann effective media study: understanding the role of csh and water saturation on the permeability of cement paste, Model. Simul. Mater. Sci. Eng. 21 (8) (2013), 085016.

[171] B. Li, D. Taniguchi, J.P. Gedara, V. Gogurancea, R. Gonzalez-Cabaleiro, J. Chen, A.S. McGough, I.D. Ofiteru, T.P. Curtis, P. Zuliani, Nufeb: a massively parallel simulator for individual-based modelling of microbial communities, PLoS Comput. Biol. 15 (12) (2019), e1007125.

[172] A. Barducci, M. Bonomi, M. Parrinello, Metadynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 1 (5) (2011) 826–843.

[173] Y. Fan, S. Yip, B. Yildiz, Autonomous basin climbing method with sampling of multiple transition pathways: application to anisotropic diffusion of point defects in hcp zr, J. Phys. Condens. Matter 26 (36) (2014), 365402.

[174] N. Mousseau, L.K. Béland, P. Brommer, J.-F. Joly, F. El-Mellouhi, E. Machado-Charry, M.-C. Marinica, P. Pochet, The activation-relaxation technique: art nouveau and kinetic art, J. At. Mol. Phys. 2012 (2012).

[175] J. Arayro, A. Dufresne, T. Zhou, K. Ioannidou, J.-F. Ulm, R. Pellenq, L.K. Béland, Thermodynamics, kinetics, and mechanics of cesium sorption in cement paste: a multiscale assessment, Phys. Rev. Mater. 2 (5) (2018), 053608.

[176] A. Dufresne, J. Arayro, T. Zhou, K. Ioannidou, F.-J. Ulm, R. Pellenq, L.K. Béland, Atomistic and mesoscale simulation of sodium and potassium adsorption in cement paste, J. Chem. Phys. 149 (7) (2018), 074705.

[177] M. Shishehbor, D. Sakaniwa, D. Stefaniuk, K. Krakowiak, M.Abdolhosseini Qomi, On the significance of interfacial chemistry on the strength of fly ash-cement composites, Cem. Concr. Res. 151 (2022), 106619.