# Cellulose consolidated with polyethylene glycol: The nanoscale mechanisms revealed by hybrid Monte Carlo/molecular dynamics modeling

Ali Shomali, Chi Zhang, Benoit Coasne, Eleanor J Schofield, Brigitte Chabbert, Dominique Derome, Jan Carmeliet

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

## To cite this version:

Ali Shomali, Chi Zhang, Benoit Coasne, Eleanor J Schofield, Brigitte Chabbert, et al., Cellulose consolidated with polyethylene glycol: The nanoscale mechanisms revealed by hybrid Monte Carlo/molecular dynamics modeling. International Journal of Biological Macromolecules, 2025, 285, pp.137661. 10.1016/j.ijbiomac.2024.137661. hal-05096478

HAL Id: hal-05096478

https://hal.inrae.fr/hal-05096478v1

Submitted on 4 Jun 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_2.png" /></figure>

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

# Cellulose consolidated with polyethylene glycol: The nanoscale mechanisms revealed by hybrid Monte Carlo/molecular dynamics modeling

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

Check for
updates

Ali Shomali<sup>a</sup>, Chi Zhang<sup>a,b,*</sup>, Benoit Coasne<sup>c</sup>, Eleanor J. Schofield<sup>d</sup>, Brigitte Chabbert<sup>e</sup>, Dominique Derome<sup>f</sup>, Jan Carmeliet<sup>a</sup>

## ABSTRACT

Polyethylene glycol (PEG) consolidation treatment is a widely used conservation strategy for wooden culture relics. However, the consolidation mechanism of PEG is still open to interpretation. PEG-cellulose, the representative component of wood cell wall, interactions are governed by various coupled multi-scale mechanisms which require nano-scale investigation. In this study, a hybrid molecular dynamics and grand canonical Monte Carlo (MD/GCMC) simulation combined with rule of mixture (RoM) analyses are employed to reveal the underlying mechanisms of PEG-induced consolidation. We found that PEG200 reduces moisture adsorption and swelling at museological conditions, confirming its consolidation effect. At high PEG content, a crossover behavior is identified at humid conditions (RH > 80) where excessive sorption and swelling are observed surpassing the untreated sample. The molecular modeling results are found to be consistent with experimental observations. Furthermore, the structural and mechanical properties of the hydrated samples are assessed by examining the porosity distribution, mechanical properties, and hydrogen bonding network. Results indicate mechanical softening induced by PEG treatment. A modified mixture model is proposed based on molecular modeling results that incorporate sorption and swelling coupling, porosity filling and mechanical softening behaviors. Two key mechanisms are identified explaining the consolidation effect of PEG: first, the PEG fills the porosities of amorphous structure thus diminishing sorption sites; second, the polymer structure prohibits PEG from further swelling thus constraining water sorption. The model and theoretical framework can serve as a guide for the design of novel consolidant materials by identifying the key molecular features of an ideal consolidant.

## 1. Introduction

During centuries of burial in undersea sediments, waterlogged archaeological wood of shipwrecks becomes severely degraded as a result of chemical deterioration and material loss [53,78]. The fragile residual structure of archaeological wood is characterized by its reduced strength and density [43,66], increased porosity and permeability [4], and increased hygroscopicity [7,66]. Without proper consolidation, there is a high risk of irreversible damage due to drastic shrinkage, distortion and collapse of these materials during drying [4]. With the discovery of more marine artifacts such as the recently discovered wreck of the Applet, Vasa's sister warship, advancing our current understanding of treatment processes and consolidation agents is deemed essential.

Among the variety of consolidation practices in the arsenal of archaeologists (see Christensen et al. [18]), the treatment of wood cell walls with polyethylene glycol (PEG) is the most common technique to consolidate waterlogged archaeological wooden artifacts. PEG treatment has been successfully applied in the preservation of the Swedish warship Vasa [30] as well as Henry VIII's warship the Mary Rose [67] for long-term study and display. PEG molecules are introduced into wood structures either by full immersion of smaller objects through a one-step impregnation procedure or by spraying a PEG solution on the surface of the larger objects through a one or two-step procedure [33]. During the decades-long consolidation procedure, PEG molecules gradually diffuse into the wood cell wall and replace the water molecules stabilizing the structure, and making the wood less susceptible to changes in humidity.

<footnote><sup>a</sup> Chair of Building Physics, Department of Mechanical and Process Engineering, ETH Zurich, 8093 Zurich, Switzerland</footnote>

<footnote><sup>b</sup> School of Flexible Electronics, Sun Yat-sen University, 518000 Shenzhen, China</footnote>

<footnote><sup>c</sup> Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France</footnote>

<footnote><sup>d</sup> The Mary Rose Trust, HM Naval Base, Portsmouth PO1 3LX, UK</footnote>

<footnote><sup>e</sup> Université de Reims Champagne Ardenne, INRAE, FARE, 51097 Reims, France</footnote>

<footnote><sup>f</sup> Department of Civil and Building Engineering, Sherbrooke University, Sherbrooke, Canada</footnote>

<footnote>* Corresponding author at: Chair of Building Physics, Department of Mechanical and Process Engineering, ETH Zurich, 8093 Zurich, Switzerland.</footnote>

<footnote>E-mail addresses: ashomali@ethz.ch (A. Shomali), outlook.zhangchi@gmail.com (C. Zhang), benoit.coasne@univ-grenoble-alpes.fr (B. Coasne), E.Schofield@maryrose.org (E.J. Schofield), brigitte.chabbert@inrae.fr (B. Chabbert), Dominique.Derome@usherbrooke.ca (D. Derome), cajan@ethz.ch (J. Carmeliet).</footnote>
PEG-cell wall interactions are governed by various entangled multi-scale mechanisms that are yet to be elucidated and distinguished. Low molecular weight PEG (200–800) is used to penetrate the cell wall to suppress drying shrinkage while high molecular weight PEG (2000–4000) is used to fill the lumens and bigger cracks to prevent collapse and seal the surfaces to prevent bleeding [8,10,31,81]. Tracing PEG and monitoring its distribution within the wood structure is still the subject of active research [55]. Added benefits of increased structural integrity and controlled shrinkage of PEG-treated samples come at the expense of increased weight, change of colour, and altered sorption behavior. PEG treatment can also lead to cell wall softening and, in some cases, acceleration of creep behavior [76]. For instance, recent oak impregnated with PEG 600 exhibits a strong decrease in mechanical properties particularly in the radial direction [43]. The lower stiffness of PEG-treated wood can be attributed to both interactions of PEG with the cell wall and increased moisture content resulting from PEG-induced hygroscopicity [77]. Moreover, recent studies have identified new challenges concerning the slow depolymerization of PEG due to Fenton-type reactions in the chemically active environment of degraded wood [3].

There are several experimental studies of the impact of PEG treatment on hygroscopic properties of various model materials such as recent oak and archaeological oak from Vasa and Mary Rose warships (Quercus robur L.) [7,62,75,77], archaeological oak from a 10th-century stronghold in Poland (Quercus spp.) [58], waterlogged archaeological and subfossil Scots pine (Pinus sylvestris L.) [22,23,45], Southern pine (Pinus spp.), sapwood [32], black willow (Salix sp.), aspen poplar (Populus sp.), white cedar (Thuja occidentalis) and elm (Ulmus sp.) [25,26]. All studies confirm a significant change in the hygroscopic response of treated wood compared to untreated or archaeological samples. Sorption isotherms of PEG-treated archaeological oak show reduced moisture content below a relative humidity (RH) of approx. 80 %, followed by excessive adsorption in humid conditions. Similar behavior has been observed in PEG-treated samples of archaeological Scots pine [47], and elm [26]. The increased hygroscopicity in humid conditions is mostly attributed to PEG present in the lumens and larger pores of degraded wood whereas the PEG in the cell wall is often assumed to be less hygroscopic compared to pure PEG due to confinement by the relatively stiffer surrounding structure of cell wall polymers [9,43]. Several mechanisms have been hypothesized to explain the hygromechanical properties of treated wood cell walls. Grattan [26] assumed that the strong interaction of PEG and wood cell wall polymers prevents the full adsorption of water into the cell wall. By implementing the Guggenheim-Anderson-de Boer (GAB) model, Olek et al.[58] and Majka et al.[45] suggested the blocking of primary sorption sites by PEG molecules to explain the lower moisture content below 80 % RH. Explaining the high relative humidity case, they suggested intensive polymolecular water adsorption due to PEG's hygroscopic nature. The suggested explanations are ambiguous in describing the underlying molecular physics and rely mainly on macroscale experiments which do not provide a complete and detailed picture of the nanoscale phenomenon.

MD simulations serve as an alternative to experiments to assess PEG-cellulose interactions in detail. Molecular modeling has been extensively used to capture the structural and hygromechanical properties of wood cell wall in particular and biopolymers in general. Some examples include the molecular modeling of cellulose allomorphs such as amorphous cellulose [41,52], paracrystalline cellulose [36,41], crystalline cellulose [50,79], and the formation and mechanical properties of cellulose microfibril aggregates [16,57]. Cell wall matrix components such as various lignins [73,82,86] and hemicelluloses [29,85] and their interaction with cellulose fibrils have been extensively studied [37,51]. By assembling cellulose fibrils, hemicellulose, and lignin different molecular models are then developed for cell wall S2 layer composite [11,20,39,84]. Derome et al. [21] used molecular dynamics to study water-cell wall interactions. They examined the microscopic origin of

the coupling between sorption and sorption-induced swelling and softening of the cell wall. During adsorption, water molecules initially adsorb sparsely in the available porosity. With further increasing the moisture content above 0.1, the material undergoes a transition by the formation of a water network through percolation [40]. Through this process, water molecules are found to rearrange polymer chains and form hydrogen bonds with available sorption sites. Water sorption induced volumetric swelling and breakage of the existing hydrogen bonds leading to mechanical softening. This molecular mechanism of moisture-induced softening through the disruption of the polymer-polymer hydrogen bonding network was observed for cell wall matrix, matrix-cellulose microfibril interphase [83], and secondary cell wall (S2) composite [84]. By examining the distribution of water molecules in the matrix, Zhang et al. [85] identified the formation of single and double layers of water around the polymer chains associated with the changes in mechanical properties. The mixture rule was shown to successfully predict the mechanical and thermodynamic properties of the composite obtained from molecular modeling [85]. Further, Zhang et al. [84] used a similar mixture rule analysis to reveal the interphase properties of the S2 composite. A similar treatment is applied in the current study to analyze the sorption and swelling isotherms.

This research aims to study the impact of PEG consolidation on the hygromechanics of cellulose, the representative component of wood cell wall, and to provide a nanoscale insight into its underlying physics. The purpose is to bridge the knowledge gap between the molecular-level mechanisms governing the interactions between PEG and cellulose and the macroscale behavior with a particular emphasis on the conservation and the preservation of artifacts such as the Swedish warship Vasa and Henry VIII's warship the Mary Rose. In the treatment of waterlogged archaeological wood, polyethylene glycol (PEG) consolidant interacts with cellulose at the molecular level, leading to altered hygromechanical properties, the mechanisms of which will be explained in the manuscript. The main objectives of this paper are to develop a robust computational modeling framework that successfully models moisture adsorption and sorption-induced deformation in mixtures of AC and PEG, and to integrate the nanoscale information from molecular modeling within a mixture model. To this aim, we adopt a hybrid MD-GCMC method where molecular dynamics is combined with grand canonical Monte Carlo simulation to study sorption and sorption-induced deformation on a molecular scale [15,19,42,65]. The hybrid MD-GCMC method has been revealed to be a powerful tool in capturing different aspects of the adsorption and deformation phenomena in wood bio-polymers. Chen et al. [12,13] were able to investigate the origin of hysteresis in deformation and hydrogen bonding network interplay upon adsorption and desorption in amorphous cellulose and cellulose nanocomposites. Following a similar framework, we investigate the sorption isotherms, swelling curves, and mechanical properties of mixtures of amorphous cellulose and PEG200 as model material representing the interaction of PEG200 with cellulose. The study is conducted on mixtures with varying ratios of PEG200 and sorption range from fully dry to RH = 0.9, approaching saturation point. Furthermore, the structural and mechanical properties of the equilibrated samples are analyzed by examining the porosity distribution and performing mechanical tests on hydrated samples. The simulation methodology and the choice of simulation parameters such as system size are validated through comparison with available simulation and experimental data [54,56]. A model is then proposed to describe the impact of PEG treatment on porosity evolution. The data attained from molecular modeling is then complemented by comparison with the rule of mixture analysis. Finally, a new modified mixture model is proposed to incorporate sorption/swelling, porosity filling, and mechanical softening unified in a single model for AC + PEG mixtures revealing the microscopic mechanism of PEG treatment. By characterizing the molecular mechanisms governing consolidation at the nanoscale, our research may guide the development of novel consolidant materials for treating and preserving cellulose-based materials, especially archaeological wood. The
nanoscale information obtained from our MD simulations can lay a crucial foundation for developing more complex atomistic models and be potentially integrated into finite element frameworks to be valuable for both material science and wood conservation fields.

## 2. Methodology

In this study, six material systems are used to explore the adsorption of water in amorphous cellulose (AC - (C₆H₁₀O₅)₂₀), PEG200 (from this point referred to as PEG - H-(O-CH₂-CH₂)₄-OH) and their mixtures. Studied systems boxes include a pure AC system containing 7 cellulose molecules each made of 10 cellobiose units, a pure PEG box containing 60 PEG molecules and four mixtures simulation boxes containing the same 7 cellulose molecules and 19, 37, 74 or 123 PEG molecules, achieving AC + 15%PEG, AC + 30%PEG, AC + 60%PEG and AC + 100%PEG mixtures respectively. The PEG percentage in mixture names is determined by the mass of added PEG per initial untreated AC mass (mPEG/mAC). The dry system size ranges from 1890 atoms for the pure PEG system to 6800 atoms for the AC + 100%PEG system. It has been shown previously that an AC system of such size and degree of polymerization is sufficient for obtaining adequate results for the determination of the hygromechanical properties of AC [14,41]. The computational costs scale with the system size by power law which is also an important factor to consider.

Molecular topologies for cellulose and PEG are constructed using Cellulose-builder [24] and Materials Studio 8.0 package respectively and molecules are inserted randomly in the box using PACKMOL [49]. The initial simulation box is then imported into GROMACS using the Intermol [69] package. It is known that non-bonded interactions, especially coulombic interactions, play a key role in simulations predicting the adsorption of polar molecules like water and carbon dioxide [19,80]. In atomistic simulations, the assignment of charges is generally done by choosing the force field. Force fields are a collection of interaction potential terms which describe the intra- and inter-molecular interactions of atoms in the simulation box. In our study, an all-atom OPLS force field, i.e., OPLS-AA, is employed to describe the cellulose and PEG molecules. OPLS is a widely used force field for the simulation of biopolymers, which is successful in capturing molecular properties of different cellulose allomorphs, such as AC [2,59], paracrystalline cellulose [63], crystalline cellulose [35,60], and their interactions with water [61]. Water molecules are modeled using the extended single point charge model (SPC/E) which has been successfully used to obtain the expected vapor pressure of water using Monte Carlo methods [6,74].

Once built, the systems are relaxed as follows. Starting with energy optimization using the steepest descent and conjugate gradient, the initial structures are then subjected to an NVT simulation at 600 K temperature for 5 ns. Then, by performing an NPT simulation at 0 bar, the system temperature is linearly decreased to 300 K over 10 ns of simulation. The NPT simulation is continued for another 5 ns at 1 bar and 300 K to ensure equilibration. For all simulations, the equations of motion are integrated using a Varlet algorithm with a 1 fs time step. The cutoff radius of coulombic and van der Waals interactions are set to 12 Å and the long-range coulombic interactions are accounted for by particle mesh Ewald summation with a precision of 1e-5 (kCal mol⁻¹ Å⁻¹). Thermostating and barostating are done by implementing the Nosé-Hoover thermostat and barostat. Periodic boundary conditions (PBC) are applied to all directions to mimic an infinite material and hence avoid finite-size effects. The resulting densities of dry AC and PEG after the equilibration are 1.41 and 1.10 (g cm⁻³) respectively, in very good agreement with both experimental and simulation data. Table 1 shows some of the physical properties of AC and PEG obtained by MD compared to values obtained in previous computational and experimental studies, showing very good agreement. In the present work, we use a hybrid MD/GCMC method to predict the moisture adsorption isotherms of AC and PEG mixtures. This method, which effectively implements the osmotic ensemble (µPT), allows to encapsulate the most

<table><caption>**Table 1**

Physical properties of AC and PEG200 of this study compared to previous studies.</caption><thead><tr><th>Component</th><th>Properties</th><th>Current simulation</th><th>Reported</th><th>Ref</th><th>Description</th></tr></thead><tbody><tr><td rowspan="4">AC</td><td>Density (g/cm<sup>3</sup>)</td><td>1.41</td><td>1.48<br>1.409<br>1.285<br>1.385 ± 0.032</td><td>[34]<br>[14]<br>[41]<br>[17]</td><td>Exp.<br>MD<br>MD<br>MD</td></tr><tr><td>Young's modulus (GPa)</td><td>8.0</td><td>8.5<br>5.6 ± 2.8</td><td>[48]<br>[41]</td><td>Exp.<br>MD</td></tr><tr><td>Bulk modulus (GPa)</td><td>7.1</td><td>6.5</td><td>[39]</td><td>MD</td></tr><tr><td>Poisson ratio</td><td>0.3</td><td>0.23 ± 0.03</td><td>[44]</td><td>MD</td></tr><tr><td>PEG</td><td>Density (g/cm<sup>3</sup>)</td><td>1.10</td><td>1.12</td><td>[87]</td><td>Exp.</td></tr></tbody></table>

prominent aspects of a mechanosensitive phenomenon, namely sorption and swelling. Adsorption or desorption of moisture can be modeled with an imposed chemical potential during the GCMC simulation, while deformation, either swelling or shrinkage, can be modeled under constant pressure in MD. In our study, GCMC is implemented using the Cassandra [68] package and MD is performed using GROMACS [1]. Following the equilibration process, the dry samples are subjected to GCMC simulation based on the µVT ensemble, which performs 10,000 Monte Carlo moves consisting of insertion and deletion of water molecules with equal probabilities while the polymer molecules are treated as rigid solid. These random configurations are then accepted/rejected based on the Boltzmann probability. MD NPT is then carried out for 10 ps to allow sorption-induced swelling and subsequent increase in available porosity that results from the occupancy of water molecules into the pore structure. Fig. 1 shows a schematic representation of the hybrid MD/GCMC simulation procedure. Depending on the chemical properties of adsorbate and adsorbent, this cycle should be iteratively repeated to achieve thermodynamic equilibrium. Due to the dense

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Schematic representation of the hybrid MD/GCMC simulation procedure.</figcaption></figure>
structure and small porosity <2.5 Å of the matrix, the success ratio of insertion and deletion moves is very low (<0.1 %) and, depending on relative humidity, the cycle needs to be repeated from 5000 cycles for low relative humidity up to 80,000 cycles for high relative humidity. Examples of the evolution of the number of adsorbed water molecules with MD/GCMC cycles can be found in Fig. S1 in the supplementary information. In short, we used a combined MD and GCMC approach to simulate the interaction of water with AC and PEG mixtures. MD allows us to model how the material swells when water is absorbed, while GCMC simulates the actual adsorption or removal of water molecules at different humidity levels. By alternating between these two methods, we can accurately capture both the absorption and swelling behaviors, repeating the process until the system reaches a stable equilibrium. Overall, the data in Fig. S1 shows that sampling in our hybrid MD/GCMC simulations is efficient and that moisture content reaches its final equilibrium value provided simulations are long enough. Sorption isotherms are reported up to vapor pressure values close to saturation pressure ($p/p_0 \le 0.9$). With further increase of the vapor pressure, uncontrolled swelling and solvation of the solid frame are observed. Fig. 2 shows a snapshot of AC + 30%PEG mixture with 0.1 moisture content. The moisture adsorption and swelling results correspond to the average value over three replicas of simulation systems.

Simulations provide the number of adsorbed water molecules as a function of chemical potential. Moisture content (MC) at equilibrium is then determined as:

$$MC = \frac{mass_{water}}{mass_{dry}} \qquad (1)$$

where mass$_{water}$ is the mass of adsorbed water and mass$_{dry}$ is the mass of the dry system. The volumetric swelling strain $\epsilon$ is calculated in reference to the initial dry state:

$$\epsilon_{(m)} = \frac{V_{(m)} - V_{(m=0)}}{V_{(m=0)}} \qquad (2)$$

where V is the volume of the box at a certain moisture content.

The ratio between the volume of PEG, $V_{PEG}$, and the volume of AC, $V_{AC}$, is defined as:

$$\alpha = V_{PEG}/V_{AC} \qquad (3)$$

<figure><img src="image_6.png" /><figcaption>Fig. 2. Snapshot of AC + 30%PEG mixture with 5 % moisture content, 3 nm × 3 nm approx. Line shows the cellulose, orange spheres show PEG molecules and red and white spheres are water molecules.</figcaption></figure>

Experimental data often report moisture content as a function of vapor pressure or relative humidity. Therefore, assuming water vapor to be an ideal gas, the chemical potential is related to vapor pressure using the equation of state:

$$\mu - \mu_0 = k_B T \ln RH \qquad (4)$$

where $\mu$ denotes the chemical potential, $\mu_0$ is the chemical potential at saturation, T is temperature and RH is relative humidity:

$$RH = \frac{p}{p_0} \qquad (5)$$

p is the vapor pressure of water and $p_0$ is the saturation vapor pressure, which is 1017 Pa for the SPC/E water model at 300 K.

Once the system is equilibrated, different hydrated samples at different relative humidity can be further examined through hydrogen bonding, mechanical tests, and porosity measurement, and with the aid of mixture rule analysis. In order to obtain the porosity $\phi$ and pore size distribution (PSD), the Pore Blazer v4.0 code [64] is used with the probe diameter set to 2.75 Å, thus approximately equal to the water molecule diameter. The available pore space is inspected by randomly inserting the probe particle into the matrix. The term “available” porosity is used to describe the porosity that can be accessed by the probe. In contrast, “unavailable” porosity corresponds to the porosity that cannot be accessed by the probe but may become available through chain deformation and dislocation. One can refer to [64] for a detailed description of the method. Although degraded wood displays a complex pore network with a significantly varied pore structure ranging from nanopores to mesopores, reaching even macropores [27], in a molecular simulation, we mainly examine the nanoscale polymer structure and capture the distribution of nanopores at that scale. Hydrogen bonds are recognized using a pair of geometric criteria. Considering a donor and an acceptor, a distance smaller than 0.35 nm and a hydrogen-donor-acceptor angle smaller than 30°, the pair is assumed to share a hydrogen bond. Furthermore, the mechanical moduli, namely bulk and Young's moduli, are determined from the slope of the linear regime of the stress-strain curve calculated from volumetric and uniaxial straining tests, respectively.

The rule of mixture analysis has proved to be an insightful practice in revealing the interphase properties and synergic effects in a mixture of polymers [84]. Here, employing the mixture rule, the volumetric density $\rho_{mix}$ of the mixture can be calculated by:

$$\rho_{mix} = \sum f_{vi} \rho_i \qquad (6)$$

where $f_{vi}$ is the volume fraction of component i.

Moisture content $MC_{mix}$ and swelling strain $\epsilon_{mix}$ of the mixture can then be calculated based on the moisture content and swelling strain of components.

$$MC_{mix} = \sum w_i . MC_i \qquad (7)$$

$$\epsilon_{mix} = \sum f_{vi} \cdot \epsilon_i \qquad (8)$$

where $w_i$ and $\epsilon_i$ are mass fraction and swelling strain of component i.

## 3. Results

### 3.1. Sorption and swelling isotherms, swelling curve of AC, PEG and AC-PEG mixtures

The sorption isotherms, representing moisture content (MC) versus relative humidity (RH), the swelling isotherms, represent volumetric swelling strain ($\epsilon$) versus RH, and the swelling curve, representing $\epsilon$ versus MC, for amorphous cellulose (AC) and PEG are shown in Fig. 4. The sorption isotherms are compared with experimental data and a
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** MD/GCMC simulation results (colored dots) and parametric functions inspired by GAB model in Eq. (9) (line) for (a) sorption isotherm, (b) swelling isotherm, and (c) swelling curve for AC; (d) sorption isotherm, (e) swelling isotherm and (f)swelling curve for PEG and experimental results [54,56] (black dots) for the sorption isotherms (a-d). Error bars added in (c-f).</figcaption></figure>

reasonable agreement is observed. It is important to note that the experimental data is obtained for low crystallinity cellulose [54], and not pure AC. Moreover, it is difficult to achieve a fully dehydrated AC sample in practice, which can explain why experimental results may differ from MD/GCMC results at low relative humidity due to the presence of residual non-desorbable water [14]. The small system size in atomistic simulations (< 100 nm³) restricts the maximal pore size and may lead to smaller moisture contents at high RH compared to experimental data.

The sorption and swelling isotherms are approximated by a parametric function inspired by the GAB model [72].

$$f_s(\text{RH}) = c_s \frac{a_s b_s \text{RH}}{(1 - b_s \text{RH})(1 - b_s \text{RH} + a_s b_s \text{RH})} \quad (9)$$

where $a_s$, $b_s$ and $c_s$ are parameters to be determined. This function provides suitable flexibility for describing different types of sorption isotherms and is also used to describe experimental sorption isotherms of PEG-treated wood [58]. The values of the parameters are given in Table S1.

The sorption isotherm of AC in Fig. 4a reflects a type II isotherm according to the classification of IUPAC (International Union of Pure and Applied Chemistry) [71] which is in agreement with previous molecular modeling studies on AC [13]. At low relative humidity, moisture content increases with RH as water molecules fill the porosity present despite the chains of AC being tightly packed, referred to as the initial porosity-filling stage [84]. Adsorption proceeds progressively and, from RH > 0.7, moisture content increases somewhat more importantly due to the swelling of the material. The swelling curve (ε versus MC) in Fig. 4c is described by a bilinear function, in which two slopes $β_A$ and $β_B$, referred to as swelling coefficients are reported in Table 2. Swelling is lower in the MC range from 0 to 0.05 with a swelling coefficient $β_A < 1$, where the initial porosity of AC is being filled. With further increasing moisture content, the swelling strain increases, and the polymeric matrix opens up accommodating additional water molecules. In this range,

<table><caption>**Table 2**
Swelling coefficients of AC, PEG, and AC-PEG mixtures.</caption><thead><tr><th>Material</th><th>β<sub>A</sub></th><th>β<sub>B</sub></th></tr></thead><tbody><tr><td>AC</td><td>2.89e-04</td><td>6.32e-04</td></tr><tr><td>AC + 15%PEG</td><td>2.77e-04</td><td>6.58e-04</td></tr><tr><td>AC + 30%PEG</td><td>2.29e-04</td><td>7.06e-04</td></tr><tr><td>AC + 60%PEG</td><td>2.35e-04</td><td>8.12e-04</td></tr><tr><td>AC + 100%PEG</td><td>4.87e-05</td><td>8.52e-04</td></tr><tr><td>PEG</td><td>-</td><td>9.39e-04</td></tr></tbody></table>

the swelling coefficient $β_B$ is almost equal to 1 and the volume change due to swelling equals the volume of the added water molecules. The separation of polymer chains due to swelling leads also to a mechanical weakening, which will be discussed later (see Fig. 6), showing a decrease of the elastic modulus with increasing MC.

PEG is a highly water-soluble liquid polymer [56] and shows a sorption isotherm (see Fig. 1d) significantly different from AC. Initially, at low RH, moisture content and swelling strain are small followed by a significant increase in MC and swelling at higher RH, resulting in a solvation of PEG. Similar behavior is also reported for high molecular weight (MW = 1450–8000) PEGs where at a critical relative humidity the semi-crystalline polymer undergoes a solid-to-solution phase transformation which is also regarded as a deliquescence process [5]. Experimental results suggest elevated sorption and lower deliquescence relative humidity for low molecular weight PEGs [28]. The swelling curve shows a single linear behavior with a swelling coefficient $β_B$ close to 1.

The sorption and swelling isotherms, and swelling curves for AC, PEG and all AC-PEG mixtures are compared in Fig. 5a-c, respectively. We remark that the moisture content of PEG is lower at low RH (<55 %) than the MC of AC, but that at high RH, the MC of PEG surpasses considerably MC of AC.

Compared to AC, AC treated by PEG shows a reduction in moisture content and swelling strain at low RH. This reduction is most evident in
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Simulation results (colored dots) and parametric functions (lines) for (a) sorption isotherms, (b) swelling isotherms, (c) swelling curves with error bars on dots, for AC, PEG, and AC-PEG mixtures. The mean deviation for sorption and swelling results over 3 replicates of simulation system is <3 %.</figcaption></figure>

mixtures with lower PEG content, such as AC + 15%PEG and AC + 30% PEG, showing compared to AC, an average reduction of swelling of 25% at RH = 0.45-0.55 (museum conditions). However, further increasing the PEG content adversely affects the sorption and swelling of the mixture at high RH. Hence, AC + 60% PEG and AC + 100% PEG mixtures show a low moisture content at low RH, but MC increases significantly at high RH and their sorption and swelling isotherms crossover with the curves for AC. The high swelling strain of AC + 60% PEG and AC + 100% PEG mixtures at high RH indicates that the swelling of PEG at high relative humidity determines the mixture behavior, since AC, which is stiffer than PEG, does not hinder the swelling of PEG anymore. In short, the overall trend of sorption/swelling results suggests that PEG treatment results in a reduction of sorption and swelling at lower relative humidity followed by an increase in sorption and swelling at high relative humidity, the magnitude of which depends on PEG content.

Interestingly, AC + 15%PEG, AC + 30%PEG, and AC + 60% PEG mixtures show at high RH a moisture content lower than that of both components. This observation indicates that these mixtures do not follow the classical mixing rule, where the resulting MC is the average of both components according to their respective volume fractions. This observation is also seen in experimental studies on PEG-treated waterlogged Scots pine [46,47], PEG-treated waterlogged Elm [26], archaeological oak [58] and recent oak [7].

As shown in Fig. 5c and Table 2, the components and their mixtures show that the swelling coefficient $\beta_A$ in the first part is smaller than the

swelling coefficient $\beta_B$ in the second part. Moreover, the swelling curves of the mixtures show a consistent shift to lower moisture contents, indicating an initial pore-filling phase of less capacity compared to the one of pure AC. The swelling coefficients, $\beta_B$, of all mixtures are close to 1, with a slight increase in $\beta_B$ with increasing PEG content.

Sorption tests at constant volume ($\mu$VT), referred to as restrained

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Sorption isotherms for AC, PEG and AC-PEG mixtures in restrained cases.</figcaption></figure>
cases, are performed to investigate the effect of no-swelling on sorption. Fig. 6 shows the sorption isotherms for the restrained case for AC, PEG and AC-PEG mixtures. The data are approximated by an exponential function (solid line). In these restrained cases, the maximum MC is substantially lower compared to the free swelling case (shown in Fig. 5a), which indicates that swelling and sorption are directly related to each other in the free swelling case. In the restrained case, sorption is related to the filling of the initial available porosity and the sorption of additional water molecules in the polymer matrix due to some chain rearrangements. The restrained sorption isotherms show a gradual reduction in MC with increasing PEG content. This is due to the fact that, at high PEG content, more of the initial porosity of AC is filled by PEG. Interestingly at high relative humidity, mixtures show a lower MC than for pure PEG does. A reason might be that the short PEG chains show more ability to dislocate and deform in order to create space for additional water molecules in pure PEG than when confined by AC.

## 3.2. Mechanical properties of AC, PEG and AC-PEG mixtures

Fig. 7a-b show the bulk and Young's moduli for AC, PEG and AC-PEG mixtures as a function of moisture content. In a first approximation, we approximate the data by a power law function following

$$K = K_0 \exp(-a_K \varepsilon), E = E_0 \exp(-a_E \varepsilon), \quad (10)$$

where $\varepsilon$ is the strain, $K_0$ and $E_0$ are the elastic properties of the dry material and $A_K$ and $A_E$ are the power exponents reported in Table S2 of

the supplementary information. The results show that AC, PEG and AC-PEG mixtures exhibit a moisture-induced weakening. An additional weakening effect is due to the PEG treatment. For instance, mixtures with 30 % PEG content show a 21 % reduction in elastic modulus compared to pure AC, while an increase of MC of 30 % leads to a 50 % drop in elastic modulus concerning AC. The weakening induced by adding PEG to AC can be attributed to the very low mechanical modulus of PEG compared to the stiffer AC. We found that the mechanical moduli of AC-PEG mixtures at dry state decrease linearly with PEG content, following a rule of mixture (see Fig. S2 in supporting information). Fig. 7c shows the breakage of AC-AC hydrogen bonds with increasing moisture and PEG content. As PEG molecules are mixed with AC, they push away the AC chains and AC-AC hydrogen bonds are partially replaced with AC-PEG hydrogen bonds (Fig. 7d). Moisture-induced weakening is shown to be associated with the breakage of hydrogen bonds between AC molecules through the separation of polymer chains upon moisture adsorption and swelling [38]. We note that the MD/GCMC data follow the weakening trend of cell walls as observed in PEG-treated archaeological oak from the Vasa warship [75].

For AC, the moisture-induced weakening effect is not observed at low moisture content (MC < 0.05), and even a small stiffening effect is present. This stiffening effect is related to water molecules filling the initial porosity and then acting as bridges binding the polymer chains by forming hydrogen bonds [84]. This phenomenon is less noticeable for AC-PEG mixtures for which the initial porosity of AC is already partly filled by PEG.

<figure><img src="image_10.png" /><figcaption>Fig. 7. Simulation results of mechanical properties and hydrogen bond density of AC, PEG, and AC-PEG mixtures versus moisture content: (a) bulk modulus, with fitted functions (b) Young's modulus, with fitted functions (c) AC-AC hydrogen bond density (d) AC-PEG hydrogen bond density.</figcaption></figure>
### 3.3. Porosity and pore size distributions of AC, PEG and AC-PEG mixtures

Fig. 8a shows the porosity accessible for water in the dry materials as a function of the PEG content. Porosity rapidly decreases up to approximately 15 % PEG content, followed by a less steep decrease. The first rapid decrease in porosity is attributed to the filling of the initial porosity of AC with PEG. The less steep decrease for PEG content >15 % is attributed to lower porosity of the added PEG decreasing the average. The decrease of porosity at increasing PEG content is also reflected in the volumetric density of AC-PEG mixtures. As it is shown in Fig. S3 in supplementary material, the density of mixtures is reduced by increasing PEG content. The density of the mixture $ρ_{mix}$ can simply be determined by the rule of mixture (RoM) where $f_{vi}$ and $ρ_i$ are the volume fraction and density of component $i$, respectively. However, at low PEG content, RoM underpredicts the density of the mixture since the mixture model does not account for the decrease in the initial porosity of AC by its filling with PEG.

The pore size distribution of the different AC-PEG mixtures is shown in Fig. 8b and is described by a three-parameter Weibull distribution function. PEG is mostly populated by pores in a range from 2.5 Å to 4 Å in diameter. In comparison, AC shows a wider pore size spectrum with larger diameters ranging from 2.6 Å to 5 Å. By increasing the PEG content in a mixture, the pore size spectrum shrinks, removing large pores of 4 Å to 5 Å. This means that PEG reduces the porosity of AC and most likely occupies the larger pores of AC. Similar to what was seen with sorption isotherms, the most pronounced impact is seen at around

15 % PEG content, where most of the initial large pores are filled with PEG. Porosity and pore size distribution are also reported for hydrated samples in Figs. S4 and S5 in the supplementary material. It is found that the porosity increases due to sorption-induced swelling indicating the appearance of larger pores, revealed by a slight shift of pore size distribution to higher diameters.

Fig. 8c shows the porosity accessible to water molecules for the restrained case, where the sample volume remains constant (no swelling strain). At RH = 0, the porosity is equal to around 4.8 %. Remarkably, the porosity increases even when the sample is restrained, which is explained by the fact that the polymer chains deform and rearrange upon hydration forming extra pore space accessible to water molecules.

## 4. New mixture model and discussion

### 4.1. Rule of mixture (RoM) model

Fig. 9 shows the RoM results for sorption isotherms of AC-PEG mixtures in restrained case. Note that, since swelling is restrained, moisture adsorption is only due to porosity filling. We observe that RoM markedly overpredicts the MC of the mixture. This overprediction in MC stems from the RoM model not taking into account the reduction in porosity when the polymers are mixed and PEG partially fills the initial porosity of AC.

Figs. 10 and 11 show the RoM results for sorption and swelling isotherms of AC-PEG mixtures in the free swelling case. For all mixtures, RoM predicts higher MC and swelling compared to MD/GCMC data

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Simulation results of (a) porosity of AC-PEG mixture as a function of PEG volume fraction (0 = pure AC, 1 = pure PEG), fitted with function (Eq. (19)); (b) pore size distribution of AC, PEG, and AC-PEG mixtures, (c) total accessible porosity as a function of relative humidity in the restrained swelling case, and (d) schematic representation of porosity of mixture. Each material mass is shown as a filled area and its pore space is represented by the hatched area.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 9. Modeling the sorption isotherms in the restrained case for different AC-PEG mixtures: (a) AC + 15%PEG, (b) AC + 30%PEG, (c) AC + 60PEG, (d) AC + 100% PEG. MD/GCMC results (symbols), RoM prediction (dashed line), new mixture model (solid line). The curve for AC is in solid blue, and PEG in solid red.</figcaption></figure>

suggesting that AC and/or PEG molecules in the mixture absorb and swell less compared to their pure form. At low RH, similarly to what is seen in the restrained case, the overprediction in MC is due to disregarding the filling of the porosity of AC by PEG. At higher relative humidity, where sorption is mainly due to swelling, the mixture swelling behavior approaches that of cellulose. This suggests that AC effectively inhibits PEG swelling capacity within the mixture.

In the following, we will develop a new mixture model which takes into account that PEG fills part of the pore space of AC and that there is a coupling between sorption and swelling, the so-called sorption-induced swelling. The approach is presented in three parts. The filling of the AC porosity by PEG for dry mixtures is described in Section 4.2 New porosity model. The coupling between adsorption and swelling for pure AC and PEG is described in Section 4.3 New hygromechanical model. Finally, the New mixture model is presented in Section 4.4.

## 4.2. New porosity model

We distinguish between total porosity and the porosity accessible by water molecules. The water-accessible porosity is the pore space that can accommodate a sphere with a diameter equal to the equivalent diameter of a water molecule, i.e. 2.75 Å. In the following, we will refer to water-accessible porosity simply as porosity. We remark that part of the total porosity, which is initially inaccessible to water, will become accessible to water molecules, as swelling or rearrangement of polymer chains

occurs.

In Fig. 8a, we found that the available porosity for water adsorption in the initial dry state reduces as more PEG is mixed with AC. PEG is assumed to fill a part of the pore volume of AC up to a certain limit $\epsilon_{0,AC} V_{pore,AC}$, where $\epsilon_{0,AC}$ is referred to as filling limit varying between 0 and 1. For $\epsilon_{0,AC} = 0$, PEG does not occupy the pore space of AC, while for $\epsilon_{0,AC} = 1$ the total initial pore volume of AC can be filled with PEG. The volume of PEG that does not fill the initial porosity of AC is assumed to retain the porosity of its pure form. The part of the volume of PEG that fills the pore space of AC is equal to $(1 - \delta)V_{0,PEG}$, where $\delta$ is referred as the filling ratio. Fig. 8d shows a schematic representation of a mixture of AC and PEG and the filling of AC porosity by PEG.

We now determine the porosity of the mixture as a function of PEG content. When PEG does not totally occupy the available pore volume of AC (superscript av), or $(1 - \delta)V_{0,PEG} < \epsilon_{0,AC}V_{0,AC}^{av}$, the total volume of the mixture equals

$$V_{0,M} = V_{0,AC} + \delta V_{0,PEG} \quad (11)$$

where $V_{0,M}$, $V_{0,AC}$ and $V_{0,PEG}$ are the dry (subscript 0) volumes of the mixture, AC and PEG respectively. The pore volume of the dry mixture $V_{0,M}^{av}$ is

$$V_{0,M}^{av} = \phi_{0,AC}^{av} V_{0,AC} - (1 - \delta) V_{0,PEG} + \phi_{0,PEG}^{av} V_{0,PEG} \quad (12)$$
<figure><img src="image_13.png" /><figcaption>Fig. 10. Modeling the sorption isotherms in the free swelling case for different AC-PEG mixtures: (a) AC + 15%PEG, (b) AC + 30%PEG, (c) AC + 60PEG, (d) AC + 100%PEG. MD/GCMC results (symbols), higher and lower bounds of RoM prediction (dashed line), new mixture model (solid line). The curve for AC is in solid blue, and PEG in solid red.</figcaption></figure>

where $\phi_{0,AC}^{av}$ and $\phi_{0,PEG}^{av}$ are the available porosity of AC and PEG, respectively. Using Eqs. (3), (11) and (12), the porosity of the mixture (subscript M) $\phi_{0,M}^{av}$ becomes

$$ \phi_{0,M}^{av} = \frac{V_{0,M}^{av}}{V_{0,M}} = \frac{\phi_{0,AC}^{av} + \alpha(\delta + \phi_{0,PEG}^{av} - 1)}{1 + \delta\alpha} \quad (13) $$

When the initial pore space in AC is totally filled by PEG, $(1 - \delta)V_{0,PEG} > \epsilon_{0,AC}V_{pore,AC}$, the total volume, pore volume and porosity of the mixture are, respectively:

$$ V_{0,M} = V_{0,AC} - \epsilon_{0,AC}\phi_{0,AC}^{av}V_{0,AC} + V_{0,PEG} \quad (14) $$

$$ V_{0,M}^{av} = (1 - \epsilon_{0,AC})\phi_{0,AC}^{av}V_{0,AC} + \phi_{0,PEG}^{av}V_{0,PEG} \quad (15) $$

$$ \phi_{0,M}^{av} = \frac{V_{0,M}^{av}}{V_{0,M}} = \frac{(1 - \epsilon_{0,AC})\phi_{0,AC}^{av} + \alpha\phi_{0,PEG}^{av}}{1 - \epsilon_{0,AC}\phi_{0,AC}^{av} + \alpha} \quad (16) $$

Using this approach, the prediction of the mixture porosity as a function of the PEG content for $\delta = 0.85$ and $\epsilon_{0,AC} = 0.6$ is shown in Fig. 8a. A good agreement between the new pore model predictions and MD/GCMC results is observed. The value of $\delta$ indicates that 15 % of the volume of PEG fills the pore space of AC. The value of $\epsilon_{0,AC}$ indicates that PEG is filling the initial pore space of AC up to 60 %. The model predicts the sharp decrease of porosity up to 15 % PEG as well as the relatively small change in porosity hereafter. Note that, in a more general approach, one can assume that AC also fills the porosity of PEG by

considering different $\delta$ and $\epsilon$ for each component in Eqs. (3)–(8). However, doing so shows negligible filling limit and filling ratio ($\delta$ and $\epsilon$) for amorphous cellulose filling the PEG. The filling limit $\epsilon_{0,PEG} \simeq 0$ means that the AC is not filling the initial porosity of PEG, which is attributed to the low molecular size and porosity of PEG.

### 4.3. New hygromechanical model for AC and PEG considering coupling between sorption and swelling

We now develop a hygromechanical model for the coupling between sorption and swelling for the pure components, namely AC and PEG. We assume that, in the restrained case (subscript R), the total moisture content (MC) as a function of RH can be described by

$$ MC_R(RH) = MC_R(max)f(RH, c) \quad (17) $$

where $MC_R(max)$ is the maximum moisture content at RH = 1, in the restrained case. The function $f(RH, c)$ is chosen as

$$ f(RH, c) = \frac{1 - \exp(-c.RH)}{1 - \exp(-c)} \quad (18) $$

with $c$ a parameter to be determined from simulation data, and the function $f$ reaches 1 at RH = 1. In Fig. 8c, we found that the porosity in the restrained case increases with increasing RH, which is attributed to the additional pore space $\phi_R^{ad}$ that becomes available during adsorption as water molecules reconfigure and deform the polymer chains, or:
<figure><img src="image_14.png" /><figcaption>Fig. 11. Modeling the swelling isotherms in the free swelling case for different AC-PEG mixtures: (a) AC + 15%PEG, (b) AC + 30%PEG, (c) AC + 60PEG, (d) AC + 100%PEG. MD/GCMC results (symbols), RoM prediction (dashed line), new mixture model (solid line). The curve for AC is in solid blue, and PEG in solid red.</figcaption></figure>

$$ \phi_R^{av}(RH) = \phi_0^{av} + \phi_R^{ad}f(RH, d) \quad (19) $$

where $\phi_R^{av}$ is the total available porosity (superscript av) to water molecules in the restrained case, $\phi_0^{av}$ is the initial available porosity in the dry state and $\phi_R^{ad}$ is the additional porosity that becomes available (superscript ad) during adsorption in the restrained case and d is a parameter to be determined from simulation data. The prediction of porosity change (dashed line in Fig. 8c) shows a reasonable agreement considering the scatter in the MD/GCMC results.

The maximum moisture content at RH = 1 in the restrained case equals:

$$ MC_R(RH) = \rho_R \phi_R^{av}(RH) \quad (20) $$

where $\rho_R$ is the density of water in the restrained case.

In the free swelling or unrestrained case, the total porosity $\phi$ consists of the total water accessible porosity in the restrained case and the new porosity created during sorption-induced swelling, which equals the volumetric swelling strain $\varepsilon(RH)$, or

$$ \phi(RH) = \phi_R^{av}(RH) + \varepsilon(RH) \quad (21) $$

The moisture content in the unrestrained case equals

$$ MC(RH) = \rho_U \phi(RH) \quad (22) $$

with $\rho_U$ the density of water during unrestrained adsorption. Using Eqs.

(20), (21) and (22), the swelling strain becomes:

$$ \varepsilon(RH) = \frac{MC(RH)}{\rho_U} - \frac{MC_R(RH)}{\rho_R} \quad (23) $$

This equation describes the coupling between swelling strain and moisture content. It shows that swelling is related to the moisture which is adsorbed in the material during swelling, where the inverse of the density of water can be interpreted as the swelling coefficient.

The moisture content in the unrestrained case, $MC(RH)$, is described by the function $f_S$ in Eq. (9), whose parameters for AC and PEG are reported in Table S1. The determined values of the parameters that describe the restrained sorption curve, c and $MC_R(max)$, the density of water in the restrained and unrestrained case $\rho_R$ and $\rho_U$ and parameter d are reported in Table S3. Results indicate that the density of water in the restrained case is higher than in the unrestrained free swelling for both AC and PEG, which is attributed to the high confinement of the water molecules in the restrained case. Figs. 12 and 13 a-d compare the predictions by the new hygrothermal model and the MD/GCMG results for AC and PEG, respectively. A very good agreement is observed showing that the proposed model is able to describe the coupling between sorption and swelling very well. The model predicts a nonlinear swelling curve for AC (Fig. 12c), while the swelling curve for PEG is found to be linear (Fig. 13c), which is reaffirmed by the MD/GCMG results.

We now introduce a model to describe the moisture dependence of
<figure><img src="image_15.png" /><figcaption>Fig. 12. Hygromechanical properties of AC: (a) sorption isotherm, (b) swelling isotherm, (c) swelling curve, (d) restrained sorption isotherm, (e-f) Bulk and Young's moduli as a function of moisture content. MD/GCMC results (symbols), hygromechanical model (dashed line).</figcaption></figure>

the mechanical properties based on the proposed hygromechanical model for sorption-induced swelling. The general weakening of the material is described by an exponential function dependent on the change in porosity which equals the swelling strain, or $K = K_0 \exp(-a_K \epsilon)$, with $a_K > 0$. It is found that the mechanical modulus of AC shows a stiffening effect at low MC, which is attributed to the attraction forces exerted by water molecules that fill the initial porosity $\phi_0$. To describe both the stiffening and weakening effects, we correct for the initial porosity and the newly created available porosity $\phi_R^{av}$ in the restrained case, leading to the following equation.

$$K = K_0 \exp(-a_K \epsilon - b_K \phi_0^{av} + b_K \phi_R^{av} (RH)) \quad (24)$$

The values of the parameters that describe the stiffening/weakening, $K_0$, $a_K$, $b_K$ and $d$ are reported in Table S3. Figs. 12e-f and 13e-f compare the mechanical properties predicted by the new hygromechanical model

and the MD/GCMG results for AC and PEG, respectively. The model is observed to describe the trend of the data considering its large scatter. We see that the model predicts a stiffening at low MC for AC followed by a weakening effect, while the mechanical properties for PEG remain almost constant.

In conclusion, the proposed hygrothermal model allows us to describe the nonlinear coupling of the hygrothermal properties of AC and PEG adequately.

### 4.4. New hygromechanical mixture model (HMM)

In the new hygromechanical mixture model (HMM) proposed here, we consider that PEG will partially invade the initial porosity of AC, as was found in the new porosity model, while AC does not invade the pore space of PEG ($\epsilon_{PEG} \simeq 0$). The derivations for the available porosity of the
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Hygromechanical properties of PEG: (a) sorption isotherm, (b) swelling isotherm, (c) swelling curve, (d) restrained sorption isotherm, (e-f) Bulk and Young's moduli as a function of moisture content. MD/GCMC results (symbols), hygromechanical model (dashed line).</figcaption></figure>

mixture can be found in the supplementary material. The available porosity of the mixture (subscript M) reads:

$$ \phi_{R,M}^{av} = \frac{(1 - \epsilon_{AC})}{1 - \epsilon_{AC}\phi_{0,AC}^{T} + \alpha}\phi_{R,AC}^{av} + \frac{\alpha}{1 - \epsilon_{AC}\phi_{0,AC}^{T} + \alpha}\phi_{R,PEG}^{av} \quad (25) $$

considering the density of water $\rho_R$ in the restricted case for each component in Eq. (20) we can calculate the moisture content of the mixture as follows.

$$ MC_{R,M}(RH) = \frac{(1 - \epsilon_{AC})}{1 - \epsilon_{AC}\phi_{0,AC}^{T} + \alpha}MC_{R,AC}(RH) + \frac{\alpha}{1 - \epsilon_{AC}\phi_{0,AC}^{T} + \alpha}MC_{R,PEG}(RH) \quad (26) $$

where $\phi_{0,AC}^T$ is the total initial porosity of AC. We note that for $\epsilon_{AC} = 0$ the new mixture model reduces to the RoM. Fig. 9 shows a good agreement

between MD/GCMC data and the new model predictions. The obtained values for $\epsilon_{AC}$ are given in Table 4. As an average for all mixtures, we find $\epsilon_{AC} \approx 0.57$ which is close to the value we found in the new porosity model ($\epsilon_{AC} = 0.6$).

For the free swelling case, we assume the components to be in a parallel configuration, where the component with higher stiffness (AC)

<table><caption>**Table 4** Fractions of AC filled by PEG for different AC-PEG mixtures.</caption><thead><tr><th>Mixture</th><th>ε<sub>AC</sub></th></tr></thead><tbody><tr><td>AC + 15%PEG</td><td>0.4573</td></tr><tr><td>AC + 30%PEG</td><td>0.5371</td></tr><tr><td>AC + 60 %PEG</td><td>0.6158</td></tr><tr><td>AC + 100%PEG</td><td>0.6846</td></tr></tbody></table>
can hinder the swelling of the softer component (PEG). As a result, AC will reduce the swelling, especially at high relative humidity. This hindering effect by AC is referred to as confinement on PEG by AC. The swelling strain of the mixture in parallel composition is [70]:

$$ \epsilon_M = \frac{f_{v,AC}^* K_{AC} \epsilon_{AC} + f_{v,PEG}^* K_{PEG} \epsilon_{PEG}}{f_{v,AC}^* K_{AC} + f_{v,PEG}^* K_{PEG}} \quad (27) $$

where the modified volume fractions $f_{v,AC}^*$ and $f_{v,PEG}^*$ take into account the pore filling in the mixture.

$$ f_{v,AC}^* = \frac{1 - \epsilon_{AC} \phi_{0,AC}^{av}}{1 - \epsilon_{AC} \phi_{0,AC}^T + \alpha} \quad (28) $$

$$ f_{v,PEG}^* = \frac{\alpha}{1 - \epsilon_{AC} \phi_{0,AC}^T + \alpha} \quad (29) $$

The derivations for modified volume fractions $f_{v,AC}^*$ and $f_{v,PEG}^*$ can be found in the supplementary material. The total porosity $\phi$ consists of the total water accessible porosity in the restrained case and the new porosity created during sorption-induced swelling, which equals the volumetric swelling strain $\epsilon(RH)$

$$ \phi_M^{av} = \epsilon_M + \phi_{RM}^{av} \quad (30) $$

Using Eq. (25) we can write

$$ \phi_M^{av} = \epsilon_M + \frac{(1 - \epsilon_{AC})}{1 - \epsilon_{AC} \phi_{0,AC}^T + \alpha} \phi_{R,AC}^{av} + \frac{\alpha}{1 - \epsilon_{AC} \phi_{0,AC}^T + \alpha} \phi_{R,PEG}^{av} \quad (31) $$

The total moisture content in the mixture $MC_M$ equals the sum of moisture content $MC_{U,M}$ in the unrestrained case (superscript U) due to swelling and the moisture content $MC_{R,M}$ in the restrained mode (subscript R):

$$ MC_M(RH) = MC_{U,M}(RH) + MC_{R,M}(RH) \quad (32) $$

$MC_{R,M}$ is calculated using Eq. (26). $MC_{U,M}$ is determined based on the swelling of the mixture and the modified volume fractions given in Eqs. (28)–(29)

$$ MC_{U,M}(RH) = f_{v,AC}^* \rho_{U,AC} \epsilon_M(RH) + f_{v,PEG}^* \rho_{U,PEG} \epsilon_M(RH) \quad (33) $$

Figs. 9, 10 and 11 compare the sorption and swelling isotherms determined by the new hygromechanical mixture model (HHM) with MD/GCMC data and the classic RoM. The HMM sorption isotherms show a reduced moisture content and swelling compared to the RoM and are in better agreement with the MD/GCMC data. HMM especially improves the RoM predictions for the case of restrained swelling, where the role of porosity filling is prominent. For the case of free swelling, the enhancement of HMM relies on incorporating both the porosity reduction of AC due to filling by PEG and the improved estimation of swelling behavior. HMM improves the prediction of swelling isotherms by considering the restriction of PEG swelling by the stiff surrounding AC chains assuming a parallel configuration of AC and PEG which leads to reduced hygroscopicity of PEG in the mixture.

One can consider a simplified HMM where AC and PEG are considered in a parallel configuration, but pore filling is neglected ($\epsilon_{AC} = \epsilon_{PEG} = 0$). Consequently, as shown in Figs. 10 and 11, the moisture content becomes higher and displays a greater deviation from MD/GCMC results at low relative humidity. However, both HMM and simplified HMM work well at high PEG content (e.g. AC + 100%PEG), where the difference in volume fractions of the two models becomes smaller. Both mixture models still overestimate the moisture adsorption and swelling, which could hint at the role an AC-PEG interphase may play in the mixture. Such an interphase model is presented in [84] by the authors for the modeling of the hygromechanical behavior of the S2 layer of wood. To determine the influence of such an interphase, the number of hydrogen bonds between AC and PEG as a function of PEG content has to be determined. The development of such an interphase model is part of

future research.

## 5. Discussion

Our findings reveal reduced moisture adsorption and sorption-induced swelling at museological conditions. This underpins the consolidating properties of PEG treatment at the molecular scale, which is consistent with experimental observations. We identified a crossover behavior at high relative humidity, where moisture content and swelling are drastically enhanced due to the high hygroscopicity of PEG in humid conditions. Moreover, we introduced a new hygromechanical model for our pure components and a modified mixture model to explain the porosity and hygromechanical behavior of AC + PEG mixtures.

The new mixture model confirms the role of PEG not only as a bulking agent, as suggested by previous experimental studies [25], but also as an active hygroscopic agent that interacts and intertwines with the polymeric system through a complex mechanism that influences its coupled hygroscopic, mechanical, and structural properties altogether. The proposed models incorporate sorption, sorption-induced swelling, porosity filling, and mechanical softening based and building upon a RoM. However, fully accounting for PEG-cell wall interactions would require a more elaborate model, such as a poromechanical model, which could be a focus of future research.

In the current study, AC was considered as the model material representing wood holocelluloses to simplify and refine the chemical heterogeneity inherent to cell wall composition. A logical next step would be to investigate the impact of PEG on the complex chemical structure of wood holocelluloses by incorporating more advanced model systems composed of a complex composite of hemicelluloses and lignins as these components show drastically different hygromechanical properties and functionality in the wood cell wall and interact differently with PEG. Moreover, in this study, we focused solely on the adsorption branch of sorption isotherms (no desorption). PEG is known to substantially diminish the sorption hysteresis in wood which, for the moment, lacks a proper nanoscale explanation. The current MD/GCMC methodology can be used to study the desorption or scanning branches in PEG-treated mixtures by starting from hydrated systems reducing the RH in GCMC simulations and equilibrating the simulation box in MD to allow for shrinkage. This study also does not account for the effects of temperature and external mechanical stresses on polymer-PEG interactions, which would provide valuable insights in future research. Furthermore, crystalline cellulose (CC) plays a significant role in the hygromechanical properties of wood. It acts as a mechanical reinforcement that stiffens the composite and limits the swelling along the fiber direction. Additionally, the CC-matrix interphase enhances sorption and swelling due to higher porosity at the interface [12]. One can expect a non-homogenous distribution of PEG with higher concentration at the matrix-fiber interphase, leading to a weakening of the CC reinforcement effect, which can be addressed in future research. The modeling framework established in this paper can be also applied to other potential conservation materials. Comparing the atomistic properties of different consolidants and their respective consolidation mechanism and interactions with wood cell wall can inspire novel conservation material with elevated hygromechanical properties.

## 6. Conclusions

This study provides valuable molecular-scale insights into the hygromechanical properties of cellulose treated with the polyethylene glycol (PEG) consolidant by employing an iterative hybrid all-atom molecular dynamics and grand canonical Monte-Carlo (GCMC) simulation framework. The molecular mixture model consisting of AC and PEG hydrated at different RH levels informs the role of PEG treatment as a consolidant agent, reducing the moisture content and swelling of AC at museological conditions, followed by a crossover behavior at high PEG contents and humid conditions (RH > 80) where the mixture undergoes
excessive sorption and swelling beyond that of untreated samples. Examining the mechanical properties and hydrogen bonding network uncovers mechanical softening induced by PEG treatment. The simulation results for pure material and mixtures agree with experimental observation.

The molecular modeling results are then complemented by comparison with the rule of mixture (RoM) analysis. RoM is shown to highly overestimate the moisture content and swelling of the treated samples missing the key features by which PEG consolidates and improves the hygroscopicity of wood. Therefore, accounting for the missing features, a new hygromechanical model for pure components and a modified hygromechanical mixture model are proposed which incorporate sorption and swelling coupling, porosity filling and mechanical softening unified in a single new model for AC + PEG mixtures to reveal the role of AC-PEG interaction. The two key mechanisms identified are: 1) the filling of the AC porosity by PEG and by which PEG reduces the available sorption sites for water molecules 2) the restriction of PEG polymers by the surrounding structure, which prevents further swelling and reduces the hygroscopicity of PEG. At high RH, the consolidation properties are shown to unfavorably diminish as PEG shows increased hygroscopicity. The models are shown to accurately predict both the restricted and free swelling cases, with significant improvement over RoM.

Our findings on the physics of polymeric mixtures and polymer-water interactions will contribute to the fields of material science and biopolymers. Despite the widespread application of PEG treatment for the conservation and preservation of waterlogged archaeological wooden artifacts, the molecular scale physics of PEG treatment lacks a proper explanation. Moreover, polymeric blends (mixtures) and their diverse properties are an active area of research in polymer science. Our hygromechanical mixture model and our observation of interaction mechanisms of PEG and AC can be insightful in developing new polymeric mixture models with specified mechanical and hygroscopic properties. The potential applications of the identified molecular mechanisms and the provided hygromechanical and modified mixture models extend beyond PEG treatment, providing a foundation for understanding the nanoscale physics of alternative consolidants and the hygromechanical behavior of natural biopolymers. This knowledge can inform the conservation and display of wooden artifacts in museums, guide decision-making on the best methods for treating and preserving these invaluable objects, and improve the overall practice of preserving archaeological wood and wooden heritage. Needless to say, the scope of impregnation treatment application goes far beyond waterlogged archaeological wood and expands to dry archaeological wood, cork, fossil wood, other plant-based material, and many other invaluable cultural artifacts such as paintings, wooden statues, ancient buildings, bridges, etc.

### CRediT authorship contribution statement

**Ali Shomali:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Chi Zhang:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Benoit Coasne:** Writing – review & editing, Validation, Methodology, Formal analysis.
**Eleanor J. Schofield:** Writing – review & editing, Validation, Formal analysis.
**Brigitte Chabbert:** Writing – review & editing, Validation, Formal analysis.
**Dominique Derome:** Writing – review & editing, Supervision, Resources, Project administration, Investigation, Funding acquisition, Formal analysis, Conceptualization.
**Jan Carmeliet:** Writing – review & editing, Supervision, Resources, Project administration, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare no conflict of interests.

### Acknowledgments

Acknowledgment of funding by project 'intoS2' no. 182512 funded by Swiss National Science Foundation together with French National Research Agency, PIs Derome Chabbert. Zhang acknowledges the Fundamental Research Funds for the Central Universities, Sun Yat-sen University, under Grant No. 24pny247.

### Appendix A. Supplementary data

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

### References

[1] M.J. Abraham, T. Murtola, R. Schulz, et al., GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, Software X 1 (2015) 19–25.

[2] A.Y. Al-Maharma, F. Bamer, B. Markert, Molecular dynamics investigation on the effect of MWCNTs on thermo-mechanical properties of amorphous cellulose reinforced with silicon carbide nanotubes, IOP Conf. Ser. Mater. Sci. Eng. 1190 (2021) 012026, https://doi.org/10.1088/1757-899x/1190/1/012026.

[3] G. Almkvist, I. Persson, Degradation of polyethylene glycol and hemicellulose in the vasca, Holzforschung 62 (2008) 64–70, https://doi.org/10.1515/HF.2008.009.

[4] F. Antonelli, G. Galotta, G. Sidoti, et al., Cellulose and lignin nano-scale Consolidants for waterlogged archaeological wood, Front. Chem. 8 (2020), https://doi.org/10.3389/fchem.2020.00032.

[5] J.A. Baird, R. Olayo-Valles, C. Rinaldi, L.S. Taylor, Effect of molecular weight, temperature, and additives on the moisture sorption properties of polyethylene glycol, J. Pharm. Sci. 99 (2010) 154–168, https://doi.org/10.1002/jps.21808.

[6] L. Belloni, Finite-size corrections in numerical simulation of liquid water, J. Chem. Phys. 149 (2018), https://doi.org/10.1063/1.5046835.

[7] I. Bjurhager, J. Ljungdahl, L. Wallström, et al., Towards improved understanding of PEG-impregnated waterlogged archaeological wood: a model study on recent oak, Holzforschung 64 (2010) 243–250, https://doi.org/10.1515/HF.2010.024.

[8] M. Broda, C.A.S. Hill, Conservation of waterlogged wood—past, present and future perspectives, Forests 12 (2021), https://doi.org/10.3390/f12091193.

[9] M. Broda, J. Majka, W. Olek, B. Mazela, Dimensional stability and hygroscopic properties of waterlogged archaeological wood treated with alkoxysilanes, Int. Biodeterior. Biodegrad. 133 (2018) 34–41, https://doi.org/10.1016/j.ibiod.2018.06.007.

[10] B. Chabbert, C. Terryn, M. Herbaut, et al., Fluorescence techniques can reveal cell wall organization and predict saccharification in pretreated wood biomass, Ind. Crop. Prod. 123 (2018) 84–92, https://doi.org/10.1016/j.indcrop.2018.06.058.

[11] L. Charlier, K. Mazeau, Molecular modeling of the structural and dynamical properties of secondary plant cell walls: influence of lignin chemistry, J. Phys. Chem. B 116 (2012) 4163–4174.

[12] M. Chen, B. Coasne, D. Derome, J. Carmeliet, Role of cellulose nanocrystals on hysteretic sorption and deformation of nanocomposites, Cellulose 27 (2020) 6945–6960, https://doi.org/10.1007/s10570-020-03247-x.

[13] M. Chen, B. Coasne, R. Guyer, et al., Role of hydrogen bonding in hysteresis observed in sorption-induced swelling of soft nanoporous polymers, Nat. Commun. 9 (2018) 3507, https://doi.org/10.1038/s41467-018-05897-9.

[14] M. Chen, B. Coasne, R. Guyer, et al., Role of hydrogen bonding in hysteresis observed in sorption-induced swelling of soft nanoporous polymers, Nat. Commun. 9 (2018), https://doi.org/10.1038/s41467-018-05897-9.

[15] M. Chen, C. Zhang, A. Shomali, et al., Wood-moisture relationships studied with molecular simulations: methodological guidelines, Forests 10 (2019), https://doi.org/10.3390/f10080628.

[16] P. Chen, J. Wohlert, L. Berglund, I. Furó, Water as an intrinsic structural element in cellulose fibril aggregates, J. Phys. Chem. Lett. 13 (2022) 5424–5430, https://doi.org/10.1021/acs.jpclett.2c00781.

[17] W. Chen, G.C. Lickfield, C.Q. Yang, Molecular modeling of cellulose in amorphous state part II: effects of rigid and flexible crosslinks on cellulose, Polymer (Guildf) 45 (2004) 7357–7365, https://doi.org/10.1016/j.polymer.2004.08.023.

[18] M. Christensen, H. Kutzke, F.K. Hansen, New materials used for the consolidation of archaeological wood-past attempts, present struggles, and future requirements, J. Cult. Herit. 13 (2012) S183–S190, https://doi.org/10.1016/j.culher.2012.02.013.

[19] F.X. Coudert, A.H. Fuchs, Computational characterization and prediction of metal-organic framework properties, Coord. Chem. Rev. 307 (2016) 211–236, https://doi.org/10.1016/j.ccr.2015.08.001.

[20] Q. Deng, S. Li, Y. Chen, Mechanical properties and failure mechanism of wood cell wall layers, Comput. Mater. Sci. 62 (2012) 221–226, https://doi.org/10.1016/j.commatsci.2012.05.050.
[21] Derome D, Kulasinski K, Zhang C, et al (2018) using modeling to understand the hygromechanical and hysteretic behavior of the S2 cell wall layer of wood. In: Plant Biomechanics. Springer, pp. 247–269.

[22] M. Fejfer, J. Majka, M. Zborowska, Dimensional stability of waterlogged scots pine wood treated with PEG and dried using an alternative approach, Forests 11 (2020) 1–13, https://doi.org/10.3390/f1121254.

[23] M. Fejfer, M. Zborowska, O. Adamek, et al., Properties and dimensional stability of 12 500-year-old subfossil pine wood, Drewno 57 (2014) 81–95, https://doi.org/10.12841/wood.1644-3985.075.05.

[24] T.C.F. Gomes, M.S. Skaf, Cellulose-builder: a toolkit for building crystalline structures of cellulose, J. Comput. Chem. 33 (2012) 1338–1346.

[25] D.W. Grattan, A practical comparative study of several treatments for waterlogged wood, Stud. Conserv. 27 (1982) 124, https://doi.org/10.2307/1506148.

[26] D.W. Grattan, A practical comparative study of treatments for waterlogged wood. Part II: the effect of humidity on treated wood, in: Proceedings of the Icom Waterlogged Wood Working Group Conference, Ottawa, 15-18 September 1981 Treated Wood, 1981, pp. 243–252.

[27] J. Guo, J. Chen, R. Li, et al., Thermoporometry of waterlogged archaeological wood: insights into the change of pore traits after the water-removal by supercritical drying, Thermochim. Acta 715 (2022) 179297, https://doi.org/10.1016/j.tca.2022.179297.

[28] Håfors B (2010) Conservation of the wood of the swedish warship vasa of a.d. 1628.

[29] H. Hao, Tam L. Ho, Y. Lu, D. Lau, An atomistic study on the mechanical behavior of bamboo cell wall constituents, Compos. Part B Eng. 151 (2018) 222–231, https://doi.org/10.1016/j.compositesb.2018.05.046.

[30] E. Hocker, G. Almkvist, M. Sahlstedt, The vasa experience with polyethylene glycol: a conservator's perspective, J. Cult. Herit. 13 (2012) S175–S182, https://doi.org/10.1016/j.culher.2012.01.017.

[31] E. Hocker, G. Almkvist, M. Sahlstedt, The vasa experience with polyethylene glycol: a conservator's perspective, J. Cult. Herit. 13 (2012) S175–S182, https://doi.org/10.1016/j.culher.2012.01.017.

[32] J. Jiang, Y. Zhou, C. Mei, J. Cao, Polyethylene glycol and silica sol penetration improves hydrophobicity and dimensional stability of wood after a short-time treatment, Eur. J. Wood Wood Prod. 79 (2021) 1395–1404, https://doi.org/10.1007/s00107-021-01710-5.

[33] S.P.P. Jones, N.K.H. Slater, M. Jones, et al., Investigating the processes necessary for satisfactory freeze-drying of waterlogged archaeological wood, J. Archaeol. Sci. 36 (2009) 2177–2183, https://doi.org/10.1016/j.jas.2009.05.028.

[34] W. Kern, Encyclopedia of Polymer Science and Technology. Plastics, Resins, Rubbers, Fibers. Vol. 10: Phenolic Resins to Polyelectrolytes. Vol. 11: Polymesters to Rayons. Vol. 12: Reinforced Plastics to Starch. Herausgeg. von HF Mark, NG Gaylord und NM Bikales, 1972.

[35] J. Ketoja, S. Paavilainen, J.L. McWhirter, et al., Mechanical properties of cellulose nanofibrils determined through atomistic molecular dynamics simulations, Nord. Pulp Pap. Res. J. 27 (2012) 282–286, https://doi.org/10.3183/npprj-2012-27-02-p282-286.

[36] A. Khodayari, U. Hirn, S. Spirk, et al., Recrystallization and size distribution of dislocated segments in cellulose microfibers—a molecular dynamics perspective, Cellulose 28 (2021) 6007–6022, https://doi.org/10.1007/s10570-021-03906-7.

[37] Y. Kong, L. Li, S. Fu, Insights from molecular dynamics simulations for interaction between cellulose microfibers and hemicellulose, J. Mater. Chem. A 10 (2022) 14451–14459, https://doi.org/10.1039/d2ta03164g.

[38] K. Kulasinski, R. Guyer, D. Derome, J. Carmeliet, Water adsorption in wood microfibril-hemicellulose system: role of the crystalline-amorphous interface, Biomacromolecules 16 (2015) 2972–2978, https://doi.org/10.1021/acs.biomac.5b00878.

[39] K. Kulasinski, R. Guyer, S. Keten, et al., Impact of moisture adsorption on structure and physical properties of amorphous biopolymers, Macromolecules 48 (2015) 2793–2800, https://doi.org/10.1021/acs.macromol.5b00248.

[40] K. Kulasinski, S. Keten, S.V. Churakov, et al., Molecular mechanism of moisture-induced transition in amorphous cellulose, ACS Macro Lett. 3 (2014) 1037–1040, https://doi.org/10.1021/mz500528m.

[41] K. Kulasinski, S. Keten, S.V. Churakov, et al., A comparative molecular dynamics study of crystalline, paracrystalline and amorphous states of cellulose, Cellulose 21 (2014) 1103–1116, https://doi.org/10.1007/s10570-014-0213-7.

[42] G. Kupgan, A.G. Demidov, C.M. Colina, Plasticization behavior in polymers of intrinsic microporosity (PIM-1): a simulation study from combined Monte Carlo and molecular dynamics, J. Membr. Sci. 565 (2018) 95–103, https://doi.org/10.1016/j.memsci.2018.08.004.

[43] J. Ljungdahl, L.A. Berglund, Transverse mechanical behaviour and moisture absorption of waterlogged archaeological wood from the vasa ship, Holzforschung 61 (2007) 279–284, https://doi.org/10.1515/hf.2007.048.

[44] E. Mahmoudinezhad, A. Marquardt, G. Eggeler, et al., Molecular modeling of cellulose in amorphous state. Part I: model building and plastic deformation study, Polymer (Guildf) 45 (2004) 1063–1071, https://doi.org/10.1021/jp0219395.

[45] J. Majka, L. Babinski, W. Olek, Sorption isotherms of waterlogged subfossil scots pine wood impregnated with a lactitol and trehalose mixture, Holzforschung 71 (2017) 813–819, https://doi.org/10.1515/hf-2017-0006.

[46] J. Majka, L. Czajkowski, W. Olek, Effects of cyclic changes in relative humidity on the sorption hysteresis of thermally modified spruce wood, BioResources 11 (2016) 5265–5275, https://doi.org/10.15376/biores.11.2.5265-5275.

[47] J. Majka, M. Zborowska, M. Fejfer, B. Waliszewska, Dimensional stability and hygroscopic properties of PEG treated irregularly degraded waterlogged scots pine wood, J. Cult. Herit. 31 (2018) 133–140, https://doi.org/10.1016/j.culher.2017.12.002.

[48] H.F. Mark, Encyclopedia of Polymer Science and Technology: Phenolic Resins to Polyelectrolytes, Interscience Publishers, 1969.

[49] L. Martínez, R. Andrade, E.G. Birgin, J.M. Martínez, PACKMOL: a package for building initial configurations for molecular dynamics simulations, J. Comput. Chem. 30 (2009) 2157–2164.

[50] J.F. Matthews, G.T. Beckham, M. Bergenstråhle-Wohltert, et al., Comparison of cellulose Iβ simulations with three carbohydrate force fields, J. Chem. Theory Comput. 8 (2012) 735–748.

[51] K. Mazeau, L. Charlier, The molecular basis of the adsorption of xylans on cellulose surface, Cellulose 19 (2012) 337–349, https://doi.org/10.1007/s10570-011-9643-7.

[52] K. Mazeau, L. Heux, Molecular dynamics simulations of bulk native crystalline and amorphous structures of cellulose, J. Phys. Chem. B 107 (2003) 2394–2403, https://doi.org/10.1021/jp0219395.

[53] G. McConnachie, Air Drying Behaviour of Waterlogged Archaeological Woods From the Tudor Warship Mary Rose, 2005.

[54] A. Mihranyan, A.P. Llagostera, R. Karmhag, et al., Moisture sorption by cellulose powders of varying crystallinity, Int. J. Pharm. 269 (2004) 433–442, https://doi.org/10.1016/j.ijpharm.2003.09.030.

[55] M. Mori, S. Kuhara, K. Kobayashi, S. Suzuki, Nondestructive visualization of polyethylene glycol impregnation in wood using ultrashort echo time 3D imaging, J. Cult. Herit. 50 (2012) 43–48, https://doi.org/10.1016/j.culher.2021.05.008.

[56] L. Ninni, M.S. Camargo, A.J.A. Meirelles, Water activity in poly(ethylene glycol) aqueous solutions, Thermochim. Acta 328 (1999) 169–176, https://doi.org/10.1007/s0040-6031(98)00638-8.

[57] Y. Ogawa, Y. Nishiyama, K. Mazeau, Drying-induced bending deformation of cellulose nanocrystals studied by molecular dynamics simulations, Cellulose 27 (2020) 9779–9786, https://doi.org/10.1007/s10570-020-03451-9.

[58] W. Olek, J. Majka, A. Stempin, et al., Hygroscopic properties of PEG treated archaeological wood from the rampart of the 10th century stronghold as exposed in the Archaeological Reserve Genius loci in Poznań (Poland), J. Cult. Herit. 18 (2016) 299–305, https://doi.org/10.1016/j.culher.2015.08.002.

[59] A. Paajanen, A. Rinta-Paaola, J. Vaari, High-temperature decomposition of amorphous and crystalline cellulose: reactive molecular simulations, Cellulose 28 (2021) 9897–9005, https://doi.org/10.1007/s10570-021-04084-2.

[60] A. Paajanen, Y. Sonavane, D. Ignasiak, et al., Atomistic molecular dynamics simulations on the interaction of TEMPO-oxidized cellulose nanofibrils in water, Cellulose 23 (2016) 3449–3462, https://doi.org/10.1007/s10570-016-1076-x.

[61] S. Paavilainen, T. Rög, I. Vattälainen, Analysis of twisting of cellulose nanofibrils in atomistic molecular dynamics simulations, J. Phys. Chem. B 115 (2011) 3747–3755, https://doi.org/10.1021/jp111459b.

[62] E. Piva, Conservation of a Tudor Warship: Investigating the Timbers of the Mary Rose, 2017.

[63] B.D. Rabideau, A. Agarwal, A.E. Ismail, Observed mechanism for the breakup of small bundles of cellulose Iα and Iβ in ionic liquids from molecular dynamics simulations, J. Phys. Chem. B 117 (2013) 3469–3479, https://doi.org/10.1021/jp310225t.

[64] L. Sarkisov, R. Bueno-Perez, M. Sutharson, D. Fairen-Jimenez, Materials informatics with PoreBlazer v4.0 and the CSD MOF database, Chem. Mater. 32 (2020) 9849–9867, https://doi.org/10.1021/acs.chemmater.0c03575.

[65] Schlaich A, Kowalik B, Kanduč M, et al. Simulation Techniques for Solvation-induced Surface-interactions at Prescribed Water Chemical Potential. 159–189.

[66] A.P. Schniewind, Physical and mechanical properties of archaeological wood, Adv. Chem. Ser. 87–109 (1990), https://doi.org/10.1021/ba-1990-0225.cho04.

[67] E.J. Schofield, R. Sarangi, A. Mehta, et al., Of the Mary Rose the preservation of waterlogged archaeological wooden finds, such as, Mater. Today 14 (2011) 354–358, https://doi.org/10.1016/S1369-7021(11)07166-3.

[68] J.K. Shah, E. Marin-Rimoldi, R.G. Mullen, et al., Cassandra: an open source Monte Carlo package for molecular simulation, J. Comput. Chem. 38 (2017) 1727–1739, https://doi.org/10.1002/jcc.24807.

[69] M.R. Shirts, C. Klein, J.M. Swails, et al., Lessons learned from comparing molecular dynamics engines on the SAMPL5 dataset, J. Comput. Aided Mol. Des. 31 (2017) 147–161, https://doi.org/10.1007/s10822-016-9977-1.

[70] N. Srisuk, A Micromechanics Model of Thermal Expansion Coefficient in Fiber Reinforced Composites, The University of Texas at Arlington, 2010.

[71] M. Thommes, K. Kaneko, A.V. Neimark, et al., Physisorption of gases, with special reference to the evaluation of surface area and pore size distribution (IUPAC technical report), Pure Appl. Chem. 87 (2015) 1051–1069, https://doi.org/10.1515/pac-2014-1117.

[72] Van den Berg C, Bruin S (1978) Water activity and its estimation in food systems. In: Proceedings Int. Symp. Properties of Water in Relation to Food Quality and Stability, Osaka, 1978.

[73] J.V. Vermaas, L.D. Dellon, L.J. Broadbelt, et al., Automated transformation of lignin topologies into atomic structures with LigninBuilder, ACS Sustain. Chem. Eng. 7 (2019) 3443–3453, https://doi.org/10.1021/acssuschemeng.8b05665.

[74] J. Vorholz, V.I. Harismiadis, B. Rumpf, et al., Vapor + liquid equilibrium of water, carbon dioxide, and the binary system, water + carbon dioxide, from molecular simulation, Phase-Phase Equilib. 170 (2000) 203–234, https://doi.org/10.1016/S0378-3812(00)00315-0.

[75] A. Vorobyev, G. Almkvist, Van Dijk NP, E.K. Gamstedt, Relations of density, polyethylene glycol treatment and moisture content with stiffness properties of Vasa oak samples, Holzforschung 71 (2017) 327–335, https://doi.org/10.1515/hf-2016-0202.

[76] A. Vorobyev, N.P. van Dijk, E. Kristofer Gamstedt, Orthotropic creep in polyethylene glycol impregnated archaeological oak from the Vasa ship, Mech.
Time-Dependent Mater. 23 (2019) 35–52, https://doi.org/10.1007/s11043-018-9382-3.

[77] L. Wagner, G. Almkvist, T.K. Bader, et al., The influence of chemical degradation and polyethylene glycol on moisture-dependent cell wall properties of archeological wooden objects: a case study of the Vasa shipwreck, Wood Sci. Technol. 50 (2016) 1103–1123, https://doi.org/10.1007/s00226-016-0861-x.

[78] Z. Walsh-Korbs, L. Avérous, Recent developments in the conservation of materials properties of historical wood, Prog. Mater. Sci. 102 (2019) 167–221, https://doi.org/10.1016/j.pmatsci.2018.12.001.

[79] J. Wohlert, M. Bergenstråle-Wohlert, L.A. Berglund, Deformation of cellulose nanocrystals: entropy, internal energy and temperature dependence, Cellulose 19 (2012) 1821–1836, https://doi.org/10.1007/s10570-012-9774-5.

[80] Q. Yang, C. Zhong, Electrostatic-field-induced enhancement of gas mixture separation in metal-organic frameworks: a computational study, ChemPhysChem 7 (2006) 1417–1421, https://doi.org/10.1002/cphc.200600191.

[81] G.S. Young, R. Sims, Microscopical determination of polyethylene glycol in treated wood—the effect of distribution on dimensional stabilization, in: Conservation of Wet Wood and Metal: Proceedings of the ICOM Conservation Working Groups on Wet Organic Archaeological Materials and Metals, 1987, pp. 109–140.

[82] S. Youssefian, J.E. Jakes, N. Rahbar, Variation of nanostructures, molecular interactions, and anisotropic elastic moduli of lignocellulosic cell walls with moisture, Sci. Rep. 7 (2017) 2054, https://doi.org/10.1038/s41598-017-02288-w.

[83] C. Zhang, M. Chen, B. Coasne, et al., Hygromechanics of softwood cellulosic nanocomposite with intermolecular interactions at fiber-matrix interface investigated with molecular dynamics, Compos. Part B Eng. 228 (2022) 109449, https://doi.org/10.1016/j.compositesb.2021.109449.

[84] C. Zhang, M. Chen, S. Keten, et al., Hygromechanical mechanisms of wood cell wall revealed by molecular modeling and mixture rule analysis, Sci. Adv. 7 (2021) 1–18, https://doi.org/10.1126/sciadv.abi8919.

[85] C. Zhang, B. Coasne, R. Guyer, et al., Moisture-induced crossover in the thermodynamic and mechanical response of hydrophilic biopolymer, Cellulose 27 (2020) 89–99, https://doi.org/10.1007/s10570-019-02808-z.

[86] C. Zhang, A. Shomali, R. Guyer, et al., Disentangling heat and moisture effects on biopolymer mechanics, Macromolecules 53 (2020) 1527–1535, https://doi.org/10.1021/acs.macromol.9b01988.

[87] Poly(ethylene glycol) | CAS#:25322-68-3 | Chemsrc. https://www.chemsrc.com/en/cas/25322-68-3_766012.html. Accessed 17 Oct 2022.