A multi-scale model for diffusion of large molecules in steam-exploded wood

In this paper, multi-scale modeling was used to resolve diffusion in steam-exploded wood at tracheid scales including sub-micrometer features of bordered pits. Simulations were performed using the lattice Boltzmann method with high-resolution X-ray tomography image data as the input for the microstructure. The results show an effective method for utilizing a variable diffusion coefficient to implement two length scales. This circumvents the need to resolve the bordered pits in detail, which requires massive computing power. Instead, the effective diffusion coefficient for one bordered pit is used as input for this model. Results based on the present model are comparable to experimental data. This methodology can be extended to more structural features at the microscale of wood, such as latewood and the cell wall. Obtaining a map of different diffusion coefficients based on features and scale gives a better overall understanding of diffusion and the importance of mass transport with regard to the pretreatment of wood.


Introduction
A challenge for society today is to revert from and find alternatives to the fossilbased industry and aim for sustainable solutions. At the center of this is the biorefinery based on the developing bio-based economy with biomass as a renewable raw material. The utilization of biomass, especially wood, as a renewable feedstock for bio-energy and value-added bio-products is a central step in reducing dependency on fossil fuels (Gomes et al. 2014). In a material biorefinery (Mattsson et al. 2017), it is desirable to obtain high molecular weight fractions of the main constituents of the wood cell wall: cellulose, hemicellulose, and lignin. The larger fractions have inherent valorization; however, to obtain these fractions from the highly recalcitrant cell wall in wood is quite difficult. Often, pretreatment steps are required to make the cell wall more accessible for enzymes (Azhar et al. 2011); however, the treatment must not be too harsh in order for the biopolymers not to degrade.
Softwoods, which are the most common type of lignocellulosic material in the northern hemisphere (Galbe and Zacchi 2002), are the focus of the present paper. The microstructure of softwoods is built up by cells known as tracheids, parenchyma, and epithelial cells. Longitudinal tracheids make up approximately 90% (Siau 1984) of the total volume of cells. The main function of earlywood tracheids is the conduction of fluids, while the function of latewood tracheids is structural support, which is seen by the much thicker cell walls in latewood. The main transport of fluid occurs between tracheids through bordered pits. Inside the chamber of a bordered pit is a centered impermeable torus surrounded by a membrane comprised of thin strands of cellulose microfibrils (Sjöström 1993). The membrane and part of the torus are covered by a cell wall that arches over the chamber with a small circular aperture leading into the chamber. The mass transport, especially for larger molecules, such as enzymes, is dependent on the pathway through the tracheids and the interconnected bordered pits. This, as well as the cell wall, gives rise to the main mass transport limitations in the microstructure of wood and is often alleviated to varying degrees by employing a pretreatment step to make the wood cell wall more accessible (Agbor et al. 2011).
The most common physicochemical pretreatment step is steam explosion (Alvira et al. 2010). In steam explosion, biomass is treated with high-pressure steam at temperatures ranging from 160 to 240 °C for residence times up to 10 min (Agbor et al. 2011). The pressure is quickly released to atmospheric conditions, and the pressure difference causes an "explosion" that can alter the microstructure of wood. Understanding how the microstructure is altered by steam explosion is important and has been studied by means of scanning electron microscopy (Donaldson et al. 1988;Zhang and Cai 2006;Muzamal et al. 2015). Further investigation of the internal structure of steam-exploded wood was performed by Muzamal et al. (2016), who used X-ray tomography. In that study, 3D data were visualized, and the degree of treatment was examined with subsequent enzyme activity experiments. The experiments showed that activity increased significantly for samples that had undergone steam explosion. Cell wall diffusion was investigated (Paës et al. 2017) for several other pretreatment methods by means of fluorescent recovery after photobleaching (FRAP) to assess accessibility. However, there is a lack of studies investigating how the change in structure impacts diffusion at the microscale in steam-exploded wood.
In the present study, the dependence of the diffusion coefficient on multi-scale implementation of bordered pits in the microstructure of wood was investigated by utilizing lattice Boltzmann simulations. The lattice Boltzmann method (LBM) has become increasingly popular for complex heterogeneous structures due to both the ease of implementing boundary conditions and its parallelization capabilities (Zhou 2004;Bernsdorf 2008;Gebäck et al. 2015). The study also aims to develop a simple method to implement two length scales, the microscale of tracheids and the submicrometer features of the bordered pit. For mild steam explosion, effects can be observed at several scales and can occur simultaneously. The main effect that causes the structure to break in steam explosion is the collision of the wood with walls in the equipment, and the type and amount of damage caused seem to depend highly on velocity and the impact angle of the wood (Muzamal and Rasmuson 2017). However, a significant part of the wood remains intact and is mainly affected by the expansion of water vapor during the release of pressure. This pressure release has been shown to increase the porosity of wood in the pore size range of 1-4 µm (Muzamal et al. 2015), which corresponds to the size range of bordered pits. Cracks in the cell wall close to pits have been seen using steam explosion (Zhang and Cai 2006;Muzamal et al. 2015) and micro-explosion (Ma et al. 2016), where high-pressure air is used in place of steam.
The microscale in wood was resolved, and diffusion was simulated in a previous work by the authors (Kvist et al. 2019). However, it was difficult to resolve the small-scale features of the bordered pit. It would not be computationally feasible to solve the diffusion equation with such a large difference in length scale. In the present work, two different length scales were utilized by using previous simulation results by the authors (Kvist et al. 2017) of diffusion through one bordered pit, as input for the local diffusion coefficient of the pits. The pits were found in the 3D structure generated by X-ray tomography images using a watershed segmentation algorithm (Meyer 1994).
It should be noted that in this paper only diffusion in the microstructure is considered, although other phenomena as adsorption/desorption phenomena and enzyme activity are relevant factors to fully understand the effect of pretreatment severity and the action of lignocellulolytic enzymes.

Modeling
To model the microstructure of wood, raw data from X-ray tomography experiments performed in a previous study (Muzamal et al. 2016) were used. The experiments were performed at the Ångström Laboratory at Uppsala University, Sweden. Wood samples of Norway spruce (Picea abies) were used and were prepared to a cubical size of 1.4 mm 3 . The X-ray detector used was an 11-megapixel 12-bit dynamic range-cooled charge-coupled CCD camera. It utilized an X-ray source voltage of 31 kV and a source current of 187 µA. The exposure time was 4000 ms, and the sample was rotated 192° in the X-ray beam at increments of 0.2°, yielding a resolution with an isotropic voxel size of 0.81 µm. A possible improvement in the images may be to adjust the threshold.
The following sections describe the use of the subsequent images obtained from the X-ray tomography projection data. A smaller volume of interest was found suitable for study. In this paper, only interconnected tracheids in earlywood were considered. Features, such as rays in the radial direction, were not considered. The diffusion equation was solved for the generated wood microstructure using the lattice Boltzmann method (LBM), and the effective diffusion coefficient was computed. A variable diffusion coefficient was used in the simulation to differentiate the bordered pits from the rest of the open space in the lumina. The diffusion coefficient inside the space of the pits was based on a previous study by the authors (Kvist et al. 2017), where the effect of steam explosion on a bordered pit was investigated. Figure 1 shows a cross section and a volume of the images from the tomography data. Ray tracheids can be seen between every fourth to fifth longitudinal tracheid. However, the rays are not resolved well in 3D and appear more as noisy abruptions. Gaussian smoothing was applied (σ = 1) to reduce background noise and create a smooth surface. All image post-processing was performed in ImageJ (version 1.51j) (Schneider et al. 2012) andMATLAB (2017, MathWorks, Inc.).

Wood microstructure
After post-processing, the images were converted into 3D TIFF stacks of 200 images, as shown in Fig. 1b. An 8-bit TIFF stack was used as input data for the surface generation of the geometry based on the threshold of the grayscale in the images (0 to 255). This determines what is considered a wall and what is open pore space. A region of interest covering 8 full tracheid pairs without ray cells was selected and generated, and the corresponding inner surface of the tracheid is shown in Fig. 2. This region was used as the base case for variation in number and size of pits. A triangulated isosurface was generated using a marching tetrahedra method (GTS library, version 0.7.6, http://gts.sourc eforg e.net).
Three more surfaces were created in the selected region of interest by increasing the threshold level of the image input data, thus increasing the number and size of the bordered pits. This enabled a comparison of the importance of porosity and the size of the pits for the simulated results, while keeping the microstructure of the tracheids relatively constant. An example of a cross section of the surface showing the pits is presented in Fig. 3. To compare the size of the pits, an equivalent diameter of a sphere with the same volume as each pit was calculated. The volume for each pit was compared with an average value based on published data (Petty 1972;Siau 1984;Hacke et al. 2004;Trtik et al. 2007;Schulte 2012;Schulte et al. 2015), similar to what was done in a previous study by the authors (Kvist et al. 2017). The volume of the pit was calculated by assuming a spheroid shape with an average margo diameter of 16.10 μm and a length from aperture to aperture of 5.27 μm. A general view and the nomenclature of the features of a bordered pit are shown in Fig. 4. The volume of the torus and margo was not considered for the computed average volume. The equivalent diameter of the pit was calculated to be 7.6 μm based on published data.
The equivalent diameter of the pits in the generated geometry was found with an algorithm used to locate bordered pits (the algorithm is described below). The total volume of the voxels that encompassed the pit was recorded along with the number of pits. The number of pits was taken as the total number of continuous individual connections between tracheids.

Lattice Boltzmann simulations and multi-scale implementation
The lattice Boltzmann method (LBM) was used to solve the diffusion equation until steady state with a variable diffusion coefficient D: where C is the concentration. D will vary in space depending on whether diffusion occurs in the lumen or the interconnected space of the bordered pit. A zero-flux condition, i.e., no diffusion through cell walls, was implemented at the boundaries of the generated surface according to Eq. (2): (1) where C is the concentration, D is the variable diffusion coefficient, and n is the outward normal. The boundary conditions were applied using a ghost-node scheme, according to the methodology by Gebäck and Heintz (2014). The scheme utilizes mirror points across a boundary to interpolate macroscopic variables to obtain the correct boundary condition. A concentration difference was applied in the radial direction to create a concentration gradient across the tracheids. A mirror boundary condition was applied in the tangential and longitudinal directions. After solving the diffusion equation to steady state, the effective diffusion coefficient was calculated based on the average flux in the direction of the concentration gradient: where Δx is the thickness of the simulation box and ΔC is the applied concentration difference. Simulations were performed on a uniform grid with 315 × 158 × 300 voxels on a D3Q19 lattice scheme and used a two-relaxation-time method for diffusion (Ginzburg 2005). All computations were performed on an Intel Core i7-8700 64-bit workstation with 32 GB of memory with simulation times up to 24 h. Diffusivities were supplied to the simulations through a VTK file based on the voxel data. An algorithm written in MATLAB was used based on watershed segmentation (Meyer 1994). The watershed transform is normally applied to separate objects that touch in images. In the present study, it was used to separate the connecting tracheids and, thus, effectively find the location of a bordered pit. A cross section of the voxel data was converted to a binary image, as shown in Fig. 5a. The Euclidean distance transform of the binary image was computed with the bwdist function in MATLAB and combined with the extended-minima transform, as shown in Fig. 5b. These two transforms are useful for objects that are circular in shape and touch, and it ensures that the watershed transform will not oversegment the image. Lastly, the watershed transform was applied, and all the tracheids were segmented, as shown in Fig. 5c.
The diffusion coefficient was then inserted into the voxels surrounding the area found using the algorithm, and as such, the total space occupied by the pit was obtained and varied, as illustrated in Fig. 6. After this had been implemented, diffusion simulations were run with the variable diffusion coefficients.
To compare the different cases and to represent steam-exploded wood, the ratios of the coefficients were set to 1:1, 1:2, 1:4, and 1:10. Thus, in the open pore space in the lumina, the coefficient will be the same, or 2 times, 4 times, and 10 times larger than that of the pit.

Results and discussion
The results are presented first with a look at the structural features of geometries, such as cell wall thickness, size of pits, and porosity. This is followed by results from the lattice Boltzmann simulations performed on said geometries based on the computed effective diffusion coefficient and the flux field. Grid-size Finding the specific location of connectivity between the larger pores by using a watershed algorithm. Depicted in a is a cross section of the larger data set, which was converted to binary. In (b), Euclidean distance and extended-minima transforms were applied to the binary image. Finally, the watershed transform was applied in (c) Fig. 6 Spatial locations of the variable diffusion coefficients supplied to the simulation. The cell wall (black) is zero, while the lumens (white) and bordered pits (gray) are supplied as ratio where the lowest is 1 (for the pit) resolution was investigated to ensure results with grid independence. The change in the solution with additional refinement was less than 0.2%.

Structure of the generated geometry
Four different representative geometries of the earlywood microstructure were simulated (Table 1). Case 1 is a base case related to the structure shown in Fig. 2, and Cases 2-4 are variations with higher porosity and number of pits. The cases compare different sizes and frequencies of pits and analyze how this in turn effects the effective diffusion coefficient. These four cases are then all simulated at different severities of the steam explosion given by 1:1, 1:2; 1:4, and 1:10. The different ratios represent different degrees of opening or rupturing of the bordered pit caused by steam explosion, where 1:1 would represent a bordered pit where the borders had been completely removed and 1:10 is only mildly affected. The lower limit for an undisturbed pit in untreated native wood is 1:25 (Kvist et al. 2017); however, such a large difference could not be simulated due to numerical stability issues in the simulations. The cell wall thickness in the radial and tangential directions was on average 4.6 μm with an increase in the cell corner due to higher intensity in the X-ray images. Published data on Norway spruce cell wall thickness range from 2.10 to 3.52 μm (Brändström 2001;Havimo et al. 2008). This discrepancy led to a somewhat lower porosity than what was expected for earlywood (about 0.7 for bulk including cell wall porosity (Plötze and Niemz 2010;Zauer et al. 2013)). However, increasing the image threshold further to decrease cell wall thickness would compromise the structure around the pits because of variations in intensity. Case 4 seems to be the limit for these image data for pushing porosity while maintaining the diameter of the pits within the calculated average for a softwood pit of 7.6 μm.
According to Sjöström (1993), the average length of a tracheid is 2-4 mm and it contains about 200 bordered pits. Assuming an average and not accounting for increased pit density toward the end of the tracheid (Brändström 2001), we can expect around 13 pits per 200 μm of length. The geometry in the present study was 162 μm long, covering about 12 tracheids, giving 6 pairs, and a reasonable number of pits of around 50. Both the size and number of bordered pits are comparable with the average found in the literature.

Free diffusion coefficient
The results from the LBM simulations were scaled with the free diffusion coefficient in aqueous solution without obstructions. Two different diffusion probes were used, taken from Kvist et al. (2017), and are shown in Table 2. The first probe was a fluorescently marked dextran with a molecular weight of 40 kDa measured using FRAP. The second was a protein, also with a molecular weight of 40 kDa. The protein was used in this study to simulate an enzyme of similar molecular weight that could be used in the biorefinery setting. The free diffusion coefficient for the protein was calculated by using the Stokes-Einstein equation. Even though the probes were of the same molecular weight, the coefficients were different because dextran is a longchain polymer compared to a protein of approximately spherical globular structure. This was also observed in the hydrodynamic radius of both probes.

Lattice Boltzmann simulations
The computed effective diffusion coefficients based on LBM simulations for the four cases at different severities of steam explosion are shown in Figs. 7 and 8. The variation based on geometry and multi-scale implementation of the variable diffusion coefficient for a bordered pit is shown on the x-and z-axes, respectively. The variable diffusion coefficient is represented by ratios, and in the first 1:1, the same diffusion coefficient was used for both the lumina and pits, i.e., no multiscale implementation. For the other three simulations, the open space within the pits was scaled with a ratio according to 1:2, 1:4, and 1:10. Effectively, this means that the coefficient in the lumina was 2, 4, and 10 times larger than that of the bordered pits. As the results have been scaled by the free diffusion coefficient, a higher effective coefficient for the protein is seen than for dextran. Although the molecular weights of the probes were the same, the configuration would be vastly different for a polysaccharide compared to the globular structure of an enzyme. As expected, when the number and size of pits were increased, the coefficient increased, as seen from Case 1 to Case 4. As the ratio between the diffusion coefficients increased, the effective diffusivity decreased, as expected. The differences between the cases also became smaller. Recall that the diffusivity ratio for native wood is 1:25.
Measurements taken with 40 kDa dextran diffusion probes using FRAP by Kvist et al. (2018) gave diffusion coefficients in the range of 7 μm 2 /s. These measurements were taken over one single earlywood tracheid. In another study by Fukuyama and Hiroyuki (1986), using polyethylene glycols (PEG) of varying molecular weight, a diffusion coefficient of 0.7 μm 2 /s was found for the largest PEG (1.5 kDa). These experiments were performed using a diffusion cell and a refractometer to measure concentration (Fukuyama and Hiroyuki 1980).
The predicted results in the present simulations were in the middle compared to both studies (Fukuyama and Hiroyuki 1986;Kvist et al. 2018). The reason for this is most likely due to the difference in the scale between the simulations and the experiments performed. As the scale increases, so does the tortuous path throughout the structure, especially when performed at the centimeter scale as with PEG. On the other hand, the path is as short as possible for measurements over one single Fig. 8 Effective diffusion coefficients for all geometry cases for a 40 kDa dextran where the free diffusion coefficient was determined using fluorescent recovery after photobleaching. The variable diffusion coefficient ratio is shown on the y-axis where the left most is the original without multi-scale implementation 1 3 tracheid. Observing the flux field from the simulations, one can study the general path through the structure, as shown in Fig. 9.
The increase in the number of pits, of course, increases the number of pathways the diffusing species traverse through the structure, as seen in Case 4 in Fig. 10. It is important to accurately portray the connections between tracheids to get a correct view of the mass transport in the wood microstructure. Even though the number of pits increased between the case geometries, the change in the effective diffusion coefficient was small, with a factor of 2. The change in the coefficient became even smaller when the pits experience greater resistance with variable diffusion coefficients. To accurately portray the actual resistance for an undisturbed pit, an even Fig. 9 The simulation box surrounding the surface of Case 1 with streamlines generated from a point source following the flux field. The radial and tangential directions are y and x, respectively, and z is longitudinal. The scale bar in the lower right corner is 25 μm. The color of the streamlines indicates different pathways from tracheid to tracheid through the bordered pits Fig. 10 Simulation box surrounding the surface of Case 4 with streamlines generated from a point source following the flux field. The radial and tangential directions are y and x, respectively, and z is longitudinal. The scale bar in the lower right corner is 25 μm. The color of the streamlines indicates different pathways from tracheid to tracheid through the bordered pits higher ratio of close to 1:25 should be used. However, such a large difference would be difficult to account for due to stability issues in the simulation.
In the model based on the image data, it was assumed that each pit was open, i.e., no pit aspiration. Based on the X-ray images, it is difficult to judge whether a pit has aspirated or not, even though several authors have argued that steam explosion may influence aspiration and have shown cracks in the cell wall close to the pit openings (Zhang and Cai 2006;Muzamal et al. 2015). In an attempt to model a more critical condition of the pit, Kvist et al. (2017) made simulations with a larger diameter of aperture. This significantly increased the effective diffusion coefficient over the pit and would correspond to a ratio in the present simulations of 1:4.

Conclusion
Diffusion of large molecules in the microstructure of wood was investigated with multi-scale implementation of bordered pits using a variable diffusion coefficient. A simple and efficient model was presented that incorporated high-resolution data on a single pit in a larger microscale 3D model. By adjusting the local diffusion coefficient of small features in the microstructure, the effect of pretreatment on different scales could be investigated. Furthermore, the methodology can be used to incorporate diffusion coefficients in different parts of the wood structure, such as latewood and cell wall, and tie scales together for a better overall understanding of the interaction of structure and mass transport properties, with regard to the pretreatment of wood. The method is essentially limited by available computer resources and numerical stability issues.