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

HAL
open science

# Drained and undrained heat capacity of swelling clays

Tulio Honorio, Laurent Brochard

► To cite this version:

Tulio Honorio, Laurent Brochard. Drained and undrained heat capacity of swelling clays. Physical Chemistry Chemical Physics, 2022, 24 (24), pp.15003-15014. 10.1039/D2CP01419J. hal-03676559

HAL Id: hal-03676559

https://hal.science/hal-03676559

Submitted on 6 Mar 2023

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
Cite this: DOI: 00.0000/xxxxxxxxxx

Received Date
Accepted Date

DOI: 00.0000/xxxxxxxxxx

# Drained and undrained heat capacity of swelling clays†

Tulio Honorio*ᵃ and Laurent Brochardᵇ

In microporous materials, poromechanical drained (open pore) and undrained (close pore) conditions play a role in the thermo-mechanical properties, which can also be affected by specific ion effects. Heat capacity is a key property to understanding the thermo-hydro-mechanical behavior of clays. However, there is a lack of bottom-up approaches aiming at capturing the influence of confined water and ion effects on the heat capacity of clays. We perform molecular dynamics simulations to obtain the heat capacity of Na- and Ca-montmorillonite using classical and semi-classical approaches. The heat of adsorption is computed from Grand Canonical Monte Carlo (GCMC) simulations at ambient conditions. We provide the theoretical framework to compute the *drained* heat capacity using the undrained heat capacity and the heat of adsorption. The semi-classical approach is required to capture the heat capacity of the solid mineral but cannot properly capture that of the inter-layer electrolyte. A hybrid approach is proposed, which combines the semi-classical approach for the mineral and the classical approach for the electrolyte. The undrained and drained heat capacities of clays evolve linearly with the water content, which leads to the conclusion that the effects of the mixing of hydration states are not significant. Undrained heat capacity is almost independent of the counterion nature. In contrast, the overall trends of drained heat capacity plotted against the water content exhibit larger values for Na-Mmt than for Ca-Mmt. Our results are in agreement with experimental data for a variety of smectites.

## 1 Introduction

Clays are critical for a variety of environmental processes (including some natural hazards) as well as industrial applications, notably the ones related to energy production (e.g. oil and gas exploration and extraction, geothermal fields, and nuclear waste disposal)¹. Other industrial uses of clays cover pharmaceutical, cosmetic, building materials, agriculture and food production¹⁻⁵. More recently, synthetic organic-inorganic composites made with clays are being used to enhance thermal and mechanical properties of polymers⁶,⁷. A fundamental understanding of clay thermal properties is needed to either predict the behavior of clay-rich rocks and soils or to design the properties of clayey materials.

The thermal behavior of clays is a subject of interest at least

since the works of Le Chatelier⁸. Heat transfer is related to a variety of problems including thermal stimulation and artificial freezing of soils, as well as temperature regimes affecting micro-climatological and agriculture research⁵. Thermo-hydro-mechanical (THM) analysis of the stability of clay-rich soils and rocks requires precise knowledge of heat capacities of clays⁹. Watanabe et al.¹⁰ point out that the main factors to which the thermal response in THM analysis are more sensible are the permeability and heat capacity of the porous media. The heat capacity is the fundamental property governing the increase of temperature upon heating (e.g.¹¹,¹²) and playing a role in the thermal pressurization of clayey materials⁹,¹³. Similarly, THM analysis of other geomaterials such as concrete structures (calcium silicate hydrates, which are the main product of cement hydration, are phyllosilicates as clays) reveals a strong dependency of temperature and stresses development on the variations of heat capacities¹⁴.

However, despite the importance of heat capacity, few fundamental studies investigated this property. Those that do so are often restricted to macroscopic measurements. Skauge et al.¹⁵ measured the heat capacities of Na- and Ca-kaolinites, Na- and Ca-montmorillonite, illite, and attapulgite from 300 to 700 K using a differential scanning calorimeter. These authors found heat capacities at a constant pressure of 811 and 779 for J/(kg.K) for Na-Mmt and Ca-Mmt, respectively, and quite significant varia-

<footnote>ᵃ Université Paris-Saclay, CentraleSupélec, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France. Fax: +33(0)1 81 87 51 02; Tel: +33(0)1 81 87 51 02; E-mail: tulio.honorio-de-faria@ens-paris-saclay.fr</footnote>

<footnote>ᵇ Laboratoire Navier, Ecole des Ponts ParisTech, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France.</footnote>

<footnote>† Electronic Supplementary Information (ESI) available: "Derivation of the vibrational density of states and Thermoelastic properties of Na- and Ca-montmorillonite" accompany this article. See DOI: 10.1039/cXCP00000x/</footnote>

<footnote>‡ Additional footnotes to the title and authors can be included e.g. 'Present address:' or 'These authors contributed equally to this work' as above using the symbols: ‡, §, and ¶. Please place the appropriate symbol next to the author's name and include a \footnotetext entry in the text.</footnote>
tions of the heat capacity with the temperature. Their measure-
ments are also reported to be sensitive to moisture content. Gar-
rite et al.<sup>16</sup> report that changes in water content strongly affect
the heat capacity of porous media and, in particular, bentonites.
An increase of the specific heat capacity with water content is ex-
pected<sup>5,17–19</sup> since the specific heat capacity of (bulk) water is
more than three to four times the heat capacity of solid minerals
in clays. Other authors provide pressure-constant heat capacity of
bentonites and sepiolites (both with different contents of smec-
tite, kaolinites, and illite)<sup>3,17</sup> or of purified samples of smectites,
beidellites, and illites<sup>20</sup> measured from calorimetry.

Atomistic simulations have proven to be a powerful tool in the elucidation of the physical origins of the mechanical and transport properties of clays (e.g.²¹). Studies on other phyllosilicates using molecular simulations show both that the heat capacity of water decreases²²,²³ or increases²³ with confinement in micro pores. The heat capacity of confined water in micropores of silica, as reported by Etzler and co-workers²⁴⁻²⁶, displays a non-monotonic evolution when plotted against pore size and temperature. These authors indicated structural transitions in adsorbed water as the cause of the pore size and temperature dependence of the heat capacity of water in micro- and meso-pores. Water confinement effect is generally neglected in macroscopic models which extrapolate calorimetric data as a function of water content, instead, assuming free water only¹⁷. Indeed, the usual description of porous medium as a bulk fluid/solid mixture yields a simple mixing law for the total heat capacity²⁷. Moreover, confined water is reported to show anomalous thermal expansion in clay micropores⁹. Also, molecular simulations allow assessing the heat capacity as a function of the composition of clay layers as well as counterion nature. Important insights regarding the behavior of clays can be obtained from bottom-up studies. Adsorption phenomena, taking place at molecular scales, is proven to be fundamental in the thermo-mechanical behavior of clays²⁸,²⁹. To the authors' knowledge, no specific studies relating adsorption to the development of thermal properties of clays such as heat capacities have been proposed so far in the literature.

In poromechanics, *drained* and *undrained* conditions, refer-
ring respectively to open pore (the pore fluid can exchange
with the environment) and close pore (the pore fluid cannot ex-
change with the environment) conditions, are known to affect
the thermo-mechanical response of porous materials and can be
linked to a timescale related to long-term and instantaneous be-
haviors, respectively (e.g.³⁰). A theoretical background linking
adsorption phenomena at the molecular scale and estimates of
the heat capacity accounting for the chemical composition is still
to be provided in the case of clays.

The motivation of this work relies on the lack of a bottom-up
model capturing the influence of adsorption and counterion na-
ture in the development of heat capacity of clays. Here, the goal
is to provide the heat capacity of swelling clays at the nanoscale
as a function of hydration state and counterion nature. To do
so, we perform molecular simulations of the clay layer mimicking
drained geomechanical conditions. Na- and Ca- montmorillonites
are chosen due to their abundance and industrial interest. This
aspect is highlighted by the fact that montmorillonite is a well-

studied system in molecular simulations (e.g.<sup>31–36</sup>). The heat capacity is computed using classical and semi-classical methods. We discuss the obtained heat capacity of hydrated clays and confined water within clay micropores. From the GCMC simulation, we also compute the heat of adsorption. Both drained and undrained conditions are considered. The theoretical framework required to compute the *drained* heat capacity as a function of the undrained heat capacity and the heat of adsorption is provided. Based on previous works<sup>28,29</sup> indicating the propensity of clays to mix hydration states at a mesoscale, we also discuss the implications of the hydration transitions on the observed heat capacity of clay particles.

2 Models and Methods

**2.1 Undrained heat capacity from molecular simulations**

At constant volume V, the heat capacity C<sub>v</sub> quantifies how internal energy E changes with temperature T; whilst, at constant pressure P, heat capacity C<sub>p</sub> quantifies how the enthalpy H = E + PV changes with temperature:

$$
C_v = \left. \frac{\partial E}{\partial T} \right|_V ; \quad C_p = \left. \frac{\partial H}{\partial T} \right|_P \qquad (1)
$$

In general, the experimental determination of C<sub>p</sub> is preferred
since it is difficult to perform isochoric experiments.

In a full classical approach, both C<sub>v</sub> and C<sub>p</sub> can be determined
from molecular simulations in canonical (NVT) and isothermal-
isobaric ensemble (NPT), respectively. A basic approach consists
in performing a series of simulations at various T and estimating
the temperature derivatives by finite differences. A more direct
way of obtaining C<sub>v</sub> and C<sub>p</sub> is through fluctuations formula<sup>37</sup>.
The fluctuation of a thermodynamic quantity A in a given ensem-
ble Ens is defined as: ⟨δA<sup>2</sup>⟩<sub>Ens</sub> = ⟨A<sup>2</sup>⟩<sub>Ens</sub> - ⟨A⟩<sup>2</sup><sub>Ens</sub>. In this context,
C<sub>v</sub> can be computed from a simulation in the NVT ensemble from
the fluctuations of the internal energy E<sup>37</sup>:

$$
C_v = \frac{\langle \delta E^2 \rangle_{NVT}}{kT^2} \qquad (2)
$$

where *k* is the Boltzmann constant and *N* is the number of atoms.
Similarly, *C*<sub>*p*</sub> can be computed from the fluctuations of the en-
thalpy *H* in the NPT ensemble<sup>37</sup>:

$$
C_p = \frac{\langle \delta H^2 \rangle_{NPT}}{kT^2} \qquad (3)
$$

Constant volume c_v and constant pressure c_p specific heat capacities (i.e., heat capacities per unit mass) are related according to (e.g. ³⁸):

$$
c_p - c_v = \frac{T}{\rho} \alpha : \mathbf{C} : \alpha \qquad (4)
$$

where ρ is the density, α is the second-rank tensor of thermal
expansion, C is the stiffness forth-rank tensor and ":" stands for
tensor double contraction. Details on the computation of α and
C are provided in the next section.

The above estimates are based on classical physics and neglect
the quantification of energy that may become significant for ma-
terials with high phonon frequencies v (hv ~ kT, with h being
the Planck's constant). When quantum effects become significant,
the heat capacity becomes lower than the classical estimate, and quantum corrections are required to compute the heat capacity³⁹. In this respect, a semi-classical approach has been adopted by several authors to compute the heat capacity of minerals with high phonon frequencies²²,³⁹,⁴⁰. This approach considers only the vibrational contribution to the heat capacity and neglects anharmonic effects. The heat capacity is computed as the superposed contributions of the various vibrational modes weighted by their respective quantum correction⁴⁰:

$$C_v = \frac{dE_{vib}}{dT} = 3Nk \int_{0}^{\infty} \left( \frac{\hbar\omega n(\omega, T)}{kT} \right)^2 \exp\left[\frac{\hbar\omega}{kT}\right] g(\omega)d\omega \quad (5)$$

where $\omega = 2\pi v$ is the angular frequency, $\hbar = h/(2\pi)$ is the reduced Planck constant, $n(\omega, T) = (\exp(\hbar\omega/kT) - 1)^{-1}$ is Bose-Einstein distribution, and $g(\omega)$ is the density of the states. The approach is semi-classical because it involves classical evaluations of the density of states $g(\omega)$ from the Fourier transform of velocity autocorrelation function (see the Supporting Material information for details on the computation):

$$g(\omega) = \frac{1}{\pi} \int_{-\infty}^{\infty} \sum_{j} \frac{\langle m_j \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle}{\langle m_j \mathbf{v}_j^2(0) \rangle} \exp[i\omega t] dt \quad (6)$$

where $\mathbf{v}_j(t) \cdot \mathbf{v}_j(0)$ is the velocity autocorrelation function to be computed via classical simulations, and $m_j$ is the mass of particle j.

## 2.2 Thermodynamics of swelling in drained conditions

In drained conditions, bound water within clay micropores is in osmotic equilibrium with an external reservoir of free water. As discussed in a previous work²⁸,²⁹,⁴¹, drained conditions correspond to a mixed thermodynamic ensemble (NVT/ µVT), for which equilibrium corresponds to the minimization of the Helmholtz free energy of the solid skeleton and counterions plus the grand potential of the water. This hybrid thermodynamic potential $\lambda$, reduced per unit area A, can be computed from the confining pressure isotherms $P(d)$, with d the inter-layer distance of the clay layer:

$$\lambda(d, T, \mu_w) = \frac{F - \mu_w N_w}{A} = \lambda(d_0, T, \mu_w) - \int_{d_0}^{d} P(e) de \quad (7)$$

where $\mu_w$ and $N_w$ are the chemical potential and number of water molecules, respectively and $d_0$ is a reference inter-layer distance.

Due to the ordering of water adjacent to the layers in clay micropores, $\lambda$ exhibits multiple wells each corresponding to a hydration state. Clay layers are generally organized in stacks (often called a clay particle). At the scale of a particle, thermodynamic equilibrium corresponds to the minimization of the convex hull of $\lambda$ for an average deformation (which is, in layered materials, proportional to the average basal spacing $\langle d \rangle$) imposed on the boundary of the particle. An inter-layer distance d is stable if $\lambda(d)$ it is part of the convex hull of $\lambda$, so that it is homogeneous under a given load, i.e. $\langle d \rangle = d$. These stable spacings are associated with a specific hydration state. In Figure 1 (left) we identify the stable hydration states of Na- and Ca-Mmt in the confin-

ing pressure isotherms and the $\lambda$ energy profiles (plotted versus the inter-layer distance) according to the Grand Canonical Monte Carlo (GCMC) simulation detailed in Section 2.3. The hydration states corresponding to 0W, 1W, and 2W are clearly defined on Na-Mmt; whereas only the hydration states corresponding to 0W and 1W are clearly defined on Ca-Mmt.

For imposed deformations in-between two stable states, the system responds by linearly combining fractions of stacks in each state²⁹,⁴². This phenomenon is a proper characteristic of materials exhibiting multi-well energy profiles⁴³ and is thermodynamically analogous to a phase transition. Figure 1 (right) depicts a scheme of the mixing stable states in response to a given imposed displacement. Under confining pressure control, a stack of hydrated clay layers follows the stable (reversible) path indicated in Figure 1 (left).

According to the theory presented in this section, we retain that under drained conditions and confining pressure control:

* Some inter-layer distances (or correspondingly, deformations) of hydrated clay layers are thermodynamically unstable. So, the corresponding heat capacities at the layer scale would not be observed.

* When submitted to an average "unstable" deformation, the mechanical response of a stack of layers is to mix hydration states, as in a phase transition. The corresponding heat capacity at the scale of a stack of layers will be the linear combination of the heat capacities of the stable hydration states.

Transitions between hydration states involve large volumes of drained water which is expected to have a significant impact on the observed heat capacity of clay because of the heat of adsorption.

### 2.2.1 Heat of adsorption and drained heat capacities

The heat associated with adsorption may play a major role in the total heat capacity of adsorbing materials under drained conditions. The isosteric heat (or enthalpy) of adsorption is defined as:

$$q_{st} = h_{bulk} - \left. \frac{\partial H}{\partial N} \right|_{T, P_w} = h_{bulk} - h_a \quad (8)$$

where $h_{bulk}$ is the molar enthalpy of the bulk fluid and $h_a$ is the molar incremental enthalpy of the adsorbed phase (H is the enthalpy of the system upon adsorption of N molecules of fluid at a temperature T and external fluid pressure $P_w$). The isosteric heat of adsorption is the heat released by the system during the adsorption of a fluid molecule. It is generally calculated in the GCMC simulations via the cross fluctuations between the internal energy E and the number of fluid molecules $N_w$ (e.g.⁴⁵):

$$h_a = \frac{\langle EN_w \rangle - \langle E \rangle \langle N_w \rangle}{\langle N_w^2 \rangle - \langle N_w \rangle^2} \quad (9)$$

It should be noted that the application of this formula requires two independent simulations: one to obtain the enthalpy of the bulk fluid (e.g., liquid water for saturated clays) and one for the adsorbed fluid (in the grand canonical ensemble).
<figure><img src="image_3.png" /><figcaption>Fig. 1 Stable path in pressure isotherms (bottom) according to free energy λ (top) for Na- and Ca-Mmt at ambient conditions (300 K and Pw of 1 atm). The results of pressure were filtered with Savitzky-Golay algorithm$^{44}$. The mixture of hydration states in response to a thermo-mechanical loading as discussed by Brochard et al.$^{28}$. The hydration states corresponding to 0W, 1W, and 2W are clearly defined on Na-Mmt isotherms (light gray rectangles); whereas only the hydration states corresponding to 0W and 1W are clearly defined on Ca-Mmt isotherms (light blue rectangles). The limit of 1W state in Ca-Mmt curves can be approximately identified considering the average diameter of a water molecule (circa 3.5 Å).</figcaption></figure>

The usual definition of the isosteric heat of adsorption is for a system maintained at pressure Pw, which corresponds to the normal experimental conditions in adsorption experiments. However, in GCMC simulations the volume is fixed instead. Therefore, the heat of adsorption obtained from the fluctuation formula in Eq. 9 is the heat of adsorption for isochoric conditions. As in the case of the heat capacities at constant pressure or constant volume, in condensed systems, isosteric and isochoric heats of adsorption are also expected to be similar.

In this context, the drained heat capacity at constant pressure can be determined according to:

$$C_v^{\text{drained}} \approx C_v + h_a \left(\frac{\partial N_w}{\partial T}\right)_{V,P_w} \quad (10)$$

where the derivative $(\partial N_w / \partial T)_{V,P_w}$ is computed from series of grand canonical simulations under the same imposed bulk fluid pressure Pw at different temperatures.

## 2.3 Molecular models and methods

In this work, we investigate montmorillonite as a particular case of swelling clay. Montmorillonite atomic structure is built from pyrophyllite's through the following substitutions$^{46,47}$: one out of each eight octahedral aluminum atoms is substituted by a magnesium atom, and one out of each 32 tetrahedral silicon atoms is substituted by an aluminum atom. The substitution sites are randomly chosen but are not allowed to share the same bridging oxygen. The resulting structural formula is M.[Si$_{62}$Al$_2$]$_{tet}$[Mg$_4$Al$_{28}$]$_{oct}$O$_{160}$(OH)$_{32}$.nH$_2$O, where M stands for

<figure><img src="image_4.png" /><figcaption>Fig. 2 Equilibrium snapshot of Na- and Ca-Mmt at 300 K and RH=100%. The triclinic periodic super-cell is shown according to ac plane. Atoms are colored as follows: Na, blue; Ca, gray; Mg, pink; Si, yellow; O, red; H, white; Al, green. The in-plane dimensions obtained after equilibration are L$_{a}$ = 20.72 Å and L$_{b}$ = 17.96 Å. Various interlayer distances were tested.</figcaption></figure>

the counterions in the inter-layer space. We consider sodium (Na-Mmt, M = 6Na$^+$) and calcium (Ca-Mmt, M = 3Ca$^{2+}$) saturated montmorillonite. The total surface charge density of the layers (-0.124 C.m$^2$) is exactly balanced by the counterions. The simulation supercell is obtained by replicating the layers in each direction of the unit cell. The resulting triclinic supercell has two slit pores as shown in Figure 2.

Non-bonded interactions between atoms are described by Lennard-Jones and Coulomb potentials. Lennard-Jones parameters and partial charges from ClayFF force field$^{48}$ are used together with the extended simple point charge (SPC/E) model for water$^{49}$. Lennard-Jones cross-parameters follow Lorentz-
Berthelot mixing rules. Long-range electrostatic interactions are treated with the Ewald sum method with accuracy in forces of $10^{-5}$. We adopt a cutoff distance of 9.0 Å for short-range dispersive interactions. Simulations were run with periodic boundary conditions. Following ClayFF, harmonic bonds are defined for hydroxyls and harmonic angles for angles between octahedral Al and hydroxyls.

Regarding water, rigid and flexible SPC/E models are considered. ClayFF provides angle and covalent bond rigidities to make water flexible, and both SPC and SPC/E models have been used with ClayFF⁵⁰. Rigid three-site water models such as SPC, SPC/E, and TIP3P are reported to best capture the experimental heat capacities of water⁵¹. This discrepancy will be analyzed in detail in the next section by comparing rigid and flexible SPC/E waters.

Simulations in NVT and NPT ensembles are performed with LAMMPS⁵². In these simulations, we use Nosé-Hoover thermostat and barostat with damping parameters of 100 and 1000 timesteps, respectively. After relaxation at 300 K and 1 atm the simulation super-cell have the ab-plane dimensions $L_a = 20.72$ Å and $L_b = 17.96$ Å. To compute undrained heat capacity, NVT simulations are run for 0.5 ns for classical estimates and NVE simulations (after proper equilibration under NVT conditions during 0.1 ns followed by 0.1 NVE equilibration) for semi-classical estimates. Water molecules are constrained with SHAKE algorithm.

The number of water molecules is a function of the interlayer distance, temperature, and chemical potential of water (which is related to the pressure of the fluid reservoir). Here, we consider inter-layer distances ranging from 9 to 20 Å, which comprises the dehydrated up to capillary pore water states in clays¹. The water content is based on hybrid GCMC-NVT simulations detailed in the Supporting Material. We adopt a chemical potential corresponding to an RH of 100% (i.e. the water pressure $P_w$ to be imposed is the liquid-vapor coexistence pressure, which is 0.01004 atm for SPC/E water at 300 K⁵³). Figure 3 displays the number of water molecules $N_w$ per unit cell ([Si<sub>7.75</sub>Al<sub>0.25</sub>]<sub>tot</sub> [Mg<sub>0.5</sub>Al<sub>3.5</sub>]<sub>oct</sub>O<sub>20</sub>(OH)<sub>4</sub>) as a function of interlayer distance obtained from GCMC simulations. The water content for Na-Mmt is in agreement with previous studies using configurational bias GCMC²⁹.

The drained heat capacity is computed from GCMC simulations. First, the isosteric heat of adsorption is obtained from the fluctuations of potential energy and the number of particles. And then, the derivative ($\partial N_w/\partial T$)$_{V,P_w}$ is computed from a series of grand canonical simulations at three different temperatures (290, 300, and 310 K) under the same imposed water pressure $P_w = 1$ atm.

# 3 Results and discussion

## 3.1 Heat capacity of bulk water

The analysis of bulk water behavior is crucial to determine the molecular model and approach better suited to capture the heat capacity of hydrated clays.

Figure 4 shows the specific heat capacity $c_v$ of bulk water as obtained from molecular simulations and experiments. Specific $c_v$ and total heat capacity $C_v$ are related according to $c_v = C_v / (N_w m_w)$, with $N_w$ the number of water molecules and $m_w$ the

<figure><img src="image_5.png" /><figcaption>Fig. 3 Water content in terms of the number of water molecules $N_w$ per unit cell ([Si<sub>7.75</sub>Al<sub>0.25</sub>]<sub>tot</sub> [Mg<sub>0.5</sub>Al<sub>3.5</sub>]<sub>oct</sub>O<sub>20</sub>(OH)<sub>4</sub>) as a function of the inter-layer distance obtained from GCMC simulation of Na-Mmt and Ca-Mmt. The standard deviations of $N_w$ in Ca-Mmt curve have the order of magnitude of the circles.</figcaption></figure>

mass of a water molecule. Experimental data on vapor show that $C_v$ goes from (6/2)$N_w k$ at 300 K to $6N_w k$ at 5000 K with the transition taking place around 2000 K. The first case ((6/2)$N_w k$) corresponds to Dulong-Petit law when only kinetic contributions are accounted for, i.e., when each one of the six degrees of freedom, related to the three possible translations and three possible rotations of each water molecule, contributes with $k/2$. The second case ($6N_w k$) corresponds to the Dulong-Petit law when both kinetic degrees of freedom and internal vibrations are accounted for, the latter include bending, symmetric, and anti-symmetric stretching modes for water molecules. This transition occurring at high temperatures shows that internal vibrations in water are not contributing to the heat capacity at ambient temperature.

In this context, considering internal vibrations inactive, $C_v$ of the condensed phase is expected to be $6N_w k$, excluding anharmonicity. Experimental evidence shows, however, that the $C_v$ of liquid water is much higher than $6N_w k$, which is expected since anharmonicity is major in liquids. The semi-classical approach assumes harmonicity and provides results that cannot exceed Dulong-Petit law. For this reason, the semi-classical approach is not applicable to liquid water.

As shown in Figure 4, the application of the full classical approach to the rigid SPC/E water model predicts reasonably well the $c_v$ of liquid water, being, therefore, the preferred method to study this property by molecular simulation. In contrast, the classical approach applied to the flexible SPC/E water model yields estimates that far exceeds the experimental heat capacity of liquid water (Fig. 4). Flexible SPC/E overestimates significantly $c_v$ because it accounts for internal vibrations, which should not contribute at ambient temperature.

In summary, these observations show that Dulong-Petit law does not apply to liquid water, and, that the appropriate molecular simulation method to estimate the heat capacity of liquid water is the classical approach applied to rigid models.
<figure><img src="image_6.png" /><figcaption>Fig. 4 Specific heat capacity $c_v$ of bulk water: comparisons of simulations results using rigid and flexible SPC/E at constant density as a function of the temperature with experimental data from NIST database. Dulong-Petit estimates regarding only kinetics contribution $C_v = (6/2)N_w k$, and both kinetics and internal vibrations ($C_v = 6N_w k$) are shown together with experimental measurements of water vapor as a function of the temperature.</figcaption></figure>

## 3.2 Heat capacity of dry clay

When only the dry mineral is considered, the classical approach yields an estimated specific heat capacity of $c_v = 1339 \pm 27$ J/(kg.K) = $24.2 \pm 0.5$ J/(mol.K) for Na-Mmt and $c_v = 1354 \pm 55$ J/(kg.K) = $24.4 \pm 0.9$ J/(mol.K) for Ca-Mmt. This compares well with the Dulong-Petit law $C_v = 3N_w k$, where $N_w$ is the number of atoms in the dry mineral, corresponding to 3 vibrations per atom, that is $c_v = 3R = 24.9$ J/(mol.K). The fact that Dulong-Petit law is recovered by the classical estimate shows that anharmonicity is negligible for the dry mineral, and the energy is mostly vibrational. Therefore, the semi-classical approach is valid for the dry mineral and can be used to evaluate whether some high-frequency vibrations should be excluded from the heat capacity estimate. The semi-classical approach applied to the dry clay yields an estimated specific heat capacity of $c_v = 763 \pm 13$ J/(kg.K) = $13.8 \pm 0.2$ J/(mol.K) for Na-Mmt and $c_v = 812 \pm 2$ J/(kg.K) = $14.7 \pm 0.1$ J/(mol.K) for Ca-Mmt (or, when the counterions are not accounted for $c_v = 758 \pm 10$ J/(kg.K) = $13.6 \pm 0.2$ J/(mol.K) for Na-Mmt and $c_v = 820 \pm 1$ J/(kg.K) = $14.6 \pm 0.1$ J/(mol.K) for Ca-Mmt), which is significantly lower than the estimation with the classical approach. Therefore a significant fraction of the vibrational modes are inactive. Figure 5 shows the vibrational density of states (VDOS) for dry Na- and Ca-Mmt and the contribution per atom type. According to the semi-classical formula (Eq. 5), the vibrational modes with high frequencies compared to the thermal agitation energy (i.e. $\hbar\omega$ compared to $kT$) do not contribute to the heat capacity. At ambient temperature, vibrations with $\omega$ larger than a hundred of Thz do not contribute significantly to the heat capacity. This is the case of a significant fraction of modes in the dry clay, in particular the hydroxyl groups (Hh with $\omega \approx 700$ Thz). The temperature associated with this energy level (so that $\hbar\omega \approx kT$) exceeds 5000 K. Accordingly, a significant part of the difference between the Dulong-Petit prediction and the

semi-classical estimate for the mineral can be attributed to the vibration of the hydroxyl groups. Interestingly, the counter-ion (Na or Ca) has little impact of the heat capacity and both are involved in low frequency vibrations.

Since the difference between $C_p$ and $C_v$ are minor for solids, we can directly compare experimental data with the $C_v$ obtained from simulations. Most of the experimental studies on dehydrated swelling clays report a specific heat capacity around 800 J/(kg.K) for the dry mineral: $c_p = 741 - 873$ J/(kg.K)) for Na-, Ca- and K-bentonites$^{54}$, $c_p = 811$ and 779 J/(kg.K)) for Na- and Ca-montmorillonite, respectively$^{15}$, and 837 J/(kg.K)) for bentonites as reported by Siebel et al. (as cited in$^2$, this study also report other authors obtaining almost twice the values above: $c_p = 1438$ and 1675 J/(kg.K) for dry bentonites computed from the extrapolation of the results of Charm et al. and Dieckerson et al., respectively). Most available experimental results, in particular recent ones, are consistent with the molecular simulation estimates based on the semi-classical approach. Accordingly, the appropriate methodology for molecular simulation estimates is the semi-classical approach.

## 3.3 Undrained heat capacity of hydrated clay

The discussion of the last two sections leads to the conclusion that the good methodology to estimate the heat capacity of hydrated clays should be a mix of (i) the full classical approach combined with a rigid water model for the inter-layer fluid, and (ii) a semi-classical approach for the mineral. Note that, in this hybrid approach, the inter-layer fluid includes the water and the counter-ions, because anharmonicity is expected to become important for the counter-ions once they become hydrated. As a result, the counter-ions are excluded from the minerals. How to implement such hybrid estimation is not so trivial. A possible strategy to obtain the inter-layer contribution is to compute the classical approach for the entire system ($C_{v-tot}^{class}$). Since the classical approach applied to the dry mineral yields Dulong-Petit law ($C_{v-min}^{class} = 3N_w k$), the contribution of the inter-layer is obtained by subtracting this quantity from the total heat classical capacity:

$$C_{v-w}^{class} = C_{v-tot}^{class} - 3N_w k \quad (11)$$

Regarding the mineral contribution, one can compute the vibrational density of state for the mineral atoms only, by considering the velocity auto-correlation specifically of those atoms. This yields the semi-classical heat capacity of the mineral ($C_{v-min}^{semi-class}$). The total heat capacity is then obtained has the sum of the inter-layer and mineral contributions:

$$C_{v-tot} = C_{v-w}^{class} + C_{v-min}^{semi-class} = C_{v-tot}^{class} - 3N_w k + C_{v-min}^{semi-class} \quad (12)$$

In terms of the specific heat capacities Eq. 12 becomes:
<figure><img src="image_7.png" /><figcaption>Fig. 5 VDOS of dry mineral Na-Mmt (left) and Ca-Mmt (right): (A) and (B) total VDOS of the dry mineral; (C) and (D) VDOS per species.</figcaption></figure>

$$c_{v-tot} = f_w c_{v-w}^{class} + (1 - f_w) c_{v-min}^{semi-class} = c_{v-tot}^{class} + (1 - f_w) \left( c_{v-min}^{semi-class} - \frac{3N_s k}{m_s} \right) \quad (13)$$

where $f_w = m_w/m_{total}$ and $(1 - f_w = m_s/m_{total})$ are the mass fractions of water and mineral, respectively.

Heat capacities at constant volume of Na- and Ca-Mmt layers in undrained conditions as a function of the inter-layer distance at 300 K are shown in Figure 6(A) and (B). Heat capacities were computed following the hybrid approach detailed above (see the Supporting Material for the details of the semi-classical computations). Heat capacities of Na- and Ca-Mmt appear similar. This is in agreement with experimental observations⁵⁵ in which the heat capacity of clays is reported to be independent of the nature of cation.

Figure 6 (C) and (D) show the heat capacity $C_{v-min}^{semi-class}$ of the mineral (not accounting for counterions) computed using the semi-classical approach for all configurations of Na- and Ca-Mmt studied here. The $C_{v-min}^{semi-class}$ of Na-Mmt shows a dependence to the inter-layer distance for values below circa 11 Å; whereas $C_{v-min}^{semi-class}$ remains approximately constant for Ca-Mmt irrespective of the inter-layer distance. The average values of $C_{v-min}^{semi-class}$ are 821.4 ± 22.4 J/(kg.K) for Na-Mmt, and 824.4 ± 11.3 J/(kg.K) for Ca-Mmt.

Figures 6(E) and (F) show the heat capacity $C_v$ of confined electrolyte (water and counter-ions) for both Na- and Ca-Mmt as a function of $d$. The heat capacity of the confined fluid reaches a plateau for $d$ above 13 Å, corresponding to the heat capacity

of bulk water. This suggests that the effects of confinement become negligible for larger inter-layer distances. Due to confinement, fluid molecules order themselves by the surface of the clay. In other microporous materials, adsorbed water was reported to exhibit a different heat capacity compared to bulk water²²,²⁵. Na and Ca cations are considered kosmotropes, i.e. order- or structure-makers with respect to water, in Hofmeister series⁵⁶. In both cases, since the mobility of water molecules is constrained, the expected overall behavior would be a decrease in the heat capacity of the fluid as observed in our results.

Because of the almost linear evolution of the heat capacity with the inter-layer distance, even for small inter-layer distances, the effect of the mixing of hydration states under pressure control (as discussed in Section 2.2) is expected to be minor. A comparison of the heat capacity directly computed from MD with estimates assuming a mix of hydration states in the unstable zones is shown in the Supporting Material (Figure S7).

## 3.4 Heat of adsorption of clays

Figures 7 (A) and (B) show the molar incremental enthalpy $h_a$ of water, which rapidly decreases with the inter-layer distance reaching the asymptotic value corresponding to the excess enthalpy of bulk liquid SPC/E water. We computed the enthalpy of bulk SPC/E and obtained -39.6±1.5 kJ/mol.

Keren and Shainberg⁵⁸ reported differences between the heat of adsorption of water ($h_{bulk} - h_a$) in Na- and Ca-Mmt. With Na-Mmt, the first layer of water adsorbed (1W state) releases circa 7.5 kJ/mol, whilst the subsequent layers release about 1.5 kJ/mol. On the other hand, for Ca-Mmt, the heat released varied roughly from 14 to 1.5 kJ/mol upon hydration. In Figure
<figure><img src="image_8.png" /><figcaption>Fig. 6 At the top, specific heat capacity $c_v$ of (A) Na- and (B) Ca-Mmt as a function of the inter-layer distance $d$ in undrained conditions at 300 K. At the center, specific heat capacity $c_v$ of the solid mineral, obtained by the semi-classical approach, of (C) Na- and (D) Ca-Mmt as a function of $d$. The red line depicts the classical approach applied to the dry mineral, that is the Dulong-Petit estimate ($c_v^{class} = 3N_s k$). At the bottom, the specific heat capacity of the confined electrolyte (water and counterions) in (E) Na- and (F) Ca-Mmt as a function of $d$: computation using the classical approach and rigid SPC/E model. The green line shows the case of bulk liquid (rigid) SPC/E water for comparison.</figcaption></figure>

7, the heat of adsorption ($h_{bulk} - h_a$) appears positive precisely in the domains of stability of the hydration states ($d \sim 13$ for the 1W state, $d \sim 16$ for 2W, and $d \sim 19$ for 3W). In these stability domains, the magnitude of $h_{bulk} - h_a$ is typically between 0 and 10 kJ/mol, which is consistent with the magnitudes measured by Keren and Shainberg⁵⁸. The heat of adsorption for Ca-Mmt does not seem larger than that of Na-Mmt, unlike what is observed by the same authors. But, the accuracy of the molecular simulation results (a few kJ/mol) is not good enough to provide a truly quantitative comparison. At this level of accuracy, the molecular simulation results are reasonably consistent with experiments.

## 3.5 Drained heat capacity of clays at constant volume

To compute the drained heat capacity of clays, we also provide estimates of the partial derivative of the number of water molecules $N_w$ with respect to the temperature under a fixed water pressure $P_w$ (Fig. 7 (C) and (D)) for both Na- and Ca-Mmt. With these results and the ones of the last section, we can compute the drained heat capacity using Eq. 10. In Figure 8, the drained and undrained specific heat capacities at constant volume $c_v$ of Na- and Ca-Mmt are plotted against the water content in terms of the

mass of water divided by the total mass. The drained heat capacity is calculated either using the results directly computed from MD (i.e. following the oscillatory profiles in Figure 1) or assuming a mix of hydration states in the unstable zones (i.e. following the stable path as depicted by the dashed lines representing the convex hull in Figure 1). As expected, drained $c_v$ is much larger than undrained $c_v$. The evolution of the drained $c_v$ as a function of the water content deviates from the clear linear trend observed for the undrained $c_v$. For both Na- and Ca-Mmt, some differences can be observed between the two estimates of the drained $c_v$, but these remain moderate with respect to the accuracy of the results. A more pronounced difference between both $c_v$ is observed for Na-Mmt. The mixing of hydration state yields higher values, especially for Na-mmt, because mixing hydration state provides lower estimates of the incremental enthalpy $h_a$.

## 3.6 Drained and undrained heat capacities of clays at constant pressure

In Figure 9, we compare our results of specific heat capacities at constant pressure $c_p$ of Na- and Ca-Mmt with the data found in the literature for a variety of montmorillonites and ben-
<figure><img src="image_9.png" /></figure>

Fig. 7 At the top, molar incremental enthalpy $h_a$ of water in (A) Na- and (B) Ca-Mmt, obtained with Eq. 9 from GCMC simulations, as a function of the inter-layer distance. For comparison, we display the enthalpy of bulk liquid SPC/E water computed in this work (-39.6±1.5 kJ/mol, the average values are represented by the dashed line and the variability by the dotted line). At the bottom, the partial derivative of the number of water molecules $N_w$ with respect to the temperature $T$ under a fixed water pressure $P_w$ for both (C) Na- and (D) Ca-Mmt. The error bars are computed through the formula of propagation of errors presented in⁵⁷ for the coefficients of linear regression. In all cases, the results were filtered with the Savitzky-Golay filter.

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

Fig. 8 Drained (circles) and undrained (squares) specific heat capacity $c_v$ of (A) Na- and (B) Ca-Mmt as a function of water mass content $m_w/m_{total}$ with $m_{total} = m_w + m_s$ the total mass. The undrained heat capacities at constant volume are computed using the semi-classical approach for the mineral and the classical approach (associated with rigid SPC/E) for the electrolyte (water and counterion). The drained heat capacity is calculated from the undrained one following Eq. 10 using the results directly computed from MD (i.e. following the oscillatory profiles in Figure 1) and with estimates assuming a mix of hydration states in the unstable zones (i.e. following the stable path at depicted by the dashed lines representing the convex hull in Figure 1). In the least-squares fitted line $y = Ax + B$, the coefficient $A$ corresponds to the heat capacities of the mineral, which are similar for Na- and Ca-Mmt (717 and 722 J/(kg.K), respectively).

The $c_p$ values are computed from $c_v$ using Mayer's relations (Eq. 4) using the full stiffness C and ther-
the Supporting Material). The drained $c_p$ is computed from the assumption that the isostatic heat of adsorption $q'$ is similar to $q_{st}$ as well as $\left.\frac{\partial N}{\partial T}\right|_{\varepsilon,P_w}$ to $\left.\frac{\partial N}{\partial T}\right|_{\sigma,P_w}$ in condensed systems. So a similar "drained" contribution term $q' \left.\frac{\partial N}{\partial T}\right|_{\sigma,P_w}$, as in Eq. 10, was added to the undrained $c_p$ to obtain the drained $c_p$.

Both drained and undrained specific heat capacities at constant pressure $c_p$ increase with water content as expected$^{5,16-19}$. Our results cover a range of water content little explored by experimental studies. In Ransom et al.$^{55}$ a variety of Na-smectites with different isomorphic substitutions are tested considering low water contents (about 0.20 g of water per total mass of clays). These authors obtained very similar heat capacities as the ones reported here. This might indicate that the heat capacity of smectites is poorly affected by isomorphic substitutions. As pointed out by the same authors$^{55}$, we also observe that the effects of ions are negligible with respect to the heat capacity of *hydrated clays*.

At inter-layer distances greater than 20 Å, the water behaves almost as free water (capillary domain), so that one can reasonably assume that volume and energy become extensive thermodynamic quantities. Therefore, heat capacity can be extrapolated to larger water contents of interest by a simple rule of mixture. Following a rule of mixture, $c_p = \sum_{i=1}^{N} f_i c_p^i$, where $f_i$ and $c_p^i$ are the mass fraction and heat capacity associated to a constituent phase $i$ in the mixture. $c_v$ and $c_p$ exhibit linear evolution with the water content even at inter-layer distances lower than 20 Å, which suggests that the rule of mixture remains approximately correct even at very small inter-layer distances. In Figure 9, we display the extrapolations of our results to the range usually explored in experiments. The extrapolation of the drained heat capacity returns results above the experimental data, whereas the extrapolation of the undrained heat capacity shows a good agreement with the experimental data of various authors. This observation can be understood when observing that the scale of the clay layer does not capture the presence of free water at larger scales. When bound water is drained from the layer scale, it is converted into free water in larger pores, so that it is not actually drained from the material. As we have seen, the incremental molar enthalpy of bound water is quite close to the molar enthalpy of free water, so that converting bound water into free water has only a minor effect on the heat capacity.

## 4 Conclusions

In this paper, we determined the drained and undrained heat capacities of sodium and calcium saturated montmorillonite at the mineral layer scale from molecular simulations. The inter-layer pores of clays were filled with water by means of grand canonical simulations mimicking geomechanical drained conditions. The water chemical potential imposed falls in the liquid domain, i.e. the material is fully saturated with liquid or adsorbed water. The main findings of this work are:

* The heat capacity obtained through a semi-classical approach leads to estimates of the heat capacity of the dehydrated clay systems that are consistent with experimental data. The quantum corrections reduce the estimates of heat

capacity to 50-60% of the value obtained in the full classical approach. Regarding water, the proper estimation of the heat capacity requires using the classical approach with a rigid water model (SPC/E here). Flexible SPC/E would lead to an overestimation of the heat capacity of clays by ~50% because of the high frequencies of internal vibrations. And the semi-classical approach is not applicable to liquid water because anharmonic effects cannot be neglected. An appropriate methodology to estimate the heat capacity of hydrated clay is to combine the classical approach for the electrolyte with a rigid water model and the semi-classical approach for the mineral.

* The computed undrained heat capacity of clays evolves linearly with the water content as expected$^{5,17-19}$, and its linear extrapolation to larger water content compares well with experimental results for smectites, as well as different compositions of bentonites$^{2,3,15,54,59,60}$. This suggests that the simple mixing rule valid for free water still applies reasonably in the domain of bound water. It is noteworthy that with the consideration of a specific kind of montmorillonite (with fixed number and location of isomorphic substitutions), it is possible to reasonably estimate the heat capacity of a large range of smectites and other clay minerals. Also, as in$^{55}$, we observe that the effects of ions are negligible with respect to the heat capacity of *hydrated clays*. In swelling clays the heat capacity of water is reduced by confinement and only slight counterion specificity is observed. In Na- and Ca-Mmt, confined water retrieves its bulk heat capacity for small inter-layer distances (above circa 13 Å).

* At the layer scale, the drained heat capacity is much larger than the undrained one, because it is dominated by the drainage of bound water upon heating. Nevertheless, experimental measurements are consistent with the undrained estimates. A reasonable interpretation is that the inter-layer drained bound water is converted into free water in larger pores, but is not actually drained from the material.

* Multistability and mixing of hydration states are key elements to understand the thermo-mechanical behavior of clays at the layer scale$^{28,29,61}$. Under confining pressure control, only some values of heat capacity of clays when plotted against the inter-layer distance are allowed. Due to the almost linear evolution of the heat capacity of clays with the inter-layer distance, these effects of mixing of hydration states have little impact on the undrained heat capacity evolution with the pore size or water content. However, these effects are more significant to understanding the drained heat capacity since the enthalpic contributions under drained conditions are directly linked to the hydration states and associated transitions. Our results show that the drained heat capacity is ion specific.

This work leaves room for a better understanding of factors affecting the thermo-mechanical behavior of clays and other phyllosilicates (e.g. calcium silicate hydrates, chrysotile). Future work
<figure><img src="image_11.png" /><figcaption>Fig. 9 Drained and undrained specific heat capacity at constant pressure $c_p$ as a function of water content ($m_w/m_{total}$). The results obtained in this article from molecular simulations concerning Na-Mmt and Ca-Mmt are shown inside the red boxes. These results are extrapolated to larger water content following the mixing rule (dashed lines, Ca-Mmt and Na-Mmt results are superposed). Drained $c_p$ is computed considering the mixing of hydration states.</figcaption></figure>

can cope with the effects of temperature on the heat capacity, provided that thermal stimulation is an important activity related to the exploration of clay-rich soils and rocks in the energy production industry.

## Author Contributions

**Tulio Honorio:** Conceptualization, Formal Analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft Writing – review and editing.
**Laurent Brochard:** Funding acquisition, Formal Analysis, Investigation, Methodology, Writing – original draft Writing – review and editing

## Conflicts of interest

"There are no conflicts to declare".

## Acknowledgements

The financial support of the French National Research Agency (ANR) through the project TEAM2ClayDesicc (ANR-14-CE05-0023-01) is gratefully acknowledged.

## Notes and references

1. A. Meunier, *Clays*, Springer Berlin Heidelberg, 2005, pp. 231–293.

2. A. Niekamp, K. Unklesbay, N. Unklesbay and M. Ellersieck, *Journal of Food Science*, 1984, **49**, 28–31.

3. L. M. Casas, J. L. Legido, M. Pozo, L. Mourelle, F. Plantier and D. Bessières, *Thermochimica Acta*, 2011, **524**, 68–73.

4. M. I. Carretero and M. Pozo, *Applied Clay Science*, 2009, **46**, 73–80.

5. O. T. Farouki, *Thermal Properties of Soils*, Defense technical report, 1981.

6. C. Ortiz and M. C. Boyce, *Science*, 2008, **319**, 1053–1054.

7. G. Hantal, L. Brochard, R. J.-M. Pellenq, F.-J. Ulm and B. Coasne, *Langmuir*, 2017.

8. H. Le Chatelier, *Bulletin de la Société Française de Minéralogie*, 1887, **10**, 204–207.

9. L. Brochard and T. Honório, *International Journal of Engineering Science*, 2020, **152**, 103296.

10. N. Watanabe, W. Wang, C. I. McDermott, T. Taniguchi and O. Kolditz, *Computational Mechanics*, 2010, **45**, 263.

11. L. Laloui, M. Nuth and L. Vulliet, *International Journal for Numerical and Analytical Methods in Geomechanics*, 2006, **30**, 763–781.

12. J. Rutqvist, *Developments in Clay Science*, Elsevier, 2015, vol. 6, pp. 329–355.

13. L. Brochard and T. Honório, *Acta Geotechnica*, 2021.

14. M. Briffaut, F. Benboudjema, J.-M. Torrenti and G. Nahas, *European Journal of Environmental and Civil Engineering*, 2012, **16**, 589–605.

15. A. Skauge, N. Fuller and L. G. Hepler, *Thermochimica Acta*, 1983, **61**, 139–145.

16. B. Garitte, H. Shao, X. R. Wang, T. S. Nguyen, Z. Li, J. Rutqvist, J. Birkholzer, W. Q. Wang, O. Kolditz, P. Z. Pan, X. T. Feng, C. Lee, B. J. Graupner, K. Maekawa, C. Manepally, B. Dasgupta, S. Stothoff, G. Ofoegbu, R. Fedors and J. D. Barnichon, *Environmental Earth Sciences*, 2017, **76**, 64.

17. J. L. Legido, C. Medina, M. Lourdes Mourelle, M. I. Carretero and M. Pozo, *Applied Clay Science*, 2007, **36**, 148–160.
18 J. D. Oster and P. F. Low, *Soil Science Society of America Journal*, 1964, **28**, 605-609.

19 P. F. Low, P. Hoekstra and D. M. Anderson, *Water Resources Research*, 1968, **4**, 541–544.

20 H. Gailhanou, P. Blanc, J. Rogez, G. Mikaelian, H. Kawaji, J. Olives, M. Amouric, R. Denoyel, S. Bourrelly, V. Montouill-out, P. Vieillard, C. I. Fialips, N. Michau and E. C. Gaucher, *Geochimica et Cosmochimica Acta*, 2012, **89**, 279–301.

21 C. Tournassat, C. I. Steefel, I. C. Bourg and F. Bergaya, *Natural and Engineered Clay Barriers*, Elsevier, 2015.

22 M. J. A. Qomi, F.-J. Ulm and R. J.-M. Pellenq, *Physical Review Applied*, 2015, **3**, 064010.

23 T. Honorio, F. Masara and F. Benboudjema, *Cement*, 2021, 100015.

24 F. M. Etzler and P. J. White, *Journal of Colloid and Interface Science*, 1987, **120**, 94–99.

25 F. M. Etzler and J. J. Conners, *Langmuir*, 1991, **7**, 2293–2297.

26 F. M. Etzler and J. J. Conners, *Langmuir*, 1990, **6**, 1250–1253.

27 O. Coussy, *Mechanics and Physics of Porous Solids*, John Wiley & Sons, Ltd, Chichester, UK, 2010.

28 L. Brochard, T. Honorio, M. Vandamme, M. Bornert and M. Peigney, *Acta Geotechnica*, 2017, 1–19.

29 T. Honorio, L. Brochard and M. Vandamme, *Langmuir*, 2017, **33**, 12766–12776.

30 M. J. Abdolhosseini Qomi, L. Brochard, T. Honorio, I. Maruyama and M. Vandamme, *Cement and Concrete Research*, 2021, **148**, 106536.

31 E. S. Boek, P. V. Coveney and N. T. Skipper, *Langmuir*, 1995, **11**, 4629–4631.

32 M. Holmboe and I. C. Bourg, *The Journal of Physical Chemistry C*, 2014, **118**, 1001–1013.

33 S. L. Teich-McGoldrick, J. A. Greathouse, C. F. Jové-Colón and R. T. Cygan, *The Journal of Physical Chemistry C*, 2015, **119**, 20880–20891.

34 M. de Lourdes Chávez, L. de Pablo and J. J. de Pablo, *Langmuir*, 2004, **20**, 10764–10770.

35 R. M. Shroll and D. E. Smith, *The Journal of Chemical Physics*, 1999, **111**, 9025–9033.

36 T. Underwood, V. Erastova and H. C. Greenwell, *Clays and Clay Minerals*, 2016, **64**, 472–487.

37 M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids*, Oxford University Press, New York, 1989.

38 D. C. Wallace, *Thermodynamics of Crystals*, Courier Corporation, 1998.

39 B. Winkler and M. T. Dove, *Physics and Chemistry of Minerals*, 1992, **18**, 407–415.

40 A. R. Oganov, J. P. Brodholt and G. David Price, *Physics of the Earth and Planetary Interiors*, 2000, **122**, 277–288.

41 T. Honorio, *Langmuir*, 2019, **35**, 3907–3916.

42 R. D. James, *Journal of the Mechanics and Physics of Solids*, 1986, **34**, 359–394.

43 K. Bhattacharya, *Theory of martensitic microstructure and the shape-memory effect*, 2003.

44 A. Roth, *Savitzky-Golay Smoothing Filters*, 2000, http://library.wolfram.com/infocenter/MathSource/789/.

45 D. Nicholson, *Computer simulation and the statistical mechanics of adsorption*, Academic Press, 1982.

46 B. Carrier, *PhD Thesis*, Université Paris-Est, 2013.

47 B. Carrier, M. Vandamme, R. J.-M. Pellenq and H. Van Damme, *The Journal of Physical Chemistry C*, 2014, **118**, 8933–8943.

48 R. T. Cygan, J.-J. Liang and A. G. Kalinichev, *The Journal of Physical Chemistry B*, 2004, **108**, 1255–1266.

49 H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, *The Journal of Physical Chemistry*, 1987, **91**, 6269–6271.

50 R. T. Cygan, J. A. Greathouse and A. G. Kalinichev, *The Journal of Physical Chemistry C*, 2021.

51 Y. Mao and Y. Zhang, *Chemical Physics Letters*, 2012, **542**, 37–41.

52 S. Plimpton, *Journal of Computational Physics*, 1995, **117**, 1–19.

53 M. Fugel and V. C. Weiss, *The Journal of Chemical Physics*, 2017, **146**, 064505.

54 T. Kozlowski, A. Rusin and E. Nartowska, *Applied Clay Science*, 2014, **101**, 82–89.

55 B. Ransom and H. C. Helgeson, *Geochimica et Cosmochimica Acta*, 1994, **58**, 4537–4547.

56 D. J. Tobias and J. C. Hemminger, *Science*, 2008, **319**, 1197–1198.

57 J. Taylor, *Published by University Science Books, 648 Broadway, Suite 902, New York, NY 10012, 1997.*, 1997.

58 R. Keren, *Clays and Clay Minerals*, 1975, **23**, 193–200.

59 T. Ferrand and J. Yvon, *Applied Clay Science*, 1991, **6**, 21–38.

60 S. Cara, G. Carcangiu, G. Padalino, M. Palomba and M. Tamanini, *Applied Clay Science*, 2000, **16**, 125–132.

61 L. Brochard, *The Journal of Physical Chemistry C*, 2021.