Lattice Boltzmann simulations of diffusion in steam-exploded wood

Diffusion of large molecules throughout the porous microstructure of wood pretreated with steam explosion was investigated by using the lattice Boltzmann method for simulations. Wood samples were investigated with high-resolution X-ray tomography to effectively reconstruct an accurate geometry of the structural changes that ensue after pretreatment. Samples of approximately 1 mm3 with voxel sizes from 0.5 to 1 μm were examined with X-ray imaging. These large volumes, relative to what reasonably can be simulated, were divided into sub-volumes and were further reconstructed into geometries suited for the LBM simulations. The transient development of the concentration was investigated, and the effective diffusion coefficient at steady state was computed. Diffusion rates were found to increase significantly in the transversal direction due to the steam explosion pretreatment. The increase was observed both in the time needed for solutes to diffuse throughout the pores and in the effective diffusion coefficient. A shorter diffusion pathway and a higher connectivity between pores were found for the pretreated samples, even though the porosity was similar and the pore size distribution narrower than the native sample. These results show that local mass transport depends not only on porosity but also, in a complex manner, on pore structure. Thus, a more detailed analysis of pore space structure using tomography data, in combination with simulations, enables a more general understanding of the diffusional process.


Introduction
Diffusion of large molecules in the porous microstructure of wood is important for understanding the limiting steps of an efficient material-driven biorefinery. In a material biorefinery (Mattsson et al. 2017), it is desirable to obtain high molecular weight fractions of the three main constituents of the wood cell wall, namely cellulose, hemicellulose, and lignin, that are used as precursors for value-added bioproducts (Gomes et al. 2014). To separate these constituents from the recalcitrant cell wall while retaining high molecular weight requires milder process conditions in the upstream pretreatment steps to reduce hydrolysis and still facilitate accessibility. Cell wall separation may be further enhanced by the use of enzymatic hydrolysis of insoluble polysaccharides (Hasunuma et al. 2013). To increase the yield of an enzyme treatment, it is necessary to open the cell wall structure with pretreatment of the wood (Azhar et al. 2011). The most common physicochemical pretreatment method for wood is steam explosion (Alvira et al. 2010). In steam explosion, pressurized steam is introduced to wood chips for a short duration of time, after which the pressure is rapidly decreased to atmospheric conditions, allowing the steam to expand inside the wood and a pressure gradient to develop between the wood and its surroundings.
The porous microstructure of steam-exploded wood has been studied with the use of scanning electron microscopy (Donaldson et al. 1988;Zhang and Cai 2006;Muzamal et al. 2015). This technique gives information on the effects of the pretreatment on the surface of the wood sample. An important tool for observing the internal structure of many porous materials is X-ray tomography where subsequent 2D slices are scanned in the material and reconstructed into a 3D volume. A general study on the applicability of synchrotron radiation tomography on wood using different types of wood samples has been performed by Mannes et al. (2010). In a more specific study, tomography was used to study the decay rate of wood caused by fungi on the micron scale (Sedighi Gilani et al. 2014). The data obtained were used to investigate the 3D structure and the change in density of samples incubated with fungi over time. The 3D structure of steam-exploded wood has previously been studied by Muzamal et al. (2016) using X-ray computed tomography. A comparison of the 3D microstructure of a steam-exploded wood sample with an untreated wood sample revealed that several cracks were created during the explosion step, and the microstructure of the wood sample was vigorously destroyed. In the present work, these images were used as input for numerical lattice Boltzmann method (LBM) simulations of diffusion.
LBM has, in recent years, become increasingly utilized as a simulation platform when tomographic data are used to reconstruct the geometry. This is due to the ease of implementation of boundary conditions in a heterogeneous sample. To assess how pore structure affects the mass transport through polymer film coatings intended for controlled drug release, Gebäck et al. (2015) have combined tomography data obtained using confocal microscopy and numeric simulations. The film was reconstructed from computed tomography images, and LBM was used to simulate diffusion through the film. Using this method, the effective diffusion coefficient that mainly controls the average drug release rate was found. Permeability is another important property of porous materials and can numerically be obtained for specific materials by simulating laminar flow in real geometries obtained from tomography data (Eshghinejadfard et al. 2016). Several studies (Guo et al. 2016;Liu and Wu 2016;Xia 2016;Zhang et al. 2016a, b) have used LBM for simulating various transport phenomena combined with a variety of tomographic imaging techniques. Thus, LBM shows great promise as a simulation platform for modeling pore-scale transport phenomena in heterogeneous structures.
Models of pure diffusion in wood are limited to older studies (Cady and Williams 1934;Burr and Stamm 1947;Christensen and Williams 1951;Behr et al. 1952;Stamm 1953;Petty 1973) in which simplified models have been used to estimate the impact of micro-scale structural elements in the porous structure of wood. Other studies regarding transport phenomena in general have focused on the moisture content of wood. Moisture content has an important role for the physical properties of wood (Siau 1984) as well as for drying wood (Fyhr and Rasmuson 1996).
The effect of steam explosion on mass transport rates has only been indirectly studied through the use of enzyme hydrolysis (Muzamal et al. 2016), which has shown a significant increase in activity caused by pretreatment. There is a lack of knowledge in the area of how steam explosion affects diffusion in the porous structure of wood.
In the present study, diffusion in the porous structure of steam-exploded and native wood was simulated using LBM in geometries reconstructed with the use of high-resolution X-ray tomography images. The effect of the structural modifications caused by the pretreatment method on mass transport was investigated and compared to the native wood structure. Local variations within the samples are inherent due to the nature of wood, and, for this reason, diffusion was investigated in similar localized volumes. This was done by dividing the 1 mm 3 samples into smaller subvolumes, as the resolution required to resolve the whole volume would require massive computing power.

Materials and methods
Wood samples and preparation of samples for high-resolution X-ray tomography were performed in Muzamal et al. (2016). The subsequent images obtained from the X-ray tomography were used in this paper. Wood samples were prepared with the dimensions of 20 × 20 × 4 mm 3 (longitudinal × tangential × radial) from a stem of Norway spruce. All samples were cut from mature wood and the same annual ring. One sample was saved as a reference, while the others were used in the steam explosion treatment. The pretreatment was carried out at 14 bar with saturated steam for 10 min after which the pressure was released to atmospheric conditions by opening a release valve rapidly. Upon release, the samples were forced down into a so-called blow tank where pretreated samples are collected. Wood chips collided with vessel walls during the release, causing a strong impact which resulted in disintegration of various degrees (Muzamal and Rasmuson 2017). The pretreatment created large structural changes of different size fragments, some completely defibrillated, while others remained relatively intact. To investigate the difference between these fragments, two samples were selected for X-ray tomography. Sample SE-1 was chosen with relatively little damage compared to sample SE-2 which sustained large damage. These chosen samples were further cut into small cubes, approximately 1.4 mm, to obtain X-ray scans at a high resolution. The same procedure was done with the native reference sample. Mainly earlywood was sampled because of the larger size of both lumina and bordered pits compared to latewood.
High-resolution X-ray tomography was performed at the Ångström Laboratory of Uppsala University (Uppsala, Sweden). The projections obtained from the X-ray scan were reconstructed into stacks of 2D 16-bit BMP images using the software nRecon 1.6.10.1 (Bruker, Sweden). The stacks contained up to 1700 images of cross sections perpendicular to the longitudinal tracheids, as seen in Fig. 1. The native and SE-1 sample had an isotropic voxel size of 0.81 μm. The SE-2 sample had a slightly higher resolution and an isotropic voxel size of 0.54 μm. The images were cropped and converted into 3D TIFF stack of 200 images where each image slice represents a thickness of 0.81 and 0.54 μm, respectively.

Fig. 1
Image obtained from the X-ray data from the native wood. The highlighted area shows the image in full quality before any image processing. The size of the cropped image is 393 × 361 µm

Region of interest and sub-volumes
To use the entire domain of the samples analyzed for simulation would be immensely computationally heavy, and, for this reason, smaller sub-volumes were investigated. To avoid any surface effects caused by sample preparation, only the center of the volume (if possible) was considered, as illustrated in Fig. 1.
The reconstructed 3D volume and a slice of the triangulated surface of the cell wall for the native wood sample can be seen in Fig. 2. Only one sub-volume was used for the native sample since its characteristics were relatively constant over the volume with only minor local variations. The pit openings range in size from 3 µm up to 15 µm. The thread-like margo and the torus inside the chamber of the pit could not be resolved, which is not surprising since they are both below 0.5 µm in thickness. A general average value for earlywood pits in softwoods is a pit opening diameter (the aperture) of 3.7 µm with an inner diameter (the outer diameter of the margo) of 16.1 µm (Kvist et al. 2017). The size would indicate that the X-ray tomography has difficulties to resolve the small chamber depth, approximately 2 µm, of the pit.
The first steam-exploded sample (SE-1) is seen in Fig. 3. Two different regions of interest were chosen (SV-1 and SV-2) to investigate local differences within the Fig. 2 Native reference wood sample from X-ray tomography. The top of the image shows the reconstructed 3D volume used to construct the surface for the geometry used for the simulations. The size of the 3D volume is: 191 × 188 × 162 μm. The bottom of the image shows a cross-sectional slice of the triangulated isosurface of two cell walls with pit openings sample. SV-1 was chosen as a region that was not noticeably affected by the treatment, and SV-2 showed some deformation and separation of the cell wall.
The second steam-exploded sample that sustained heavy damage (SE-2) during pretreatment is shown in Fig. 4. The sample was highly variable in all spatial dimensions. To choose comparable volumes in the irregular structure is difficult if the open pore space between the separated cell walls is too large. As such, the sub-volumes were selected to encompass as much of the cell wall structure as possible while still including some of the separation, deeper into the sample. The sub-volumes and their corresponding region of interest can be seen in Fig. 5.

Image processing and 3D reconstruction
It is necessary to process X-ray tomographic images before reconstructing a suitable geometry for performing numerical simulations. The main programs used for processing were ImageJ (version 1.51j) (Schneider et al. 2012) and MATLAB (The Math-Works, Inc., Natick, MA). The main processing concerns were removing background noise and artifacts while preserving the integrity of the sample structure. A Gaussian filter with σ = 1 was used to reduce jagged edges and to obtain a smooth surface of the The size of SV-1 is: 136 × 111 × 108, SV-2: 139 × 120 × 108 μm, and SV-3: 148 × 104 × 108 μm geometry suitable for lattice Boltzmann simulations. A marching tetrahedral method (implemented in the GTS library, version 0.7.6, http://gts.sourc eforg e.net) was used to create a triangulated isosurface at a given threshold level. The isosurface is obtained in the form of an STL file.

Porosity, pore size distribution, and tortuosity
Porosity was taken as the available fluid volume fraction in the simulation box. The simulation box spans the entire generated 3D geometry based on the STL file, defining the boundaries of the cell wall, and thus, the available voxels in the open pore space are known. Pore size distribution was calculated by fitting spheres in the image stack from the X-ray data with the Thickness plugin in the BoneJ (version 1.4.2) (Doube et al. 2010) extension for ImageJ. Multiple spheres were fitted into the open pore space until all voxels are covered. The smallest sphere needed was 2 µm. Gamma probability distributions were fitted to the histogram of sphere diameters for the entire image stack. An empirical relation to quantify how a porous material retards diffusion is tortuosity and is defined in two ways (Shen and Chen 2007) where the first is: scaling the free diffusion coefficient D 0 (m 2 /s) with porosity and tortuosity to find the effective diffusion coefficient. The other class defines the tortuosity as what it physically represents, as the ratio of the distance traversed ( ΔL ) by a solute per the length of the medium ( Δx): In the present work, the second definition was used. Using the obtained diffusive flux field from the LBM simulation, the distance traversed can be calculated as streamlines of the flux vector field. To obtain an average value for the distance, diffusive flux lines were randomly sampled throughout the simulation domain for each sample. Calculations of tortuosity were performed in MATLAB.

Lattice Boltzmann method
The lattice Boltzmann method (LBM) is a numerical method used to solve partial differential equations. The computer code used was based on the Palabos open-source code and implemented in C++. Computations were performed on an Intel Core i7-8700 64-bit workstation with 32 GB of memory. The D3Q19 model was used in the simulations with a two-relaxation-time collision model (Ginzburg 2005). Using a multiple relaxation model for the collision operator increases stability while being equal in terms of computational time and simplicity compared to single relaxation approximations. In this study, LBM was used to solve the diffusion equation: where C (mol/m 3 ) is the concentration and D 0 (m 2 /s) is the free diffusion coefficient in the open pore space. A zero normal flux boundary condition was implemented at all walls in the wood cell structure, i.e., the cell wall is considered to be impermeable. The boundary was found based on the reconstructed 3D geometry from the tomography images. The boundary condition was implemented with the use of ghost nodes, according to the methodology presented by Gebäck and Heintz (2014). Compared to the classical bounce-back rule applied to boundaries, the ghost nodes produced an accurate tangential flux near the surface by interpolating macroscopic values with a mirror point inside the domain based on the normal direction of the surface geometry. The zero flux condition was then applied to curved boundaries according to: where n is the outward unit normal. With an applied concentration difference, the diffusion equation was solved until steady-state conditions were reached. The initial concentration was specified as zero in the simulation domain to investigate the progression of the concentration front over time. If steady-state results were desired for the calculation of the effective diffusion coefficient, the initial concentration was specified as linear throughout the domain, which sped up simulation time significantly. The effective diffusion coefficient, D eff (m 2 /s), was then obtained from the average flux according to: where Δx is the thickness of the simulation box, ΔC is the applied concentration difference, and J x (mol/s m 2 ) is the average flux in the direction of the concentration gradient. For time-dependent results, data were collected at specified time intervals for a specified number of total iterations. Simulations were performed on a uniform grid in the range of 300 × 300 × 200 voxels, depending on the size of the geometry.

Results and discussion
Tortuosity and porosity for the corresponding pore size distribution of the native sample can be seen in Figs. 6 and 7. The choice of sub-volume region and the size chosen are of importance and should be consistent between samples to obtain comparable values for porosity and tortuosity. The native samples had the highest porosity of all the samples. Pore sizes are seen to decrease for the pretreated samples, especially for SE-2, as the choice of sub-volume for SE-2 mainly encompasses the compacted microstructure due to collisions of the wood after pressure release. If the entire samples were analyzed, the large cracks and ruptures of the cell walls would increase pore size significantly. In a study by Zauer et al. (2014), they observed a decrease in pore size for thermally treated wood. Similarly, Van den Bulcke et al. (2013) showed that aspen wood thermally treated at 160 °C shrinks, but the overall microstructure remains alike. In general, the pore size presented in published data is larger. The average value for an earlywood tracheid is in the range 33-39 μm in the radial direction and 31-33 μm in the tangential direction for Norway spruce (Brändström 2001;Havimo et al. 2008). The average presented here for native wood is slightly lower than 20 μm. One explanation for this discrepancy is due to the shape of the cells and that the algorithm used to estimate pore size is based on fitting spheres inside the open pore space. This will underestimate the size because several spheres will be required to cover the entire space of a lumen and as such will produce lower diameters. Although the porosity reported here is lower compared to published values close to 0.7 (Plötze and Niemz 2010;Zauer et al. 2013), it has to be taken into account that this is a local volume compared to a macroscopic average. Plötze and Niemz (2010) also reported the large amount of nano-pores present in spruce wood that cannot be resolved here, although it is a smaller volume of pores. The general trend observed is that for native wood we have a larger porosity and tortuosity, indicating that although more pore space is available for diffusion, the path is generally longer.
Another effect lowering both pore size distribution and porosity is the delineation of the cell wall. Reported average cell wall thickness values of Norway spruce are 3.52 µm and 2.90 µm in radial and tangential direction, respectively (Brändström 2002), while even values as low as 2.1-2.2 µm are reported by Havimo et al. (2008) for cell wall thickness in tangential direction. Measurements on micrographs taken with a confocal microscope of a native wood sample yielded a cell wall thickness of 3-4 μm. Measurements on the tomography scans of the native samples result in an average cell wall thickness of 4.80-6.48 μm for radial and tangential direction, respectively. The high value in tangential direction is due to ray cells being resolved as cell wall. This in combination with smoothing by Gaussian filters overestimates the cell wall thickness. One may use filters that preserve edges, such as bilateral or anisotropic diffusion, to avoid loss in sharpness. However, smoothing was necessary for both noise removal and to obtain a smooth surface such that the boundary condition could be applied accurately.
The only noticeable effect of the pretreatment on the SE-1 sub-volumes is the lower tortuosity for SV-2, seen in Fig. 7. The lower tortuosity for SV-2 seemingly has more of an influence than porosity on the effective diffusion coefficient, as the porosity is lower as well for SV-2. In a previous study by Muzamal et al. (2015), an increase in small pores in the range of 1-5 µm for steam-exploded wood was observed using mercury porosimetry. The increase in smaller pore sizes and the effect of thermal treatment may explain the decrease in porosity for the steamexploded samples, especially since it is a local volume and thus not representative of the entire sample obtained from steam explosion.
Tortuosity was found to drop significantly for the steam-exploded samples, which sustained heavy damage from the pretreatment, as shown in Fig. 7, while porosity for SE-2 is harder to judge in a comparative way because some of the regions have open channels wherein the fibers have completely split and the cell wall structure has been compressed due to collisions. As such, the sub-volumes were chosen to mainly represent the cell wall structure which is also seen in similar pore size distributions. A general trend of increasing diffusion coefficients and lower tortuosity was found for all the SE-2 sub-volumes.
The simulations were solved until steady state was reached after which the effective diffusion coefficient was calculated. The steady-state results are presented as the ratio of the effective-to the free diffusion coefficient. This yields a geometry factor, i.e., a measure of how the structure retards diffusion. For specific molecules, the results can be scaled by the free diffusion coefficient (coefficient in solvent without obstructions) of mainly large molecules not diffusing through the cell wall. The evolution of the concentration front was studied over time to obtain information on diffusion pathways and how long it takes for diffusing solutes to traverse the porous structure being studied. For the transient study, the free diffusion coefficient of 40 kDa FITC-dextran was used at 42.6 (µm 2 /s), calculated based on fluorescence recovery after photobleaching (FRAP) measurements in aqueous solution (Kvist et al. 2018).
Diffusion occurs in three directions in wood, longitudinal, radial, and tangential. Tangential diffusion of ions has been reported (Christensen 1951;Ra et al. 2001) to be about 2-4 times lower compared to that of radial diffusion. Although ions diffuse through the cell wall, the difference observed between radial and tangential diffusion is most likely because of pits being more common in the radial direction, as well as the occurrence of ray cells. Due to difficulty to properly resolve the ray cells as well as the variability of the bordered pits based on the tomography data, it was assumed that radial and tangential diffusivities were similar in this case. As such, only transversal diffusion will be considered here. A similar assumption was made by Kazi et al. (1996) in a combined reaction-diffusion model where experimental results showed fairly good agreement.
The concentration front for the native wood moved significantly slower than the front for wood pretreated with steam explosion, as illustrated in Figs. 8 and 9. In the native wood, diffusion between tracheids was through bordered pits, and Time evolution of the concentration front for a steam-exploded sample. The increase in connectivity affects the path of diffusion and, thus, the time it takes for the concentration to progress over the sample can clearly be seen the transport depended heavily on connectivity. Additional connectivity was seen for steam-exploded samples due to the degree of disintegration. If a sample that has undergone steam explosion in any way ruptures and opens up, an increase in the connectivity between tracheids will be seen, and diffusion will significantly increase.
The transversal diffusion in different regions in the native sample behaved similar to what is shown in Fig. 8. The concentration front moved slowly from tracheid to tracheid through the bordered pit openings. A similar result was observed for SE-1. However, in certain regions of the SE-1 sample, where the connectivity was greater, the front progressed faster and spread out more quickly.
Observations of the concentration over time in the SE-2 samples reveal an increase in the diffusion rate throughout the structure, as was expected. Though pore sizes can be somewhat small due to the structure being compressed by collisions, the increase in connectivity was significant, which leads to a much smoother concentration profile in the exploded sample compared to the native.
Diffusion in wood is highly direction-dependent, and in the longitudinal direction, along the fibers, diffusion takes place many times faster. In the transversal direction, diffusion occurs from lumen to lumen via bordered pits. The pit openings are significantly smaller than the lumina and act as a large obstacle to diffusion, which is clearly seen in Fig. 10. If a pit is resolved in the geometry, it is completely open, i.e., no pit aspiration. Due to the small size of the torus, it is unlikely that aspiration of pits can be properly considered with this method. Due to this and the fact that the aperture of the pits is of variable size, up to that of the inner diameter of the margo, the diffusion coefficient should be somewhat inflated. Another effect may be filtering of large molecules by the fine pores in the margo. These effects are not considered in the present model.
In a study on diffusion of polyethylene glycol in native Sitka spruce, Fukuyama and Hiroyuki (1986)  In earlier work performed by the authors (Kvist et al. 2018) on steam-exploded wood, FRAP (fluorescent recovery after photobleaching) was used to determine local diffusion coefficients of dextran. The measurements were taken over a single tracheid in the transversal direction, and values of D eff /D 0 for a 3 kDa FITCdextran probe were for native: 0.067, steam-treated: 0.20, and steam-exploded: 0.33. Although the pretreated samples are vastly different, it highly depends on the local tracheid chosen and how it was affected. However, the values for the native sample are of the same order of magnitude as in the simulations.
The release of pressure during steam explosion causes an expansion of vapor within the structure of a sample and is known to increase porosity (Muzamal et al. 2015) in pores having diameters of less than two micrometers. Due to the limitations of tomography imaging, it is hard to capture these small changes, which is reflected in the similar results for the native and SE-1 samples in Fig. 11. It is difficult to see a micro-scale difference in diffusion due to this; however, for nanoscale phenomena, such as enzyme binding and penetration inside the cell wall, this is of great importance as has been shown in Azhar et al. (2011), Jedvert et al. (2012 and Muzamal et al. 2016).
For the impacted sample (SE-2), the wood disintegrated and the structural changes were on a scale that could be captured with X-ray tomography. The wood disintegration was found to have quite a large influence on transversal diffusion. This was highly localized in the samples and heavily dependent on the choice of region of interest. As such, these results are not representative for the entire steam-exploded samples, but rather give indications of the local diffusion coefficients.  Fig. 11 The ratio of effective diffusivity over free diffusivity for different regions of interest for native wood and two steam-exploded wood samples. The native sample only has one sub-volume (SV-1), and the steam-exploded samples (SE-1 and SE-2) have two and three sub-volumes (SV-1 to SV-3), respectively

Conclusion
In a material-driven biorefinery setting, the mass transport of large molecules, such as enzymes and polysaccharides, in wood structure is of fundamental importance for understanding limiting steps of future biorefineries. To this end, simulations of diffusion were performed on wood samples pretreated with steam explosion. The wood structure was reconstructed based on X-ray tomography images with a resolution down to 0.54 µm. Information obtained from the simulations gives insight into the connection between structure and diffusion pathway, and this information can easily be scaled to find the effective diffusion coefficient based on the diffusion coefficient in a free solution.
It was found that transversal diffusion rates increased after the pretreatment method. The increase in diffusion rates varied locally depending on the degree of disintegration of the cell wall structure and the local volume investigated. A lower tortuosity, or shorter diffusion pathway, was found for the pretreated samples. Despite that pore size distribution was low for these samples due to compaction during pretreatment, diffusion rates nevertheless increased.
X-ray tomography was limited in that sub-micrometer features could not be resolved. However, several studies have shown that steam explosion effects were found to occur in the cell wall on a much smaller scale, for example leading to an increase in enzymatic hydrolysis. Future studies should consider cell wall diffusion on the nanometer scale and the influence of pretreatment as the next stage in understanding mass transport phenomena in wood.