Macroscopic X-ray computed tomography aided numerical modelling of moisture flow in sawn timber

Mathematical models are essential for the development of schedules for the air-circulation drying of timber in Swedish sawmills, but earlier models have been shown to be conservative leading to longer drying times than necessary. In the current study, macroscopic (macro) X-ray computed tomography (CT) has been used in both the development and validation of a finite element (FE) model, to enable the macro-CT aided FE modelling of the nonlinear transient moisture flow in wood. The model uses more advanced theory than has previously been used in Swedish sawmills, by incorporating a surface emission coefficient to simulate the surface resistance to moisture flow. A single piece of Norway spruce [Picea abies (L.) Karst.] timber was subjected to that part of a traditional kiln-drying schedule, which is associated with diffusion-driven moisture transport. The incorporation of macro-CT data into the FE model resulted in a more realistic representation of the board’s geometry, the initial moisture state, and the definition of material parameters. It also led to a better simulation of flow speed and moisture gradient, especially the asymmetric MC development within the cross section throughout the drying process.


Introduction
The industrial X-ray computed tomography (CT) scanners available in large Swedish sawmills enable non-destructive tests that can recreate the three-dimensional internal inhomogeneous structure of timber on the macroscopic (macro) material level. The spatial resolution of these macro-CT scanners is ca. 0.1 mm for a 100 mm wide specimen. Macro CT was first used to investigate wood in the mid-1980s, when Japanese researchers used a medical CT scanner to determine the moisture content (MC) distribution in wood Kanagawa and Hattori 1985). In the 1990s, this technique was also adopted in Sweden (Lindgren 1991), where a medical CT scanner was used to determine the MC distribution in Norway spruce. The use of macro CT in the field of wood has been summarised by Bucur (2003).
Macro-CT-aided finite element (FE) modelling uses macro-CT data as a non-destructive method to determine material parameters in various stages of model development and in the experimental calibration and validation of such models. The history of macro-CT-aided FE modelling of wood is brief. Pyrkosz et al. (2010) first used the method to obtain the most dominant vibration mode of a Stradivari violin, while Coaldrake (2020) identified the different vibration modes of a traditional Japanese harp, and Hartig et al. (2021) and Huber et al. (2022) used the method to validate the elastic behaviour of moulded wood products and timber, respectively. Pang and Wiberg (1998), Eriksson et al. (2007) and Li et al. (2018) have also actively used macro-CT data in combination with FE modelling. However, the CT data were only used to perform an experimental validation to see whether the model could accurately describe moisture flow.
The aim of the current study was to use macro-CT aided FE modelling to study the moisture flow in wood during the air-circulation drying of timber ( < 100 • C ). Mathematical models play an important role in understanding drying processes in sawn timber and in the development of new drying schedules for air-circulation drying (Pang and Haslett 1995;Perré and Passard 2004;Salin and Wamming 2008;Zhao et al. 2016). The mathematical model currently used in the development of new drying schedules in Swedish sawmills is however conservative in nature, and relies on a 1 3 moisture-and temperature-dependent diffusion coefficient (Salin 1999). In the current study, a more recent FE model has been adopted. This model has been shown to give a better prediction of moisture gradient and flow speed during drying by the introduction of a temperature-and moisturedependent surface emission coefficient (Florisson et al. 2020), which describes the resistance to moisture leaving the wood at the exchange surface. Moisture gradient and flow speed are two important characteristics that influence the hygro-mechanical behaviour of wood and the end-quality of wood products (Simpson 1983a, b).
In the present study, CT data were used to create the model's geometry using a simple manual segmentation procedure. The data were also used to estimate the location of the pith to define the annual ring pattern in the cross-section of the board, to define the MC state to initiate the numerical simulation, to define the dry-density state to estimate the flow parameters, and to perform an experimental validation of the numerical model. The validation considered only the drying regimes that make up the drying schedule associated with the diffusion of water in wood.

Material and methods
The research methodology is illustrated in Fig. 1. It consists of (a) an experiment to image the kiln drying of timber using macro X-ray CT (Couceiro et al. 2020), (b) an image processing step to determine the MC of timber from CT data (Hansson and Cherepanova 2012), (c) a numerical step to develop a CT-aided FE model to simulate moisture flow within timber during kiln drying, and (d) an experimental validation of the numerical simulations based on the outcomes of steps b and c. The figure also illustrates the required input and output for each step. The experiment is described in subsect. 2.1, the image processing step in Fig. 1 Research methodology, where T w,a is the recorded wet-bulb temperature, T d,a is the recorded dry-bulb temperature, T A 1 is the recorded temperature inside the timber, w is the wet density, 0 is the dry density, w v is the moisture content estimated for each voxel, w n is the moisture content estimated for each FE mesh node, w n,i is the moisture content assigned to each mesh node at the start of the simulation, D r,t is the diffusion coefficient related to the radial and tangential material directions, and s is the surface emission coefficient. The subscript '45' refers to a drying time of 45 h subsect. 2.2, the numerical model in subsect. 2.3, and the experimental validation in subsect. 2.4.

Experiment
A sawn Norway spruce timber specimen in the green condition, labelled A 1 , was obtained from the Stenvalls Trä Lövholmen sawmill in Northern Sweden. The specimen had nominal cross-sectional dimensions of 47 × 150 mm (thickness × width) and a length of 700 mm. The cross sections at the ends of the board were sealed with a heat-resistant sealant named Sikafex 221 (Sika, New Jersey, USA) and kept in plastic under frozen condition to avoid uncontrolled drying before the start of the experiment. A medical Siemens Somatom Emotion Duo CT-scanner (Siemens, Munich, Germany) was used in combination with a custom-made integrated kiln (Hansson and Cherepanova 2012) to obtain cross-sectional density images of the specimen during drying. A detailed description of the experimental setup has been reported earlier by Couceiro et al. (2020).

Kiln drying
The drying schedule (wet-bulb temperature, T w,s , and dry-bulb temperature, T d,s ) used in the experiment is presented in Fig. 2a. This schedule was provided by a sawmill outside Skellefteå in Northern Sweden. The figure also presents the recorded wet-bulb temperature, T w,a , and drybulb temperature, T d,a , provided by the integrated sensors in the climate chamber, and it gives an indication of the onset of each drying regime (1-5). The schedule consisted of a sequence of heating, capillary, transition, diffusion, and conditioning regimes.
Within the heating regime, the T d,s and the T w,s were raised to 66 °C and 63 °C, respectively, in 14 h. In the capillary regime, T d,s was slowly increased to 67.5 °C and the T w,s was kept constant for 32 h to remove capillary (free) water from the lumen above FSP. The transition regime was used to bring T d,s up to 70 °C and T w,s up to 64.5 °C in a period of 14 h to achieve a smooth transition between capillary-driven and diffusion-driven drying. In the diffusion regime, T d,s was kept constant and T w,s was gradually lowered to 60 °C over a period of 28 h to lower the relative humidity inside the kiln and to remove bound water from the cell walls below FSP. The conditioning regime was characterized by a constant T d,s and an increase in T w,s to 67.5 °C to raise the relative humidity in the kiln in 7.65 h, in order to elevate the MC directly beneath the board surface. In the experimental validation presented in this paper, the data collected from the transition regime onwards are used. Figure 2a also shows the recorded temperature inside the timber board, T A 1 , obtained with a TGP-4017 Tinytag plus 2 (Intab, Stenkullen, Sweden) in combination with a PB-5006-1m5 Allround Tinytag sensor (− 40 °C to +125 °C, ± 0.2 °C between 0 °C and +70 °C). It can be seen that T A 1 was generally lower than T d,a . The position of the sensor within the cross section of the timber board is shown in Fig. 2b. The air velocity, v a , within the chamber was set to a constant value of 5 m s −1 .

Computed tomography scanning
The medical CT scanner used to scan the timber consisted of an X-ray tube and four rows of detectors, both of which were aligned at opposite positions and could rotate around the gantry where the board was located to produce a fan beam. For each full rotation, ca. 1000 projections were collected and used to reconstruct a two-dimensional image through a back-projection algorithm (Hsieh 2015). The scan parameters were set to 110 kVp and 88 mAs, and the total scanning time was 0.8 s. The reconstruction process resulted in images with a pixel resolution of 0.586 × 0.586 mm 2 and a 5 mm voxel depth. Since only single scans were considered, the reconstruction was two-dimensional. The location of this image was at the centre of the board in the lengthwise direction.
The image collection was done at regular time intervals. A one-hour (h) time interval was adopted during the heating, capillary, transition, and diffusion regimes, and a 15 min time interval during the conditioning regime. Each image was stored in the medical image format DICOM. The raw images contained a CT number (CTN) in Hounsfield units at each voxel point. The CTN was determined according to:  (5) conditioning regimes, where a illustrates the scheduled wet-bulb, T w,s , and dry-bulb, T d,s , temperature, the recorded wet-bulb, T w,a , and dry-bulb, T d,a , temperature, as well as the measured temperature within the specimen T A 1 , and b shows the location of the temperature sensor within the cross section of the specimen to measure T A 1 where is the linear attenuation coefficient of the material, and w and a are the attenuation coefficients of water and air, respectively (Kalender 2011;Hsieh 2015). There may be further processing of the CT images by the CT scanner that can affect the CTNs, such as filtering, but this is part of the proprietary software and unknown to the user of the scanner. The CTNs were converted to density in kg m −3 by the metadata of the DICOM files based on an internal law that gives the linear relationship between these two variables as: where is the density of the scanned wood. Equation (2) is applicable to all ranges of MC. The accuracy of the scans is ensured by a calibration of the machine using phantoms of known density, such as water and air. If necessary, a corrected calibration function can be included in the density calculation. The calibration for the current CT scans was shown to be well within the range of accepted values. The existence of a linear relationship between the CTN and wood density was proven by, among others, Lindgren (1991) and has been widely applied since then.

Image-processing algorithm
The image-processing algorithm presented by Hansson and Cherepanova (2012) was adopted to estimate the MC at each voxel point, w v , using the wet-density, w , and dry-density, 0 , images and the expressions: where max is a coefficient indicating the dimensional change due to shrinkage, V w is the wet volume, and V 0 is the dry volume (Hansson and Cherepanova 2012). The proportionality between the grey scale values of a CT image and the MC of wood was established in the early 1990s by for example Lindgren (1991). The temperature-dependent fibre saturation point, w f , was estimated with Eq. (5) established by (1995). The Hough transform for circles was used to determine the pith location. A challenge in the estimation of the MC is the shrinkage of the timber during drying. This distorts the crosssection of the timber and prevents the one-on-one comparison between density images. The bidirectional elastic registration algorithm published by Arganda-Carreras et al. (2006) was therefore integrated into the image-processing algorithm and used to transform the shape of the dry-density image to the shape of the wet-density image. The mapped density values at each pixel were corrected by using polygon clipping (Hansson and Fjellner 2013).

Surface detection
A manual threshold segmentation method was adopted to define the geometric boundary of the specimen. This segmentation resulted in a two-dimensional binary image, B 0 (x, y) , using the dry-density image, P 0 , and a threshold value, B . The threshold value was based on the two density modes (air and wood) found in the histogram of P 0 (Prewitt and Mendelsohn 1966). The value used for B was 250 kg m −3 . The segmentation procedure was implemented in the MATLAB ® software (MathWorks, Natick, USA).

Data smoothing
The dense latewood of Norway spruce generally has a higher MC than the less dense earlywood (Fromm et al. 2001). These differences in density and MC between latewood and earlywood were also detected in both the density images obtained through CT scanning ( Fig. 3) and in the MC images created with the image-processing algorithm (Fig. 4). However, the numerical model assumed a homogenized cross section which made no distinction between earlywood and latewood. To be able to use the density and MC data in the development of the numerical model and in the execution of the experimental validation, the density and MC were therefore smoothed. This smoothing also improved the numerical convergence and meant that parts of the data were emphasized. For example, the cross-sectional variation in density between sapwood and heartwood became more distinct, as well as the moisture gradient within the cross section of the board. An averaging scheme was chosen that used a two-dimensional Gaussian kernel, G , which represents a normal distribution based on a standard deviation, (Russ 1995). To deal with smoothing at the boundary of the board, the surface detection (see subsection 2.2.2) preceded the smoothing. The smoothing was executed in the MATLAB ® software using statement 'imgaussfilt'.

Dry-density data
The oven-dry density image, P 0 , of specimen A 1 presented in Fig. 3a was determined by exposing the specimen to a temperature of 103 ± 2 °C until the difference in mass between two successive weighings (interval of 2 h) was less than 0.1% (CEN 2003). Figure 3a also shows the actual thickness (154.2 mm) and width (52.2 mm) of the specimen at the beginning of the transition regime, and the location of the pith (80.4 mm, 54.6 mm) determined with the imageprocessing algorithm.
In Fig. 3b and c, the oven-dry density is given along a row of pixels centered along the width and thickness of the specimen's cross section, respectively. The plots give both the original data (exp) and the smoothed data (gauss). Based on the smoothed data, the mean dry-density of the specimen was determined to be 349 kg m −3 with minimum and maximum dry-density values of 256 kg m −3 and 391 kg m −3 , respectively. Figure 3d gives the green state (wet) density image obtained at the beginning of drying ( t = 0). The border between heartwood and sapwood was derived from the difference in color between heartwood and sapwood in the green state. This difference is due to the rapid increase in MC between the border and the surface of the board. The percentage of heartwood in the specimen was 54%.  Siau (1995). In the case of the smoothed data, most of the cross section was below FSP at the start of the transition regime, indicating that the moisture flow was driven only by diffusion. Figure 4c and d present smoothed moisture vs distance plots at different stages of the drying process starting at the transition regime. The onset of the transition regime is indicated as t = 0. The plots show that the MC distribution was unsymmetric and differed in magnitude between the two edges.

Numerical model
The non-linear transient flow of moisture in wood is here described by a single-Fickian model and a Neumann boundary condition dependent on a surface-emission coefficient. The numerical model was created in the Abaqus FEA ® (Dassault Systèmes 2017b) FE software. The material model and the boundary condition were implemented in the user-subroutines UMATHT and FILM, respectively. The material orientation needed to describe the fiber orientation of wood was defined in the user-subroutine ORIENT. The non-isothermal diffusion of water and sorption hysteresis are not considered by the model. The sub-routines mentioned are all supported by Abaqus FEA ® (Dassault Systèmes2017a). A full description of the theory used to simulate the moisture flow and its implementation into the user-subroutines is given in Florisson et al. (2020). The corresponding theory Fig. 4 Moisture content determined at the voxel level, w v : a w v centered along the width of the board, b , b w v centered along the thickness of the board, h , c w v centered along the width of the board, b , covering different times during the transition, diffusion and conditioning regimes, and d w v centered along the thickness of the board, h , also covering different times during the last three regimes. a and b include both the original data, the smoothed data using a filter with a of 10, and an indication of the FSP at the beginning of the transition regime to describe the material orientation is given by Ormarsson (1999).

Theory
The balance equation and constitutive relationship used to describe the non-linear transient moisture flow in wood are given by Eqs. 9 and 10, respectively: where c w is the moisture capacity, w is the density of water, q is the moisture flux vector, D is a moisture-, dry-densityand temperature-dependent diffusion matrix, 0 is the dry density, w is the unit-less MC, T is the temperature and ∇ defines gradients. The symbol (•) in Eq. (8) refers to the orthogonal local coordinate system for wood. The moisture capacity is an indicator of how well the wood retains moisture. For simplicity, it was assumed that the balance equation is divided by c w and w and that the inverses of c w and w are implicitly included in the diffusion coefficient. In general, w is dependent on the concentration of water inside the wood, c , and the dry density in the sense that w = c∕ 0 . In the current model this is not incorporated, which means instead that w is assumed to be based on a constant dry density.
The Neumann boundary condition that describes the flux normal to the exchange surface q n is defined as: where, s is a density-, moisture-, and temperature-dependent surface emission coefficient and w ∞ is the unit-less MC of the ambient air. Further information concerning the derivation of the FE formulation and how to iteratively solve the non-linear system of equations is given by Florisson et al. (2020) and Johannesson (2019).

Model geometry
The meshed geometry of board A 1 is shown in Fig. 5a and b. Figure 5a shows the initial moisture content distribution, while Fig. 5b shows the dry-density distribution. These are discussed in further detail in subsection 2.3.4. The model is assumed to be 1 mm thick, creating a more-or-less twodimensional scenario. Figure 5c presents the dry-density CT image with an indication of the absolute thickness (154.2 mm) and width (52.2 mm) of the specimen at the beginning of the transition regime. The figure also shows the local (l, r, t) coordinate system associated with the wood's fiber direction, and the global (x, y, z) coordinate systems. The estimated pith location (80.4 mm, 54.6 mm in the global domain) obtained with the image-processing algorithm (see subsection 2.2.1) was used to define the annual ring pattern in the cross-section of the board using user-subroutine ORIENT.
To mesh the board, a geometric meshing technique was used. The surface geometry was first obtained from the CT image corresponding to the beginning of the transition regime. The threshold segmentation method presented in subsection 2.2.2 was used to identify the voxels located at the surface of the board in relation to the global coordinate system. These data were then imported into the FE program as individual points and interpolated linearly. A partition, indicated by the black solid line in Fig. 5a, was made to create a structured 1 mm hexahedral mesh of type DC3D8 in the central region of the specimen, and an unstructured hexahedral mesh at the edges, in order to benefit the numerical Fig. 5 Model geometry based on CT data: a meshed model geometry, including a representation of the MC at the start of the simulation, b the dry-density, 0,n , profile used by the FE program, and c the dry-density image from CT-scanning with an indication of the global (x, y, z) and local (l, r, t) coordinate systems analysis and to enable the more efficient extraction of results at certain locations.

Material parameters
The expressions used to determine the diffusion coefficient, D , and the surface emission coefficient, s , are: where r and t indicate the radial and tangential material directions, respectively. Both material parameters were assumed to be dependent on oven-dry density, 0 , moisture content, w , and temperature, T (Avramidis and Siau 1987;Sehlstedt-Persson 2001;Yeo et al. 2008). The coefficients D w , D T , s w and s T were determined on the basis of experimental data previously collected by Florisson et al. (2021) for different types of pine spp. from for example Avramidis and Siau (1987), Rosenkilde (2002) and Yeo et al. (2008). The data are listed in Table 1. The linear relationship between the diffusion coefficient and dry density, which results in larger coefficients for smaller densities, is based on experimental data published by Sehlstedt-Persson (2001). Due to a lack of data, the same relationship was used to describe the dependence of the surface emission coefficient on the density.
In Fig. 6a and c, the diffusion and surface emission coefficient, respectively, based on the coefficients presented in Table 1 and used in the numerical validation, are given for the temperature ranges and a mean density of 349 kg m −3 obtained through the experiment. The assumed relationship between the material parameters and density is shown in Fig. 6b and d for a relative humidity of 85% and a temperature of 60 °C. The experimental data (exp) for heartwood and sapwood specimens of Norway spruce in Fig. 6b are taken from Sehlstedt-Persson (2001). The values lie between 3.24e-6 and 1.08e-6 m 2 h −1 given by Eq. (10). Figure 5a presents the MC distribution at the beginning of the simulation, i.e. at the beginning of the transition regime. Figure 5b shows the dry-density distribution needed (10) D r,t 0 , w, T = D ρ1 0 + D w1 e D w2 w + D w3 D T1 e D T2 T (11) s 0 , w, T = s ρ1 0 + s w1 e s w2 w + s w3 s T1 e s T2 T to define the material parameters. These initial conditions were generated by the FE program based on the smoothed MC and dry-density data described in subsections 2.2.4 and 2.2.5, respectively. The FE program used a distance weighting algorithm as a deterministic method for multivariable interpolation of the imported data to the nodes that make up the FE mesh. The default settings for mapping set by the FE program were used. The change in temperature inside the sawn timber was included in the simulation by assigning the measured temperature T A 1 (see Fig. 2) to the mesh nodes that lie within the partition that encloses the inner region of the cross-section. This created a constant temperature field within each time increment used to determine the size of the material parameters. Table 1 Variables needed to describe the diffusion coefficient, D , related to the radial, r , and tangential, t , material direction and the surface emission coefficient, s

Initial conditions
The value of D related to the longitudinal material direction is assumed to be zero

Fig. 6
Material parameters: a relation between the diffusion coefficient, D , and moisture content, w at different temperatures, b relation between D and dry density, 0 , c relation between the surface emission coefficient, s , and w at different temperatures, and d relation between s and 0 .

Boundary condition
The climate surrounding the board during drying was simulated as the MC of the ambient air (see Fig. 7), w ∞ , presented in Eq. (9). The value of w ∞ within each time increment was determined in two steps, where the first step was to determine the RH, , based on the recorded T d,a and T w,a , and the method presented in ASHRAE Handbook ASHRAE (2017). This method was implemented in MATLAB ® . The equations used to determine were: where p v,w∕d is the vapor pressure determined with the Antoine formula based on either the recorded dry-bulb temperature inside the climate chamber, T d,a , or on the recorded wet-bulb temperature, T w,a . p a is the atmospheric pressure and had a value of 101.325 kPa. is the degree of saturation and w is the degree of saturation corresponding to the wet-bulb temperature. T d,a and T w,a in Eq. (14) are expressed in Kelvin. The second step was to estimate w ∞ using Simpson's formula (Simpson 1971), and T d,a . The procedure was executed inside user-subroutine FILM. The recorded temperature T d,a was assigned to the mesh nodes that lie outside the partition, close to the timber surface (see Fig. 5). Simpson's formula was originally created for Sitka spruce (Picea sitchensis) and is based on an average of the adsorption and desorption curves.

Experimental validation
The experimental validation was based on the relationship between the simulated moisture content, w n , and the moisture content calculated with the image-processing algorithm, w v , during the transition, diffusion, and conditioning regimes of the drying schedule. For this purpose, moisture vs distance plots were used to validate the moisture gradients and drying speed, and moisture vs time plots were used to validate the drying speed. The validated model was used in a parametric study to see the effect of temperature and dry density on the moisture gradients and drying speed during kiln drying.

Results and discussion
The experimental validation of the numerical model involved a single board of Norway spruce dried with a conventional kiln drying schedule. Since results for a single specimen cannot be assigned any statistical probability, tests on more specimens are recommended in the future.
It is well known that the Gibbs phenomenon, beam hardening, and ring artifacts influence the results obtained with the medical CT scanner. The Fourier series used to create the two-dimensional image from the CT projections can lead to oscillatory behavior between the air and wood in the CT image, called the Gibbs phenomenon, and this can lead to locally higher wood density values beneath the specimen surface. The affected length is usually about 5 voxels. This phenomenon, however, was not observed in the CT images used in this study. Beam hardening, which is caused by the selective attenuation of low energy photons, results in lower density values in the center of the specimen and higher density values at the edge of the specimen (Murphy and Glick 2022). For a low-density material such as wood, this effect is small, especially below FSP. Ring artifacts due to defects in the detector Fig. 7 Climate schedule needed for the numerical model, including data used as input: is the relative humidity, T d,a is the dry-bulb temperature inside the climate chamber, T w,a is the wet-bulb temperature inside the climate chamber, and w ∞ is the estimated MC of the ambient air. The three drying regimes considered are the transition regime (3), the diffusion regime (4), and the conditioning regime (5) are easily detected in the CT images and were not present in this case (Fig. 7). Figure 8 shows the validation of the numerical model in the time domain at two points within the cross section of the simulated piece of timber. One point is located in the middle of the timber piece in the heartwood region and the other is located 6 mm beneath the exchange surface in the sapwood region. The two points represent two extremes within the cross section. Figure 8 shows a good fit between the experimental (w v ) and simulated (w n ) data. In Fig. 8a, the largest difference between w v and w n occurs in the diffusion regime and is 1.9% MC. The largest difference in Fig. 8b is found in the conditioning regime and is 1.2% MC.

Validation in the time domain
At the start of the analysis, the MC differed by only 1.9% MC between the selected points in the heartwood (24.8% MC) and sapwood (22.9% MC) regions. The diffusion regime led to a MC slightly higher than the target MC. This led to moisture gradients within the cross section of the sawn timber and to an area below the timber surface that is subjected to case hardening (McMillen 1958;Lazarescu et al. 2010). Case hardening is a term used to refer to the compression stress state just below the board surface at the end of kiln drying due to stress reversal. After about 40 h of drying, the point closest to the surface of the timber has reached a lower MC value (12.8% MC) than the point in the center (17.3% MC), i.e. a difference of 4.5% MC.
To undo case hardening and minimize the occurrence of cracks and distortion, the timber is subjected to a conditioning regime characterized by an increase in the relative humidity in the kiln. Figure 8 shows that during the conditioning regime, the center of the timber continues to dry, whereas the region close to the surface becomes more moist leading to MCs of 16.2% and 15.0%, respectively, an absolute difference of only 1.2% MC.

Validation in the space domain
The MC vs distance plots obtained over the thickness and width of the cross section are shown in Fig. 9. To create the initial MC, the FE program used an automatic linear interpolation procedure to overcome the difference between the voxel positions in relation to the mesh nodes. This resulted in a slight overestimation of the MC at the exchange surface. Figure 9a shows a maximum overestimation of 0.1% at the right-hand side of the MC vs thickness plot, and a maximum overestimation of 1.1% at the right-hand side of the MC vs width plot.
At the beginning of the transition regime, in Fig. 9a, the MC of almost the entire cross section was below the estimated FSP, but although the asymmetric distribution of the MC over the thickness of the specimen decreased during the drying process, it remained until the end of the conditioning regime. A comparison of Fig. 9e and f, respectively, shows that the increase in MC of the ambient air during the conditioning regime increased the MC beneath the surface, but that the inner region of the board became drier, reducing the MC gradient within the cross section.
A reasonable correlation was found between drying speed and moisture gradient. The MC vs thickness plots in Fig. 9 based on the MC obtained from the image-processing algorithm (exp) show that the region at the left-hand side dried slightly slower than the region at the right-hand side, despite the slightly lower dry density (see Fig. 3). The model was not able to simulate this larger drying speed on the righthand side. After 25 h of drying (Fig. 9f) this led to a 3.3% difference between the experimental and simulated MC (∆w) values. The moisture gradient simulated over the width of the timber was generally smaller than the experimental gradient. The maximum difference in MC (∆w) of the plots Fig. 8 Validation of the MC change with time at points a in the center of the timber and b 6 mm below the surface, sharing a pixel in relation to the experimental data (wv) or a nodal point in relation to the simulated data (wn) 1 3 Fig. 9 Moisture content, w , over the thickness, h , and width, b , of the timber cross-section obtained with the image-processing algorithm (wv) and with the numerical model (wn) at different times during the drying process, together with an indication of the FSP at the beginning of the transition regime and a color wn-plot of the entire crosssection over the width of the board was 3.5% MC after 10 h of drying at the right-hand surface (Fig. 9c).
When interpreting the MC data, it should be kept in mind that the gradient obtained with the Gaussian filter is dependent on the chosen size of the standard deviation, σ. The difference between experimental and simulated results can also be due to the absence of a realistic temperature gradient, which leads to non-isothermal diffusion, and to the modelling assumption that the MC is estimated with a constant dry density.
The specimen was selected at the sawmill in the green state based on its heartwood to sapwood ratio. Its mass and volume were taken from CT data and not directly from the specimen. The average dry density after smoothing was 349 kg m −3 and the minimum and maximum dry densities were 256 kg m −3 and 391 kg m −3 , respectively. These values are relatively low. In Niemz and Sonderegger (2021) (page 148), the dry density of spruce is said to vary between 370 and 540 kg m −3 , with a mean value of 430 kg m − 3 . Boutelje and Rydell (1995) (page 27) mention values between 370 and 440 kg m −3 for Norway spruce.
The diffusion coefficient used to perform the simulation was based on experimental data for different pine species. The dry density of pine is generally higher than that of spruce, although the diffusion coefficient is generally lower than that of spruce (Sehlstedt-Persson 2001). In Niemz and Sonderegger (2017) (page 144), the dry density of pine is reported to be between 300 and 860 kg m −3 , with a mean value of 490 kg m −3 . Boutelje and Rydell (1995) (page 22) mention values between 450 and 500 kg m −3 . In Niemz and Sonderegger (2021) (page 86), the diffusion coefficient of spruce ranges between 2.52e-7 and 7.56e-7 m 2 /h for a MC range between 5 and 20% MC. The temperature is 60 °C but the density is unknown. According to Sehlstedt-Persson (2001) the diffusion coefficient for Norway spruce for 15% MC and the same temperature is between 1.4e-06 and 3.0e-06 m 2 /h for a density range between 550 and 375 kg/m 3 . The diffusion coefficient value used by the numerical model at 20% MC and 60 °C was 8.5e-06 m 2 /h, which is 91% larger than the value found by Niemz and Sonderegger (2021) for a similar MC and temperature. At an average density of 349 kg m −3 , however, the diffusion coefficient in the simulation at 15% MC and 60 °C was 6.4e-6 m 2 /h, which is 53% larger than the value found by Sehlstedt-Persson (2001) for a similar dry density.
Experimental data for the surface emission coefficient of spruce are sparse and this prevents an active comparison with the data used in the numerical model. In a recent publication, however, Yeo et al. (2008) reported a value for pine of 0.25e-3 m h −1 at 20% MC and 65 °C. The dry density was not mentioned. The value used in the numerical model for similar climatic conditions and a dry density of 349 kg m −3 resulted in a value of about 0.59e-3 m h −1 . This value is of the same order as that of pine.
The simulations made in Florisson et al. (2020) show how important the correlation is between the size of the diffusion coefficient and the size of the surface emission coefficient on the size of moisture gradients within the cross section of timber and the flow speed of moisture. The diffusion coefficient used in the current numerical model was higher than the experimental values given for spruce in the literature, but the experimental values are based on a theory that does not take into account the surface resistance. Avramidis and Siau (1987) reported a considerably higher diffusion coefficient for pine when the surface emission coefficient was taken into consideration when the diffusion coefficient was estimated.

Parametric study
A parametric study was undertaken to investigate the effect of dry density and temperature on the development of moisture gradients and the drying rate within a piece of timber. The first part of the parametric study included two simulations, one of which considered the dry density of the timber to be constant at 300 kg m −3 and one that considered the dry density to be constant at 700 kg m −3 , instead of the variation in Fig. 5. The chosen densities cover dry-density values that can be expected for both Norway spruce (370 and 550 kg m −3 ) and Scots pine (between 300 and 850 kg m −3 ).
The second study also consisted of two simulations, one that assumed the temperature of the cross section to be constant at 50 °C and one that assumed the temperature to be constant at 80 °C, instead of the temperatures suggested in Sect. 2.1.1. This change in temperature meant that not only material properties but also the MC of the ambient air, w ∞ , was affected. The moisture vs distance plots over the thickness of the timber obtained after 15 h of drying are presented in Fig. 10, together with the results of the validated model. In addition, an indication is given of the difference in MC, Δw n , between the two simulations in each study.
The density range (300-700 kg m −3 ) used to generate the plots in Fig. 10a is much larger than the density variation (256-391 kg m −3 ) assumed in the numerical validation. At 20% MC and 65 °C, the range resulted in a diffusion coefficient between 1.1e-5 and 0.81e-5 m 2 h −1 and a surface emission coefficient between 6.2e-4 and 3.8e-4 m h −1 , leading to a maximum MC difference of 1% close to the exchange surface. The temperature range (50-80 °C) used for Fig. 10b was chosen to be slightly larger than the temperature range assumed in the numerical validation (60-70 °C). At 20% MC and 65 °C, this temperature range resulted in a diffusion coefficient between 0.57e-5 and 1.7e-5 m 2 h −1 and a surface emission coefficient between 2.3e-4 and 12.0e-4 m h −1 . The simulations resulted in a maximum difference in MC of 5.7% close to the exchange surface.
The lines indicating Δw n in Fig. 10a and b show the difference in MC obtained between the two simulations in each parametric study, where a difference in either density or temperature affects both the speed of flow and especially the moisture-content gradient. A lower density or a higher temperature leads to a steeper gradient and more rapid drying. The comparison between Fig. 10a and b shows that the effect of density on the MC development is smaller than that of temperature. The speed and gradient are influenced not only by the size of the diffusion coefficient and the size of the surface emission coefficient but also by the ratio of these two coefficients. The size and range of the flow parameters considered are smaller in the density study than in the temperature study. Especially, at higher temperatures the surface resistance becomes low, and this increases both the flow and the gradient.

Conclusion
The incorporation of CT data into the FE modelling resulted in a more realistic representation of the board geometry, the initial moisture state, and the material parameters. It also led to a more realistic representation of flow speed and moisture gradients during drying since the asymmetric development of MC within the cross section of the board could then be simulated.
The incorporation of an actual drying schedule used in kiln drying showed that the model can simulate the flow behavior during the transition, diffusion, and conditioning regimes, and the simulation confirmed that the moisture gradients within the cross section were reduced during the conditioning regime. A reasonable fit was found between the simulated data and the experimental data obtained from CT studies. The fit could probably be improved by incorporating non-isothermal diffusion into the numerical model and the dependency between concentration, dry density, and moisture content in the theoretical definition.
The dry density of the timber cross section had a smaller effect than the temperature on the MC development during drying, and it was observed that a lower dry density or a higher temperature created a larger MC gradient and led to faster drying. Funding Open access funding provided by Lulea University of Technology. This work was supported by the Swedish Kempe and Gabrielssons Foundation. Support of the CT WOOD, a center of excellence at Luleå University of Technology supported by the Swedish wood industry, is also gratefully acknowledged.

Data availability Not applicable.
Code availability Not applicable.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not Fig. 10 Moisture content, w n , over the thickness, h , of the simulated timber piece obtained after 15 h of drying for a the validated model, with two constant values of dry density, permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.