# Thermoporoelastic Moduli from Molecular Fluctuations and Application to Moisture Effect on Amorphous Cellulose

Laurent Brochard

► To cite this version:

Laurent Brochard. Thermoporoelastic Moduli from Molecular Fluctuations and Application to Moisture Effect on Amorphous Cellulose. Journal of Engineering Mechanics - ASCE, 2023, 149 (11), 10.1061/JENMDT.EMENG-7020. [https://doi.org/10.1061/JENMDT.EMENG-7020](https://doi.org/10.1061/JENMDT.EMENG-7020). [hal-04206781](https://HAL.inistre.fr/hal-04206781)

HAL Id: hal-04206781

https://hal.science/hal-04206781v1

Submitted on 4 Sep 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Thermo-poro-elastic moduli from molecular fluctuations and application to moisture effect on amorphous cellulose

Laurent Brochard¹

¹Laboratoire Navier (UMR 8205), CNRS, ENPC, Univ. Gustave Eiffel, 6 & 8 avenue Blaise Pascal, 77455 Marne-la-Vallée, France. Email: laurent.brochard@enpc.fr

## ABSTRACT

Although imperceptible at the macroscale, fluctuations of physical properties are ubiquitous at the nanoscale because of thermal agitation. While fluctuations are detrimental to the accuracy of molecular simulations, their magnitude is actually a great source of information since it quantifies the elasticity of the system with respect to any evolution of the control variables. Fluctuation formulas are well known for fluids or solids, but have not been derived yet for porous media. In this article, we derive fluctuation formulas for all the moduli involved in thermo-poro-elasticity. These fluctuation formulas are applicable even for adsorbing media, and can be used for poroelasticity extended to adsorption, whereas usual poroelasticity assumes a bulk fluid in the pores. Moreover, the practical application of the fluctuation formulas in molecular simulation is straightforward as one only needs to collect extensive quantities readily accessible and at almost no additional computational cost. The application of these fluctuation formulas is illustrated by considering a typical example of mecano-sorption coupling, namely amorphous cellulose submitted to moisture at various relative humidities. The thermo-poro-elastic behavior is fully characterized from the fluctuations computed during molecular simulation of amorphous cellulose. Some peculiarities of moisture-induced behavior of cellulose are discussed, such as the high Biot coefficient, the negative Biot modulus, the negative thermal expansion, or the high heat capacity. This work investigates only the case of a pure fluid in the pore with volumetric deformations, but can be easily extended to fluid mixture and shearing effects.
## INTRODUCTION

Statistical physics bridges the atomistic and continuum descriptions of matter (Landau and Lifshitz 1980). The atomistic description of matter is probabilistic and therefore physical observables are fluctuating around their equilibrium values. Although, the relative fluctuations vanish in the limit of large macroscopic systems (thermodynamic limit), they become significant at the micron-scale and are major at the nanoscale. Historically, these fluctuations were first evidenced by the well-known brownian motion (Einstein 1905). Nanoscale systems are nowadays commonly studied, and fluctuations are ubiquitous in these studies. In this work, we focus on the use of molecular simulations applied to porous media. The presence of fluctuations is detrimental to the accuracy of molecular simulation results, since the larger the fluctuations $\Delta A$ of a quantity $A$ are, the less accurate the estimation of the statistical average $\langle A \rangle$ is and the longer it takes to estimate $\langle A \rangle$ with a given accuracy (central limit theorem). But the fluctuations can also be turned to an advantage, to measure additional thermodynamic properties. If one focuses on the fluctuations associated to the thermodynamic state parameters describing a system, statistical physics shows that the magnitude of these fluctuations is related to second order derivatives of the thermodynamic potential minimum at equilibrium (Landau and Lifshitz 1980). For example, in the case of a pure fluid in the isobaric-isothermal ensemble, the second order derivatives of the Gibbs free enthalpy quantify the compressibility (second derivative with respect to pressure), the thermal expansion (cross derivative with respect to pressure and temperature) and the heat capacity (second derivative with respect to temperature). The knowledge of these second order derivatives provides a complete description of the elasticity and coupling with respect to any possible evolution of the system (complete constitutive behavior). Here, 'elasticity' refers to all the thermodynamic derivatives between conjugated and state variables, e.g., the stiffness tensor as the derivative of stresses with respect to strains for a solid. Accordingly, combining the statistical averages (first order derivatives) and the magnitude of statistical fluctuations (second order derivatives) offers a route to fully characterize the thermodynamics of a system in a single molecular simulation. In practice, more accurate estimates of the second order derivatives are often obtained by performing
finite difference of first order derivatives, i.e., by repeating the molecular simulation with a small change of a control parameter. But, the approach by fluctuations becomes interesting when many state parameters are needed to describe the system, since the number of simulations to repeat for the finite difference approach rapidly becomes prohibitive. Moreover, the approach by fluctuations is also interesting when first order derivatives are poorly estimated by molecular simulation (e.g., chemical potential), or when the statistical ensemble of interest is not suited to the estimation by finite difference (see the case of the heat of adsorption hereafter).

Many fluctuation formulas are well-established for standard systems (Allen and Tildesley 1989), e.g., for pure fluids or non-porous solids. In this paper, we focus on drained porous media. For undrained porous media, fluctuation formulas are the same as for non-porous solids, and are frequently encountered in the literature. For drained porous media however, fluctuation formulas are seldom used, only the isosteric heat of adsorption is frequently estimated from fluctuations. Following the IUPAC, the isosteric heat of adsorption $\Delta q$ is defined as the heat released per adsorbed molecule at a given temperature $T$, and (reservoir) fluid pressure $P^b$ (Nič et al. 2009): $\Delta q = q^b - q$ with $q = \frac{\partial H}{\partial N}|_{T,P^b}$ the incremental molar enthalpy ($q$ and $q^b$ refer to the fluid adsorbed (inside the pores) and to the fluid bulk (reservoir), respectively). The derivative defining $q$ is at constant pressure, which is unambiguous for a bulk fluid, but may appear unclear for a porous medium since the mechanical stress $\sigma$ on the medium can differ from the fluid pressure. The appropriate thermodynamic condition is the so-called unjacketed condition ($\sigma = -P^b$) which corresponds to a porous medium simply immersed in the fluid reservoir (with no other mechanical constraint). Alternatively, $q$ can also be related to the change of chemical potential $\mu$ with temperature, thanks to a Maxwell relation of the Gibbs free enthalpy: $q = \frac{\partial(\mu/T)}{\partial(1/T)}|_{N,P^b}$. The quantity $q$ for the adsorbed fluid is generally obtained by the cross fluctuations of internal energy $E$ and number of fluid molecules $N$ in the grand canonical (GC) ensemble (Nicholson et al. 1982):

$$q = \frac{\langle \Delta E \Delta N \rangle_{GC}}{\langle \Delta N^2 \rangle_{GC}} \qquad (1)$$

It is well recognized that this fluctuation formula is valid for a rigid adsorbent only, so that
the incremental molar enthalpy at constant pressure coincides with the incremental internal energy at constant volume: $q \approx \frac{\partial E}{\partial N}|_{V,T}$ and one actually measures the differential heat of adsorption (Karavias and Myers 1991). Accounting for the deformability of the solid requires a consistent thermodynamics description of all thermo-poro-mechanical couplings between the fluid and the solid. Such an approach is proposed in this paper, and, in particular, we shall see that one can generalize the fluctuation formula (1) to account for the solid deformability. Interestingly, in this example of the heat of adsorption, the use of the fluctuation formula (1) is often preferred to a finite difference approach, because one generally prefers studying adsorption in the grand canonical ensemble whereas the finite difference approach would require a canonical simulation.

Porous media are fluid-solid mixtures that are submitted to changes both in fluid state parameters (reservoir bulk pressure $P^b$, temperature $T$) and in solid state parameters (deformations, temperature). Excluding the deviatoric deformations, which are often considered uncoupled with the other state parameters for isotropic media, the description of the deformation reduces to the volumetric strain $\varepsilon$ (Coussy 2003). A theory that encompasses any possible evolution is the Biot-Coussy theory of poromechanics (Biot 1941; Coussy 2010). For a system with 3 control parameters $\varepsilon$, $P^b$, and $T$, one needs a total of 6 independent moduli to fully describe the constitutive behavior. Following the Biot-Coussy theory of poromechanics, the framework of non-linear poroelasticity defines the 6 independent moduli as the moduli relating the changes in the control parameters $\varepsilon$, $P^b$, and $T$ to the changes of their conjugated quantities $\sigma$ (spherical stress), $\phi$ (porosity), and $s$ (entropy):

$$ \left\{ \begin{aligned} d\sigma &= K d\varepsilon - b dP^b - K\alpha dT \\ d\phi &= b d\varepsilon + \frac{dP^b}{N} - \alpha_{\phi} dT \\ ds &= K\alpha d\varepsilon - \alpha_{\phi} dP^b + \frac{c_v}{T} dT \end{aligned} \right. \qquad (2) $$

where $K$ is the bulk modulus, $b$ is the Biot coefficient, $N$ is the Biot modulus, $\alpha$ is the volumetric thermal expansion, $\alpha_{\phi}$ is the pore thermal expansion, and $c_v$ is the volumetric heat capacity at constant volume.
The purpose of this paper is to derive fluctuation formulas to estimate all the poromechan-
ical moduli. Doing so, we will extend the poromechanical moduli to adsorbing media (usual
poromechanics neglects adsorption), and the fluctuation formulas will apply to those generalized
moduli. This will make it possible to capture the mechano-sorption couplings in addition to the
usual thermo-poro-mechanical couplings. In particular, we will discuss the adsorption-induced
swelling in terms of fluctuations, and we will revisit the concept of heat of adsorption. Finally,
we illustrate the use of these fluctuation formulas by considering the case of amorphous cellulose
exposed to water adsorption, a typical example of material exhibiting mechano-sorption couplings.
The present work is limited to volumetric deformations, but extension to deviatoric deformations
can be achieved likewise.

In the first section, we recall how fluctuation formulas are derived in statistical physics and we
illustrate with the case of a bulk fluid. In particular, we formulate the constitutive behavior of a
bulk fluid in the grand canonical ensemble, which is the appropriate ensemble for a pore fluid in
a drained porous medium. In the second section, we generalize the usual constitutive behavior of
porochemanics (Eqs. 2) to adsorbing media by combining the descriptions of a fluid phase in the
grand canonical ensemble and of a solid phase in the canonical ensemble. We then derive fluctuation
formulas for all the poromechanical moduli. The implications of these formulas is discussed, and
we revisit the concept of heat of adsorption for flexible adsorbent. Finally, the last section is
dedicated to an example of application to amorphous cellulose exposed to water adsorption. The
cellulose is simulated in the osmotic ensemble and its full thermo-hydro-mechanical constitutive
behavior is derived from appropriate fluctuations. This example illustrates the peculiarity of
adsorption-induced couplings, and shows that it can be captured by the fluctuation formulas.

FLUID PHASE

In this section, we first recall the derivation of fluctuation formulas for a pure bulk fluid. While
the most usual formulation considers the fluid in the isobaric-isothermal ensemble, this ensemble
is not adapted to the case of pore fluid in a drained porous medium. The appropriate ensemble that
will be needed to address the drained porous medium is the grand canonical ensemble. Therefore,
in this section, we reformulate the fluid constitutive behavior in the grand canonical ensemble, and
derive the corresponding fluctuation formulas.

Fluctuations in statistical physics

In this section, we recall how to relate the fluctuations and second-order derivatives in statistical physics (Allen and Tildesley 1989). Let us consider an extensive state variable *F* of the system described in a statistical ensemble described by a partition function *Z*<sub>*F*</sub>. The intensive variable *f* conjugated to *F* is defined by *f* = − ∂ψ<sub>*F*</sub>/∂*F*|<sub>*X*</sub>, where ψ<sub>*F*</sub> = −ln (*Z*<sub>*F*</sub>) is the thermodynamic potential minimum at equilibrium, and *X* stands for the other control variables of the system. One can consider the system in an alternative statistical ensemble in which one controls *f* instead of *F* and the associated partition function is *Z*<sub>*f*</sub> = Σ<sub>*F*</sub> *Z*<sub>*F*</sub> exp(−*Ff*), where the sum is over all the possible values of *F*. In this new ensemble, *F* is no more constant but fluctuates, and its average value can be expressed as the derivative *F* = ∂ψ<sub>*f*</sub>/∂*f*|<sub>*X*</sub>, with ψ<sub>*f*</sub> = −ln (*Z*<sub>*f*</sub>) = min<sub>*F*</sub> (ψ<sub>*F*</sub> + *fF*) the thermodynamic potential minimum at equilibrium in the new ensemble. The statistical fluctuations of any physical observable depend on the statistical ensemble we consider. For two physical quantities of interest ℜ<sub>1</sub> and ℜ<sub>2</sub>, the general relation between the cross fluctuations in both ensembles is given by (Lebowitz et al. 1967):

$$
\langle \Delta \mathcal{A}_1 \Delta \mathcal{A}_2 \rangle_{F,X} = \langle \Delta \mathcal{A}_1 \Delta \mathcal{A}_2 \rangle_{f,X} + \left. \frac{\partial f}{\partial F} \right|_X \frac{\partial}{\partial f} \left[ \langle \mathcal{A}_1 \rangle_{f,X} \right]_X \frac{\partial}{\partial f} \left[ \langle \mathcal{A}_2 \rangle_{f,X} \right]_X \quad (3)
$$

where 〈·〉<sub>F,X</sub> and 〈·〉<sub>f,X</sub> denotes the ensemble averages in the ensembles at controlled *F* and *f*, respectively, and Δ*Y* = *Y* − 〈*Y*〉 represents the deviation from the average. A particular case of interest is when the physical quantities *A*<sub>*i*</sub> are *F*. Indeed, since Δ*F* is zero in the ensemble at controlled *F*, equation (3) reduces to:

$$
\langle \Delta F^2 \rangle_{f,X} = - \left. \frac{\partial F}{\partial f} \right|_X = \left. \frac{\partial^2 \psi_f}{\partial f^2} \right|_X \qquad (4)
$$

Therefore, the fluctuations of *F* in the ensemble at controlled *f* quantify the elasticity of the
system with respect to *f*. One can easily extend this relation to the case of two extensive state
parameters $F_1$ and $F_2$, conjugated to $f_1$ and $f_2$. The cross fluctuations of $F_1$ and $F_2$ in the ensemble at controlled $f_1$ and $f_2$ quantify the coupling moduli between the two parameters.

$$ \langle \Delta F_1 \Delta F_2 \rangle_{f_1, f_2, X} = - \left. \frac{\partial F_1}{\partial f_2} \right|_{f_1, X} = - \left. \frac{\partial F_2}{\partial f_1} \right|_{f_2, X} = \left. \frac{\partial^2 \psi_{f_1, f_2}}{\partial f_1 \partial f_2} \right|_X \quad (5) $$

where the thermodynamic potential $\psi_{f_1, f_2}$ minimum at equilibrium at controlled $f_1$ and $f_2$ is obtained as a Legendre transform of the potential $\psi_{F_1, F_2}$ minimum at controlled $F_1$ and $F_2$: $\psi_{f_1, f_2} = \min_{F_1, F_2} (\psi_{F_1, F_2} + f_1 F_1 + f_2 F_2)$.

## Thermodynamic description of a bulk fluid

Before considering the case of a porous medium, one needs first to formulate the constitutive behavior in a suitable thermodynamic ensemble. Let us consider a pure fluid. A common formulation of the constitutive behavior of a pure fluid is obtained by considering the isobaric-isothermal ensemble in which the control parameters are the number of molecules $N$, the bulk pressure $P$, and the temperature $T$. The thermodynamic potential minimum at equilibrium in this ensemble is the Gibbs free enthalpy $G = E - TS + PV$, where $E$ is the internal energy, $S$ the entropy, and $V$ the volume. And the associated energy balance is $dG = VdP + \mu dN - SdT$. Because of extensivity, the molar Gibbs free enthalpy is equal to the chemical potential $G/N = \mu$, so that the energy balance reduces to only two control parameters:

$$ d\left(\frac{G}{N}\right) = d\mu = \frac{V}{N}dP - \frac{S}{N}dT \quad (\text{Gibbs-Duhem equation}) \quad (6) $$

Accordingly, the complete constitutive behavior can be formulated with only two constitutive equations relating the changes of molar volume $V/N$ and molar entropy $S/N$ to the changes in $P$ and $T$. Accounting for the Maxwell relation $\left. \frac{\partial(V/N)}{\partial T} \right|_P = - \left. \frac{\partial(S/N)}{\partial P} \right|_T$, 3 thermo-mechanical moduli are enough to fully describe the constitutive behavior. Those moduli are the compressibility $\kappa = - \left. \frac{\partial(V/N)}{\partial P} \right|_T$, the thermal expansion $\alpha = \left. \frac{N}{V} \frac{\partial(V/N)}{\partial T} \right|_P = - \left. \frac{\partial(S/N)}{\partial P} \right|_T$, and the volumetric heat capacity at constant pressure $c_p = T \left. \frac{N}{V} \frac{\partial(S/N)}{\partial T} \right|_P$. The constitutive behavior of a pure fluid takes the form:
$$ \begin{cases} d\left(\frac{V}{N}\right) = -\frac{V}{N}\kappa dP + \frac{V}{N}\alpha dT \\ d\left(\frac{S}{N}\right) = -\frac{V}{N}\alpha dP + \frac{V}{N}\frac{c_p}{T} dT \end{cases} \qquad (7) $$

This formulation is not adapted to derive fluctuation formulas. Indeed, one should consider the thermodynamic potential corresponding to the isobaric-isothermal ensemble as it is simulated in molecular simulation (ensemble where the fluctuations will be computed eventually). The appropriate potential is the one corresponding to the logarithm of the partition function ($\psi = -\ln(Z)$). This potential is obtained as a Legendre transform of the entropy $S(-S)$ is minimum for isolated systems in the micro-canonical ensemble): $\psi_{N,\beta P,\beta} = \min_{E,V} \left(-\frac{S}{k} + \beta E + \beta PV\right) = \beta G$, with $k$ the Boltzmann constant and $\beta = \frac{1}{kT}$ the statistical temperature. The energy balance of this thermodynamic potential takes the form: $d\psi_{N,\beta P,\beta} = Ed\beta + Vd(\beta P) - (-\beta\mu) dN$ so that the pairs of conjugated quantities are now $(E, \beta)$, $(V, \beta P)$, and $(N, -\beta\mu)$. As for $G$, in the particular case of a fluid, the potential per molecule reduces to $\psi_{N,\beta P,\beta}/N = \beta\mu$ so that two control parameters are enough to describe the evolution of the system: $d\left(\frac{\psi_{N,\beta P,\beta}}{N}\right) = \frac{E}{N}d\beta + \frac{V}{N}d(\beta P)$. Rearranging equations (7), we obtain the constitutive behavior in this ensemble:

$$ \begin{cases} d\left(\frac{V}{N}\right) = -\frac{\kappa V}{\beta N} d(\beta P) + \frac{V}{\beta N} \left(\kappa P - \frac{\alpha}{k\beta}\right) d\beta \\ d\left(\frac{E}{N}\right) = \frac{V}{\beta N} \left(\kappa P - \frac{\alpha}{k\beta}\right) d(\beta P) - \frac{V}{\beta N} \left(\frac{c_p}{k\beta} - \frac{1}{\kappa} \left(\frac{\alpha}{k\beta}\right)^2 + \frac{1}{\kappa} \left(\kappa P - \frac{\alpha}{k\beta}\right)^2\right) d\beta \end{cases} \qquad (8) $$

In the isobaric-isothermal ensemble, the extensive variables that are fluctuating are $V$ and $E$, and their conjugated quantities are $\beta P$ and $\beta$. Applying the fluctuation formulas (4) and (5), we have:

$$ \left\{ \begin{aligned} \langle \Delta V^2 \rangle_{N,\beta P,\beta} &= - \left. \frac{\partial V}{\partial(\beta P)} \right|_{N,\beta} = \left. \frac{\partial^2 \psi_{N,\beta P,\beta}}{\partial(\beta P)^2} \right|_{N,\beta} \\ \langle \Delta E^2 \rangle_{N,\beta P,\beta} &= - \left. \frac{\partial E}{\partial \beta} \right|_{N,\beta P} = \left. \frac{\partial^2 \psi_{N,\beta P,\beta}}{\partial \beta^2} \right|_{N,\beta} \\ \langle \Delta E \Delta V \rangle_{N,\beta P,\beta} &= - \left. \frac{\partial E}{\partial(\beta P)} \right|_{\beta} = - \left. \frac{\partial V}{\partial \beta} \right|_{\beta P} = \left. \frac{\partial^2 \psi_{N,\beta P,\beta}}{\partial(\beta P) \partial \beta} \right|_{N,\beta} \end{aligned} \right. \qquad (9) $$

Developing $d\left(\frac{V}{N}\right)$ and $d\left(\frac{E}{N}\right)$ in equations (8), we readily identify these derivatives in terms of
the fluid thermo-mechanical moduli $\kappa$, $\alpha$, and $c_p$:

$$ \left\{ \begin{aligned} \langle \Delta V^2 \rangle_{N,\beta P,\beta} &= \frac{V\kappa}{\beta} \\ \langle \Delta E^2 \rangle_{N,\beta P,\beta} &= \frac{V}{\beta} \left( \frac{c_p}{k\beta} - \frac{1}{\kappa} \left( \frac{\alpha}{k\beta} \right)^2 + \frac{1}{\kappa} \left( \kappa P - \frac{\alpha}{k\beta} \right)^2 \right) \\ \langle \Delta E \Delta V \rangle_{N,\beta P,\beta} &= -\frac{V}{\beta} \left( \kappa P - \frac{\alpha}{k\beta} \right) \end{aligned} \right. \qquad (10) $$

Replacing the energy $E$ by the enthalpy $H = E + PV$, one can decouple the contributions of $\kappa$, $\alpha$, and $c_p$, and we find that:

* The fluctuations of volume quantify the compressibility: $\langle \Delta V^2 \rangle_{N,\beta P,\beta} = \frac{V}{\beta} \kappa = kT\kappa$

* The fluctuations of enthalpy quantify the heat capacity at constant pressure: $\langle \Delta H^2 \rangle_{N,\beta P,\beta} = \frac{V}{\beta} \frac{c_p}{k\beta} = kT^2 V c_p$

* The cross fluctuations of enthalpy and volume quantify the thermal expansion: $\langle \Delta H \Delta V \rangle_{N,\beta P,\beta} = \frac{V}{\beta} \frac{\alpha}{k\beta} = kT^2 V \alpha$

## Bulk fluid in the grand canonical ensemble

The description of the fluid in the isobaric-isothermal ensemble is not adapted for the study of a porous medium. Indeed, one cannot control the pressure of the fluid inside the pores. One can only control the pressure of an external reservoir in osmotic equilibrium (same $\mu$ and $T$) with the fluid inside the porous medium, which corresponds to drained conditions. Therefore, the appropriate control variables to describe the fluid phase in a drained porous medium are the chemical potential $\mu$ and the temperature $T$, which correspond to the grand canonical ensemble. The thermodynamic potential minimum at equilibrium in the grand canonical ensemble is the grand potential $\Lambda = E - TS - \mu N$ of energy balance $d\Lambda = -PdV - Nd\mu - SdT$. Again, thanks to the extensivity of the behavior of a fluid, the grand potential takes the simple form $\Lambda/V = -P$ so that the energy balance reduces to two control parameters $d(\frac{\Delta V}{V}) = d(-P) = -\frac{N}{V}d\mu - \frac{S}{V}dT$. Therefore, the constitutive behavior takes the form:
$$ \begin{cases} d\left(\frac{N}{V}\right) = \frac{N}{V}\kappa^{\mu}d\mu + \frac{N}{V}\alpha^{\mu}dT \\ d\left(\frac{S}{V}\right) = \frac{N}{V}\alpha^{\mu}d\mu + \frac{N}{V}\frac{c^{\mu}}{T}dT \end{cases} \qquad (11) $$

where the relevant moduli $\kappa^\mu$, $\alpha^\mu$, and $c^\mu$ in the grand canonical ensemble are readily identified with respect to the usual ones $\kappa$, $\alpha$, and $c^p$ by comparing equations (11) and (7):

* $\kappa_\mu = \frac{V}{N} \left. \frac{\partial(N/V)}{\partial \mu} \right|_T = \frac{N\kappa}{V}$ is the chemical compressibility

* $\alpha_\mu = \frac{V}{N} \left. \frac{\partial(N/V)}{\partial T} \right|_\mu = \frac{V}{N} \left. \frac{\partial(S/V)}{\partial \mu} \right|_T = \kappa_\mu \left( \frac{S}{N} - \frac{\alpha}{\kappa_\mu} \right)$ is the thermo-chemical coupling coefficient

* $c_\mu = \frac{TV}{N} \left. \frac{\partial(S/V)}{\partial T} \right|_\mu = \frac{V}{N} \left( c_p - \frac{T\alpha^2}{\kappa} \right) + \frac{T\alpha_\mu^2}{\kappa_\mu}$ is the molar heat capacity at constant chemical potential

As in the isobaric-isothermal case, the grand potential is not the appropriate potential to derive fluctuation formulas, since it does not correspond to the actual statistical ensemble that is sampled by molecular simulations. The appropriate potential associated with the grand canonical partition function is $\psi_{-\beta\mu,V,\beta} = \min_{E,N} \left(-\frac{S}{k} + \beta E - \beta\mu N\right) = \beta\Lambda$ and the associated energy balance is: $d\psi_{-\beta\mu,V,\beta} = Ed\beta - \beta P dV + Nd(-\beta\mu)$. Again, because of extensivity, $\psi_{-\beta\mu,V,\beta}/V = -\beta P \Rightarrow d\left(\frac{\psi_{-\beta\mu,V,\beta}}{V}\right) = \frac{E}{V}d\beta + \frac{N}{V}d(-\beta\mu)$, and the constitutive equations takes the form:

$$ \left\{ \begin{aligned} d\left(\frac{N}{V}\right) &= -\frac{\kappa_\mu N}{\beta V} d(-\beta\mu) - \frac{N}{\beta V} \left(\kappa_\mu\mu + \frac{\alpha_\mu}{k\beta}\right) d\beta \\ d\left(\frac{E}{V}\right) &= -\frac{N}{\beta V} \left(\kappa_\mu\mu + \frac{\alpha_\mu}{k\beta}\right) d(-\beta\mu) - \frac{N}{\beta V} \left(\frac{c_\mu}{k\beta} - \frac{1}{\kappa_\mu} \left(\frac{\alpha_\mu^2}{k\beta}\right)^2 + \frac{1}{\kappa_\mu} \left(\kappa_\mu\mu + \frac{\alpha_\mu}{k\beta}\right)^2\right) d\beta \end{aligned} \right. \qquad (12) $$

In the grand canonical ensemble, the fluctuating extensive variables are $N$ and $E$, and their conjugated quantities are $-\beta\mu$ and $\beta$, respectively. Applying the fluctuation formulas (4) and (5), we have:
234

$$ \left\{ \begin{aligned} \langle \Delta N^2 \rangle_{-\beta\mu, V, \beta} &= - \left. \frac{\partial N}{\partial(-\beta\mu)} \right|_{V, \beta} \\ \langle \Delta E^2 \rangle_{-\beta\mu, V, \beta} &= - \left. \frac{\partial E}{\partial\beta} \right|_{-\beta\mu, V} \\ \langle \Delta E \Delta N \rangle_{-\beta\mu, V, \beta} &= - \left. \frac{\partial E}{\partial(-\beta\mu)} \right|_{\beta} = - \left. \frac{\partial N}{\partial\beta} \right|_{-\beta\mu} \end{aligned} \right. \qquad (13) $$

235

236

Developing $d(\frac{N}{V})$ and $d(\frac{E}{V})$ in equations (12), we readily identify these derivatives in terms of the moduli $\kappa^\mu$, $\alpha^\mu$, and $c^\mu$:

237

$$ \left\{ \begin{aligned} \langle \Delta N^2 \rangle_{-\beta\mu, V, \beta} &= \frac{N\kappa_\mu}{\beta} \\ \langle \Delta E^2 \rangle_{-\beta\mu, V, \beta} &= \frac{N}{\beta} \left( \frac{c_\mu}{k\beta} - \frac{1}{\kappa_\mu} \left( \frac{\alpha_\mu^2}{k\beta} \right)^2 + \frac{1}{\kappa_\mu} \left( \kappa_\mu\mu + \frac{\alpha_\mu}{k\beta} \right)^2 \right) \\ \langle \Delta E \Delta N \rangle_{-\beta\mu, V, \beta} &= \frac{N}{\beta} \left( \kappa_\mu\mu + \frac{\alpha_\mu}{k\beta} \right) \end{aligned} \right. \qquad (14) $$

238

239

Replacing the energy $E$ by the Hill potential $L = E - \mu N$, one can decouple the contributions of $\kappa^\mu$, $\alpha^\mu$, and $c^\mu$, and we find that:

240

241

242

243

244

245

246

247

248

249

250

251

252

* The fluctuations of number of molecules quantify the chemical compressibility: $\langle \Delta N^2 \rangle_{-\beta\mu, V, \beta} = \frac{N}{\beta}\kappa_\mu = kTN\kappa_\mu$

* The fluctuations of Hill potential quantify the heat capacity at constant chemical potential: $\langle \Delta L^2 \rangle_{-\beta\mu, V, \beta} = \frac{N}{\beta} \frac{c_\mu}{k\beta} = kT^2 N c_\mu$

* The cross fluctuations of Hill potential and number of particle quantify the thermo-chemical coupling coefficient: $\langle \Delta L \Delta N \rangle_{-\beta\mu, V, \beta} = \frac{N}{\beta} \frac{\alpha_\mu}{k\beta} = kT^2 \alpha_\mu$

## ADSORBING POROUS MEDIUM

In this section, we now consider the case of a porous solid filled with a pure fluid. The usual theory of poromechanics assumes that the pore fluid is bulk, but this is not valid at the molecular scale since fluid adsorption is ubiquitous when the pore size is comparable to the range of inter-molecular interaction. Therefore, we first generalize the thermodynamic description of usual poromechanics to account for adsorption. Doing so, we make use of the constitutive behavior of the fluid phase derived here before. Next, we derive the fluctuation formulas for all the poromechanical moduli.
And we finish with a discussion, in which we revisit the concept of heat of adsorption when the solid is not assumed rigid.

## Thermodynamic description

A porous medium in drained conditions is a fluid-solid mixture for which the fluid is in osmotic equilibrium with an external reservoir, i.e., in the grand canonical ensemble. Under controlled deformation and temperature, the thermodynamic state parameters are the undeformed volume $V_0$, the deformation $\epsilon$ (here limited to the volumetric deformation), the fluid chemical potential $\mu$, and the temperature $T$. Since the behavior is extensive with respect to the undeformed volume $V_0$, one can adopt a thermodynamic description per unit undeformed volume. The thermodynamic potential minimum at equilibrium is a mix between the Helmholtz free energy of the solid and the grand potential of the fluid: $\omega = \frac{F-\mu N}{V_0}$, with $N$ the number of fluid molecules. The associated energy balance is:

$$d\omega = \sigma d\epsilon - nd\mu - sdT \qquad (15)$$

with $n = N/V_0$ and $s = S/V_0$. In usual poromechanics, the pore fluid is assumed to be in its bulk state (superscript $\square^b$ in what follows), so that, using Gibbs-Duhem equation (6), the energy balance can be rewritten:

$$d\omega = \sigma d\epsilon - \phi dP^b - s_s dT \qquad (16)$$

with $\phi$ the Lagrangian porosity (ratio of pore volume to $V_0$), $P^b$ the bulk fluid pressure, and $s_s = s - \phi s^b$ the entropy per unit undeformed volume of the solid only ($s^b$ is the entropy per unit volume of the bulk fluid). According to this energy balance, the pairs of conjugated state variables are $(\epsilon, \sigma)$, $(\phi, P^b)$, and $(s_s, T)$. Therefore, one can describe any evolution of the system by the usual constitutive equations of poromechanics (Eqs. 2). Interestingly, the energy balance (16) is identical to that of the solid only, submitted to the deformation $\epsilon$ at the external boundaries and to the pressure $P^b$ at the internal boundaries (pore walls). As a consequence, the potential $\omega$ depends
on the solid only, and all the drained poromechanical moduli ($K$, $b$, $N$, $\alpha$, $\alpha_{\phi}$, and $c$) depend on the solid only.

A major interest in using molecular simulation to study porous media is precisely to address the issue of fluid adsorption, since adsorption is the consequence of the molecular interactions between the fluid and the solid. Therefore, it is important to consider a thermodynamic framework that does not assume the pore fluid in its bulk state. In that case, the energy balance (16) does not hold. Coming back to the initial energy balance (15), one can still reformulate it in terms of the fluid bulk pressure with the Gibbs Duhem equation (6), but this pressure does not represent the pore pressure anymore, it represents the pressure of an external reservoir of fluid in osmotic equilibrium with the pore fluid. We have:

$$d\omega = \sigma d\varepsilon - \frac{n}{n^b} dP^b - \left(s - \frac{ns^b}{n^b}\right) dT \quad (17)$$

with $n^b$ the number density of the bulk fluid. According to this new energy balance, the bulk pressure of the fluid $P^b$ is conjugated to the quantity $n/n^b$. One can readily see that $n/n^b$ is the porosity $\phi$ if the density of the pore fluid is equal to the bulk density. Therefore, in the general case, the bulk pressure $P^b$ is conjugated to an effective porosity defined by $\phi^{eff} = n/n^b$, which coincides with the actual porosity only if the adsorbed and bulk densities are equal. Likewise, the temperature $T$ is conjugated to an effective entropy $s^{eff} = s - \frac{ns^b}{n^b}$, which coincides with the solid entropy $s_s = s - \phi s^b$ in the limit of a bulk fluid in the pores. Following this new energy balance, we can generalize the constitutive equations (2) of drained poroelasticity to adsorbing media as follows:

$$\left\{
\begin{aligned}
d\sigma &= K d\varepsilon - b dP^b - K\alpha dT \\
d\left(\frac{n}{n^b}\right) &= b d\varepsilon + \frac{dP^b}{N} - \alpha_{\phi} dT \\
d\left(s - \frac{ns^b}{n^b}\right) &= K\alpha d\varepsilon - \alpha_{\phi} dP^b + \frac{c_v}{T} dT
\end{aligned}
\right.
\quad (18)$$

In these generalized constitutive equations, the definition of the usual poromechanical moduli
is extended to the generalized pairs of conjugated variables $(\epsilon, \sigma)$, $(\phi^{eff}, P^b)$, and $(s^{eff}, T)$:

$$ K = \left. \frac{\partial \sigma}{\partial \epsilon} \right|_{P^b, T} ; b = - \left. \frac{\partial \sigma}{\partial P^b} \right|_{\epsilon, T} = \left. \frac{\partial (n/n^b)}{\partial \epsilon} \right|_{P^b, T} ; \alpha = - \left. \frac{1}{K} \frac{\partial \sigma}{\partial T} \right|_{\epsilon, P^b} = \left. \frac{1}{K} \frac{\partial (s-ns^b/n^b)}{\partial \epsilon} \right|_{P^b, T} \quad (19) \\ \frac{1}{N} = \left. \frac{\partial (n/n^b)}{\partial P^b} \right|_{\epsilon, T} ; \alpha_{\phi} = - \left. \frac{\partial (n/n^b)}{\partial T} \right|_{\epsilon, P^b} = - \left. \frac{\partial (s-ns^b/n^b)}{\partial P^b} \right|_{\epsilon, T} ; c_v = T \left. \frac{\partial (s-ns^b/n^b)}{\partial T} \right|_{\epsilon, P^b} $$

All these generalized moduli coincide with the usual poromechanical moduli in the limit of a bulk fluid in the pores. In the general case, the generalized moduli are not properties of the solid only and do depend on the nature of the fluid. In what follows, we derive fluctuation formulas to estimate those generalized moduli.

## Fluctuation formulas

To establish the fluctuation formulas for the porous medium, one has to rewrite the constitutive equations (18) with the appropriate pairs of conjugated variables corresponding to the statistical ensemble that is actually sampled by molecular simulation. Although one can derive fluctuation formulas in any statistical ensemble, a convenient choice to simulate a drained porous medium is the so-called osmotic ensemble in which one controls the undeformed volume $V_0$, the spherical stress $\sigma$, the fluid chemical potential $\mu$ and the temperature $T$ (note that the solid is at constant number of atoms and $V_0$ thus constitutes an extensive control parameter, which ensures that the osmotic ensemble is not pathological (Allen and Tildesley 1989)). The interest of this ensemble is that the evaluation of the fluctuating quantities (deformation, amount of fluid, and energy) is direct at almost no additional computational cost. Controlling the number of fluid molecules and measuring the chemical potential relies on the Widom insertion method, which is computationally quite inefficient for dense phases such as adsorbed fluids (Nezbeda and Kolafa 1991). Controlling the deformation and measuring the stress, is of rather moderate computational cost, but the fluctuation formulas for stress require to evaluate the hypervirial (Allen and Tildesley 1989) which involves the second derivative of the inter-atomic potential, and is most often not implemented in standard molecular simulation codes. In this respect, the osmotic ensemble appears as the most convenient choice to apply fluctuation formulas. The corresponding thermodynamic potential minimum at equilibrium
is $\omega_{os} = \omega - \sigma\epsilon$ of energy balance $d\omega_{os} = -\epsilon d\sigma - nd\mu - sdT$. Actual molecular simulations in the osmotic ensemble mix the isobaric-isothermal ensemble for the solid and the grand canonical ensemble for the fluid. The appropriate thermodynamic potential corresponding to the logarithm of the partition function and its energy balance, are given by:

$$ \psi_{-\beta\sigma, -\beta\mu, \beta} = \min_{V,N,E} \left( -\frac{S}{k} - \beta\sigma V - \beta\mu N + \beta E \right) = V_0\beta\omega_{os} \quad (20) $$

$$ d\psi_{-\beta\sigma, -\beta\mu, \beta} = Vd(-\beta\sigma) + Nd(-\beta\mu) + Ed\beta \quad (21) $$

Accordingly, the fluctuations observed in molecular simulations of the osmotic ensemble are fluctuations of volume $V = V_0\epsilon$, number of fluid molecules $N = V_0n$, and energy $E = V_0e$, which are conjugated to the statistical stress $-\beta\sigma$, the statistical fluid chemical potential $-\beta\mu$, and the statistical temperature $\beta$, respectively. Application of the fluctuation formulas (4) and (5) gives:

$$ \left\{ \begin{aligned} \langle \Delta V^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \left. \frac{\partial V}{\partial(-\beta\sigma)} \right|_{-\beta\mu, \beta} \\ \langle \Delta N^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \left. \frac{\partial N}{\partial(-\beta\mu)} \right|_{-\beta\sigma, \beta} \\ \langle \Delta E^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \left. \frac{\partial E}{\partial\beta} \right|_{-\beta\sigma, -\beta\mu} \\ \langle \Delta V \Delta N \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \left. \frac{\partial V}{\partial(-\beta\mu)} \right|_{-\beta\sigma, \beta} = - \left. \frac{\partial N}{\partial(-\beta\sigma)} \right|_{-\beta\mu, \beta} \\ \langle \Delta V \Delta E \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \left. \frac{\partial V}{\partial\beta} \right|_{-\beta\sigma, -\beta\mu} = - \left. \frac{\partial E}{\partial(-\beta\sigma)} \right|_{-\beta\mu, \beta} \\ \langle \Delta N \Delta E \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \left. \frac{\partial N}{\partial\beta} \right|_{-\beta\sigma, -\beta\mu} = - \left. \frac{\partial E}{\partial(-\beta\mu)} \right|_{-\beta\sigma, \beta} \end{aligned} \right. \quad (22) $$

To derive the expressions of these derivatives in terms of the poromechanical moduli, one must rewrite the constitutive equations of poromechanics (18) with the set of conjugated variable $(V, -\beta\sigma)$, $(N, -\beta\mu)$, and $(E, \beta)$. To do so, one needs to express the changes in bulk fluid properties $(P^b, 1/n^b$, and $s^b/n^b$) in function of $-\beta\mu$ and $\beta$. Using the constitutive behavior of the bulk fluid in the grand canonical ensemble (Eqs. 12), one obtains:
$$ egin{cases} dP^b = -\frac{n^b}{\beta}d(-\beta\mu) - \frac{h^b}{\beta}d\beta \\ d\left(\frac{1}{n^b}\right) = \frac{1}{K^b\beta}d(-\beta\mu) + \frac{1}{\beta K^b n^b}\left(h^b - \frac{K^b\alpha^b}{k\beta}\right)d\beta \\ d\left(\frac{s^b}{n^b}\right) = \frac{\alpha^b}{\beta}d(-\beta\mu) - \frac{1}{\beta n^b}\left(c_v^b - \alpha^b\left(h^b - \frac{K^b\alpha^b}{k\beta}\right)\right)d\beta \end{cases} \quad (23) $$

where $K^b = 1/\kappa^b$ is the fluid bulk modulus, $h^b = e^b + P^b = \frac{s^b}{k\beta} + n^b\mu$ is the volumetric bulk enthalpy, and $c_v^b = c_p^b - \frac{T(\alpha^b)^2}{\kappa^b}$ is the volumetric heat capacity at constant volume. Combining these expression with the constitutive equations (18) of poromechanics and the energy balance of internal energy $de = \sigma d\epsilon + \frac{1}{k\beta}ds + \mu dn$, we obtain:

$$ egin{cases} dV = -\frac{V_0}{\beta K}d(-\beta\sigma) - \frac{V_0 b n^b}{\beta K}d(-\beta\mu) - \frac{V_0}{\beta K}\left(\frac{K\alpha}{k\beta} + bh^b + \sigma\right)d\beta \\ dN = -\frac{V_0 b n^b}{\beta K}d(-\beta\sigma) - \frac{V_0 (n^b)^2}{\beta}\left(\frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b}\right)d(-\beta\mu) \\ + \frac{V_0 n^b}{\beta}\left(\frac{1}{k\beta}\left(\alpha_\phi - b\alpha + \frac{n}{n^b}\alpha^b\right) - \left(\frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b}\right)h^b - \frac{b}{K}\sigma\right)d\beta \\ dE = -\frac{V_0}{\beta K}\left(\frac{K\alpha}{k\beta} + bh^b + \sigma\right)d(-\beta\sigma) \\ + \frac{V_0 n^b}{\beta}\left(\frac{1}{k\beta}\left(\alpha_\phi - b\alpha + \frac{n}{n^b}\alpha^b\right) - \left(\frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b}\right)h^b - \frac{b}{K}\sigma\right)d(-\beta\mu) \\ - \frac{V_0}{\beta}\left(\frac{1}{k\beta}\left(c_v + \frac{n}{n^b}c_v^b\right) - N\left(\frac{\alpha_\phi}{k\beta}\right)^2 + N\left(\frac{\alpha_\phi}{k\beta} - \frac{h^b}{N}\right)^2\right) \\ + \frac{1}{K}\left(\frac{K\alpha}{k\beta} + bh^b + \sigma\right)^2 + \frac{n}{K^b n^b}\left(\frac{K^b\alpha^b}{k\beta} - h^b\right)^2 \end{cases} \quad (24) $$

One can then readily identify the expressions of the derivatives in the fluctuation formulas (22):
$$
\left\{
\begin{aligned}
\langle \Delta V^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} &= \frac{V_0}{\beta K} \\
\langle \Delta N^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} &= \frac{V_0 (n^b)^2}{\beta} \left( \frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b} \right) \\
\langle \Delta E^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} &= \frac{V_0}{\beta} \left( \frac{1}{k\beta} \left( c_v + \frac{n}{n^b} c_v^b \right) - \mathcal{N} \left( \frac{\alpha_\phi}{k\beta} \right)^2 + \mathcal{N} \left( \frac{\alpha_\phi}{k\beta} - \frac{h^b}{N} \right)^2 \right. \\
&\qquad \left. + \frac{1}{K} \left( \frac{K\alpha}{k\beta} + b h^b + \sigma \right)^2 + \frac{n}{K^b n^b} \left( \frac{K^b \alpha^b}{k\beta} - h^b \right)^2 \right) \\
\langle \Delta V \Delta N \rangle_{-\beta\sigma, -\beta\mu, \beta} &= \frac{V_0 b n^b}{\beta K} \\
\langle \Delta V \Delta E \rangle_{-\beta\sigma, -\beta\mu, \beta} &= \frac{V_0}{\beta K} \left( \frac{K\alpha}{k\beta} + b h^b + \sigma \right) \\
\langle \Delta N \Delta E \rangle_{-\beta\sigma, -\beta\mu, \beta} &= - \frac{V_0 n^b}{\beta} \left( \frac{1}{k\beta} \left( \alpha_\phi - b\alpha + \frac{n}{n^b} \alpha^b \right) - \left( \frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b} \right) h^b - \frac{b}{K} \sigma \right)
\end{aligned}
\right.
\quad (25)
$$

As we did with the enthalpy and the Hill potential for a bulk fluid, we can decouple the different poromechanical moduli, by introducing the fluctuations of suitable quantities. We find that:

• The bulk modulus $K$ is quantified by the fluctuations of $V$:

$$ \frac{1}{K} = \frac{\beta}{V_0} \langle \Delta V^2 \rangle_{-\beta\sigma, -\beta\mu, \beta} \quad (26) $$

• The unjacketed compressibility $\frac{1-b}{K}$ is quantified by the cross fluctuations of $V$ and $V_s = V - \frac{N}{n^b}$ (solid volume in the limit of a bulk fluid in the pores):

$$ \frac{1-b}{K} = \frac{\beta}{V_0} \left\langle \Delta V \Delta \left( V - \frac{N}{n^b} \right) \right\rangle_{-\beta\sigma, -\beta\mu, \beta} \quad (27) $$

• The thermal expansion $\alpha$ is quantified by the cross fluctuations of $V$ and $H_s = E - \sigma V - \frac{h^b}{n^b} N$ (solid enthalpy in the limit of a bulk fluid in the pores)

$$ \frac{\alpha}{k\beta} = \frac{\beta}{V_0} \left\langle \Delta V \Delta \left( E - \sigma V - \frac{h^b}{n^b} N \right) \right\rangle_{-\beta\sigma, -\beta\mu, \beta} \quad (28) $$
• The undrained Biot modulus $\frac{1}{M} = \frac{1}{N} + \frac{1}{K^b n^{b/2}}$ is quantified by the fluctuations of $V_p = \frac{N}{n^b} - bV$ (in usual poromechanics this would correspond to $V_0 (\phi - \phi_0 - b\epsilon)$ to an additive constant which does not contribute to the fluctuations. This quantity corresponds to the change of pore volume corrected by the change of volume due to deformation only):

$$ \frac{1}{M} = \frac{\beta}{V_0} \left\langle \Delta \left( \frac{N}{n^b} - bV \right)^2 \right\rangle_{-\beta\sigma, -\beta\mu, \beta} \qquad (29) $$

• The sum of the pore thermal expansion $\alpha_\phi$ and of the fluid thermal expansion $\alpha^b$ is quantified by the cross fluctuations of $V_p = bV - \frac{N}{n^b}$ and $H_s = E - \sigma V - \frac{h^b}{n^b}N$

$$ \frac{\alpha_{\phi}}{k\beta} + \frac{n}{n^b}\frac{\alpha^b}{k\beta} = \frac{\beta}{V_0} \left\langle \Delta \left( bV - \frac{N}{n^b} \right) \Delta \left( E - \sigma V - \frac{h^b}{n^b}N \right) \right\rangle_{-\beta\sigma, -\beta\mu, \beta} \qquad (30) $$

• The sum of the volumetric heat capacities at constant pressure of the porous solid $c_p = c_v + \frac{K\alpha^2}{k\beta}$ and of the fluid $c_p^b = c_v^b + \frac{K^b(\alpha^b)^2}{k\beta}$ is quantified by the fluctuations of $H_s = E - \sigma V - \frac{h^b}{n^b}N$

$$ \frac{c_p + \frac{n}{n^b} c_p^b}{k\beta} = \frac{\beta}{V_0} \left\langle \Delta \left( E - \sigma V - \frac{h^b}{n^b} N \right)^2 \right\rangle_{-\beta\sigma, -\beta\mu, \beta} \qquad (31) $$

## Discussion

In summary, the fluctuation formulas (26-31) make it possible to determine all the poromechanical moduli from fluctuations of appropriate quantities during a simulation in the osmotic ensemble. It is assumed that the properties of the bulk fluid (density $n^b$, molar enthalpy $h^b/n^b$, thermomechanical moduli $(K^b, \alpha^b, c_p^b)$) are known a priori, either from analytic equation of state, or from separate molecular simulations of the bulk fluid. All the quantities involved in the fluctuations $(N, V, E)$ can be evaluated unambiguously at a modest additional computational cost.

The derivation of the fluctuation formulas (26-31) involves no particular assumption about the state of the fluid in the medium (adsorbed or non-adsorbed) or about the medium itself. Therefore, it can be used to analyze unusual poromechanical couplings due to adsorption, such as the well-known adsorption-induced swelling in unjacketed conditions. Unjacketed conditions correspond
to the adsorbent immersed in the bulk fluid so that $\sigma = -P^b$. Under isothermal conditions, the deformation of the medium in unjacketed conditions is given by: $d\varepsilon = -\frac{1-b}{K}dP^b$. In usual poromechanics (bulk fluid in the pores), one can formally show that $0 < b < 1$ (Coussy 2010). Accordingly, unjacketed conditions are expected to lead to a contraction as the fluid pressure $P^b$ is increased ($d\varepsilon < 0$). Yet, swellings have long been observed for adsorbing media in unjacketed conditions (Brochard et al. 2012b; Gor et al. 2017; Vandamme 2019). In presence of adsorption, the restriction $0 < b < 1$ does not hold anymore. According to the fluctuation formula (27), the unjacketed modulus $-\frac{1-b}{K}$ can well turn positive (i.e., swelling and $b > 1$), because adsorption of a vapor generally leads to much higher density in the adsorbent than in the bulk phase (liquid-like density): the quantity $N/n^b$, which is the pore volume for a bulk fluid, can well dominate over the total volume $V$ in the fluctuation formula (27).

We can also use this approach to revisit the concept of isosteric heat of adsorption in the framework of a deformable porous medium. As mentioned in the introduction, the heat of adsorption is the difference of incremental molar enthalpy between the bulk fluid and the adsorbed fluid: $\Delta q = q^b - q$, with $q = \frac{\partial H}{\partial N}|_{T,\sigma=-P^b} = \frac{\partial(\frac{\mu}{T})}{\partial(\frac{1}{T})}|_{N,\sigma=-P^b}$. For the bulk fluid satisfying Gibbs-Duhem equation ($dP^b = -\frac{n^b}{\beta}d(-\beta\mu) - \frac{h^b}{\beta}d\beta$), one readily identifies $q^b$ as the bulk molar enthalpy: $q^b = -\frac{\partial(-\beta\mu)}{\partial\beta}|_{N,P^b} = \frac{h^b}{n^b}$. For the adsorbed fluid, $q$ is usually estimated from the cross fluctuations of internal energy and fluid amount (Eq. 1), but this is valid only in the limit of a rigid solid. Rewriting the definition of the quantity $q$ in terms of the appropriate control variables of the osmotic ensemble, we have:

$$ \left. \frac{\partial \left(\frac{\mu}{T}\right)}{\partial \left(\frac{1}{T}\right)} \right|_{N,\sigma} = - \left. \frac{\partial (-\beta\mu)}{\partial \beta} \right|_{-\beta\sigma,N} + \sigma \left. \frac{\partial (-\beta\mu)}{\partial (-\beta\sigma)} \right|_{N,\beta} \quad (32) $$

Using the second constitutive equation (24), we obtain:

$$ q = \frac{h^b}{n^b} - \frac{\frac{1}{k\beta} \left( \alpha_\phi - b\alpha + \frac{n}{n^b}\alpha^b \right)}{n^b \left( \frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b} \right)} \quad (33) $$
The isosteric heat of adsorption is given by the difference between the bulk and adsorbed fluids, that is:

$$ \Delta q = \frac{\frac{1}{k\beta} \left( \alpha_{\phi} - b\alpha + \frac{n}{n^b} \alpha^b \right)}{n^b \left( \frac{b^2}{K} + \frac{1}{N} + \frac{n}{K^b n^b} \right)} \quad (34) $$

Using the fluctuation formulas derived for the porous medium (Eqs. 25), we can express the isosteric heat $\Delta q$ in function of the fluctuations of $N$ and $H_s(\sigma = -P^b) = E + P^b V - \frac{h^b}{n^b} N$:

$$ q = \frac{\langle \Delta N \Delta (E + P^b V) \rangle_{-\beta\sigma, -\beta\mu, \beta}}{\langle \Delta N^2 \rangle_{-\beta\sigma, -\beta\mu, \beta}} \Rightarrow \Delta q = - \frac{\langle \Delta N \Delta (E + P^b V) \rangle_{-\beta\sigma, -\beta\mu, \beta}}{\langle \Delta N^2 \rangle_{-\beta\sigma, -\beta\mu, \beta}} + \frac{h^b}{n^b} \quad (35) $$

In the limit of a rigid solid (V constant), the osmotic ensemble coincides with the grand canonical ensemble, and we recover the usual fluctuation formula of the isosteric heat (Eq. 1). The fluctuation formula (Eq. 35) generalizes the usual one for deformable solids. Interestingly, the correction due to the solid deformability is simply to consider the cross fluctuations of the medium enthalpy $E + P^b V$ instead of the internal energy $E$. The expression (34) shows that the heat of adsorption is identically null only if $\frac{n}{n^b} \alpha^b = b\alpha - \alpha_{\phi}$. For a homogeneous solid matrix and a non-adsorbing fluid ($\alpha_{\phi} = (b - \phi)\alpha$ and $\frac{n}{n^b} = \phi$), the heat of adsorption is null only if the thermal expansion of the fluid matches that of the solid ($\alpha = \alpha^b$).

## APPLICATION TO AMORPHOUS CELLULOSE EXPOSED TO WATER ADSORPTION

In this section, we illustrate the use of the fluctuation formulas for a typical case of adsorption-sensitive material, namely amorphous cellulose subjected to moisture at various relative humidity. We first present the molecular models and simulation methods considered, along with the bulk water properties needed to apply the fluctuation formulas. Then, all the thermo-hydro-mechanical moduli are estimated from the appropriate fluctuations. The quantitative results are discussed, we highlight in particular the unusual effects of adsorption, and we investigate the impact of considering capillary water as the bulk water phase instead of water vapor. Finally, the unjacketed response as a function of relative humidity is predicted from the poromechanics and compared to the actual response obtained by molecular simulation.
Simulation methods

Cellulose, one of the major constituents of biomass, is well known for its sorption-induced
swelling when exposed to humidity. It is therefore a good example of application to illustrate
the use of the fluctuation formulas. We considered the molecular model of Chen et al. (2018),
which is elaborated from the compaction of several polymeric chains made of cellobiose. Pre-
cisely, we studied the 'Sample_1' system made publicly available by Chen et al. (2018) at the
ETH Research Collection repository: https://doi.org/10.3929/ethz-b-000272221. This
molecular model reproduces well the adsorption-induced swelling and water content of cellulose
as a function of relative humidity at ambient temperature and pressure. In Fig. 1, a few snapshots
illustrate the evolution of the system with relative humidity. Chen et al. (2018) provide detailed
information about the porosity of this cellulose model: the pore size distribution contains sub-
nanometric pores of mean size around 0.15 nm in the dry state (±0.05 nm), and the distribution
evolves significantly toward larger pore sizes as the RH is increased, reaching a mean pore size
of 0.2 nm at 100% RH (±0.1 nm). Accordingly, a significant part of the pore space is formed
during the adsorption and results from the disjoining pressure exerted by the water molecules on
the polymer chains. This quite peculiar type of adsorption is of special interest in the present work
because it results in very significant poromechanical couplings.

In our simulations, we used the OPLS-aa force field (Jorgensen et al. 1996) along with the rigid
SPC/E model of water (Berendsen et al. 1987). The system was simulated at controlled temperature
*T*, confining stress *σ*, and water chemical potential *μ*, by combining Configurational Bias Grand
Canonical Monte Carlo (CB-GCMC), and isobaric-isothermal Molecular Dynamics (NPT-MD).
The CB-GCMC follows the biasing scheme of Brochard (2021) adapted to water SPC/E, and the
NPT-MD considers a velocity-Verlet integration with a Nose-Hoover thermostat and barostat. The
barostat allows independent volumetric deformations in each axial directions, and no shear was
allowed. 30 CB-GCMC insertion-deletion trials were performed every 100 MD steps, and each
trial considered 100 possible water positions. We considered a MD time step of 1 fs, and damping
times of 100 fs for the thermostat and 1000 fs for the barostat. The systems were simulated over
a total of 10⁸ MC trials including 3 · 10⁷ for equilibration. All the simulations were performed
with the LAMMPS software (Thompson et al. 2022) (https://www.lammps.org) modified to
incorporate the CB-GCMC.

We simulated the cellulose hydration following the adsorption branch, i.e., starting from the
dry cellulose. For all simulations, we considered a temperature of 300K and zero confining stress
(all directions relaxed independently). A total of 27 relative humidities (RH) were considered.
The corresponding fugacities f = RH · fsat were used to fix the water chemical potential for the
CB-GCMC. fsat = 998.70 Pa is the reference fugacity at saturation, and was obtained by molecular
simulation of bulk water SPC/E to determine the liquid-vapor coexistence point (Brochard 2021).
Bulk water simulations also provide estimates of the bulk density nᵇ, pressure Pᵇ, bulk modulus
Kᵇ, specific enthalpy hᵇ/nᵇ, thermal expansion αᵇ and specific heat capacity c_pᵇ/nᵇ that are needed
to apply the fluctuation formulas. All these properties are listed in Table 1. While water vapor is the
stable bulk phase, for RH larger than 30% bulk water can also exist in the meta-stable liquid phase.
This is the capillary water that fills the smallest pores in hydrophilic porous media. Therefore,
above 30%, one can also apply the fluctuation formulas by considering this liquid water as the bulk
reference. Since the adsorbed water is much closer to a liquid phase than to a vapor phase, this is
expected to lead to values of poromechanical moduli much more similar to usual values (i.e., values
for non-adsorbing media). In this respect, we report in Table 1 both the bulk vapor and bulk liquid
properties. Note that the bulk vapor is almost an ideal gas at the considered temperature (300K).
Note also that it is essential to consider a rigid water model because internal vibrations are almost
inactive at room temperature due to energy quantification (Honorio and Brochard 2022).

Simulation results

The molecular simulations at the different *RH* readily provide the adsorption and swelling isotherms of Fig. 2. The isotherms exhibit an initial Langmuir-like adsorption at low *RH*, then a plateau around 30-70 % *RH*, followed by an inflection and sharper increase above 80% *RH*. Such type II isotherm is typical for amorphous cellulose (Kulasinski et al. 2014b; Chen et al. 2018). Confronting carefully the two isotherms one observes that the water content and swelling are not
exactly proportional. At low *RH* the swelling per mass of water adsorbed is smaller than at high *RH*. Again, this is typical of amorphous biopolymers (Kulasinski et al. 2015b). This difference between low *RH* response and high *RH* response has been studied in detail in the literature, and is attributed to a change in the energetics of water adsorption (Chen et al. 2018): dominated by water-cellulose hydrogen bonds at low *RH*, and by water-water hydrogen bonds at large *RH*. We display in Fig. 2, the water content isotherm because this is the most straightforward with respect to experiments. One could also consider the effective porosity *n*/*n*<sup>b</sup> instead, which is the relevant quantity with respect to thermodynamics (quantity conjugated to the bulk fluid pressure *P*<sup>b</sup>, see Eq. 18). The change of adsorption from low *RH* to high *RH* would be also visible in the effective adsorption isotherm (not shown here): while the bulk density *n*<sup>b</sup> is small (vapor) and increases linearly with *RH* (almost ideal gas) the adsorbed density *n* is large (liquid-like) and follows a type II isotherm non-linear with *RH* (Fig. 2). As a consequence, the effective porosity *n*/*n*<sup>b</sup> decreases in the low *RH* regime (∝ *RH*<sup>-0.6</sup>) because the adsorbed density *n* is sublinear in *RH* (concave), and increases in the high *RH* regime because *n* becomes superlinear (convex). The minimum effective porosity *n*/*n*<sup>b</sup> is reached near 70% *RH* which corresponds to the transition between the two adsorption regimes.

In Fig. 2, the adsorption and swelling isotherms are compared with the original results of Chen
et al. (2018), who proposed the cellulose model investigated here. For the water content, it was
possible to compare with the exact same cellulose model and the results are well consistent, except
at 100% RH, Chen et al. (2018) reporting larger water contents. For the swelling, Chen et al.
(2018) provide only the average strain for three different cellulose models, so that the comparison
is not as good but the overall trend is reasonably consistent, but, again, the swelling at 100% RH
is significantly larger. The present results are not as smooth as in Chen et al. (2018), suggesting a
lower accuracy. The typical errors of the present results are ±0.02 g/g in water content and ±0.07
in volumetric strain, so that the difference with Chen et al. (2018) is within the error bars, except for
100% RH. The present study considers a CB-GCMC method with a total of 10⁸ trials, while Chen
et al. (2018) have used a GCMC with 2.4 × 10⁸ trials. The present methodology would require
longer simulations to reach the same accuracy as Chen et al. (2018). As for the discrepancy at
100% RH, careful investigation of the system evolution does not show any equilibration issue: both
water content and swelling are well equilibrated for the last 7 × 10⁷ CB-GCMC trials. Accordingly,
the reason for the discrepancy at 100% RH is non-obvious. A possible reason could be the chosen
force field for the cellulose which is not the same (OPLS-aa vs. PCFF). Another possible reason
could be a difference in the value considered for the water SPC/E chemical potential at saturation,
since water content and deformation appear quite sensitive near 100% RH: a fugacity at saturation
of 998.7 Pa is considered in the present work (from Brochard (2021), consistent with the value
of 996.7 Pa derived from the activity at saturation reported in the NIST water SPC/E database,
https://www.nist.gov), but it is not mentioned in the work of Chen et al. (2018).

Let us now investigate the poromechanical moduli. The fluctuation formulas (26-31) were used
to provide estimates of all the poromechanical moduli from the fluctuations of number of water
molecules, energy, and volume. The results are displayed in Fig. 3. Some estimates exhibit much
more uncertainty than others depending on the moduli: The Biot coefficient (b), Biot modulus
(N), and pore thermal expansion (αₕ) seem less uncertain than the bulk modulus (K), thermal
expansion (α), and heat capacity (cₚ). If all fluctuations arise from the thermal agitation and
Boltzmann statistics, the accuracy of their estimation is highly system-dependent and it relies on
how well the phase space of molecular configurations is explored. In this respect, cellulose is
not an easy system as the polymer chains can adopt many conformations separated by significant
energy barriers. Moreover, some fluctuation estimates are intrinsically more precise than others.
For instance, fluctuations of volume turn out to be quite hard to estimate accurately (hence the poor
precision on bulk modulus K), whereas the cross-correlation between volume and number of fluid
molecules is much more precise (hence the better precision on the Biot coefficient b).

To the best of our knowledge, most of those poromechanical moduli have never been investigated
in the existing literature, to the exception of the bulk modulus (Kulasinski et al. 2014b; Kulasinski
et al. 2015a; Kulasinski et al. 2015b; Chen et al. 2018). Yet, those studies report the undrained bulk
modulus K_u = ∂σ/∂ε|_{n,T}, and not the drained bulk modulus K = ∂σ/∂ε|_{p^b,T}.
Combining the constitutive
equations (18) of the porous medium and (23) of the bulk fluid, one can relate the two bulk moduli according to: $K_u = K + Mb^2$ with $\frac{1}{M} = \frac{1}{N} + \frac{1}{Kb} \frac{n}{n^b}$ the undrained Biot modulus. In Fig. 3, we display both the drained and undrained bulk moduli. The undrained modulus is systematically about 2 GPa larger than the drained modulus and decreases from about 9 GPa in the dry state to about 6 GPa at 100% RH, which is rather consistent with the results of Chen et al. (2018) who considered the same molecular model of amorphous cellulose (from 9.5 GPa to 6.5 GPa).

For the other moduli, no results can be found in the literature to confront them, as they are not usually investigated. Yet, let us mention the work of Kulasinski et al. (2015a), which developed a poromechanical formulation adapted to amorphous cellulose, but at controlled $\sigma$ and $n$, conjugated to $\varepsilon$ and $\mu$. They also provide estimates of generalized poromechanical moduli but not the same as ours. Confrontation is theoretically possible by deriving relations like the one between drained and undrained moduli. But, the biopolymer model of Kulasinski et al. (2015a) seems much softer than the one we considered (Chen et al. 2018) with an undrained bulk modulus varying from 7 GPa to almost 2 GPa. Therefore, any quantitative comparison does not seem relevant.

## Discussion

Let us comment the quantitative values of the other poromechanical moduli. The values of Biot coefficients $b$ are decreasing from about $10^6$ at low RH to about $10^5$ at high RH. Such values are extremely high compared to that expected in usual poromechanics. For a bulk fluid in the pores, the confining pressure induced by a pore pressure is expected to be a fraction of the pore pressure, i.e., $0 \le b \le 1$, the upper boundary corresponding to the case of an incompressible solid phase. Here, the values are orders of magnitude larger than this limit because the pore pressure exerted by the fluid has nothing to do with the bulk fluid pressure. The bulk fluid pressure is in the range of the saturation pressure ~ kPa, whereas the pore fluid is in a liquid-like state, so with a magnitude of pore pressure similar to the magnitude of a capillary pressure, i.e., ~ $10^2$ MPa. Conversely, one can also consider the meta-stable liquid state as the reference bulk fluid. The Biot coefficients estimated with this alternative reference are displayed in Fig. 4. This time, the values are in the usual range between 0.6 and 1, which suggests that the adsorbed water pore pressure is of similar
magnitude as the capillary pressure.

The value of Biot modulus *N* are also unusual. They are mostly negative decreasing from 0 Pa to -0.05 Pa. Negative values are unexpected in usual poromechanics, because it would correspond to a mechanically unstable system. Yet, when the pore fluid is no more bulk, the stability criterion becomes much more complex and the usual stability condition does not apply anymore. A negative Biot modulus can be interpreted as the fact that the effective porosity *n*/*n*<sup>b</sup> decrease as the bulk water vapor pressure increases at fixed deformation of the cellulose. Otherwise said, the relative increase in pore water content, is smaller than the relative increase of the bulk density. As shown by Chen et al. (2018), the water content isotherm at constant deformation is Langmuir-like (concave) whereas the bulk density is almost linear (ideal gas). Therefore, it appears consistent to obtain negative Biot modulus (the ratio *n*/*n*<sup>b</sup> is always decreasing with *RH* ∝ *P*<sup>b</sup>). Since the Langmuir isotherm is almost linear at small *RH*, it is also expected to observe a Biot modulus that converges to zero in this limit (ratio *n*/*n*<sup>b</sup> almost constant).

The thermal expansion α is also negative and decreases linearly from near 0 in the dry state to about −3 · 10<sup>-3</sup> K<sup>-1</sup> at 100% RH. Not only the values are negative but 10<sup>-3</sup> K<sup>-1</sup> is a rather large magnitude for a thermal expansion: for most solids, it is a few 10<sup>-5</sup> K<sup>-1</sup>. The thermal expansion of dry amorphous cellulose is estimated in the literature around 0.06 · 10<sup>-3</sup> K<sup>-1</sup> (Kulasinski et al. 2014a). Given the intrinsic uncertainty of our estimations, our value of thermal expansion in the dry state seems consistent with the literature. In presence of moisture, the thermal expansion becomes strongly negative. This is again a consequence of the drainage of the pore water: as temperature is increased while maintaining a constant vapor pressure, water desorbs which induces a shrinkage. Note that maintaining a constant vapor pressure is not equivalent to maintaining a constant *RH* since the saturation pressure increases with temperature. Thus, a drained heating corresponds to a net decrease in *RH*. The magnitude of drained thermal contraction is much larger than usual thermal expansion of solids, because desorption-induced contraction is much larger than usual thermal effects.

The pore thermal expansion $\alpha_{\phi}$ exhibits values that decrease from about $10^4$ K$^{-1}$ in the dry
state to $10^2 \text{ K}^{-1}$ at 100% RH. Again these values are extremely high compared to usual values. In usual poromechanics with a bulk fluid in the pores and a homogeneous solid phase of thermal expansion $\alpha_s$, the pore thermal expansion takes the form $\alpha_{\phi} = \alpha_s (b - \phi)$. It quantifies the increase of porosity upon drained heating while the overall deformation is fixed. Its magnitude is expected to be comparable or smaller than the thermal expansion of the solid phase, i.e., about $10^{-5} \text{ K}^{-1}$.

The values obtained here are orders of magnitude larger. Again, one can understand such values in light of the adsorption phenomenon. $\alpha_{\phi}$ quantifies the increase of the effective porosity $n/n^b$ as temperature is increased while maintaining a fixed deformation and vapor pressure. Such heating is expected to lead to a desorption (decrease of $n$) but also to a decrease of the bulk density ($n^b$).

Highly positive values of $\alpha_{\phi}$ mean that the relative decrease of bulk density is much faster than that of the adsorbed water content. In a domain of temperature well below the characteristic temperature of desorption, this is indeed expected.

Finally, the heat capacity is in the range of 10-20 MJ.m$^{-3}$.K$^{-1}$. Such values are rather large, and stem to the desorption upon heating which comes with a latent heat of desorption that is comparable to a vaporization enthalpy. By considering the meta-stable liquid phase as the bulk water reference, we obtain the values of Fig. 4 which are much more usual (~ 4.5 MJ.m$^{-3}$.K$^{-1}$), since the desorption from the adsorbed state to the liquid state comes with a much smaller heat of desorption. The values obtained with a bulk liquid reference are very close to what can be estimated classically with Dulong-Petit law: $c_v = 3N_{solid}k/V$ and $c_p = c_v + KT\alpha^2$ (see curve with crosses in Fig. 4). Therefore, most of the heat capacity arises from the internal vibration within the material. Yet, these values still overestimate significantly the experimental heat capacity of amorphous cellulose, which is usually measured around 2 MJ.m$^{-3}$.K$^{-1}$ (Blokhin et al. 2011). We attribute this discrepancy to the high frequency of some internal vibrations (typically hydroxyl groups), for which a quantum correction would be needed to properly estimate the heat capacity (Honorio and Brochard 2022).

Since the poromechanical moduli are known as a function of relative humidity, one can integrate the constitutive equations (18) of poromechanics to estimate the evolution with RH. Doing so, one
can recover the moisture content and swelling isotherms of Fig. 2. To do so, we considered the fits of the poromechanical moduli displayed by the dashed lines in Fig. 3 (only *K*, *b*, and *N* are needed since the temperature is unchanged). The estimated water content and swelling are shown as lines in Fig. 2. Integration for the water content is started at 5% *RH*, because any tentative integration of the second state equation (18) of poromechanics turns out to be highly inaccurate at lower *RH*. Indeed, on the left-hand side, the inverse of the bulk density $1/n^b$ diverges at 0% *RH*, while on the right-hand side, both the Biot coefficient *b* and the inverse Biot modulus $1/N$ take extremely high values (Fig. 3). 5% *RH* is the smallest *RH* for which the integration provides consistent predictions with respect to the molecular simulation results.

The integration of the constitutive equations leads to reasonable water content and swelling isotherm, except above 80% *RH*. When approaching high *RH*, the predictions underestimate significantly both isotherms. We suspect to possible reasons for this:

1. The adsorption and induced swelling at *RH* > 80% do not evolve continuously but in a stepwise manner so that any tentative continuous integration fails and would require a much finer resolution in *RH* to be captured. This can be seen as a small length scale effect, since such stepwise evolution are expected to smoothen and vanish as bigger and bigger systems are considered.

2. The molecular simulations were not sufficiently long at these *RH* to fully capture the volume and water content fluctuations. This can be seen as a small time scale effect.

One can further investigate these hypothesis by inverse analysis of the molecular simulations isotherms, i.e., by determining what values of fluctuations would be needed to explain the observed water content and swelling isotherms. In particular, assuming a smooth swelling isotherm (i.e., excluding hypothesis 1), strain increments $\frac{d\epsilon}{dP_p}|_{\sigma,T} = b/K$ (Eq. 18) provide an evolution of the cross correlation $\Delta V \Delta N$ (4th relation in Eqs. 25) with *RH* that exceeds the $\Delta V \Delta N$ computed at *RH* above 70%. The difference is significant, in the sense that it well exceeds the intrinsic variability of our results of $\Delta V \Delta N$. This tends to suggest that more accurate $\Delta V \Delta N$ estimates (hypothesis 2)
is unable to explain the swelling isotherm. Accordingly, the most likely reason of the discrepancy
between the poromechanics and the molecular simulation seems to be first hypothesis (stepwise
isotherms because of small length scale).

CONCLUSION

In this paper, we show how one can estimate all moduli involved in the constitutive laws
of poromechanics from fluctuations in molecular simulations of a drained porous medium. We
derive fluctuation formulas for each of the poromechanical moduli, applicable in the osmotic
ensemble, which are easy to implement at negligible additional computational cost. Although,
molecular simulations in the osmotic ensemble are often intensive computationally, the fluctuation
formulas makes it possible to assess in a single simulation a wide variety of potential couplings.
The fluctuation formulas characterize generalized poromechanical moduli which are defined as
extension of the usual moduli of poromechanics to adsorbing media (usual poromechanics neglects
adsorption). These definitions coincide with that of the usual moduli for non-adsorbing media
but also quantify the unusual poromechanical couplings in presence of adsorption. Therefore,
the derived fluctuation formulas make it possible to characterize any poromechanical couplings,
including the unusual ones due to adsorption. In particular, we provide the fluctuation formula for
the heat of adsorption, which generalizes the usual fluctuation formula valid for rigid porous media
only.

In this work, we restricted ourselves to volumetric deformation for sake of simplicity. Indeed,
thermo-poro-mechanical couplings for isotropic media are often assumed limited to volumetric
effects. But deviatoric effects cannot be excluded in the general case, especially for anisotropic
media. Following the same route as above and taking into account the full strain and stress tensors,
one can derive similar fluctuation formulas to address the deviatoric effects. In particular, non-zero
cross fluctuations between the amount of fluid and the deviatoric deformations would point to the
existence of couplings between the fluid and the deviatoric behavior of the solid.

We also restricted ourselves to a pure fluid. For multi-component fluids, the energy balance
(15) of a drained porous medium takes the form:
$$d\omega = \sigma d\varepsilon - \sum_{i} n_{i} d\mu_{i} - s dT \qquad (36)$$

For a pure fluid, one can substitute the chemical potential control variable by the bulk fluid pressure thanks to the Gibbs-Duhem equation (energy balance 17). But, this is no more possible for a multi-component fluid. Thus, with a multi-component fluid, one cannot apply the generalized constitutive equations (18) and the generalized poromechanical moduli (19). It is necessary to introduce separate control parameters for each fluid component, such as activity, or partial pressure. Doing so, one would have to define component-specific poromechanical moduli (e.g., component-specific Biot coefficient). And corresponding fluctuation formulas can be derived with the same procedure as followed in this paper. The resulting fluctuation formulas would characterize poromechanical couplings that depend on the composition of the fluid, for instance the differential swelling upon a change in composition at constant bulk pressure (Brochard et al. 2012a).

To illustrate the application of the fluctuation formulas to a case of interest with respect to mecano-sorption couplings, we considered the case of amorphous cellulose submitted to moisture at various relative humidity. Molecular simulation were performed with a molecular model of amorphous cellulose available in the literature. All the thermo-poro-mechanical moduli were estimated from the fluctuations formulas, and we highlighted the unusual effect of adsorption on those moduli, such as a very Biot coefficient, or a negative thermal expansion. This example illustrates a first practical implementation of the fluctuation formulas.

## DATA AVAILABILITY STATEMENT

All data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

Funding by the Labex MMCD is gratefully acknowledged. Labex MMCD is joint research project of ENPC, Univ. Gustave Eiffel, UPEC, and CNRS, in the framework of the Programme
Investissement d'Avenir supervised by the French National Research Agency (ANR-11-LABX-022-01).

REFERENCES

Allen, M. P. and Tildesley, D. J. (1989). *Computer simulation of liquids*. Oxford University Press.

Berendsen, H. J. C., Grigera, J. R., and Straatsma, T. P. (1987). "The missing term in effective pair potentials." *The Journal of Physical Chemistry*, 91, 6269-6271.

Biot, M. A. (1941). "General theory of three-dimensional consolidation." *Journal of Applied Physics*, 12, 155-164.

Blokhin, A. V., Voitkevich, O. V., Kabo, G. J., Paulechka, Y. U., Shishonok, M. V., Kabo, A. G., and Simirsky, V. V. (2011). "Thermodynamic properties of plant biomass components. heat capacity, combustion energy, and gasification equilibria of cellulose." *Journal of Chemical & Engineering Data*, 56, 3523-3531.

Brochard, L. (2021). "Swelling of montmorillonite from molecular simulations: Hydration diagram and confined water properties." *The Journal of Physical Chemistry C*, 125, 15527-15543.

Brochard, L., Vandamme, M., Pellenq, R. J. M., and Fen-Chong, T. (2012a). "Adsorption-induced deformation of microporous materials: Coal swelling induced by co 2-ch 4 competitive adsorption." *Langmuir*, 28, 2659-2670.

Brochard, L., Vandamme, M., and Pellenq, R.-M. (2012b). "Poromechanics of microporous media." *Journal of the Mechanics and Physics of Solids*, 60, 606-622.

Chen, M., Coasne, B., Guyer, R., Derome, D., and Carmeliet, J. (2018). "Role of hydrogen bonding in hysteresis observed in sorption-induced swelling of soft nanoporous polymers." *Nature Communications*, 9, 3507.

Coussy, O. (2003). *Poromechanics*. Wiley, <https://onlinelibrary.wiley.com/doi/book/10.1002/0470092718> (12).

Coussy, O. (2010). *Mechanics and Physics of Porous Solids*. Wiley, <https://onlinelibrary.wiley.com/doi/book/10.1002/9780470710388> (7).
718

Einstein, A. (1905). "Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen." *Annalen der Physik*, 322, 549–560.

719

720

Gor, G. Y., Huber, P., and Bernstein, N. (2017). "Adsorption-induced deformation of nanoporous materials-a review." *Applied Physics Reviews*, 4, 011303.

721

722

Honorio, T. and Brochard, L. (2022). "Drained and undrained heat capacity of swelling clays." *Physical Chemistry Chemical Physics*, 24, 15003–15014.

723

724

Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J. (1996). "Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids." *Journal of the American Chemical Society*, 118, 11225–11236.

725

726

727

Karavias, F. and Myers, A. L. (1991). "Isosteric heats of multicomponent adsorption: thermodynamics and computer simulations." *Langmuir*, 7, 3118–3126.

728

729

Kulasinski, K., Guyer, R., Derome, D., and Carmeliet, J. (2015a). "Poroelastic model for adsorption-induced deformation of biopolymers obtained from molecular simulations." *Physical Review E*, 92, 022605.

730

731

732

Kulasinski, K., Guyer, R., Keten, S., Derome, D., and Carmeliet, J. (2015b). "Impact of moisture adsorption on structure and physical properties of amorphous biopolymers." *Macromolecules*, 48, 2793–2800.

733

734

735

Kulasinski, K., Keten, S., Churakov, S. V., Derome, D., and Carmeliet, J. (2014a). "A comparative molecular dynamics study of crystalline, paracrystalline and amorphous states of cellulose." *Cellulose*, 21, 1103–1116.

736

737

738

Kulasinski, K., Keten, S., Churakov, S. V., Guyer, R., Carmeliet, J., and Derome, D. (2014b). "Molecular mechanism of moisture-induced transition in amorphous cellulose." *ACS Macro Letters*, 3, 1037–1040.

739

740

741

Landau, L. D. and Lifshitz, E. M. (1980). *Statistical Physics*. Pergamon Press, third edition.

742

Lebowitz, J. L., Percus, J. K., and Verlet, L. (1967). "Ensemble dependence of fluctuations with application to machine computations." *Physical Review*, 153, 250–254.

743

744

Nezbeda, I. and Kolafa, J. (1991). "A new version of the insertion particle method for determining
the chemical potential by monte carlo simulation.” *Molecular Simulation*, 5, 391–403.

Nicholson, P., Nicholson, D., Nicholson, W., and Parsonage, N. G. (1982). *Computer Simulation and the Statistical Mechanics of Adsorption*. Academic Press.

M. Nič, J. Jirát, B. Košata, A. Jenkins, and A. McNaught, eds. (2009). *IUPAC Compendium of Chemical Terminology*. IUPAC, <http://goldbook.iupac.org> (6).

Thompson, A. P., Aktulga, H. M., Berger, R., Bolintineanu, D. S., Brown, W. M., Crozier, P. S., in 't Veld, P. J., Kohlmeyer, A., Moore, S. G., Nguyen, T. D., Shan, R., Stevens, M. J., Tranchida, J., Trott, C., and Plimpton, S. J. (2022). "Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales." *Computer Physics Communications*, 271, 108171.

Vandamme, M. (2019). "Coupling between adsorption and mechanics (and vice versa)." *Current Opinion in Chemical Engineering*, 24, 12–18.
<table><caption>TABLE 1. Bulk water properties at various relative humidities and 300K. Two values are provided for each property, the first one refers to the vapor phase (stable phase) and the second one refers to the liquid (metastable phase that exists only at RH above 30%). The vapor pressure and density are almost that of an ideal gas, so that we considered the properties of an ideal gas for the other properties (IG in the table): $K^b = P^b$, $h^b/n^b = 4kT = 9.977 \text{ kJ.mol}^{-1}$, $\alpha^b = 1/T = 3.333 \cdot 10^{-3} \text{ K}^{-1}$, and $c_p^b = 4k = 33.26 \text{ J.mol}^{-1}.\text{K}^{-1}$. Note that these expressions neglect the internal vibrations of the water molecule which are almost inactive at room temperature (Honorio and Brochard 2022).</caption><thead><tr><th>f (Pa)</th><th>RH (%)</th><th>P<sup>b</sup> (kPa)</th><th>n<sup>b</sup> (mol.m<sup>-3</sup>)</th><th>K<sup>b</sup> (kPa)</th><th>h<sup>b</sup>/n<sup>b</sup> (kJ.mol<sup>-1</sup>)</th><th>α<sup>b</sup> (K<sup>-1</sup>)</th><th>c<sub>p</sub><sup>b</sup> (J.mol<sup>-1</sup>.K<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1.0</td><td>0.1</td><td>0.0010 / -</td><td>4.029 · 10<sup>-4</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>10.0</td><td>1</td><td>0.0100 / -</td><td>4.029 · 10<sup>-3</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>20.0</td><td>2</td><td>0.0200 / -</td><td>8.057 · 10<sup>-3</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>49.9</td><td>5</td><td>0.0499 / -</td><td>2.014 · 10<sup>-2</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>99.9</td><td>10</td><td>0.1007 / -</td><td>4.029 · 10<sup>-2</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>149.8</td><td>15</td><td>0.1512 / -</td><td>6.049 · 10<sup>-2</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>199.7</td><td>20</td><td>0.2016 / -</td><td>8.070 · 10<sup>-2</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>249.7</td><td>25</td><td>0.2521 / -</td><td>1.009 · 10<sup>-1</sup> / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td><td>IG / -</td></tr><tr><td>299.6</td><td>30</td><td>0.3026 / -1.573 · 10<sup>5</sup></td><td>1.211 · 10<sup>-1</sup> / 4.954 · 10<sup>4</sup></td><td>IG / 8.789 · 10<sup>5</sup></td><td>IG / -41.28</td><td>IG / 1.342 · 10<sup>-3</sup></td><td>IG / 98.778</td></tr><tr><td>349.5</td><td>35</td><td>0.3530 / -1.388 · 10<sup>5</sup></td><td>1.414 · 10<sup>-1</sup> / 5.042 · 10<sup>4</sup></td><td>IG / 1.135 · 10<sup>6</sup></td><td>IG / -41.03</td><td>IG / 1.085 · 10<sup>-3</sup></td><td>IG / 96.868</td></tr><tr><td>399.5</td><td>40</td><td>0.4035 / -1.220 · 10<sup>5</sup></td><td>1.616 · 10<sup>-1</sup> / 5.116 · 10<sup>4</sup></td><td>IG / 1.337 · 10<sup>6</sup></td><td>IG / -40.80</td><td>IG / 8.938 · 10<sup>-4</sup></td><td>IG / 94.789</td></tr><tr><td>449.4</td><td>45</td><td>0.4540 / -1.066 · 10<sup>5</sup></td><td>1.819 · 10<sup>-1</sup> / 5.172 · 10<sup>4</sup></td><td>IG / 1.463 · 10<sup>6</sup></td><td>IG / -40.58</td><td>IG / 7.945 · 10<sup>-4</sup></td><td>IG / 93.408</td></tr><tr><td>499.3</td><td>50</td><td>0.5046 / -9.267 · 10<sup>4</sup></td><td>2.021 · 10<sup>-1</sup> / 5.218 · 10<sup>4</sup></td><td>IG / 1.554 · 10<sup>6</sup></td><td>IG / -40.38</td><td>IG / 7.366 · 10<sup>-4</sup></td><td>IG / 92.380</td></tr><tr><td>549.3</td><td>55</td><td>0.5551 / -7.997 · 10<sup>4</sup></td><td>2.224 · 10<sup>-1</sup> / 5.261 · 10<sup>4</sup></td><td>IG / 1.648 · 10<sup>6</sup></td><td>IG / -40.19</td><td>IG / 6.750 · 10<sup>-4</sup></td><td>IG / 91.336</td></tr><tr><td>599.2</td><td>60</td><td>0.6057 / -6.837 · 10<sup>4</sup></td><td>2.428 · 10<sup>-1</sup> / 5.295 · 10<sup>4</sup></td><td>IG / 1.710 · 10<sup>6</sup></td><td>IG / -40.02</td><td>IG / 6.517 · 10<sup>-4</sup></td><td>IG / 90.741</td></tr><tr><td>649.2</td><td>65</td><td>0.6563 / -5.776 · 10<sup>4</sup></td><td>2.631 · 10<sup>-1</sup> / 5.327 · 10<sup>4</sup></td><td>IG / 1.771 · 10<sup>6</sup></td><td>IG / -39.86</td><td>IG / 6.269 · 10<sup>-4</sup></td><td>IG / 90.069</td></tr><tr><td>699.1</td><td>70</td><td>0.7070 / -4.800 · 10<sup>4</sup></td><td>2.835 · 10<sup>-1</sup> / 5.357 · 10<sup>4</sup></td><td>IG / 1.834 · 10<sup>6</sup></td><td>IG / -39.71</td><td>IG / 5.988 · 10<sup>-4</sup></td><td>IG / 89.381</td></tr><tr><td>749.0</td><td>75</td><td>0.7576 / -3.895 · 10<sup>4</sup></td><td>3.040 · 10<sup>-1</sup> / 5.382 · 10<sup>4</sup></td><td>IG / 1.880 · 10<sup>6</sup></td><td>IG / -39.57</td><td>IG / 5.906 · 10<sup>-4</sup></td><td>IG / 89.096</td></tr><tr><td>799.0</td><td>80</td><td>0.8084 / -3.048 · 10<sup>4</sup></td><td>3.244 · 10<sup>-1</sup> / 5.408 · 10<sup>4</sup></td><td>IG / 1.940 · 10<sup>6</sup></td><td>IG / -39.44</td><td>IG / 5.633 · 10<sup>-4</sup></td><td>IG / 88.641</td></tr><tr><td>823.9</td><td>82.5</td><td>0.8337 / -2.643 · 10<sup>4</sup></td><td>3.347 · 10<sup>-1</sup> / 5.419 · 10<sup>4</sup></td><td>IG / 1.963 · 10<sup>6</sup></td><td>IG / -39.38</td><td>IG / 5.579 · 10<sup>-4</sup></td><td>IG / 88.419</td></tr><tr><td>848.9</td><td>85</td><td>0.8591 / -2.246 · 10<sup>4</sup></td><td>3.449 · 10<sup>-1</sup> / 5.431 · 10<sup>4</sup></td><td>IG / 1.990 · 10<sup>6</sup></td><td>IG / -39.32</td><td>IG / 5.466 · 10<sup>-4</sup></td><td>IG / 88.204</td></tr><tr><td>873.9</td><td>87.5</td><td>0.8845 / -1.858 · 10<sup>4</sup></td><td>3.552 · 10<sup>-1</sup> / 5.442 · 10<sup>4</sup></td><td>IG / 2.016 · 10<sup>6</sup></td><td>IG / -39.26</td><td>IG / 5.360 · 10<sup>-4</sup></td><td>IG / 87.998</td></tr><tr><td>898.8</td><td>90</td><td>0.9100 / -1.476 · 10<sup>4</sup></td><td>3.654 · 10<sup>-1</sup> / 5.452 · 10<sup>4</sup></td><td>IG / 2.032 · 10<sup>6</sup></td><td>IG / -39.20</td><td>IG / 5.380 · 10<sup>-4</sup></td><td>IG / 87.898</td></tr><tr><td>923.8</td><td>92.5</td><td>0.9354 / -1.099 · 10<sup>4</sup></td><td>3.757 · 10<sup>-1</sup> / 5.463 · 10<sup>4</sup></td><td>IG / 2.064 · 10<sup>6</sup></td><td>IG / -39.14</td><td>IG / 5.197 · 10<sup>-4</sup></td><td>IG / 87.636</td></tr><tr><td>948.8</td><td>95</td><td>0.9608 / -7.243 · 10<sup>3</sup></td><td>3.859 · 10<sup>-1</sup> / 5.473 · 10<sup>4</sup></td><td>IG / 2.084 · 10<sup>6</sup></td><td>IG / -39.09</td><td>IG / 5.160 · 10<sup>-4</sup></td><td>IG / 87.491</td></tr><tr><td>973.7</td><td>97.5</td><td>0.9863 / -3.511 · 10<sup>3</sup></td><td>3.962 · 10<sup>-1</sup> / 5.481 · 10<sup>4</sup></td><td>IG / 2.100 · 10<sup>6</sup></td><td>IG / -39.03</td><td>IG / 5.166 · 10<sup>-4</sup></td><td>IG / 87.382</td></tr><tr><td>998.7</td><td>100</td><td>1.0117 / 1.012</td><td>4.064 · 10<sup>-1</sup> / 5.490 · 10<sup>4</sup></td><td>IG / 2.113 · 10<sup>6</sup></td><td>IG / -38.98</td><td>IG / 5.194 · 10<sup>-4</sup></td><td>IG / 87.284</td></tr></tbody></table>
List of Figures

767

1 Snapshots illustrating the evolution of the amorphous cellulose model when relative humidity is increased from 0% to 100%. The balls represent the atoms of the water molecules, the sticks represent the covalent bonds within the amorphous cellulose structure, and the colors depend on the nature of the atoms. The solid line boxes represent the system size at the given *RH* and the dashed line boxes represent the system size at 0% *RH*.

768

769

770

771

772

773

37

2 Water content and swelling isotherms of amorphous cellulose at 300K under zero confining stress. The circles are the results of the combined CB-GCMC / MD simulations, and the lines are estimates obtained by integration of the poromechanical constitutive equations. The diamonds are results of Chen et al. (2018) (note that for the water content, the comparison is for the same cellulose model, while for the swelling it is not, see text).

774

775

776

777

778

779

38

3 Poromechanical moduli of the amorphous cellulose submitted to moisture at various relative humidities. All moduli are obtained by applying the fluctuation formulas (26-31), with water vapor as the bulk reference (bulk properties in Table 1). The dashed lines are fitted trends and the associated *R*² determination coefficients are provided on the charts. For *b* and *α*<sub>φ</sub> that vary on logarithmic scales, the *R*² values are calculated by considering the logarithm of the properties.

780

781

782

783

784

785

39

4 Poromechanical moduli of the amorphous cellulose when the reference bulk fluid is the meta-stable liquid water. All moduli are obtained by applying the fluctuation formulas (26-31), with liquid water as the bulk reference (bulk properties in Table 1). Since the liquid water phase does not exist at *RH* below 30% (cavitation), the moduli are estimated for *RH* greater or equal to 30%.

786

787

788

789

790

40
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Snapshots illustrating the evolution of the amorphous cellulose model when relative humidity is increased from 0% to 100%. The balls represent the atoms of the water molecules, the sticks represent the covalent bonds within the amorphous cellulose structure, and the colors depend on the nature of the atoms. The solid line boxes represent the system size at the given *RH* and the dashed line boxes represent the system size at 0% *RH*.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Water content and swelling isotherms of amorphous cellulose at 300K under zero confining stress. The circles are the results of the combined CB-GCMC / MD simulations, and the lines are estimates obtained by integration of the poromechanical constitutive equations. The diamonds are results of Chen et al. (2018) (note that for the water content, the comparison is for the same cellulose model, while for the swelling it is not, see text).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Poromechanical moduli of the amorphous cellulose submitted to moisture at various relative humidities. All moduli are obtained by applying the fluctuation formulas (26-31), with water vapor as the bulk reference (bulk properties in Table 1). The dashed lines are fitted trends and the associated *R*² determination coefficients are provided on the charts. For *b* and $\alpha_{\phi}$ that vary on logarithmic scales, the *R*² values are calculated by considering the logarithm of the properties.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Poromechanical moduli of the amorphous cellulose when the reference bulk fluid is the meta-stable liquid water. All moduli are obtained by applying the fluctuation formulas (26-31), with liquid water as the bulk reference (bulk properties in Table 1). Since the liquid water phase does not exist at *RH* below 30% (cavitation), the moduli are estimated for *RH* greater or equal to 30%.</figcaption></figure>