{'fraction_non_alphanumeric': 0.029542173644050318, 'fraction_numerical': 0.010955867189111157, 'mean_word_length': 4.751371386212009, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Humidity Influence on Mechanics and Failure of Paper Materials: Joint Numerical and Experimental Study on Fiber and Fiber Network Scale

Paper materials are natural composite materials and well-known to be hydrophilic unless chemical and mechanical processing treatments are undertaken. The relative humidity impacts the fiber elasticity, the fiber-fiber bonds and the failure mechanism. In this work, we present a comprehensive experimental and computational study on the mechanical and failure behaviour 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. 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. On the other hand, we 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. In addition, it regards the fiber-fiber bonding and damage by using a humidity dependent cohesive zone interface model. Finite element 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 curves, wet strength, extensiability and the local fiber-fiber debonding. In meantime, fiber network failure in a locally wetted region is revealed by tracking of individually stained fibers using in-situ imaging techniques. Both the experimental data and the cohesive finite element simulations demonstrate the pull-out of fibers and imply the significant role of the fiber-fiber debonding in the failure process of the wet paper.


Introduction
Paper materials have been used since decades as packing, printing media, and now they even have become interesting for electronic, microfluidic and medical devices on small scale due to its recyclability to reduce pollution and save resources [1,2,3,4,5,6]. As these devices are often exposed to a humid environment, their mechanical reliability and durability are often constraint and need to be understood before used for massive application.
The general understanding is that the mechanical property of natural composite-like paper-based materials is sensitive to the variation of relative humidity (RH). Salmen et al. [7] reported how RH affects the overall stiffness and tensile strength of the paper sheet. Upon increasing moisture content at higher RH, paper starts to exhibit more ductile and elastic behaviour, whereas upon drying the material becomes more brittle. However, the RH induced variation of the elastic property of single fibers and the influence of this variation on the overall paper sheet is still insufficiently characterized, especially on the scale of the fiber and fiber network.
Their load bearing and failure mechanism varies also with moisture state. There are still controversial statements in the literature regarding the dependency of the fiber mechanical property on RH. Jajcinovic et al. [8] has shown that the strength of individual fibers increases upon exposure to high RH, where as others have shown that the strength decreases. Based on single fiber testing at different RH, they reported that the breaking load of individual softwood fibers and fiber contacts displayed a maximum breaking load at 50% RH, with the values showing a decreasing trend towards higher or lower RH. On the other hand, hardwoods show no change in the breaking load observed at different RH. More importantly, the fiber-fiber bond strength is very sensitive to humidity. Dry paper strength is determined mostly by both the strength of their fibers and the bond strength between fibers. When paper breaks, the failure of both fibers and bonds has been observed. In the case of wet paper with sufficiently high RH, the overall strength decreases. When it ruptures, no fibers break but mostly fiber bonds [9]. The fiber-fiber bonding strength mainly results from the hydrogen bonds and depend on the fiber/fiber cross contacts, which can be broken by wetting of the fibers [10]. The force holding the fiber-fiber bonding is very sensitive to water, and the extent of bonding decreases with increasing water content [11,12].
Besides the loss of bond strength and decrease of elastic modulus, the swelling or hygroscopic expansion, which describes the moisture uptake, also plays an important role in RH sensitivity of the paper strength and extensiability. [13,14,15,16,17] 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 were investigated by means of both experimental mechanics and numerical modelling techniques. Thereby moisture is typically characterised either in terms of RH in surrounding 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 [13]. Clearly, the characterization of the hygroscopic expansion and its impact on the mechanics and failure of paper fiber networks require a more comprehensive study.
In the present paper, we combined numerical and experimental approaches to unveil systematically the RH impact on the mechanical behavior and the failure mechanism of paper materials on the fiber network scale. Using the Atomic Force Microscopy and Confoncal Laser Scanning Microscopy, we firstly characterized the humidity dependency of the elastic modulus and the hygroscopic expansion coefficient (HEC) of a cellulose fiber. Instead of a raw fiber before the paper-making process, we studied a cellulose fiber manually extracted from a finished paper sheet. The determined parameters were then used as input in a finite element model to simulate the mechanical behavior and the failure of fiber networks. The force-elongation curves of the network gives an implication on the overall mechanical behavior. We studied particularly the mechanical features like the tensile strength, the effective stiffness and the extensiability. Additionally, by using a cohesive zone damage model, the fiber network simulations addressed the humidity dependent fiber-fiber contact bonding, which was resolved by the finite element cohesive interface elements. It allows the simulation of the local damage at individual fiber-fiber joints and its dependency on the humidity. At the same time, we evidenced single fiber pull-out during the tensile test of a locally wetted paper sheet, by labelling individual fibers and tracking their movement using in-situ imaging. We compared the cohesive zone finite element simulations and the experimental results, and confirmed the damage of the fiber network in a local wet state due to the fiber/fiber separation.

Fiber sample and bending test
For the bending experiments, the cellulose fibers were manually extracted from a cotton Linters paper sheet prepared according to DIN 54358 and ISO 5269/2 (Rapid-Köthen process). The extracted fiber was mounted on a 3D printed sample holder, which supplied a fixed trench distance of 1 mm between the two attachment points. 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 µm colloidal probe (Glass-beads,Kisker Biotech GmbH & Co. KG, Steinfurt, Germany) in diameter. 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.

Determination of humidity dependant Young's modulus
One extracted fiber sample is shown in Fig. 1(a). The fiber was glued at both ends and then force is applied on different positions on the fiber as illustrated in the sub-figure (b). The force-distance curves were recorded for 9 segments framed in the sub-figure (a). The fiber is modelled mechanically by using the beam model illustrated in Fig. 1(b). Based on the images, the fiber cross section is close to a hollow ellipse shown in Fig. 1(c). By use of the Euler beam theory, the differential equation governing the deflection w of the fiber and the corresponding solutions are given as: where w denotes the deflection, E the longitudinal Young's modulus, F the loading force, M the internal reaction moment, I = π 4 c 3 0 d 0 − c 3 i d i 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. Eqs. (2) and (3) can be obtained by calculating the internal reaction moment M and integrating Eq. 1. Then, using the boundary values w(0) = w(a + b) = 0 and w (0) = w (a + b) = 0 for the integration constants. Detailed procedure about solving statically indeterminate beam systems can be found in [18]. Using these equations, the average longitudinal Young's modulus E is obtained by minimizing the difference 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Ī, and then determine the average E along the fiber sections by further dividing the valueĪ = mean(I), which denotes the averaged area moment along the the fiber segments.

Determination of hygroscopic expansion coefficient
A VK-8710 (Keyence, Osaka, Japan) Confocal Laser Scanning Microscope (CLSM) was used to investigate the swelling behaviour of the fiber. For that, the fiber was put into a climate chamber, where the RH was varied as in the AFM measurements. The fiber was suspended 45 minutes to the RH before starting the measurements. The swelling behaviour was analysed by the VK analyser software from Keyence (Osaka, Japan), after creating cross-section images along the fiber, where the local radii were measured. Further, in order to determine the HEC, the change in volume is calculated based on the experimentally recorded change of the cross-section at different RH. The reference volume of the non-prismatic fiber with assumed geometry is Fig. 1 where is the sum of 3D anisotropic hygroscopic expansion tensor β kl with in total 9 independent components in the general anisotropic case. Hereby δ kl is the Kronecker-delta, β T is the HEC in the transverse direction, and β L 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 [15], that the hygroscopic expansion in the fiber length direction, β L , shown to be one order lower than that in the transverse directions, and therefore can be neglected. It leads to the following equation: Afterwards, the HEC in the transverse direction β T can be linearly fitted at different ∆RH for every segments along the fiber length, as explained in Sec. 4.2 later on.

Paper sheet sample and failure mechanism
For studying the fracture mechanism of the paper materials on microscopic level, lab-engineered paper samples with bleached eucalyptus sulfate pulp (median fiber length (length-weighted): 0.76 mm; curl: 15.9 %; fibrillation degree: 5.1 %; fines content: 9.1 %) were used. The paper samples with grammages of 50 ± 1 gm −2 were prepared using a Rapid-Köthen sheet former according to DIN 54358 and ISO 5269/2. The standard procedure was slightly changed, in order to incorporate fluorescent labelled fibers. By adding the dye (Pergasol Yellow F6-GZ liq. cationic dye) to the disintegrated fibers in the fiber suspension, the amount of labelled fibers in the paper samples was controlled. For the best results, 0.1 wt.-% of the fibers were labelled with the dye. Afterwards the paper samples were conditioned for at least 24h under standard conditions (23°C, 50 % RH). A Zwick Z1.0 with a 20 N load cell using the software testXpert II V3.71 (ZwickRoell GmbH & Co. Kg) in a controlled environment with 23°C and 50 % RH was then used for tensile testing with a constant strain-rate of 5 mm/min. In order to analyze the failure mechanism on a single fiber scale, a high magnification with a small field of view was chosen. This was made possible by using a small amount of distilled water to wet the paper samples in a defined region. To analyze the the failure mechanism of the paper samples, a commercially available full-frame mirrorless camera from Panasonic (DC S1) with a macro lens from Canon (MP-E 65mm f/2.8 1-5x Macro Photo) and an adapter from Novoflex (SL/EOS) was used. The camera was mounted on a manual x/y/z-stage on a table that was decoupled from vibrations of the tensile testing equipment. The aperture was set to 5.6, the shutter speed was 1/30 sec and the ISO was set to 800. The videos were recorded with a resolution of 3840 x 2160 px at a frame rate of 29.97 frames/second. A UV-lamp (365 nm) was used from the backside of the paper samples to excite the fluorophore of the labelled fibers. In order to analyze the single images of the recorded videos, they were converted from MP4 to AVI with a FFMPEG-script. Afterwards the image-stacks in AVI format are processed using the program Fiji [19]. Processing included cropping of the images, so that only the part of the failing fiber network is visible, extracting the green-channel (where the fluorescing fibers are most visible) and saving the image-stack in an uncompressed TIF-format for further analysis.

Mechanical model of the single fiber
We apply the linear elasticity to the fibers, including the stress equilibrium, the linear kinematics and the transversely istropic linear elastic material law: in which σ ij is the Cauchy stress tensor, C ijkl the stiffness tensor and ε kl the strain tensor. The strain due to the hygroscopic expansion is given as ε h kl = β kl ∆RH, where β kl and ∆RH 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 follows: with five independent parameters: E L , E T , G T the longitudinal modulus, the transverse modulus, the shear modulus and ν LT , ν T T the two Poisson's ratios. The values of the Poisson's ratios are specified in Tab. 1. Base on [20], 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 Tab. 1 are employed in the following finite element simulations. Note that the longitudinal Young's modulus E L and the HEC are experimentally determined as explained in Sec. 2.2. Further, on the single fiber level, no fiber damage is assumed. Elastic parameters

Cohesive zone interface model for fiber-fiber contact in the dry state
Similar to works [21,20], a cohesive zone-based approach is utilized to characterize the debonding behavior. See Fig. 3(d) for the illustration. In the current work, we applied a non-potential based CZM [22], which aims to give a proper behavior in mixed-mode loading scenario to avoid the fiber penetration. The traction-separation law are given as: where T n , T t are the traction components of T in their normal and tangential loading state. ∆ n , ∆ t are the displacement separation at the fiber-fiber bonds. σ max and τ max are the maximal stresses in pure normal and pure shear separation, δ c n and δ c t the critical length to the maximum stresses, respectively. The energy per surface for both separation modes are given with φ n = δ c n σ max exp (1) and φ t = δ c t τ max 0.5 exp (1). The damage variable D is defined as: ranging from 0 to 1 and represents the intact state and fully damaged state. ∆ ef f = ∆ 2 n + ∆ 2 t , ∆ c = (δ c n ) 2 + (δ c t ) 2 are the effective mixed-mode separation and the separation corresponding to the maximum stresses, and ∆ f = the final separation at complete failure. Eq. 12 and 13 are referred as the traction-separation law 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 is modified by multiplying a decrease term: This is motivated by [23], mainly based on using CZM for hydrogen embrittlement of steel structures, where the effect of hydrogen on the accelerated material damage are studied. Here, the simple assumption is made that the cohesive strength decreases linearly with increasing RH, where the separation distance remains the same. K is a softening parameter and describes the decrease and the bound of the allowable cohesive strength, which needs to be adjusted to the experiment. For a K = 0.5, the traction-separation function are plotted in Fig. 2. a) b) Figure 2: Decay of traction-separation law due to RH with K = 0.5 (a) Normal cohesive behaviour (b) Tangential cohesive behaviour

Weak formulation for cohesive FE implementation
The FE implementation is 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 [24] but with extension to include further the cohesive interface contributions: where δε, δu, δ∆ are the virtual strain tensor, the virtual displacement vector and the virtual separation vector at the interface, respectively. σ is the tensor notation of the Cauchy stress tensor σ ij , T is the traction vector at interface Γ int with the normal and tangential components T n , T i and T ext the external traction on outer boundary Γ. Eq. 16 describes the energy balance between the strain energy or the work of the internal force in the domain Ω and the work contributed by the interface elements and by external forces. The non-potential based cohesive zone model is implemented in an user material kernel in the open source FE software MOOSE [25]. A detailed step-by-step instruction of general cohesive FE implementation from Eq. 16 can be found in [24]. The Software Geodict ® was used to create synthetic fiber network samples, which permits a similar deposition procedure as reported in [26,21]. It uses voxels to represent the structure of the fiber network after the structure is generated for the given geometric parameters. Upon initialization, a voxel size of one micron was chosen, see Fig. 3(b). The fibers are 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 is solid and has 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 step 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 [27]. 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 is given in the supplementary information. For the present work, the fiber network setting and features are given in Tab. 2. The size of the simulated fiber network sample is similar to that used in [21]. For the finite element simulation, 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 M D = 0 is specified on the left boundary, u M D = u on the right boundary and u CD = 0 on the front and back boundary. These correspond to displacement boundary condition in the machine direction (MD) with the boundary cross-section kept as a plane. The out-ofplane (OD) degree of freedom and others on each boundary side were unspecified and stress-free. CD is hereby referred as the cross-machine direction. FEM simulations were performed on a high performance computer with 24 cores for around 8 hours.  Fig. 4(a), along with the polynomial fitting. For simplicity, E L will be written as E in the following.

Fiber network generation and FE simulation setups
It is noticeable that E decreases with increasing RH for different loading levels. This is in agreement with the literature. The Young's modulus reduces because water molecules act as lubricant between the cellulose molecules, as the water molecules infiltrate the network and destroy the hydrogen bonds between the cellulose molecules [28,29,30,31]. The maximum value was obtained at 2 % RH for about 270 MPa. Compared to the values available in the literature, the obtained elastic modulus is relatively low. Fernando et al. [32] obtained from an AFM-based three point bending test a Young's modulus of 24.4 GPa of softwood fibers. From tensile tests of single fibers, a range from 20-80 GPa has also been reported in [33,34,35]. There is one exception [36], which reported a Young's modulus in the MPa range. The low level of our Young' modulus can be attributed to the fibers under testing. As stated previously, we manually extracted the single fibers from a finished paper sheet. Thus, our fiber was 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, the Young's modulus seems to increase with the load. This may stems from nonlinear elasticity or nonlinear geometrical deformation which are neglected in the current mechanical simulations. We average thus the obtained values at different loading. The variation of the averaged Young's modulus w.r.t the relative humidity is then fitted with a polynomial function of 3rd order: where E max represents the averaged modulus at 2% RH. The fitted coefficientsā,b,c,d are given in Tab. 3. The fitted E-RH curve is depicted in Fig. 4(b), along with the deviation due to different loading forces.

Hygroscpic expansion coefficient
The experimentally determined volume change of different segments of a single fiber exposed to different level humidity are summarized in Fig. 5. Due to the local inhomogeneity, deviation of the data for different segments are observed. Therefore volumetric strains are averaged over the segments. The averaged volume strain demonstrates linear dependency on the relative humidity, based on the relation Eq. (7). The transverse HEC β T can be extracted as the slope of the curve, which is indicated by the thick red line, and the value turns out to be around 0.41. This is close to the value of 0.44 reported in [37] as collected by [15] from different data available in the literature, ranging from 0.2 -0.44. It should be noted that the recorded value for the segment 4 was not considered in our calculation, since this segment contained inflated microfibrils that could bias the overall results. There are several other methods in determining the HEC reported in [15], such as measuring from X-µCT or back-calculation from composite properties as in [37,38]. Since our focus is on the influence of RH on mechanical behaviour of the fiber network, the readers are referred to literature for further investigations about the dimensional stability loss in relation with HEC.

Force-elongation relation
Bulk and cohesive finite element simulations were carried out on the generated fiber network under the prescribed boundary conditions in Sec. 3.5. One fiber network sample without humidity influence is shown for demonstration purpose in Fig. 6(a), with the cohesive fiber-fiber interface highlighted in dark blue. The sample deforms under the tensile displacement along the MD. A few snapshots on the deformation and the interface damage level are sequentially presented in Fig. 6(b-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 occur at 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 processes continue 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 varies with the different kind 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 is considered, as indicated by point (c) in Fig. 6. force-elongation behaviour at post-damage state are out of our current scope, where the numerical convergence becomes difficult for implicit finite element solver as previously mentioned. Tab. 4 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 damage parameters. The cohesive parameters describes the mechanical properties of the bond and the energy stored in the fiber-fiber contacts. Fig. 7(a) demonstrates the influence of the critical separation δ c n , δ t n on the force-elongation curves. When increase δ c t , δ c n from 0.01 µm, to 0.1 µm, and then to 1 µm, (in this case the parameters δ f t , δ f n are also set 10 times larger, respectively), the calculated failure force changes correspondingly from 0.8 mN to 1.6 mN, and to 3mN. Similar increasing is observed for the failure strain, from 0.4%, to 0.8%, and finally to around 1.7%. It implies that when δ c t , δ c n is increased by a factor of 10, both the calculated elongation at failure and the failure force are doubled. The effective stiffness of the fiber network does not change noticeably, which is reasonable since the damage parameters discussed here enlarge simply the linear elastic region.
The parameters σ max , τ max effect the force-elongation behavior differently. As it is shown in Fig. 7(b), for the small value of σ max , τ max , namely 1 MPa, the interface damages take place already at an elongation of 0.2%, and the curves becomes nonlinear afterwards due to the damage development at more fiber bonds and the progressive and cooperative interactions between damaged interfaces. The sample reaches a maximal force of about 0.75 mN eventually at the elongation of 0.6%. If we increase now σ max , τ max to 10 MPa (100 MPa), the elongation at failure is improved to about 0.85% (1.05%).
Influence of Young's modulus. As discussed previously, our determined fiber Young's modulus appears to be lower than reported values in the literature. We study hence here the influence of the variation of the elastic modulus on the force-elongation behaviour by varying the elastic constants in the reported arrange. By increasing the Young's modulus by a factor of 10, the stiffness of the fiber network is enhanced by a factor of 4 to 6, as shown in Fig. 8(a). Whereas, the tensile strength of the fiber network increases by a factor of around 3. However, in the meantime, the extensiability decreases with increasing fiber Young's modulus by a factor of about 2. It can be concluded that the stiffness of the fiber network is greatly driven by the stiffness of the constituent fibers, instead of bond properties given by the cohesive parameters. Influence of relative humidity. RH influence on the force-elongation curves are presented in Fig. 8(b). Note that in our simulations two key parameters are subjected to RH dependency, namely the fiber Young's modulus and the cohesive traction-separation law. The elastic moduli based on our previous bending test of different RH state are used here accordingly. Due to the lack of available experiment data on the fiber/fiber bonding at different RH state, we assume the decaying parameter K used in the humidity dependent cohesive zone model in Eq. 15 to be 0.5. It means that in the extreme case for ∆RH = 100%, the cohesive strength is decreased by 50%, while the critical separation remains the same. As depicted in the figure, for a reference dry state as RH = 2%, with an increase of RH from 2% to 90% , the stiffness of the fiber network reduces to 20% of the dry sample. The force to failure decreases to 40% of the maximal force in the dry state. However, the extensiability is almost doubled. Given the sorption isotherm data published in [39], we compare our RH dependent simulation results with the moisture dependent experimentally data in [40]. They obtained their data from tensile test of a paper strip with dimension of 50mm x 15mm x 200 µm. To make a clear comparison, the ratios/changes of the mechanical properties for the simulated data between two different humidity states, namely at RH=75% and RH=90%, and the experimental case for the moisture content MC=10% and MC=14% at room temperature are studied, given the sorption isotherms previously mentioned. The simulated and the measured results on the effective stiffness, the failure force and the elongation at failure are compared in Tab. 5. Simulation results show that the RH change leads to decrease of the effective stiffness and of the failure strength, but increase of the elongation at failure and thus, the extensiability. The observations show qualitative agreement with the experimental data. Note that due to the difference in the microstructure, as well as the length scales, it is not reasonable to compare the exact numbers. Nevertheless, our simulations deliver important information for understanding and designing of paper materials on the fiber and fiber network scale.
Since the microscopic fibrous network manifests always a certain degree of randomness, the mechanical properties are certainly sensitive to the geometrical variations. The effect of geometrical variation and the uncertainty quantification of micro-structural features have been investigated in our previous work [41] using a data-driven and machine learning based approach. It must be mentioned that the mechanical behaviour of the fiber network is also strongly size dependent as reported in work [42]. The strength of the paper sample follows a weakest-link scaling law. The density of the fiber network naturally plays an important role, as increasing the basis weight of the paper sheet can lead to increase in its strength and stiffness as well. More details regarding these aspects can be also found in studies by [26,43].

Fiber network failure: Fiber pull-out
As described in the section 2.4, tensile test of paper sheet were carried out. Paper sheets were locally wetted, and individual fibers were fluorescently stained to allow tracking their movement and deformation during the successive rapture. To be in line with axis convention of the simulation condition in Sec. 3.5, loading direction is denoted as MD and the perpendicular direction is denoted as CD. However, it must be noted that the lab-engineered paper sheets are close to isotropic and therefore, the fibers do not show significant orientation alignment. MD and CD are rather axis convention in this setting.
In the following, the considered two exemplary fiber networks with specific single fiber alignment w.r.t MD are presented. At the same time, corresponding finite element simulations of the fiber network samples including the cohesive fiber contact model were carried out for comparison. Simulation boundary conditions as in Sec. 3.5 were slightly modified, namely u CD = 0 on the front and back boundary are discarded to allow more fiber movement in the failure analysis. The results for the tensile test with fiber alignment (The Fiber of interest is marked as region of interest (ROI)) within MD are presented in Fig. 9. The experimental images in Fig. 9(a) shows the initial positions of the fiber before the tensile test. As the subsequent snapshots Fig. 9(b-c) taken during the rupture clearly demonstrate, those fibers, which lie rather parallel to the loading direction (MD in this case), are pulled out from the neighbouring fiber networks. For better understanding, the crack surfaces are indicated in the figures. No fiber breakage was observed for all the visible illuminated fibers near the crack surfaces. This is evidenced by the facts that the illuminated single fibers retain almost their initial lengths during the whole process. It implies that the fiber-fiber debonding is the main failure mechanism in this tensile test. This is in agreement with the observation reported in [44] for the wet tensile test. Moreover, it is noticeable that the fibers oriented in the loading direction barely undergo any deformation or reorientation. These experimental observations are strongly supported by the fiber network simulation results shown in Fig. 9(d-f). In particular, one fiber aligned in the MD is highlighted in red and tracked during the tensile test. Due to the numerical convergence issue related to the rigid body movement after progressive rupture, simulation results cannot be shown for the final fully separated stage of the fiber pull-out. But, Fig. 9(e-f) show the apparent debonding of the fiber in red from the neighbouring fibers. The fiber is slightly deformed, but remains the original orientation in the MD.
The scenario changes in the case of the tensile test with fiber orientated in the CD. The initial image with the illuminated fibers is presented in Fig. 10(a). Again, we particularly studied the fiber marked as ROI, as the fiber is aligned perpendicular to the loading direction. Fig. 10(b-c) illustrate that this fiber is progressively deformed and reoriented during the test. This can be partly explained by the bending effect resulted from the relative geometry relation between the fiber orientation and the loading direction. It can also be attributed to the bond orientation and separation process, where most of the separation are orthogonal to the loading direction. Considerable deformation and reorientation of one similar fiber has also been observed in the finite element simulations Fig. 10(d-f). Both the experimental images and the simulations indicate that after the reorientation of the fiber, it is pulled out from the opposite crack surface.

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 poly-fitted the humidity dependent fiber elasticity. The obtained maximum Young's modulus at 2% relative humidity is about 270 MPa and rather lower than values published in the literature. It can be attributed to the fact that the tested fiber is rather a fiber from a finished paper sheets than a raw one without the influence of a paper-making process. Further, the obtained relationship allows calculation of fiber Young's modulus at any humidity level. Hygroscopic expansion coefficient was obtained to be 0.41 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 is applied, jointly with the humidity dependent elasticity to simulate the mechanical behaviour of fiber network under tensile test on a higher scale. On the fiber network scale: The mechanical behaviour of simulated tensile test under humidity influence was in good qualitative agreement, compared to the experimental data on the paper sheet scale found in literature. Experimental tensile tests of paper-sheet were carried out to study the local failure system. The fiber network failure in a locally wetted region was revealed by tracking individual fluorescently stained fibers using fluorescence macro videography. Compared to the simulated fiber network failure, both experimental tensile test and the cohesive finite element simulations demonstrated the pull-out of fibers and implied the significant role of the fiber-fiber debonding in the failure process of the wet paper.