Humidity influence on mechanics of paper materials: joint numerical and experimental study on fiber and fiber network scale

Paper materials are well-known to be hydrophilic unless chemical and mechanical processing treatments are undertaken. The relative humidity impacts the fiber elasticity, the interfiber joint behavior and the failure mechanism. In this work, we present a comprehensive experimental and computational study on mechanical properties of the fiber and the fiber network under humidity influence. The manually extracted cellulose fiber is exposed to different levels of humidity, and then mechanically characterized using atomic force microscopy, which delivers the humidity dependent longitudinal Young’s modulus. We describe the relation and calibrate the data into an exponential function, and the obtained relationship allows calculation of fiber elastic modulus at any humidity level. Moreover, by using confoncal laser scanning microscopy, the coefficient of hygroscopic expansion of the fibers is determined. We further present a finite element model to simulate the deformation and the failure of the fiber network. The model includes the fiber anisotropy and the hygroscopic expansion using the experimentally determined constants, and further considers interfiber behavior and debonding by using a humidity dependent cohesive zone interface model. Simulations on exemplary fiber network samples are performed to demonstrate the influence of different aspects including relative humidity and fiber-fiber bonding parameters on the mechanical features, such as force-elongation curve, strength and extensibility. Finally, we provide computational insights for interfiber bond damage pattern with respect to different humidity level as further outlook.


Introduction
Cellulose-based fiber materials have been used for decades as packing, printing media. Nowadays, they Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/ s10570-021-04355-y. even have become popular as base material for electronic and microfluidic devices on a small scale (Gong and Sinton (2017), Shen et al. (2019), Schabel and Biesalski (2019), Carstens et al. (2017)), and biocomposite (Pantaloni et al. (2021), Regazzi et al. (2019), Lee and Wang (2006)), due to their recyclability to reduce pollution and save resources. As these devices are often exposed to a humid environment, their mechanical reliability and durability are often constrained and need to be understood before being used for massive application. The general understanding is that mechanical properties of natural compositelike paper-based materials are sensitive to the variation of relative humidity (RH). Salmen and Back (1980) have reported how RH affects the overall stiffness and tensile strength of the paper material. Upon increasing moisture content at higher RH, paper material starts to exhibit a ductile and more elastic behavior, whereas upon drying the material becomes more brittle. On the single fiber scale, moistureinduced elasticity change of fibers from Douglas-Fir latewood tracheids has been firstly investigated by Kersavage (1973) in the early 1970s, then in recent works by Ganser et al. (2014) and Czibula et al. (2021), successfully applying an AFM-based characterization method. In the former work, Kersavage (1973) found the maximum elastic modulus at about 6 % moisture content with decreasing trend towards higher moisture levels. This agrees with the later work, observing a lower elastic modulus towards higher RH levels. However, most of the literature investigating the RH-dependent elastic property is based on lignocellulosic fibers from wood sources. But apart from the source of fibers, the fiber-fiber bond strength plays a crucial role for the load-bearing of paper materials. The interfiber joint strength can be contributed by different sources. While it is believed that hydrogen bounds are the main factor holding paper together, Hirn and Schennach (2015) have shown using AFMbased method that other bonding mechanisms such as interdifussion, Van der Walls forces and mechanical interlocking or capillary bridges coexist. As for the RH influence on the load bearing and failure mechanism of the fiber joint, there are still unclear observations in the literature regarding the dependency of the fiber joint mechanical properties due to RH. Jajcinovic et al. (2018) have shown that the strength of individual fiber bonds increases upon exposure to high RH, where as others have shown that the bond strength decreases.
Based on single fiber bond tests at different RH levels, they have reported that the breaking load of individual softwood fibers and fiber bonds displays a maximum breaking load at 50 %RH, with the values showing a decreasing trend towards higher or lower RH levels. On the other hand, hardwood fibers show rather a decreasing trend in breaking loads with increasing RH levels. Besides the loss of bond strength and decrease of elastic modulus, the swelling or hygroscopic expansion, which describes the moisture uptake and dimensional stability, also plays an important role in RH sensitivity of the paper strength and extensibility. Gamstedt (2016), Joffre et al. (2013), Neagu et al. (2005) and Motamedian and Kulachenko (2019) have reported that moisture adsorption due to humidity is the key feature to the dimensional stability loss of the cellulose-based paper structures. In order to improve the moisture resistance and identifying controlling parameters for engineering application, the underlying mechanisms has been investigated by means of both experimental mechanics and numerical modelling techniques. Thereby moisture is typically characterised either in terms of RH in the surrounding air or the moisture content (relative moisture mass) in the specimen itself. The two measurements are intimately linked, and the relationship is characterised by the dynamic vapour sorption isotherms. While the RH is easily characterised by hygrometers in the ambient air, the determination of moisture content needs to weight the sample and compare it to its dry weight ( Gamstedt (2016)). On the part of computational studies, a recent literature review by Simon (2020) extensively summarizes the numerical modelling approaches in study of material behavior of paper and paperboard on different scales. Fiber network simulation is usually performed by means of direct simulation techniques, where individual fibers are explicitly modelled, and interfiber joint behavior studied by a cohesive zone based modelling approach ( Borodulina et al. (2018), Li et al. (2018)). However, in the mentioned works, key RH influence factors as humidity dependent fiber elasticity and fiber debonding models have not or not sufficiently been included, which remain an open topic for future research in the modelling studies. In the present paper, we aimed to combine experimental and numerical approaches to systematically unveil the RH impact on mechanical properties of paper materials on the fiber, and fiber network scale. Using atomic force microscopy (AFM) and confoncal laser scanning (CLSM) microscopy, we firstly characterized the humidity dependency of the fiber elastic modulus and the hygroscopic expansion coefficient (HEC) of manually extracted cellulose fibers from a finished paper sheet. The determined humidity dependent elastic modulus was then used as input for a transversely isotropic fiber model for the finite element simulations to study the force-elongation behavior. Additionally, by proposing a humidity dependent cohesive zone model of an exponential form, the fiber network simulation addressed the humidity dependent fiber-fiber contact bonding, resolved by the finite element cohesive interface elements. It allowed the simulation of the local damage at individual fiber-fiber joints at different RH levels. The force-elongation curves of the network gave an implication on the overall mechanical behavior. We particularly studied mechanical properties as the tensile strength, the effective stiffness and the extensibility. The micro-mechanical model was then qualitatively compared with tensile test at paper sheet level. Interfiber debonding pattern due to RH in combination with RH-dependent fiber elasticity and cohesive strength reduction were computationally provided as further outlook.

Fiber sample and bending test
For the bending experiments, the cellulose fibers were manually extracted from a cotton Linters paper sheet (median fiber length (length-weighted): 1.13 mm; fiber width: 16.5 lm; curl: 12.87 %; fibrillation degree: 1.64 %; fines content: 19.55 %; Crystallinity: 56.31 %). The sheet was prepared according to DIN 54358 and ISO 5269/2 (Rapid-Köthen process). The extracted fiber was glued (Glue Essie Polish, Düsseldorf, Germany) on a 3D printed sample holder, which supplied a fixed trench distance of 1 mm between the two attachment points. Worth mentioning that the used glue had not been soaked into the fibers as demonstrated in our previous work (Auernhammer et al. (2021)). AFM was then performed by using a Dimension ICON (Bruker, Santa Barbara, USA) to measure static force-distance curves along the cellulose fiber with a colloidal probe. The cantilever (RTESPA 525, Bruker, Santa Barbara, USA) with the nominal spring constant 200 N/m was modified with a 50 lm colloidal probe (Glass-beads,Kisker Biotech GmbH & Co. KG, Steinfurt, Germany) in diameters. All experiments were done in a climate chamber at the AFM. Hence, it was possible to vary the RH during the experiments. The chosen RH were 2 %, 40 %, 75 % and 90 %. As the RH was adjusted, the fiber was exposed 45 minutes to the environment before starting the measurements. According to Jajcinovic et al. (2018), reaching a full adsorption equilibrium, the fiber has be exposed at least for 12 hours to the certain level RH. Therefore, the measurements were carried while the fibers were not yet fully equilibrated. But, to guarantee the consistency between CLSM and AFM experiments and to avoid drift problems in the AFM, we defined 45 minutes of ''exposure time'' as a feasible compromise to balance effects from instrumental drift and slow swelling processes. For each fiber at each RH level, four different load levels with F = 500 nN, 1000 nN, 1500 nN, 2000 nN were applied and deflections of nine different segment points were evaluated.
Determination of humidity dependent Young's modulus One extracted fiber sample is shown in Fig. 1a. The fiber was glued at both ends. Apart from standard bending experiments, we performed a scanning-like bending test, as we applied the load along the fiber length direction as illustrated in the sub-figure (e). This aims to cover local inhomogeneity induced variations and give a proper statistical evaluation later.
The fiber was modelled mechanically by using the beam model illustrated in Fig. 1b. Based on the images, fiber cross sections are close to hollow ellipses shown in Fig. 1c. By use of the Euler beam theory, the differential equation governing the deflection d of the fiber and the corresponding solutions are given as: where d denotes the deflection, E the longitudinal Young's modulus, F the loading force, M the internal reaction moment, the area moment of the assumed cross-section with c i ; d i and c o ; d o the main axes of the inner and of the outer ellipse, respectively. Equations (2) and (3) can be obtained by calculating the internal reaction moment M and integrating Eq. 1. Then, using the boundary values dð0Þ ¼ dða þ bÞ ¼ 0 and d 0 ð0Þ ¼ d 0 ða þ bÞ ¼ 0 for the integration constants. Detailed procedure about solving statically indeterminate beam systems can be found in Gross et al. (2009). Using these equations, the average longitudinal Young's modulus E is obtained by minimizing the differences between the measured and the calculated deflection from tapping different segments of the fiber via the least square approach: where the positive constant tol is the tolerance. One obtains, in fact, first the solution for E I, and then determines the average E along the fiber sections by further dividing the value I ¼ meanðIÞ, which denotes the averaged area moment along the the fiber segments. The cross-sectional dimensions along the fiber to calculate I were determined using CLSM. This procedure was performed for every tested fiber at given load and RH level. The results are presented in the later section.

Determination of hygroscopic expansion coefficient
A VK-8710 (Keyence, Osaka, Japan) confocal laser scanning microscope was used to investigate the swelling behavior of the fiber. In the climate chamber, where the RH was varied as in the AFM measurements previously, the fiber was suspended 45 minutes to the RH before starting the measurements. The swelling behavior was analysed by the VK analyser software from Keyence (Osaka, Japan), after creating crosssection images along the fiber, where the local radii were measured. Further, to determine the HEC, the change in volume was calculated based on the experimentally recorded change of the cross-section at different RH. The reference volume of the nonprismatic fiber with assumed geometry is Fig. 1c. The total length of the fiber is assumed to be constant and not vary with the changing RH, since both ends are clamped. The where is the sum of 3D anisotropic hygroscopic expansion tensor b kl with in total 9 independent components in the general anisotropic case. Hereby d kl is the Kronecker-delta, b T is the HEC in the transverse direction, and b L is the HEC along the fiber length direction. In the transversely isotropic case, the hygroscopic expansion tensor has only the orthogonal components, namely those in the fiber longitudinal direction and in its cross-section. Further, it is experimentally validated considering the hierarchical layered wall structure ( Joffre et al. (2016)), that the hygroscopic expansion in the fiber length direction b L is shown to be one order lower than that in the transverse directions, and therefore to be neglectable. Under this assumption, it leads to the following equation: Subsequently, the HEC in the transverse direction b T is calibrated with Eq.7 at different DRH for every segment along the fiber length, as explained in Sect. 4.2 later on.

Paper sheet sample and tensile test
For studying the tensile behavior of the paper materials, lab-engineered paper samples with bleached eucalyptus sulfate pulp (median fiber length (lengthweighted): 0.80 mm; curl: 16.7 %; fibrillation degree: 5.2 %; fines content: 8.0 %) were used. The paper samples with grammages of 30 AE 0:6gm À2 were prepared using a Rapid-Köthen sheet former according to DIN 54358 and ISO 5269/2. The tensile tests of the paper samples were done in accordance with DIN ISO 1924-2 with a Zwick Z1.0 with a 20 N load cell, and the speed of separation was set to 10 mm/min, using the software testXpert II V3.71 (ZwickRoell GmbH & Co.Kg) in a controlled environment with 23 C and 50 %RH. Prior to testing at 50 %RH, the paper samples were stored for at least 24 hours in the controlled environment with 23 C and 50 %RH. To perform tensile tests of paper samples at 90 %RH, the samples were stored in a plastic container with a saturated aqueous KNO3-solution (Wexler and Hasegawa (1972)), producing an atmosphere with approximately 90 %RH at 23 C indicated by a sensor inside the container, for 24 hours. After equilibrating the samples in the plastic container, each sample was separately taken out and tested immediately in the apparatus.
Computational simulations using the cohesive zone finite element model

Mechanical model of the single fiber
The mechanical governing equations including the stress equilibrium, the linear kinematics and the linear elastic material law are given as following: in which r ij is the Cauchy stress tensor, C ijkl the stiffness tensor and e kl the strain tensor. The strain due to the hygroscopic expansion is given as where b kl and DRH denote the anisotropic HEC tensor as described in the previous section and the relative humidity change, respectively. The transversely isotropic constitutive material law is applied for the fiber anisotropy along the fiber direction. Thus, the inverse of the stiffness tensor can be given as following: with five independent parameters: E L , E T , G T the longitudinal modulus, the transverse modulus, the shear modulus and m LT , m TT the two Poisson's ratios. The values of the Poisson's ratios are specified in Table 1. Based on the work by Magnusson and Ö stlund (2013), one can further reduce the number of elastic parameters, by assuming the correlations between other elastic parameters and the longitudinal Young's modulus E L with the unified property-related S2-Layer. This assumption is based on the fact that this layer represents the main constituent of the fiber. The relations between the elastic parameters shown in Table 1 (2013), a cohesive zone-based approach was utilized to characterize the debonding behavior. See Figure 3d for the illustration. In the current work, we applied a non-potential based CZM from McGarry et al. (2014), which aims to give a proper behavior in mixed-mode loading scenario to avoid the fiber penetration. The traction-separation laws are given as: where T n ; T t are the traction components of T in their normal and tangential loading state. D n ; D t are the displacement separations at the fiber-fiber bonds. r max and s max are the maximal stresses in pure normal and pure shear separations, d c n and d c t are the critical lengths to the maximum stresses, respectively. The energy terms per surface for both separation modes are given with ranging from 0 to 1 and represents the intact state and are the effective mixed-mode separation and the separation corresponding to the maximum stresses, and D f ¼ is the final separation at complete failure. Eq. 12 and 13 are referred as the traction-separation laws in the dry state T dry for the following humidity dependent cohesive zone model.

Humidity dependency of the cohesive zone model
In a subsequent step, the cohesive zone model due to the humidity influence was modified by multiplying a decrease term: Table 1 The elastic parameters and their dependency on the longitudinal Young's modulus E L based on Magnusson and Ö stlund (2013) Elastic parameters This was motivated by Jemblie et al. (2017), mainly based on using CZM for hydrogen embrittlement of steel structures, where the effect of hydrogen on the accelerated material damage were studied. The analogy in the existing context with hydrogen embrittlement is phenomenologically similar, as the water molecules infiltrate the network and destroy the bonds and as the moisture gets absorbed, the swelling and softening of the fibers lead to partial failure of the fiber bonds ( Hirn and Schennach (2017)). Whereas in case of steel, diffusible hydrogen dissolves into the metal and causes material to be brittle and thus reduces the yield stress in a tensile field at the crack front ( Barnoush (2007)). In the formulation, K is a softening parameter and describes the degree of decrease and the bound of the allowable cohesive strength, which need to be adjusted, when experimental data are available. This intuitive assumption was made for the first step, since there was barely any experimental data in the literature regarding RH dependency of mechanical properties of interfiber joints. Only Jajcinovic et al. Weak formulation for cohesive FE implementation The FE implementation was based on the governing equations describing the general bulk deformation and the cohesive zone damage model as given previously. The corresponding weak formulation can be obtained from the principle of virtual work based on the work of Park and Paulino (2012): where de; du; dD are the virtual strain tensor, the virtual displacement vector and the virtual separation vector at the interface, respectively. r is the tensor notation of the Cauchy stress tensor r ij , T is the traction vector at interface C int with the normal and tangential components T n ; T i . T ext is the external traction on the outer boundary C. Equation 16 describes the energy balance between the strain energy or the work of the internal force in the domain X and the work contributed by the interface elements and by external forces. The non-potential based cohesive zone model is implemented in a user material kernel in the open-source FE software MOOSE ( Permann et al. (2020)). A detailed step-by-step instruction of general cohesive FE implementation from Eq. 16 can be found in Park and Paulino (2012). The Software Geodict r was used to create periodic, synthetic fiber network samples, which allows a similar deposition procedure as reported in Kulachenko and Uesaka (2012), Li et al. (2018). It uses voxels to represent the fiber network after the structure is generated for the given geometric parameters. Upon initialization, a voxel size of one micron was chosen in our work, see Fig. 3b. The fibers were deposited by assuming Gaussian distributions of fiber orientation, diameter and length. For simplicity, we neglected the complex shape of the cross-section and assumed that the fiber was solid and had a circular-shaped cross section along the fiber length. They were created in a flat-plane, afterwards fell under gravitational force, and were then elastically deformed and deposited onto the previous fibers. The deposition process was completed when the specified grammage was reached. After generation of the voxel geometry, the surface mesh was created. Afterwards, coarsening and smoothing steps were performed until sufficiently fine surface mesh was obtained. The surface mesh was then exported and meshed to volume elements using the open-source meshing software Gmsh ( Geuzaine and Remacle (2009)). The modeled fiber network sample included approximately 1 million linear tetrahedron volume elements. About 40000 local 2D-elements were used to resolve the contact area in order to achieve mesh independence. Details on the choice of the mesh size are given in the supplementary information. For the present work, the fiber network setting and features are given in Table 2. The size of the simulated fiber network sample was similar to that used by Li et al. (2018) and our previous work (Lin et al. 2020). The diameter of the fiber was similar to our tested fiber. Due to the periodic fiber network structure, the length of the fibers in the simulation was  chosen to be around half of the median fiber length (length-weighted) from the tested paper sheet as mentioned previously. For the finite element simulations, we calculated the orientation tensor for each fiber in the voxel-based geometry by use of Euler angles. The orientation tensor was then used for transformation of the elasticity tensor for each fiber.
Regarding the boundary conditions, u x ¼ 0 was specified on the left boundary, u x ¼ u on the right boundary and u y ¼ 0 on the front and back boundaries. This setting corresponds to displacement boundary condition in the x-direction with the boundary crosssection in y-direction kept as a flat plane. The out-ofplane degrees of freedom and others on each boundary side were unspecified and stress-free. FEM simulations were performed on a high performance computer with 24 cores for around 8 hours.

Results and discussion
Parameterization of humidity dependent elastic modulus Prior to executing of the minimization procedure, the cross-sectional areas of the fibers were determined by the CLSM at indentation points to calculate the second area moment I in the deflection direction. The values were recorded at each RH level, thus also considering the morphological change of the fiber for a precise determination of humidity-dependent elastic modulus. In Figure 4, three representative cross-sectional images of a fiber at different RH levels are demonstrated. The fiber cross-sections vary from elliptical to circular shapes, typically with a lumen of the same shape for non-collapsed fibers. Therefore, assuming the cross-sections to be elliptical to calculate the I with the corresponding inner and outer diameters also included the extreme case of circular cross-section and helped to improve the accurate determination of I. The numbers of determined areas and second area moments are provided in Table 3 The experimentally determined longitudinal Young's moduli E L using aforementioned technique for the prescribed loads at different RH values are presented in form of boxplot in Fig. 5a. For simplicity, E L will be written as E in the following. For each box, Q1 to Q3 quartile values of the data, (corresponding to 25% and 75%) are drawn inside the box. The yellow line inside the box marks the median (Q2). The whiskers extend from the edges of box to show the range of the data. The extension of the data are limited to 1.5 * IQR with IQR = Q3 -Q1 from the edges of the box, ending at the farthest data point within that interval. Outliers are plotted as separate black framed dots. From the boxplot, it is clear that E decreases with  2011)). Further, it is noticeable that for RH = 2 %, a higher variation is observed due to distribution of determined data points. Apart from one outlier, the median is obtained to be approximately 4 GPa. The variation decreases with increasing RH level, which is given by the shrinkage of the boxes. Compared to the work by Hearle and Morton (2008) on inital elastic modulus of different cotton fibers (5-12 GPa), our values seem to be lower. The lower level of our Young' moduli can be attributed to the cotton fibers under testing. As stated previously, we manually extracted single fibers from a finished paper sheet. Thus, our fibers were a part of the completed paper after production process with beating or pressing. In other words, the fiber tested is not anymore in its natural raw state. In fact, in the paper making process, the layered wall structure of the fiber is often milled off, which could lead to softening of the fiber. Furthermore, a small degree of determined Young's moduli seems to increase with the applied load. This may stem from nonlinear elasticity or nonlinear geometrical deformation which is neglected in the current mechanical models. Compared the elastic moduli with other wood-based lignocellulosic fiber sources as provided by Lorbach et al. (2014), within a range from 15 -35 GPa, natural cotton fibers are generally softer. Further, we proposed the Young's modulus E as a function of RH and approximate the decaying due to the humidity effect in an exponential form. Afterwards, the coefficients of the exponential function were calibrated to the mean values that were determined by the least square approach from the experimental bending tests. The exponential function takes the following form: where a can be understood as E max , which represents the averaged modulus extrapolated at 0 %RH and is (a) (b) Fig. 5 a Distribution of Young's modulus data in dependency of RH. In total, 20 data points were evaluated for 5 fibers and 4 applied loads, at given RH levels. b Calibrated exponential law with corresponding mean values as star dots and its 95 % confidence interval around 4.5 GPa. b determines the degree of elastic decaying. The calibrated coefficients a, b are given in Table 4. The E-RH curve is depicted in Fig. 5b, along with the deviation as 95 % confidence interval. With the formula given in Eq. 17, Young's modulus can be determined at any RH values and can be used further for simulation of RH-dependent material behavior. In Dale et al. (1949), Hearle and Sparrow (1979), a decreasing trend in inital elastic modulus of the cotton fiber w.r.t moisture was also observed. However, an exponential decay in elastic modulus of cotton fibers were less reported. Only recent work by Czibula et al. (2021) using an AFM-based nano indentation method determined similar exponential decrease of elastic modulus for softwood fibers w.r.t RH.

Hygroscpic expansion coefficient
The experimentally determined volume changes of different segments of a fiber exposed to different RH are summarized in Fig. 6a. Together, data of five fibers, each of nine segments were collected and 45 data points were evaluated in total. Due to the high local inhomogeneity, high deviation of the data points for different segments was observed, which had larger distributions in the boxplot at each RH change. DRH is referred hereby to the difference between the current RH level and reference RH state, which is 2 %RH. Similar to the determination of E-RH relation, the HEC was determined in Fig. 6b. The volumetric strains were averaged over the data points at each DRH level. The average volume strain assumes a linear dependency on the relative humidity, based on the relation Eq. (7). The transverse HEC b T was thus extracted as the slope of the curve, which was indicated by the dashed black line, and the value was determined as 0.35. This is in good agreement with the values reported by Joffre et al. (2013) as collected in Joffre et al. (2016) from different data available in the literature, ranging from 0.2 to 0.44. It should be noted that the recorded outliers in the boxplot originated from the segments that contained inflated microfibrils, as e.g. segment number 4 in Fig. 1a, which would bias the overall results. There were several other methods in determining the HEC, such as measuring from X-lCT or back-calculation from composite properties ( Neagu and Gamstedt (2007), Almgren et al. (2009)). Since our focus was on the influence of RH on the mechanical behavior of the fiber network, the readers are referred to literature for further investigations about the dimensional stability loss in relation with HEC, see e.g. Neagu et al. (2005).

Simulation of force-elongation relation on fiber network scale
Bulk and cohesive finite element simulations were carried out on the generated fiber network under the prescribed boundary conditions in Sect. 3.5. One fiber network sample without humidity influence is shown for demonstration purpose in Fig. 7(a), with the  Fig. 6 a Distribution of measured volumetric strain data. In total 45 data points were evaluated, with data from 9 segments per fiber at each RH change. b Calibrated linear relation of volumetric strain with DRH and the corresponding mean values as star dots, and its 95% confidence band. The slope of the linear curve was then determined as the transverse HEC cohesive fiber-fiber interface highlighted in dark blue. The sample deforms under the tensile displacement along the x-direction. A few snapshots on the deformation and the interface damage level are sequentially presented in Fig. 7b-f. We simultaneously documented the normalized reaction force (w.r.t maximal value of the reaction force) and the elongation of the sample. The corresponding curve is depicted as the insert of the figure.
The points which correspond to the snapshots (a-f) are indicated on the curve. It can be seen that mechanical response from (a) to (b) is almost linear. From (b) to (c), damage initiation occurs in certain areas. The maximal force reaches at (c), and from (c) to (d), the force drops about 35%. The force at point (c) is referred as the maximal force or the failure force, while the corresponding elongation is termed as elongation at failure in the following context. From (c) on, damage process continues more extensively. In particular, from (e) and (f) almost all fiber contacts are damaged, and the sample bears hardly any load (20% of the maximal force) due to the damage percolation. The numerical convergence at this stage becomes difficult for implicit finite element solver due to the large deformation and the energy release stored in the fiber bonds. Using the established fiber network simulation framework, we firstly carried out a series of material parameter studies without humidity influence. It is necessary to understand the extent of the influence of those potentially important material parameters on the mechanical behavior, because the material parameters inherently vary with different kinds of fibers. Simulations with humidity influence are reported thereafter. It should be noticed that for the following parameter studies, force-elongation curve up to force at failure and the corresponding elongation at failure was considered, as indicated by point (c) in Fig. 7. Force-elongation behaviors in post-damage range are out of our current scope, where the numerical convergence may become difficult for implicit finite element solver as previously mentioned. Table 5 summarizes the default values and the corresponding variation sets of the material parameters considered. Unless it is otherwise stated, only one set of the material parameters is varied at once, while the other parameter sets are assumed to take their default values.

Influence of cohesive model parameters
The cohesive parameters describe the mechanical properties of the bond and the energy stored in the fiber-fiber contacts. The choice of cohesive parameters was based on fiber joint tests from available literature. The reported experimental values in a fiber-fiber shearing type study range from 0.2 to 16 MPa according to Magnusson et al. (2013) and about 3-5 MPa for softwood and hardwood joints, according to Jajcinovic et al. (2016), respectively. The normal cohesive strength has been reported to be approximately four times smaller than that in tangential direction as given in Magnusson and Ö stlund (2013), Borodulina et al. (2018), Marais et al. (2014). Worth mentioning, the determined strength values in peeling or shearing type studies are to some extent, always mixed in the loading and due to experimental setups never in its pure form ( Magnusson (2016)). Less information is available on the critical separations of fiber joints in experimental context, due to difficult experimental realizations as mentioned. As for the numerical simulation using cohesive zone models, work by Magnusson (2016) set the cohesive final separations to be 1 lm for both normal and tangential direction, and in Borodulina et al. (2018), Mansour et al. (2019), they are assumed to be 1.56 lm and 0.35 lm, respectively. Given the variation in the reported values in the literature, a range of parameters were therefore checked to study their influence on the mechanical response of the fiber network. Based on the literature, we set the default cohesive material parameter for s max ¼ 1 MPa, r max ¼ 0:25 MPa and followed Magnusson (2016) setting d f t ¼ 1lm and d f n ¼ 1lm, with corresponding critical separations determined by equating the energy terms / n ¼ d c n r max exp 1 ð Þ and / t ¼ d c t s max ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0:5 exp 1 ð Þ p from our work with bilinear cohesive energy / n ¼ 1 2 r max d f n and (2016), Borodulina et al. (2018). The variation range is given in Table 5. Generally, non-linearity in force-elongation behavior is observed. From the simulation results, after a small elastic region, the curves become nonlinear due to the inhomogeneous damage development at different fiber bonds and the progressive and cooperative interactions between damaged interfaces. It is noted that force drop occurs during the loading, which can be reasoned by reorientation and readjustments of the fibers under loading. Figure 8a demonstrates the influence of the parameters r max ; s max on the mechanical behavior, where the blue curve with corresponding (*) in the legend of the figure denotes the results with default parameters. As it is shown, reducing s max by 10 times from 1 to 0.1 MPa (respectively for r max ), the elongation reduces by a factor of around 8 (from 0.58 to 0.07 %), the maximal force reduces about 8.3 times as well (from 2 mN to 0.24 mN). Increasing s max by a factor of 10 (respectively for r max ), the elongation increases by a factor of around 2.4 (from 0.58 to 1.41 %), the reaction force increases around 4 times (from 2 to 7.5 mN). Further increase of s max from 10 to 20 MPa shows an increase of 1.4 times in elongation and of 1.5 times in maximal force. Clearly, the force and elongation have non-linear responses to the change of the cohesive strength, the factors w.r.t changing cohesive strength are rather same for both force and elongation quantities. For the effective stiffness of the sample, the secant modulus is considered (ratio of the maximal force per unit area to elongation), since not every simulated case has a clear elastic region. The effective stiffness is less sensitive against the change of the cohesive strength, which implies in turn, that the overall stiffness of the fiber network is less influenced by the cohesive strength. Figure 8b demonstrates the influence of the cohesive separation d f n ; d f t on the force-elongation curves. When reducing d f t ; d f n from 1 to 0.1 lm, maximal force decreases by a factor of 2.2 (from 2 to 0.9 mN) and elongation decreases by a factor of 3 (from 0.57 to 0.19 %). Increasing from default parameter 1 lm, to 5 lm and 10 lm, the failure force changes from 2 mN to 3 mN and 3.4 mN, which corresponds to an increase factor of 1.5, 1.13, respectively. For the elongation, factors of 3, 2.5, 1.8 are obtained for values increasing from 0.19 %, 0.57 %, 1.41 % and 2.55 %. The change in effective stiffness (as secant modulus) shows a decreasing trend with factors of 0.74, 0.6 and 0.63. Therefore, increasing the cohesive separation, the fiber network is considered to be softer, stronger and has a higher deformability, whereas changing cohesive strength does not influence the stiffness of the fiber network, but makes it stronger and more extensible as well.

Influence of Young's modulus
As discussed previously, our determined fiber Young's modulus appears to be lower than the reported values in the literature. We study hence here the influence of the variation of the elastic modulus on the force-elongation behavior, given our experimentally determined results and those reported in the literature. Default elastic modulus is set to be 4.5 GPa, extrapolated from our calibrated exponential law at RH = 0 %. Figure 8c shows that increasing the Young's modulus by a factor of 10, the maximal force increases from 2 to 3.4 mN, which corresponds to a factor of 1.7. The elongation decreases from 0.57 to 0.25 % and gives a factor of 2.3. Further increasing E from 45 to 90 GPa, only minor changes are observed in the mechanical responses. More interestingly, reducing 4.5 GPa to 0.45 GPa, the elongation increases by a factor of 3, and the maximal force decreases approx. 2 times. The effective stiffness changes significantly, with a factor of 1.09, 2.3 and 3 when changing 90 GPa to 45 GPa and to 4.5 GPa, finally to 0.45 GPa. The mechanical responses scale exponentially, and implies a strong impact on the extensiblity and the effective stiffness of the fiber network by the stiffness of constituent fibers. It is also notable, that the mechanical behavior changes from very brittle to a more ductile behavior, with decreasing Young's modulus.

Influence of relative humidity
In this section, the RH influence on the two key parameters in our simulations are investigated, namely the fiber elastic modulus and the cohesive strength parameter. For a better analysis of the parameter influence and individual parameter contributions, we firstly studied the force-elongation curves by RHprovoked elasticity decay, then by the RH-provoked bonding strength reduction, and finally combined influence of both factors. For a simple comparison, the force values were normalised with the maximal number. In the Figure 9a, for a reference dry state of RH = 0%, force-elongation curves are shown for RH levels of 2%, 40%, 75% and 90%. Accordingly, using calibrated exponential law of fiber elasticity decay from previous section, see Eq. 17 and Table 4, the (a) (b) (c) Fig. 9 Force-elongation behavior of a RH Influence on the fiber elasticity. b RH influence on the cohesive strength. c Combined influence of cohesive strength and fiber elasticity on the force-elongation curve. (*) denotes the default RH level resulted elastic moduli due to RH influence to determine the force-elongation curves were around 4000, 400, 50 and 20 MPa, respectively to the RH levels. As depicted in the figure, the stiffness and the maximal force of the fiber network gradually decrease, whereas the elongation or extensibility increases with RH levels. The trend of the curve qualitatively follows the exponential drop at single fiber level with increasing RH levels and is also comparable with the previous elastic modulus parameter study. Further, for the RH-provoked cohesive strength reduction, we assumed the decaying parameter K used in the humidity dependent cohesive zone model in Eq. 15 to be 0.3, which implied that for e.g. a DRH = 30%, the cohesive strength is decreased by approx. 10%. This amount of decrease is similar to the interfiber joint test for hardwood fibers subjected to varying RH observed in Jajcinovic et al. (2018). The set of force and elongation curves in Fig. 9b, indicates that while the stiffness remains stable in the initial range, the force and the elongation decrease rather linearly within the given parameter range by the RH impact. Further looking into the combined factors in Fig. 9c, the set of curves subjected to both RH dependent factors is similar to the set of curves of RH-provoked elastic modulus decay in Fig. 9a, which shows an exponentially decreasing trend in maximal force and stiffness, but an improvement in elongation towards higher RH levels. Coupled with the effect of cohesive strength reduction, stiffness, elongation and strength change slightly from (0.75, 0.36) to (0.68, 0.31) at 40 %RH ,(1.71,0.17) to (1.70,0.11) at 75 %RH and (2.85,0.08) to (2.66, 0.07) at 90 %RH, respectively. As seen earlier, change of cohesive strength influences the elongation and maximal force, but to a smaller extent the effective stiffness. For a cohesive strength decay parameter K = 0.3 (i.e. maximal cohesive strength decrease of 30%), comparatively minor effect on the mechanical behavior is observable. It is therefore understandable from the computation results, that based on the experimentally obtained E-RH behavior of exponential form from single fiber level, RH dependent elastic modulus makes a greater impact on the RH-dependent mechanical behavior of the underlying fiber network structure within the current setting.

Tensile test of paper sheet and comparison with micromechanical model
To provide additional information of RH influence on the mechanical behavior and qualitatively validate the model, tensile test on paper sheet samples were performed at RH = 50 % and 90 %, with the mentioned sample size in previous section. Results are shown in Fig. 10a. For RH = 50%, data from nine samples was collected and the maximal force was determined to be 7:56 AE 0:31N, with the corresponding elongation at break to be 1.45 ± 0.12 %. For RH = 90 %, data from ten samples was collected. The maximal force was 6:6 AE 0:45 N, with the corresponding elongation 2:18 AE 0:15%. As the (a) (b) Fig. 10 Experimental force-elongation curve from tensile test of paper sheet samples at different RH levels figure shows clearly, the maximal force and the effective stiffness decrease and the elongation increases with increasing RH levels. These observations have been reported in the available experimental literature. To further validate the micro-mechanical model with the experimental results, the numbers from experimental results were normalised, errorbars were determined and compared with simulated results at RH = 50 % and 90 % in Fig. 10b. As shown in the figure, the simulated results agree well with experimental results at 50 %RH level for the force and elongation, respectively. The stiffness in the initial region is slightly higher, but agrees towards the softening region. However, larger discrepancies are observed for the curve at 90 %RH level. Whereas the simulated elongation is slightly higher, the force and stiffness numbers are lower than the experimental results. The main reason assumed for this discrepancy lies in the obtained E-RH relation from tested single fibers. The procedure of extracting the single fibers from a finished paper sheet can damage the fiber wall, as mentioned in previous Sect. 4.1, thus leading to an higher elastic modulus decrease when exposed to increasing humidity (from approx. 4 GPa at 2 %RH to 20 MPa at 90 %RH), for which the decay relationship determined in our calculation to be an exponential form. For fibers in the paper sheet under actual testing, this decrease in fiber elasticity might be lower, while the fiber wall remains intact and the fibers well embedded in the fibernetwork. Thus, the curve tends to have a higher force and stiffness at sheet level from the experimental testing. In general, the fibernetwork simulation results presented here can be rather interpreted in terms of a statistical volume element as in Trias et al. (2006), Lin et al. (2020) framework, limitation in computation power arises due to the number of fibers and interfiber contacts, that needs to be resolved sufficiently to obtain meshindependent results. Therefore, further modelling techniques are required, e.g. homogenisation of the presented fiber network models as an element or starting point of a multi-scale modelling paradigm. For that purpose, our bottom-up approach considering RH influence delivers important information for understanding and designing of paper materials on the fiber, and fiber network scale and can serve as baseline for studying mechanical and fracture behavior of paper material models on paper sheet, or even higher scale.

Outlook on computation of RH dependent interfiber damage pattern
In this section, computational insights on the fiber network deformation and interfiber debonding pattern with respect to RH levels are provided. As demonstrated earlier, increasing RH levels causes the softening of fibers in terms of fiber elasticity decay and reduction of cohesive strength. These combined influences have large impact on the interfiber damage behavior, its accumulation of fiber deformation and stress state in a fiber network. In Figure 11, the von Mises stress and the damage pattern of the underlying fiber bonds are presented. As shown by the von Mises stress map, at RH = 2 %, the numbers at specific location of fibers are several times higher, due to the high stiffness, low deformability of the fibers. Therefore, more bonds are damaged to that end. Compared to higher RH level e.g. at 75 %, the fibers are more deformable and upon stretching, the fibers have high tendency to deform without having every bond necessarily broken, thus the interfiber bonds at higher RH are less damaged. Worth mentioning, that we do not consider any fiber breakage in our model, nevertheless, in a realistic case, that the higher stress concentration in fibers at lower RH or drier state can result in multiple fiber breakage instead of fiber bond failure, and then through the failure, fiber bonds remain less damaged. However, to conclude whether fiber breakage becomes more dominant in drier state, experimental tests and statistical evaluation of failure mechanism on a fibernetwork level, with respect to different RH level need to be performed, which remain as an interesting aspect for further work.

Conclusion
In summary, on different scales, starting from the fiber scale: We experimentally determined the humidity dependent longitudinal fiber elastic modulus using atomic force microscopy and proposed an exponential law for decaying of the humidity dependent fiber elasticity. The obtained Young's modulus was about 4 GPa at 2 %RH and lower than values available in the literature. It could be attributed to the fact that the tested fiber was rather a fiber manually extracted from a finished paper sheets than a raw one without the influence of a paper-making process. Hygroscopic expansion coefficient was obtained to be 0.35 using confocal laser scanning microscopy, by recording the change of fiber cross-section swelling at different relative humidity levels. On the fiber/fiber cross scale: humidity dependent cohesive zone model was proposed and applied jointly with the humidity dependent elasticity to study the force-elongation behavior of fiber network under tensile loading on a higher scale. On the fiber network scale: we carried out a series of parameter studies, and determined that key contribution to RH influence on the force-elongation behavior is due to the RH-dependent exponential decaying from single fiber elasticity, whereas the interfiber properties play a minor role within the reported material parameter range. The force-elongation behavior of simulated tensile test under humidity influence was in qualitative agreement compared to our experimental tensile tests of paper sheets. Further, to quantitatively address the RH influence on the paper sheet scale, a multi-scale, microstructure-informed modelling approach with demonstrated fibernetwork simulation would be helpful, which will remain as future perspective on the computational modelling part. Computational simulations of damage pattern at different RH levels were demonstrated, however, to draw further conclusion on the failure mechanism, experimental test and statistical evaluation on a fibernetwork scale, with respect to different RH level are needed. Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The authors declare that they have no conflicts of interest.
Ethical approval This work does not contain any studies with human participants or animals performed by any of the authors.
Informed consent No human participants were involved in this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.