Nanopore-Level Wood-Water Interactions—A Molecular Simulation Study

: The nanoscale wood-water interaction strength, accessible sorption sites, and cell wall pore sizes are important factors that drive water sorption and the hysteresis phenomenon in wood. In this work, these factors were quantitatively studied using molecular simulations based on a cell wall pore model, previously developed by the authors. Speciﬁcally, the wall-water interaction strength, the sorption sites network including their number, interaction range, strength, and spatial distributions were set at a series of theoretical values as simulation input parameters. The results revealed that most of the investigated parameters signiﬁcantly affected both sorption isotherms and hysteresis. Water monolayers and clusters were observed on the simulated pore surface when the wood-water interaction and sorption site strength were set at unrealistically high values. Furthermore, multiple linear regression models suggested that wood-water interaction and sorption site parameters were coupled in determining sorption isotherms, but not in determining hysteresis.


Introduction
Owing to bound water's pronounced influence on various physical and mechanical properties of wood and engineered wood products, great efforts have been devoted to study the wood-water fundamentals. A common, but effective way to explore this relationship is to study water sorption isotherms and the associated hysteresis phenomenon. The three major components of wood, namely, cellulose, hemicellulose, and lignin, are hydrophilic, albeit to very different extents with lignin being almost hydrophobic, and the corresponding sorption isotherm displays a sigmoid shape (Type II according to the International Union of Pure and Applied Chemistry (IUPAC) [1]), that shows hysteresis in its entire range.
For almost a century [2], sorption mechanisms and the origin of hysteresis have been extensively studied, however, there are still unanswered questions and some uncertainty [3]. Water layering and clustering are two proposed sorption mechanisms at low to medium relative humidity (H) regions [4][5][6][7], whereas capillary condensation is argued to take place at high or very high H (>98%) regions [8][9][10][11][12]. Corresponding to these three mechanisms, nano-level wood-water interactions, accessible sorption sites (mainly hydroxyl groups) and cell wall pore size are three major factors that drive the sorption process in wood. The roles of accessibility and quantity of sorption sites have been reassessed in some recent studies [13][14][15][16].
The studies on the origin of sorption hysteresis in wood (more specifically, reproducible hysteresis in [17]) are quite diverse. These include, but not limited to the roles of accessible sorption sites during sorption and their possible coupling with sorption-induced swelling [18][19][20][21], the rigidity or softening of wood cell walls [22][23][24], and the metastable states of adsorbed and desorbed water associated with the cell wall nanopore filling [25].
Experimentally, the effects of the wood-water interactions, accessible sorption sites and cell wall pore size on sorption and hysteresis are often investigated by modifying wood using various methods, e.g., thermal or chemical, and then relating the hygroscopic property change to changes of the number of accessible hydroxyl groups, wood-water interaction strength (WWIS) and cell wall pore size distributions [26][27][28][29][30][31][32][33]. However, it is difficult to quantitatively evaluate the aforementioned factors and possible couplings among them. In [25], the effect of cell wall pore size on sorption and hysteresis was explored using molecular simulation methods based on a nanoscale cell wall pore model (details can be found in the following Simulation Model and Method section), and promising results were reported. Comparing to the purely theoretical and experimental approach, the molecular simulation approach has the advantage of conducting systematic and quantitative investigations at the nanoscale using a model that is close to the real system. Molecular simulations of wood cell walls at the atomic level is still challenging due to the implicit cell wall structures, and is very time-consuming which limits its usage in parametric investigations [21].
In this work, the effects of sorption sites and wood-water interactions, their potential coupling, and possible interrelations with cell wall pore sizes were explored by molecular simulations based on a cell wall pore model. By exploring the extreme value scenarios of the investigated factors in the simulations, more insights may be gained into the formation of water monolayers and clusters, and therefore the overall water sorption mechanism in wood.

Simulation Model and Method
The simulation model was taken from our previous publication [25]. The wood cell wall-water system was modeled as a collection of infinitely long cylindrical nanopores of different sizes d j with water confined in them ( Figure 1a). These pores were independent of each other and had full access to the environment. The independence assumption was based on the systematic investigation of experimental water sorption hysteresis patterns [17], and was further supported by a modified Preisach model [34]. The simulated pores were assumed rigid, and the actual swelling of the cell wall and possible generation of new pores during sorption was addressed by assumed cell wall pore size distributions (ψ i ) evolving with H ( Figure 1c, details can be found in [35]). Consequently, pore walls were modeled as rigid ones composed of equally distributed Lennard-Jones (LJ) [36] carbon and oxygen atoms ( Figure 2a).  Figure 1. Illustration of the wood cell wall nanopore model and modelled sorption isotherms. (a) the cell wall substance and nanopores with an initial pore size distribution ψ 1 ; (b) simulated sorption isotherms for different pore sizes d j ; ρ wa represents the adsorbed water density (c) cell wall pore size distributions ψ i evolving with humidity; (d) the eventual wood water sorption isotherm simulated. Figure 1. Illustration of the wood cell wall nanopore model and modelled sorption isotherms. (a) the cell wall substance and nanopores with an initial pore size distribution ; (b) simulated sorption isotherms for different pore sizes ; represents the adsorbed water density (c) cell wall pore size distributions evolving with humidity; (d) the eventual wood water sorption isotherm simulated. WWIS's influence on sorption and hysteresis was investigated by simulating six additional pore walls (PW3-PW8). Varied strengths were set within possible ranges in practice, and at unrealistic extreme values to analyze theoretical trends. PW3 were comprised of pure oxygen atoms, and PW4, pure carbon atoms. WWIS of PW5 to PW7 was enhanced based on default PW1, whereas WWIS of PW8 was decreased based on another default PW2. Table 1 presents the LJ potential parameters of the studied eight wall types. To avoid complications from sorption sites, the number of energy pits on the additional walls was set to 0.  Two default pore walls, PW 1 and PW 2 were defined as representatives associated with amorphous holocellulose and lignin, respectively. The accessible sorption sites were modeled as negative energy pits attached to walls (Figure 2), and water was represented by the extended simple point charge (SPC/E) model [37]. Grand canonical Monte Carlo (GCMC) [38] was used to calculate the quantity of adsorbed water molecules at a series of given H under isothermal conditions for each fixed d j . Figure 1b (Figure 1c). Details can be found in [35].
The interaction energy (u ij ww ) of any two water molecules i and j was calculated as the sum of the LJ potential and the Coulomb potential [39]; the interaction energy of a water molecule with the pore wall (U(s 1 )) was calculated by integrating LJ potential over the wall atoms [40]. The LJ parameters (ε sf and σ sf , representing the well depth at the minimum interaction energy and the separation where the energy equals 0) of PW 1 and PW 2 -water interactions were taken from [41,42]. Due to the limited computational resources, only one representative cell wall pore diameter of 0.95 nm (an approximate average pore diameter) and one temperature level of 25 • C were selected for the simulations.
WWIS's influence on sorption and hysteresis was investigated by simulating six additional pore walls (PW 3 -PW 8 ). Varied strengths were set within possible ranges in practice, and at unrealistic extreme values to analyze theoretical trends. PW 3 were comprised of pure oxygen atoms, and PW 4 , pure carbon atoms. WWIS of PW 5 to PW 7 was enhanced based on default PW 1 , whereas WWIS of PW 8 was decreased based on another default PW 2 . Table 1 presents the LJ potential parameters of the studied eight wall types. To avoid complications from sorption sites, the number of energy pits on the additional walls was set to 0. The influence of accessible sites on sorption was analyzed by controlling the number, strength, interaction range, and spatial distributions of the energy pits on PW 1 and PW 2 . The number of sorption sites n H was controlled by n c and n l (n H = n c × n l ), namely, the number of energy pits in the lateral and longitudinal direction of the simulated pore,  (Figure 2). Since the perimeter of the simulated 0.95 nm pore is limited, n c was kept to be a constant of 3, but values of n l were adjusted. The strength of the pits (ε H ) was set to be 4.5 kcal mol −1 (normal values), 0.5 kcal mol −1 (extreme low values), and 22.5 kcal mol −1 (extreme high values). The interaction range of the energy pits was evaluated by the pit radius R H (Figure 2) that were set at 0.317 nm (normal values), 0.158 nm, and 0.475 nm. The spatial distributions of the energy pit were set "even, random, clustered, and helix" (Figure 3). The helix distributions were inspired by recent findings of two-and threefold helix xylan configurations when hydrogen bonding with cellulose fibrils or lignin [43][44][45]. Table 2 summarizes the assessed settings of energy pits on default PW 1 and PW 2 cylindrical pores. The influence of accessible sites on sorption was analyzed by controlling the number, strength, interaction range, and spatial distributions of the energy pits on PW1 and PW2. The number of sorption sites was controlled by and ( = ), namely, the number of energy pits in the lateral and longitudinal direction of the simulated pore, respectively ( Figure 2). Since the perimeter of the simulated 0.95 nm pore is limited, was kept to be a constant of 3, but values of were adjusted. The strength of the pits ( ) was set to be 4.5 kcal mol −1 (normal values), 0.5 kcal mol −1 (extreme low values), and 22.5 kcal mol −1 (extreme high values). The interaction range of the energy pits was evaluated by the pit radius RH (Figure 2) that were set at 0.317 nm (normal values), 0.158 nm, and 0.475 nm. The spatial distributions of the energy pit were set "even, random, clustered, and helix" (Figure 3). The helix distributions were inspired by recent findings of two-and threefold helix xylan configurations when hydrogen bonding with cellulose fibrils or lignin [43][44][45]. Table 2 summarizes the assessed settings of energy pits on default PW1 and PW2 cylindrical pores.    Figure 4a,b present the simulated water sorption isotherms for the studied eight types of wood cell wall pores at a fixed diameter of 0.95 nm. The relative density of adsorbed water ( ), calculated as the ratio of the density of simulated adsorbed and bulk liquid water at the same temperature and ambient pressure, was used to represent the quantity of adsorbed water molecules.

The Differences between the Shapes of Simulated and Experimental Sorption Isotherms
could be further connected to M of simulated   Figure 4a,b present the simulated water sorption isotherms for the studied eight types of wood cell wall pores at a fixed diameter of 0.95 nm. The relative density of adsorbed water (ρ wa ), calculated as the ratio of the density of simulated adsorbed and bulk liquid water at the same temperature and ambient pressure, was used to represent the quantity of adsorbed water molecules. ρ wa could be further connected to M of simulated wood samples by taking into account the mass and chemical composition of simulated pore walls [35]. The shapes of the simulated sorption isotherms were more like those from non-swelling sorption systems, e.g., carbon-water [46,47], carbon nanotube-water [48], carbon-N 2 /Ar [49], activated carbon-water [50], zeolites-water [51] systems, etc., rather than the ones from the experimental wood water sorption isotherm [23]. This was because only sorption isotherms for one cell wall pore size, 0.95 nm, were presented here. In the "Simulation Model and Method" section, it was pointed that since the cell wall pores were assumed rigid in the model, the actual swelling of cell walls during sorption had to be addressed by assumed ψ i at different H stages. Therefore, to obtain the final simulated sorption isotherm, sorption isotherms at a wide range of pore sizes have to be simulated and then superimposed based on ψ i . This process had been illustrated in detail in [35], and the simulated sorption isotherms could be very close to the experimental ones depending on the choice of ψ i . It is also reported in the literature [52] that a wide pore size distribution can dramatically change the shapes of the sorption isotherms of silica-water systems. The aforementioned process was not done in this work due to the high computation cost of molecular simulations at different pore sizes, especially when several model parameters were also changed. However, this would not affect the main conclusion of this work since the focus here was to evaluate the selected model parameters.

The Differences between the Shapes of Simulated and Experimental Sorption Isotherms
wood samples by taking into account the mass and chemical composition of simulated pore walls [35]. The shapes of the simulated sorption isotherms were more like those from non-swelling sorption systems, e.g., carbon-water [46,47], carbon nanotube-water [48], carbon-N2/Ar [49], activated carbon-water [50], zeolites-water [51] systems, etc., rather than the ones from the experimental wood water sorption isotherm [23]. This was because only sorption isotherms for one cell wall pore size, 0.95 nm, were presented here. In the "Simulation Model and Method" section, it was pointed that since the cell wall pores were assumed rigid in the model, the actual swelling of cell walls during sorption had to be addressed by assumed at different H stages. Therefore, to obtain the final simulated sorption isotherm, sorption isotherms at a wide range of pore sizes have to be simulated and then superimposed based on . This process had been illustrated in detail in [35], and the simulated sorption isotherms could be very close to the experimental ones depending on the choice of . It is also reported in the literature [52] that a wide pore size distribution can dramatically change the shapes of the sorption isotherms of silica-water systems. The aforementioned process was not done in this work due to the high computation cost of molecular simulations at different pore sizes, especially when several model parameters were also changed. However, this would not affect the main conclusion of this work since the focus here was to evaluate the selected model parameters.

The Validity of the Simulation Model
The developed simulation model had been qualitatively examined using the predicted water sorption hysteresis properties. The magnitude of simulated sorption hysteresis increased with the cell wall pore size, the lignin content of cell walls, and the reduced temperature, which is consistent with experimental observations. Details can be found in [25]. Furthermore, the model had been quantitatively examined by comparing the predicted cell wall pore size distributions at fully saturated states with the experimental ones derived from the solute exclusion method. The predicted distributions were relatively wide with several major peaks evolving in the hygroscopic range but were assessed to be fairly reasonable. Details can be found in [35].

The Effect of Wood-Water Interaction Strength on Simulated Sorption Isotherms
The effect of wood-water interaction strength on simulated sorption isotherms was revealed by comparing the eight sorption isotherms in Figure 4. WWIS was reflected by

The Validity of the Simulation Model
The developed simulation model had been qualitatively examined using the predicted water sorption hysteresis properties. The magnitude of simulated sorption hysteresis increased with the cell wall pore size, the lignin content of cell walls, and the reduced temperature, which is consistent with experimental observations. Details can be found in [25]. Furthermore, the model had been quantitatively examined by comparing the predicted cell wall pore size distributions at fully saturated states with the experimental ones derived from the solute exclusion method. The predicted distributions were relatively wide with several major peaks evolving in the hygroscopic range but were assessed to be fairly reasonable. Details can be found in [35].

The Effect of Wood-Water Interaction Strength on Simulated Sorption Isotherms
The effect of wood-water interaction strength on simulated sorption isotherms was revealed by comparing the eight sorption isotherms in Figure 4. WWIS was reflected by the LJ parameter ε sf in Table 1, and larger ε sf values meant stronger wood-water interactions. Hence, for the simulated pores in Figure 4, WWIS decreased in the following order: PW 7 > PW 6 > PW 5 > PW 3 > PW 1 > PW 2 > PW 4 > PW 8 . It was clear from Figure 4 that WWIS affected both sorption isotherms and hysteresis considerably. Overall, stronger wallwater interaction led to lower H micropore fillings (the abrupt increase on the simulated adsorption isotherm), and narrower hysteresis loops (Figure 4). When the interaction became very strong (PW 7 ), hysteresis almost disappeared, and when the interaction became very weak (PW 4 and PW 8 ), water could not condense in the simulated pores. The sorption isotherm from PW 3 (blue lines in Figure 4a), slightly deviated from the general trend since it was supposed to lie between PW 5 and PW 1 (dark green and black lines in Figure 4a), according to the order of WWIS. This was probably caused by the additional sorption sites on the wall of default PW 1 . As explained below, the sorption sites could also influence the isotherms and hysteresis significantly. Therefore, the additional sorption sites on PW 1 pores could shift the micropore filling of water towards an earlier H stage.

The Effect of Sorption Sites on Simulated Sorption Isotherms
Figures 5 and 6 demonstrate how sorption sites affected the simulated sorption isotherms and hysteresis for PW 1 and PW 2 , respectively. In most cases, micropore filling of water occurred in simulated pores. Stronger sorption site networks brought in by either increased n l , ε H or R H caused lower H micropore fillings and narrower hysteresis loops (Figures 5a-c and 6a-c). The isotherm trends from PW 2 ( Figure 6) were similar to that of PW 1 (Figure 5), but the formers were more sensitive to the sorption sites-related parameter variations, indicating the coupling of the sorption sites and WWIS.
PW7 > PW6 > PW5 > PW3 > PW1 > PW2 > PW4 > PW8. It was clear from Figure 4 that WWIS affected both sorption isotherms and hysteresis considerably. Overall, stronger wall-water interaction led to lower H micropore fillings (the abrupt increase on the simulated adsorption isotherm), and narrower hysteresis loops (Figure 4). When the interaction became very strong (PW7), hysteresis almost disappeared, and when the interaction became very weak (PW4 and PW8), water could not condense in the simulated pores. The sorption isotherm from PW3 (blue lines in Figure 4a), slightly deviated from the general trend since it was supposed to lie between PW5 and PW1 (dark green and black lines in Figure 4a), according to the order of WWIS. This was probably caused by the additional sorption sites on the wall of default PW1. As explained below, the sorption sites could also influence the isotherms and hysteresis significantly. Therefore, the additional sorption sites on PW1 pores could shift the micropore filling of water towards an earlier H stage.

The Effect of Sorption Sites on Simulated Sorption Isotherms
Figs. 5 and 6 demonstrate how sorption sites affected the simulated sorption isotherms and hysteresis for PW1 and PW2, respectively. In most cases, micropore filling of water occurred in simulated pores. Stronger sorption site networks brought in by either increased nl, εH or RH caused lower H micropore fillings and narrower hysteresis loops (Figures 5a-c and 6a-c). The isotherm trends from PW2 ( Figure 6) were similar to that of PW1 ( Figure 5), but the formers were more sensitive to the sorption sites-related parameter variations, indicating the coupling of the sorption sites and WWIS.

Relative density of adsorbed water
Relative humidity (%)    Unlike the aforementioned three sorption site parameters, the site spatial distributions did not significantly affect simulated sorption isotherms and hysteresis. The only exception was for the cases of clustered distributions (Figures 5d and 6d) where a stepping stage on the simulated isotherms was observed. That was caused by a two-stage micropore filling that firstly initiated in the pore region where sorption sites clustered and then other regions.

Water Monolayer
Interestingly, a water monolayer (or film) formed on the surface of the simulated pore at an extremely low H of 9.94 × 10 −10 % when the WWIS parameter εsf/kB was increased to an extremely high value of 1165.8K from 58.29K. Figure 7 gives such a snapshot observed in the simulations. Water molecules were well organized on the surface of the simulated pore walls, compared to the random structure from a normal εsf/kB setting. However, in this specific case, the pore wall-water interaction had been strengthened to an extent that was far beyond the one in a real system, and therefore not likely to exist. Moreover, the simulation here also offers some insights to the small amount of water that cannot desorb from wood even when heated at 105 °C for a long period (defined as water of constitution by Stamm [53]). This type of water is supposed to bound tightly to wood cell wall substances, thus forming a strong wood-water interaction scenario, at least on partial regions of the cell walls. Since the water monolayer could exist at extremely low H that is difficult to reproduce experimentally, it is very difficult to remove it in reality. Unlike the aforementioned three sorption site parameters, the site spatial distributions did not significantly affect simulated sorption isotherms and hysteresis. The only exception was for the cases of clustered distributions (Figures 5d and 6d) where a stepping stage on the simulated isotherms was observed. That was caused by a two-stage micropore filling that firstly initiated in the pore region where sorption sites clustered and then other regions.

Water Monolayer
Interestingly, a water monolayer (or film) formed on the surface of the simulated pore at an extremely low H of 9.94 × 10 −10 % when the WWIS parameter ε sf /k B was increased to an extremely high value of 1165.8 K from 58.29 K. Figure 7 gives such a snapshot observed in the simulations. Water molecules were well organized on the surface of the simulated pore walls, compared to the random structure from a normal ε sf /k B setting. However, in this specific case, the pore wall-water interaction had been strengthened to an extent that was far beyond the one in a real system, and therefore not likely to exist. Moreover, the simulation here also offers some insights to the small amount of water that cannot desorb from wood even when heated at 105 • C for a long period (defined as water of constitution by Stamm [53]). This type of water is supposed to bound tightly to wood cell wall substances, thus forming a strong wood-water interaction scenario, at least on partial regions of the cell walls. Since the water monolayer could exist at extremely low H that is difficult to reproduce experimentally, it is very difficult to remove it in reality.

Water Clusters
Hysteresis disappeared when ε H increased to 22.5 kcal/mol (Figures 5b and 6b), and in contrast to abrupt micropore fillings observed in earlier scenarios, a gradual increase of adsorbed water molecules driven by the formation of clusters initiated at the positions of sorption sites was observed. Figure 8 shows snapshots of simulated 0.95 nm pores at a series of given H values: 0.000036%, 0.00050%, 0.16%, 0.35% for PW 1 and 0.000036%, 0.60%, 5.57%, 6.34% for PW 2 . For both pore types, the water clusters were initiated at the sorption sites (indicated in blue balls in Figure 8), and then grew bigger, eventually coalescing as H increased. However, the water clusters only existed at very low H values (<0.35%) for PW 1 and slightly higher H values (6.34%) for PW 2 . Since PW 2 has weaker WWIS, the simulation results demonstrated that the less-hydrophilic cell wall pore surfaces benefited the stabilization of the formed water clusters. The sorption process in this special case could be used to examine the situation described by clustering theory mainly developed by Hartley and Kamke and Hartley and Avramidis [5,6]. Accordingly, one water molecule was attached to each sorption site at a low H region (0-30%), and then as H increased, water clusters with an average size of two formed and grew larger. The predicted maximum size of the water cluster was 98 [6]. Apparently, the H region where clusters existed and their maximum size predicted from the clustering theory was not supported by the simulation results here. The set sorption site energy (22.5 kcal/mol) in the simulations was extremely high compared to the normal value of 4.5 kcal/mol, so water sorption driven by cluster formation was not likely to occur in a real wood-water system.

Interaction of Simulation Parameters
The simulated sorption isotherm trends from different parameter settings were similar, and also resembled those from varied cell wall pore sizes in [25]. In summary, the investigated five simulation parameters-pore size, ε s f , n H , ε H and R H -largely affected the simulated sorption isotherms and hysteresis. Furthermore, different variation sensibility from PW 1 and PW 2 indicated interaction among these parameters. Figure 9 illustrates a typical water sorption isotherm from simulated wood cell wall pores at a given pore size. Two critical points where micropore filling and evaporation occurred, respectively, could be identified (points A and B in Figure 9). The critical Hs at A and B were denoted as H u and H l , and then the magnitude of the simulated hysteresis loop could be calculated as ∆H = H u − H l . Statistical models were employed to analyze the effects of investigated parameters on H u , H l and ∆H, and the potential interactions among these parameters based on the pooled simulation data.
increased, water clusters with an average size of two formed and grew larger. The predicted maximum size of the water cluster was 98 [6]. Apparently, the H region where clusters existed and their maximum size predicted from the clustering theory was not supported by the simulation results here. The set sorption site energy (22.5 kcal/mol) in the simulations was extremely high compared to the normal value of 4.5 kcal/mol, so water sorption driven by cluster formation was not likely to occur in a real wood-water system.

Interaction of Simulation Parameters
The simulated sorption isotherm trends from different parameter settings were similar, and also resembled those from varied cell wall pore sizes in [25]. In summary, the investigated five simulation parameters-pore size, , , and RH-largely affected the simulated sorption isotherms and hysteresis. Furthermore, different variation sensibility from PW1 and PW2 indicated interaction among these parameters. Figure 9 illustrates a typical water sorption isotherm from simulated wood cell wall pores at a given pore size. Two critical points where micropore filling and evaporation occurred, respectively, could be identified (points A and B in Figure 9). The critical Hs at A and B were denoted as Hu and Hl, and then the magnitude of the simulated hysteresis loop could be calculated as ∆ = − . Statistical models were employed to analyze the effects of investigated parameters on Hu, Hl and ∆ , and the potential interactions among these parameters based on the pooled simulation data. For the classical Kelvin's equation (Equation (1) in [25], by replacing the surface tension of water with or the combination of , and , two candidate variables, and , where (nm) is the radius of the pore, were derived for multiple linear regression analysis. Another variable was attempted considering the possible coupling effect from WWIS and sorption sites. Equations (1)-(3) were multiple linear regression models found to fit the simulation data with satisfactory accuracy.
= + + (1) Table 3 gives estimated coefficients ( , i = 0 to 2) in Equations (1)-(3) and corresponding adjusted coefficients of determination R 2 . All the coefficients in Equations (1)-(3) were found statistically significant under the null hypothesis H0: = 0 with a significance level ( ) of less than 0.001. The statistical models (1) and (2) indicated that the cell wall pore size, WWIS, and the sorption sites were coupled in affecting the critical properties of sorption isotherms, but the effect from the cell wall pore size was different in determining the critical points on adsorption and desorption isotherms (Hu and Hl, respectively). For the classical Kelvin's equation (Equation (1) in [25], by replacing the surface tension of water γ with ε s f or the combination of n H , ε H and R H , two candidate variables, ε s f r p and n H ε H R H r p , where r p (nm) is the radius of the pore, were derived for multiple linear regression analysis. Another variable ε s f n H ε H R H r p was attempted considering the possible coupling effect from WWIS and sorption sites. Equations (1)-(3) were multiple linear regression models found to fit the simulation data with satisfactory accuracy. Table 3 gives estimated coefficients (β i , i = 0 to 2) in Equations (1)-(3) and corresponding adjusted coefficients of determination R 2 . All the coefficients in Equations (1)-(3) were found statistically significant under the null hypothesis H 0 : β i = 0 with a significance level (α) of less than 0.001. The statistical models (1) and (2) indicated that the cell wall pore size, WWIS, and the sorption sites were coupled in affecting the critical properties of sorption isotherms, but the effect from the cell wall pore size was different in determining the critical points on adsorption and desorption isotherms (H u and H l , respectively). Model (3) indicated the WWIS and the sorption sites were not coupled in determining hysteresis, but both factors were coupled with cell wall pore size. Apart from illustrating the interaction of r p , ε s f , n H , ε H and R H , the statistical models in this work can be used to predict water sorption isotherms of wood without time-consuming simulations when more accurate input parameters are obtained from experiments.

Conclusions
The shapes of simulated sorption isotherms in this work were different from the experimental ones, which were caused by focusing on only one cell wall pore size in the simulations due to the computational cost. With proper cell wall pore size distributions further considered, the simulated sorption isotherms would get closer to the experimental ones. Nevertheless, the molecular simulations demonstrated that most investigated parameters including wood-water interaction strength and the number, strength, and interaction range of sorption sites affected both adsorption and desorption isotherms and hysteresis considerably. Water monolayers and clusters were observed on the surface of the simulated cell wall pore at an extremely low humidity when wood-water interaction and sorption site strength was set to extremely high values. However, since the simulation parameters were unrealistically high in these specific cases, water monolayers or clusters were not likely to occur in a real wood-water system. Further multiple linear regressions suggested that the cell wall pore size, wood-water interaction strength, and the sorption sites were coupled in affecting the critical properties of sorption isotherms, but the effect from the cell wall pore size was different in determining the critical points on adsorption and desorption isotherms. Wood-water interaction strength and the sorption sites were not coupled in determining hysteresis, but both parameters were coupled with cell wall pore size.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The representative simulation data has been included in the text. Other simulation data will be available upon request.