Fast transport diffusion of bound water in cellulose fiber network

A remarkable property of cellulose-based materials is that they can absorb huge amounts of water (25% of the dry mass) from ambient vapor, in the form of bound water confined at a nanoscale in the amorphous regions of the cellulose structure. The control of the dynamics of sorption and desorption of bound water is a major stake for the reduction of energy consumption and material or structure damages, but in the absence of direct observations this process is still poorly known. Here we present measurements of bound water transport thanks to Nuclear Magnetic Resonance relaxometry and Magnetic Resonance Imaging measurements. We show that the bound water is transported along the fibers and throughout the network of fibers in contact. For each material a single transport diffusion coefficient value allows to represent the processes over the whole range of saturation. The dependence of the transport diffusion coefficient on the fiber density and orientation is then analyzed to deduce the (elementary) transport diffusion coefficient of bound water along a cellulose fiber axis. This constitutes fundamental physical data which may be compared with molecular simulations, and opens the way to the prediction and control of sorption dynamics of all cellulosic materials or other hygroscopic materials.


Introduction
Cellulose, the world's most abundant natural, renewable, biodegradable polymer, is a major component of plants, wood, algae, bacteria and is widely used in food, industrial, pharmaceutical, paper, textile production, or in wastewater treatment applications. A remarkable property of cellulose-based materials is that they are hygroscopic: they can absorb huge amounts of water (typically 25% of the dry mass) from ambient vapor, in their solid structure, thanks to the high sorption affinity of the cellulose chains for water molecules in particular near their polar OH groups (Berthold et al. 1996;O'Neill et al. 2017;Kulasinski et al. 2015). This "bound water", which is at the origin of swelling or shrinkage of such materials, is confined at the nanoscale in the amorphous regions between cellulose microfibrils (Kulasinski et al. 2015;Paajanen et al. 2019). Note that in wood, lignin and hemicellulose, which are predominantly amorphous (Hill et al. 2009), do also contain bound water. Also note that, considering the high mobility of this water, which we will further demonstrate here, it is somewhat paradoxical to call it "bound".
The control of sorption and desorption of bound water is a major stake for the reduction of energy consumption and structure damages. Sorption or desorption of bound water in the structure are at the origin of the swelling or shrinkage of cellulosic materials (Kulasinski et al. 2015;Simpson 1999) which can cause wood structure damages (Simpson 1999;Jakes et al. 2019). For paper production, the extraction of the final (small) bound water fraction is the main energy consumer in the paper mill (Park et al. 2007), and the majority of the functional properties of paper are developed during this stage (Ghosh 2011). Also, gradients of bound water content caused by storage, heating or ink printing can lead to mechanical instabilities like curls or cockles (Bosco et al. 2018;Zapata et al. 2013;Dano and Bourque 2009). Bound water transport and storage is at the origin of the moisture buffering characteristics of bio-based construction materials (Osanyintola and Simonson 2006;Li et al. 2012). It also impacts the moisture transfers through textiles (Barnes and Holcombe 1996) and the heat loss due to sweating, since water-vapor adsorption is an exothermic process (Havenith et al. 2008). At last, it was recently shown that the extraction of free water in depth during wood drying relies on bound water transport diffusion towards the upper layers (Cocusse et al. 2022;Penvern et al. 2020). Thus, the control of bound water transport is a major stake for the reduction of energy consumption in a wide variety of applications involving cellulosic materials. Actually, there is a similar issue, although in general with smaller water contents, for other types of polymers containing amorphous regions (Peppas and Brannon-Peppas 1994;Litvinov 2015;Monson et al. 2008).
Despite its ubiquitous role, the transport properties of bound water within cellulosic structures, which govern the dynamics of the above-mentioned processes, are poorly known. Standard measurements (Simpson and Liu 1991) rely on the observation of mass variations of a sample placed in air at controlled humidity at some distance from the sample surface. It is then considered that the water progressively penetrates deeper into the sample according to a transport diffusion process. However, in general, the diffusion coefficient which is deduced from such tests results from an analysis which doesn't distinguish bound water from vapor transport through the porous structure and describes the boundary conditions along the sample surface with the help of unknown coefficients, i.e., the "surface emission coefficient" for wood (Siau and Avramadis 1996;Cai and Avramadis 1997;Choong and Skaar 1972), or the "clothing vapor resistance" for textiles Havenith et al. 1999;Liu et al. 2020). This precludes quantitative predictions and good control of the processes. The transport diffusion coefficient estimated through such approaches was generally observed to increase with the water content (Simpson and Liu 1991;Perkowski et al. 2017;Liu and Simpson 1999) but some works concluded that it remains constant or even decreases when the water content increases (Perkowski et al. 2017;Liu and Simpson 1999;Nelson 1986).
On the other hand, one can consider the self-diffusion of bound water, resulting solely from the thermal agitation of the molecules. In that case, in contrast with the above process, there is no net transport of matter. The self-diffusion coefficient of bound water in pure cellulose was estimated through direct internal measurements by Pulsed Field Gradient NMR (Nuclear Magnetic Resonance) (Topgaard and Söderman 2001;Li et al. 1992;Araujo et al. 1993;Perkins and Batchelor 2012) or from quasi-elastic neutron scattering (O'Neill et al. 2017), providing values typically at least order of magnitude smaller than the selfdiffusion coefficient of water (i.e. 2.3 × 10 -9 m 2 /s). Even smaller values were obtained from molecular simulations (Kulasinski et al. 2014), with a decrease when the water content decreases. However, it is worth noting that in the above experimental works, the self-diffusion observed is in fact directed along the fiber axis in a 3D fiber network with various orientations and contacts, and exchanges are possible with the surrounding vapor. All these aspects are not accounted for in the analysis. Moreover, in such complex nanoporous systems the mechanisms of selfdiffusion in a homogeneous material may differ from those of the transport diffusion associated with a net transport due to a moisture gradient (Falk et al. 2015).
Here we solve these problems by direct measurements of the (transport) diffusion coefficient associated with diffusion resulting from a gradient of moisture concentration. The pores of a cellulose piling are filled with oil, which precludes vapor transport, and the amount of bound water as well as its spatial distribution are followed over time by NMR relaxometry and MRI during drying. We show that the bound water is extracted from the structure by diffusion along the length of the fibers and via fiber contacts. This result can be extrapolated to any type of cellulosic structure (textile, wood fiberboard, etc.): simple contacts allow bound water transport throughout the structure. The phenomenon can be described by a simple diffusion process with a single transport diffusion coefficient, while properly taking into account boundary conditions through the solution of the convection-diffusion equation. Further tests with different fiber orientations make it possible to estimate the longitudinal (along fiber axis) transport diffusion coefficient of bound water in fibers. Remarkably, although the bound water molecules are confined in nanometric pores between cellulose microfibrils, we find a transport diffusion coefficient in the order of the coefficient of self-diffusion of water. This opens the way to a straightforward, quantitative prediction of the internal dynamics of humidification or drying processes resulting from variations of ambient humidity conditions, and their consequences on various cellulose-based materials such as textiles, biobased construction materials, wood, etc.

Materials
We prepare the materials from disordered natural cellulose fibers (Arbocel ® BC 200, from Kremer Pigmente) placed in a highly humid air (Relative humidity (RH) about 97% at a given temperature of 22 °C ± 1 °C) for at least two weeks. The material was evenly spread as a cellulose fiber stack, approximately 1 cm thick, over the bottom surface of a desiccator, and we could check that equilibrium was reached after a few days by monitoring the mass vs time evolution. The uncertainty on the temperature may induce some variations of the RH inside the desiccator. However, since the possible variations are slow, this is the same for the water content inside the material, and this should not affect the homogeneity of the moisture content in the final sample. This is the critical point, since for a diffusion process the results and analysis of the drying tests, followed from the saturation vs time, do not depend on the absolute MC value, as long as it is initially uniform. Due to the porous nature of cellulose fiber stack, the vapor penetrates throughout the structure, which ensures that a homogenous bound water content approaching the maximum possible value is reached at equilibrium (see sorption curve data in Supplementary Information). This was confirmed from previous MRI observations of the bound water content distribution in cellulose stacks in this initial state (see Ma et al. 2022). Some amount of such material is then compressed so as to obtain a sample of controlled porosity (void to total volume ratio, noted ) in a wide range, i.e. from 0.2 to 0.8. Note that under these conditions corresponds to the porosity in the saturated state. This compression appears to be sufficient to get a solid material, which will not break during the preparation and the test. This results from the formation of a continuous network of bonded fibers. The nature of these bondings is still debated, it has been suggested (Hirn and Schennach 2018) that either hydrogen bond, van der Waals attraction, or Coulomb interactions, could prevail, in addition to interdiffusion increasing the area of contact. Whatever its exact origin the formation of this solid network ensures that there exists some area along which the external microfibrils of two such fibers in contact are very close to each other. An example of the typical aspect of the resulting structure may be observed by microscopy (Hitachi TM4000 Tabletop SEM (Scanning Electron Microscope)) (see Fig. 1a). Since such images are a projection of all visible objects along one axis, and due to the entanglements of fibers, no differences in the aspect of the structure for different porosities and in different directions could be detected. The cellulose fibers have a strip-like shape with an average length of 300 microns, a width of 10 microns and a thickness about 1 micron (see Fig. 1a). Each fiber also exhibits some slight helical shape, which is in fact found at the different scales of its structure (Paajanen et al. 2019). In this context, we will define the longitudinal axis of a fiber as its main axis as observed at a macroscopic scale, i.e. thus including this twisting trend. The dry cellulose fiber density is considered to be s ≈ 1500 kg m −3 (Gibson and Ashby 1997). The moisture content ( s ) is the bound water to dry cellulose mass ratio, and its maximum value s m is the moisture content at saturation. Here we have s m = 0.25 . We can estimate the density ( ) of the cellulose fibers at a moisture content s assuming the density of bound water is close to that of bulk water, i.e. w ≈ 1000 kg m −3 . We thus have = (1 + s) ( −1 s + s −1 w ) . The density of the saturated fibers is then m = (s m ) . For cylindrical samples of cross-section area A and height h made of saturated fibers the initial porosity is directly deduced from = 1 − m m hA , in which m is the sample mass. The sample size was 12 mm diameter and 8 mm thick for NMR measurements, and 50 mm diameter and 17.5 mm thick for MRI measurements. The constant thickness of our samples despite porosity variations implies that the initial material mass before compression is adjusted as a function of the expected porosity. Note that, here, for these calculations, for the sake of simplicity, we take for the reference of zero water content the material state at the end of drying under our dry air flux (which is an approximation).
The solid samples are then placed again at a high humidity for some time to ensure they are in a homogeneous saturated state at the beginning of the experiment. Afterwards, still at a very high relative humidity, the sample bottom is placed in contact with a bath of olive oil which spontaneously invades this porous medium, This occurs in a few tens of minutes for high porosity samples and a few hours for low porosity samples. The volume of the bath is limited to 85% of the pore volume of the sample, but the oil does not simply fill the first 85% of the sample height, it spreads somewhat farther. This finally allows to obtain imbibition or wetting by oil over the whole sample height, except in a very thin layer at its top surface which clearly remains un-wetted with oil, as appeared from direct observations (see Fig. 2) and from MRI data (see below).

Set up
The sample is then placed inside the magnet (see Fig. 3a) of the spectrometer and subjected to a dry air flux, 2 l/min for NMR measurements and 10 l/min for MRI measurements, induced by a tube situated at distances of 1 cm and 3 cm respectively, from its free surface (see Fig. 3a). This set up essentially results in a tangential air flux along the top surface of the sample. Two types of tests were carried out: evaporation in the same direction as the (previous) compression, and evaporation in a direction perpendicular to the compression axis. In the latter case a parallelepiped sample was prepared by compression, so as to have Due to the presence of oil in the (large) pores between cellulose fibers and coating their surface, no vapor diffusion is possible. With our technique for filling the sample with oil, i.e. without imposing some pressure, there could remain some voids, i.e. regions not reached by the oil in the pore network of the fiber stack. However, since the oil spontaneously invades the sample, it necessarily wets most of the fiber stack and thus is in contact with most of the cellulose fibers, so that a continuous network of such voids which would allow some vapor transport cannot exist. Such possible voids will be essentially "neutral" with regards to the drying process. Moreover, the humidity in the air of such voids will be at equilibrium with the moisture content of the fibers, preventing the formation of free water pools.

NMR measurements
We then observe the internal processes from 1 H NMR relaxometry with a Bruker NMR minispec mq20, 0.5 T, 20 MHz. We measure the transverse relaxation time T 2 using a Carr-Purcell-Meiboom-Gill (CPMG) sequence (Carr and Purcell 1954) with the main following parameters: Number of scans: 128, Echo time: 1 ms, Number of echos: 2000, and Repetition time: 4 s (to ensure a total equilibrium recovery magnetization), for a total measurement duration of 13 min for one T 2 distribution, which is much lower than the total duration of the drying test. Note that such a standard CPMG sequence is sufficient, i.e. there is no need to use a sequence appropriate for hydrogen spin in the solid phase, as the relaxation time of bound water is still relatively high, which finds its origin in the rather high mobility of these water molecules, as will be shown below. It is also worth noticing that the choice of the echo time and the number of echos was done after a series of tests with different parameters and different samples (with or without oil) and results from a compromise between different constraints distribution over the whole sample) (b) and MRI (successive 1 D profiles of water distribution along sample height) (c) measurement techniques used to explore the water content during sample drying (data treatment, temperature increase, range of T2 covered, etc.).
After a Laplace inversion of the signal vs time curve, we finally get the pdf (probability density function) of the NMR relaxation time T2 (see Fig. 3b). Note that the Laplace inversion still suffers from some important drawbacks and artifacts (Faure and Rodts 2008;Borgia et al. 1998), and the exact shape of this distribution somewhat depends on the processing protocol (Lerouge et al. 2020). However, for given NMR sequence parameters, the relative amplitude of the peaks, as well as the evolution of their shape over time, are robust data (Lerouge et al. 2020) and may be used to follow the water distribution in porous media containing pores from the nano to the macro-scale . Different ways of analysis of the data were tested which led us to consider that the integral of the pdf below the peaks was the most reliable approach to follow the water amount evolution in the sample. Also, tests with different echo times, from 1 to 0.06 ms, but keeping all other parameters identical, were carried out to check the robustness of our pdf data. The pdf associated with water appeared to be similar in this range of echo times, which proves that the pdf obtained for an echo time of 1 ms (the value used for all our tests in this work) indeed provides the complete NMR signal of water in the sample. At last, a series of tests with different masses of copper sulfate solutions with a relaxation similar to that of water in cellulose were carried out, which confirmed that the integral below the water peak in the pdf well is simply proportional to the water mass in the system and may thus be interpreted in this way with our cellulose samples.

MRI measurements
Moreover, we can have a look at the spatial distribution of bound water along the sample axis from an original MRI approach, analogous to that used to distinguish bound from free water during wood drying (Cocusse et al. 2022), but here used to distinguish bound water and oil. These 1 H MRI measurements were carried out with a 24/80 DBX 0.5-T 1H MRI spectrometer by Bruker, operating at 20 MHz. The spatial distribution along sample axis (i.e. 1D profiles, Fig. 3c) of all liquid molecules, i.e. including oil and bound water molecules, was measured with the help of the Single Point Imaging (SPI) sequence (Emid and Creyghton 1985;Bogdan et al. 1995). The parameters of the MRI sequence used are: Number of scans: 256; Signal collecting time: 480 µs; Field of view: 40 mm; Number of pixels: 64; Flip angle: 15°; Repetition time: 200 ms, for a total measurement duration of 54 min. Secondly, with the same MRI system we measure the T 2 distribution with a multi echo (ME) sequence (Zhou et al. 2019) using the following parameters: Number of scans: 512, Echo time: 8.6 ms, Number of echos: 16, Field of view: 40 mm, Number of pixels: 64 and Repetition time: 1 s for a total measurement duration of 8.5 min. From this sequence, which is only sensitive to liquid molecules with relaxation times larger than a few tens of milliseconds we can determine the 1D profile of oil only. Note that this was done here by leaving apart the first three echos for which the bound water signal might not be fully negligible, and a mono-exponential fit was used.

Results
We find that a typical pdf for our sample exhibits mainly two peaks, one around 100 ms, and another around 2 ms (see Fig. 4). Interestingly, a sample saturated with bound water but without oil exhibits a pdf with only one peak similar to the peak at 2 ms in the presence of oil (see Fig. 4), while free water, i.e. outside these nanopores, would exhibit a much larger relaxation time (Ben Abdelouahab et al. 2021). We deduce that the peak at 2 ms for the sample with oil corresponds to bound water saturating the cellulose while the peak at 100 ms corresponds to oil, a fortiori outside the nanopores. Thus, the oil fills the porosity between the cellulose fibers, but does not penetrate the nanoporosity. These results also confirm that there is no (free) water pools in the sample, since there is no significant difference of the pdf of water in the cellulose sample whether there is oil or not, and such pools cannot exist initially in our oil-free sample as it has been prepared by simple sorption at a given RH. Water pools of free water would give some signal at higher T2, as is the case for example for a cellulose paste (cellulose fibers mixed with water) (Ben Abdelouahab et al. 2021).
Moreover, the pdf of oil filling a dry cellulose sample appears to be quite similar to the pdf part associated with oil in the bound water saturated sample filled with oil (see Fig. 4). This proves that there is no penetration of oil inside some small pores of the cellulose fibers themselves, otherwise we would in principle observe peaks at smaller T2 values. Moreover, this proves that there is no particular interaction of the oil with the water in oil + water sample: the oil is apparently in the same state, not particularly mixed with water, and distributed in the same pores, whether the cellulose initially contains water or not.
Under these conditions, since desorption in the form of vapor is prevented in most of the sample, we could expect the bound water situated below the top layer to be blocked inside the cellulose fibers. Yet, as it appears from the evolution of the pdf over time (see Fig. 4) the sample fully dries relatively rapidly, i.e. typically in less than 2 days for a 0.8 cm sample thickness: the total amount of water amount associated with the bound water peak, computed by integrating the pdf over the corresponding range of relaxation times, progressively decreases until the peak finally vanishes (see inset of Fig. 4). Note that the slight increase of the total signal for oil is due to artefacts resulting from the Laplace transform procedure.
Such a result is remarkable. Since bound water cannot be extracted from the sample by diffusing through oil, it means that bound water diffuses along the fiber network up to the sample free surface. As a consequence, our observations prove that i) Bound water is transported inside the fiber, towards the sample free surface ii) It is able to jump from one fiber to another in contact with it Moreover, we can have a look at the spatial distribution of bound water along the sample axis from MRI measurements (see Methods). The spatial distribution along sample axis of all liquid molecules, i.e. including oil and bound water molecules, was measured with the help of the SPI sequence while the 1D profile of oil only was measured with a ME sequence (see Fig. 5). From these measurements we can first see that, although the total volume of oil was less than 85% of the pore volume, it almost reaches the sample top (see Fig. 5). This results from a partial saturation of the pore in some region below the sample top, say over a 5 mm thickness. Nevertheless, the presence of some significant concentration of oil in this region ensures that the fibers are coated with oil thus preventing water vapor extraction and transport. Moreover, the oil profile remains approximately constant, which confirms that oil is not transported through the medium. In contrast, the total NMR signal (i.e., for all species) varies over time, which proves that some bound water disappears. Also, we can remark that this process starts from the sample top, where the SPI profile curves and penetrates farther in the sample. In a second stage, say, after 5 h, deep regions start desaturating and the level of the profile in depth decreases almost as a whole over time (see Fig. 5).
In order to have a more straightforward appreciation of the bound water distribution over time we need to look at the difference between these two types of profiles. However, the 1D profile for all liquid molecules (from SPI) and for oil only (from ME) cannot be directly subtracted to get the bound water filled with oil. The dashed blue curve corresponds to the initial state of a bound water saturated sample without oil, the orange circle curve to dry sample filled with oil. Note that for these latter tests the NMR signal of the pdf was rescaled to reach an amplitude similar to the initial ones for the test with both oil and bound water. Also note that no significant signal was observed out of the range of relaxation times presented here. The inset shows the integral of the signal (peak) associated with each liquid phase (bound water and oil) as a function of time for the drying test. Remark that the bound water peak for low water saturation extracted from the Laplace transform is likely influenced by the oil peak, which explains that its T 2 does not decrease exactly as observed  when decreasing the saturation in the absence of oil distribution as the signal intensity is not calibrated in the same way for two independent sequences. In principle, we could just determine the proper factor of calibration by adjusting the final SPI profile to the oil profile, but the test could not be carried out over the full duration needed for complete drying. In order to determine the factor to apply on the NMR signal of oil profile, so that it becomes consistent with the total liquid signal, i.e., the same value represents the same oil mass, we assume that, in the second stage of drying, the fraction of NMR signal due to water around the sample top surface, i.e., at the right end of the profile, should tend to a value close to zero. This assumption will be validated a posteriori through the 1D profile predictions (see inset of Fig. 5) of the diffusion model fitted to the saturation vs time curve (see below). The factor is thus determined by a trialand-error procedure allowing to identify the value for which the resulting bound water profiles, i.e. after subtraction, over most of the drying duration, converge towards zero. Note that to facilitate the above analysis, in the representation of Fig. 5 the 1D profiles for oil had been multiplied by this coefficient. We finally proceed to the subtraction between the current water + oil profile and the oil profile at the same time, with a factor, here equal to 24, on the NMR signal in the latter case. Note that in all the analyses here we neglect the sample deformations, which result from the extraction of bound water and impact the whole sample thickness. Since the apparent variations of upper position of the oil profile remain smaller than 1 mm (see Fig. 5), the sample deformation remains smaller than 6%, which has a negligible impact on the transport properties. This means that, although the fibers desaturate, and thus shrink, the overall size of the fiber network does not change significantly. This is due to the fact that the main impact of the presence of bound water between the microfibrils, whose main axis is along the fiber axis, is a variation of the fiber thickness, while the induced variations of fiber length are negligible in this context. Finally, while each fiber shrinks the network of fibers in contact remains essentially unaffected.
Although the bound water profiles thus obtained are affected by the inherent noise on NMR data and some possible material heterogeneity, the qualitative evolution of these profiles over time is reliable. The successive profiles for bound water over time resemble classical concentration distributions observed during a diffusion process (Crank 2011), for example when a saturated plane sheet is initially placed in contact with a medium which can absorb mobile species: the saturation starts to decrease abruptly around the sample top, forming a front which progresses farther in the sample and finally the level of the saturation plateau in depth starts to decrease (see inset of

Diffusion model
We can go on with this concept and assume that the bound water is transported thanks to a diffusion process in the fiber network. Under these conditions, the drying test can be fully modelled. The internal transport is described by a standard diffusion equation: in which z the sample axis, and D the transport diffusion coefficient. Assuming that the sample does not swell or shrink during drying, the corresponding The Eq. (1) can be solved as soon as one knows the boundary conditions. Obviously one of them is the no-flow condition through the sample bottom. The other boundary condition is imposed by the air flux. This flux imposes a low vapor density, described through the RH ( n ) in the air around the sample free surface. This induces desorption of the cellulose fibers in this region, since the initial bound water saturated state corresponds to equilibrium with a RH approaching 100%. The vapor resulting from this desorption then diffuses through the air flux. The vapor diffusion can vary during the drying process, as it depends on the RH along the surface, i.e., n s , which itself depends on the moisture content in the cellulose along this surface, i.e., s s , which decreases during drying. This explains the decrease of the drying rate over time generally observed in such tests.
However, a critical point to note is that the imposed dry air flux has constant characteristics, which are negligibly affected by the vapor diffusion through it. Then, for a unidirectional flow over a smooth planar interface, a dimensional analysis of the mass and momentum conservation equations shows (Cocusse et al. 2022) that the vapor mass flux can be expressed in a general way as J = 0 D 0 (n s − n ∞ ) , in which 0 = 0.02 kg m −3 is the saturation vapor density at the ambient temperature (22 °C), D 0 is the vapor diffusion coefficient in air, n ∞ is the relative humidity of the air flux (here n ∞ ≈ 0 ). In this expression, the parameter may be seen as an equivalent thickness of vapor diffusion; this is a function of the (dimensionless) Reynolds and Schmidt numbers (see Bergman et al. 2011), which depend only on D 0 and on the characteristic velocity of the air flux, the air density and viscosity, and the characteristic length of the surface. For a rough interface, would also depend on the roughness. It follows that, as long as the shape of the interface does not change, which is the case here if we neglect some slight deformations resulting from moisture content variations, is a constant for a given air flux. As a consequence, it can be determined from the initial (measured) drying characteristics, or more precisely the initial drying rate J 0 . Indeed, in the initial state we have s s ≈ s m and thus n s = 1 , so that = 0 D 0 J 0 . Note that J 0 may be deduced from the initial saturation vs time slope ( Ṡ t=0 ) through J 0 = m 0Ṡt=0 A , in which m 0 is the initial water mass in the sample, and the saturation is defined as S = s s m . For all our data based on drying tests characterized by NMR relaxometry we found = 1 ± 0.1 mm.
At last, along the air-cellulose interface we have a sorption equilibrium between the relative humidity and moisture content. We end with the following boundary condition along the sample top surface: in which f is the sorption function (i.e. s = f (n) at equilibrium, see Supplementary Information). The Eq. (1) along with the boundary condition (2) may be solved and fitted to data by adjusting the value of D.
Note that here we assumed a constant transport diffusion coefficient independent of the moisture content, which contrasts with the conclusions of various studies (See Introduction).
Our model rather well predicts the trends observed by MRI (see inset of Fig. 5), but the relatively large uncertainty on these data resulting from the successive steps of measurement and analysis does not allow a precise determination of the transport diffusion coefficient. Instead, NMR relaxometry allows to get a precise measurement of the total amount of bound water in the sample. The evolution of the corresponding bound water saturation over time may be compared with the model prediction by adjusting D . Note that we here consider the average saturation over the sample. This saturation is determined experimentally from the first NMR signal recorded during the relaxation, after a time much shorter than the relaxation time of bound water, which ensures that it encompasses all the water contained in the sample. The value of this first signal then evolves in time as the material dries and loses bound water. Typical results for different porosities and directions of evaporation with respect to the compression direction are shown in Fig. 6. For evaporation in the direction of compression the drying is slower for smaller porosities. For evaporation in a direction perpendicular to the compression direction the drying rate does not seem to depend much on the porosity in the range of observation.
The agreement between the model and such data is good for compression parallel to evaporation, and it is excellent for compression perpendicular to evaporation (see Fig. 6 and all data in Supplementary Information) even in the very last stages of drying, i.e. when the saturation tends to zero. Note that the uncertainty on the value of the transport diffusion coefficient found from this fitting is very low, typically less than 3%, whatever the fitting technique used. Here we present results from a least square method based on the Levenberg-Marquardt algorithm. From these results it is clear that in order to represent the whole desorption process of the cellulose fiber piling, a simple diffusion model with a single diffusion coefficient is sufficient. This contrasts with the expectations of a variation of the diffusion coefficient with the saturation. Moreover, this good agreement between the model prediction and the saturation vs time curve proves that there is no significant heterogeneity of porosity along the compression axis. Indeed, such heterogeneity would induce a spatial gradient of diffusivity, which would induce some discrepancy between the model and the data. Considering the complexity of the medium under consideration it seems more appropriate to refer to this transport diffusion coefficient as a diffusivity. Our experiments thus allow to find the diffusivity of a fiber network under different conditions (evaporation parallel or perpendicular to the compression axis) and different porosities (see Fig. 7). It corresponds to the transport diffusion of water in the cellulose network along various directions and contacts. Remarkably, for evaporation in the same direction as compression, the diffusivity coefficient logically increases with the porosity, i.e. it decreases for more compression, which favors the fiber orientation in the direction perpendicular to evaporation. In contrast, for evaporation perpendicular to compression axis, the fibers tend to be oriented along the evaporation axis, and the diffusivity is much larger. These observations suggest that the diffusion essentially occurs along the fiber axis. Note that there is a significant scattering of data (see Fig. 7). This likely comes from the partially controlled initial state of the material (before compression), i.e. in some cases there might be some anisotropy in the distribution of the orientation of the fibers. Also note that we could not test samples at porosity larger than 0.8, as in that case the penetration of oil induces a significant swelling of the sample which becomes too soft and cannot be handled, Fig. 6 Saturation vs time as measured by NMR relaxometry (see text) during the drying of initially saturated cellulose fiber stacks filled with oil with different porosities (see legend) and for evaporation direction either parallel (inset (a)) or perpendicular (inset (b)) to the direction of compression. The continuous lines are the model (1-2) fitted to data for tests with oil (with = 1 mm) Fig. 7 Diffusivity as determined from saturation vs time curves (see text) for different orientations of evaporation with respect to compression axis, and different porosities. The crossed (respectively open) symbols correspond to evaporation perpendicular (respectively parallel) to the compression direction. The circles correspond to tests carried out by following the sample mass. The predictions of the microstructural model (see text) for each case are represented by continuous lines. Note that the application of this model requires to assume a specific value of the transport diffusion coefficient along the fiber axis; here we used D L = 3 × 10 −9 m 2 s −1 . The intersection point between the two curves of the model corresponds to the initial state of the sample (assumed isotropic) i.e. during its insertion in the container for the drying test, without too large deformations.
Note that since we now know that the bound water is extracted thanks to a diffusion process inside the fiber network and without interactions with oil, similar experiments can be carried out under more standard conditions, i.e. by following the sample mass during its drying. These sample mass variations correspond to the transport of bound water. However, in that case the air flux on the sample induces some fluctuations of the pressure on the top surface of the sample, leading to fluctuations of the apparent mass so that we have to use a sufficiently large sample volume to reduce their impact on the measurements. As a proof of concept we used a sample diameter of 5 cm. However, with such a large diameter it is more challenging to obtain a uniform front of oil imbibition. Moreover, during drying, the sample surface can more easily swell unevenly, which can induce additional uncertainty on the results. We could however carry out two tests for which such artefacts were minimal. They show a saturation vs time curve analogous to those found from the NMR measurements (see below); note that a direct comparison of such curves with NMR data is not of great interest because the different conditions are not exactly the same: air flux, sample dimensions and porosity. Finally, the values of the transport diffusion coefficient found from these data, and which we added in Fig. 7, are close to the values found by NMR. This provides an additional validation of our NMR measurements and show that an alternative technique might be used while much more difficult to carry out without artefacts. Moreover, NMR tests remain more interesting as they allow to check that artefacts such as water pools or mixing with oil do not occur, they are not sensitive to the vertical pressure induced by the air flux and its possible fluctuations, and they allow to use samples of small diameter, which minimizes artefacts developing during imbibition and drying.

Microstructural model
Actually, although bound water diffusion along each fiber is an expected process, it is rather surprising that our systems made by piling short fibers could transport bound water throughout the sample. This necessarily means that bound water is able to "jump" from one fiber to another in contact with it. In addition, this means that these contacts were established during compression, and the introduced oil did not break them. This has important consequences in practice since it means that in any structure made by piling fibers such as in textiles, wood fiberboard, direct contacts will allow the transport of bound water throughout the structure, which may constitute a fundamental mechanism of sorption or desorption of the system.
Let us now analyze the diffusion process in the fiber network, in order to deduce the transport diffusion coefficient of bound water along the fiber axis. The knowledge of this diffusion coefficient would indeed constitute the key for the prediction of bound water transport and sorption or desorption dynamics of any type of cellulosic system as soon as its structure has been identified. The bound water molecules are transported (red arrows in Fig. 1a) along a fiber over some distance in the main flow direction. This distance is obviously larger as the fiber axis is closer to this flow direction. At some point they jump into another fiber in contact with the former (crossed blue circles in Fig. 1a, b), so as to continue their advance in the same direction. Due to the overall gradient of moisture content along the main flow direction, we expect that, at each point of contact between fibers the bound water will essentially go on moving along the fiber segment closest to the flow direction. The succession of such segments from the sample bottom to the sample top finally forms a path. Note that the above observation that water diffusion through oil is much slower than the overall extraction of water from our sample implies that the bound water path will preferentially follow segments in direct contacts between cellulose fibers, rather than attempting to jump from one fiber to another nearby through a thin film of oil. From the point of view of these paths this network may be seen as a set of fiber segments in contact and of various lengths and orientations with respect to the main flow direction ( z ) (see Fig. 1b). The total length ( L effective ) (see Fig. 1b) of the molecule path essentially corresponds to the sum of lengths of N p steps of elementary length a i in these various orientations (defined by the angles i and i in the spherical coordinate system of polar axis z ) (see Fig. 8a,b). The projected path length along the flow direction may then be written as.
1<i<N p a i cos i , so that the averaged path length writes in which P is the total number of paths. L direct is in fact equal to the sample thickness. Assuming a large number of paths and an isotropic medium so that all orientations are equally represented and independent of the step lengths (i.e. fully disordered state), this expression may be transformed as Note that this approach, through the various fiber orientations and steps, should remain valid for slight fiber curvature and as long as the length of the steps associated with jumps between fibers in all directions are sufficiently small. In this expression we recognize , the average effective length of the molecule path, i.e., through the tortuous structure, and we finally have ⟨L direct ⟩ = 1 2 � L effective � . From this relationship we deduce the tortuosity of the medium T = ⟨Leffective⟩ ⟨L direct ⟩ = 2 and the diffusivity in the direction z may be written D z = D L ∕ T , in which D L is the transport diffusion coefficient along the fiber axis. (3) This relationship describes the tortuosity in the case of a random isotropic medium. When a uniaxial compression is imposed to the medium so as to reduce the sample thickness along the z direction by a factor k , the distances along this axis are necessarily reduced on average by a factor k . This implies that on average the fiber axis will be bent so that their apparent height along z be reduced by a factor k (final to initial thickness ratio) (see Fig. 8a,b). This simply means that in the expression for the direct length ⟨L direct ⟩ , i.e., Eq. (4), each cos term must be replaced by cos � = k cos , while the total length for all the steps along the fiber axes is unchanged. Thus, we now have ⟨L direct ⟩ = k 2 � L effective � so that T = 2∕ k , and the diffusivity is D = k 2 D L . On the other hand, we can establish the relation between the factor of compression and the porosity. We assume that in the initial state, i.e., from gentle pouring and/or mixing, the fibers are in fully disordered state. The initial porosity 0 in that case is about 91%. Due to compression along a given axis we reach a porosity such that k = (1 − 0 ) (1 − ) . Then the expected diffusivity for a compression along the evaporation axis is For a compression perpendicular to the evaporation axis, the fiber vector orientation ( , ) changes to ( ′ , ′ ) as a result of compression while its tip moves along the compression axis. The distance to the central point along the compression axis is sin cos , thus sin � cos � = k sin cos . On the other hand, from geometrical considerations, the motion of the vector tip Fig. 8 Changes of orientation of a fiber, from its initial (dark green) to final position (light green) in a spherical coordinate frame resulting from a compression parallel (a) or perpendicular (b) to the evaporation axis ( z) in the plane of the two dotted red lines implies that cos ∕ sin sin = cos � ∕ sin � sin � . We deduce cos ′ = cos 1−k 2 sin 2 cos 2 1−sin 2 cos 2 1∕ 2 , so that (k) ranges from 1 to 1.27, when k ranges from 1 to 0.09 (theoretical minimum value), and the diffusivity coefficient is The predictions of this model can be compared with the values of the diffusivity deduced from the NMR tests. They appear in very good agreement, for data in the different directions of evaporation / compression and for porosity in a wide range (see Fig. 7). In this approach, the only fitting parameter of the model is the transport diffusion coefficient of bound water in the longitudinal direction of the cellulose fibers. Here a global model/ data agreement was obtained for D L = 3 × 10 −9 m 2 s −1 . Thus, the longitudinal transport diffusion coefficient of bound water between cellulose microfibrils appears to be close to the self-diffusion coefficient of (bulk) water (i.e. D L = 2.3 × 10 −9 m 2 s −1 ).
This value of the (transport) diffusion coefficient along the fiber axis obtained from this complete analysis of transport and fiber stack structure constitutes a new result. The rather large value obtained here, of magnitude similar to that of the self-diffusion coefficient of water, suggests that there exists a significant fraction of bound water molecules exhibiting high mobility. The constancy of this coefficient over a wide range of saturation values suggests that this high mobility concerns most of the water molecules absorbed in the nanostructure of the cellulose. The large value of the transport diffusion coefficient of water along the cellulose fiber axis obtained here is rather surprising, but as far as we know no data of experiments or simulations for bound water in cellulose exist in literature which could be discussed for comparison. Instead, various works (Gubbins et al. 2011;Kärger and Ruthven 2016;Jobic et al. 1999;Falk et al. 2015) suggest that for other liquids (k) = ∬ 0< < ∕ 2 0< <2 1−k 2 sin 2 cos 2 1−sin 2 cos 2 1∕ 2 cos sin d d cos sin d d in confined systems the self-and transport diffusion coefficients can be significantly different, by up to 2 orders of magnitude in some cases, often (but not always (Han et al. 2015)) in favor of a larger transport diffusion coefficient, but one must also consider the limits of the molecular simulations (Tsimpanogiannis et al. 2018) used for such predictions. Actually, Chmelik et al. (2010) considers that one should expect in general a larger transport diffusion coefficient as the finiteness of pore space accelerates transport diffusion in comparison with self-diffusion. To date, this effect of confinement is confirmed by all microscopic measurements of transport diffusivities in nanoporous solids which are found to be larger than or, at least, equal to the self-diffusivities (citation Chmelik et al. (2010)). At last note that the elasticity of the medium, i.e. the cellulose molecules, could play a role in the transport by diffusion, as was suggested for the self-diffusion (Zhang et al. 2023). On another note, anomalously (or surprisingly) fast "hydrodynamic" transport of water in nanopores is well documented in literature (Bocquet and Charlaix 2009;Thomas et al. 2014), a phenomenon which might be related to the effect of diffusion enhancement due to hydrodynamic fluctuations (Detcheverry and Bocquet 2012), but the theoretical frame differs from that of diffusion. These different observations do not provide straightforward explanations of the physical origin of the large transport diffusion coefficient of bound water in cellulose, but they show that a value larger than the self-diffusion coefficient of (bulk) water of much larger than the self-diffusion coefficient of bound water in cellulose is not fundamentally in contradiction with theoretical approaches.

Conclusion
From test preventing vapor transport we have been able to directly observe bound water transport by diffusion inside cellulose fibers, which made it possible to quantify unequivocally the transport diffusion coefficient of bound water. This made it possible to show a new effect: bound water is transported along a network of fibers in contact, even when immersed in a liquid. This effect has important practical consequences, since it shows that in various materials such as paper, textiles, wood, bio-based insulating materials used in construction, in which the cellulose fibers are in direct contact the sorption and desorption processes may rely on a transport of bound water throughout the fiber network. We could observe that this effect starts occurring once critical compression is produced (here associated with a porosity of 0.8) to get a sufficiently rigid sample, i.e., which can be handled without significant deformation. This suggests that a sufficient compression places the fibers in contact with each other, thus inducing sufficient bonding between fibers in contact, which gives its rigidity to the network. These close contacts will then also allow for the transport of bound water from one fiber to another in contact with it. In this context, note that changing the surface properties of the cellulose might induce changes in the diffusion throughout the network, if this can affect the characteristics of the contacts between fibers.