Experimental investigation of microscale mechanisms during compressive loading of paperboard

Compression of paperboard is a common procedure during industrial package forming and better knowledge of the material response is needed to avoid defective packages and waste. To go beyond current modelling approaches, experimental identification of mechanisms underlying the macroscopic stress–strain responses is needed. In this study, in-situ uniaxial compression of paperboard is studied through synchrotron tomography at high spatiotemporal resolutions. Both the microstructural evolution of the fibre network and the actual boundary conditions of the loading were quantified and analysed. At the microscale, the loading equipment plates were not perfectly flat resulting in an increasing sample-equipment contact area with loading. This is, however, shown to only have a small effect on the form of the macroscopic stress–strain curves. The evolution of 3D strain fields showed that strain accumulated close to the sample surfaces in the early part of the compression process, whereafter the main deformation zone shifted to the out-of-plane centre. Both fibre walls and pore volumes were observed to decrease during loading (and recover partly after unloading). Regarding the pore volume, the main reduction mechanism was seen to be closure of layers between fibres. Even if the total pore volume reduction was seen to be the dominant deformation mechanism in a second stage of compression, the volumetric change of fibre walls was non-negligible. Fibre wall compression is not commonly considered in theoretical treatments of paperboard compression, but this work suggests that the stored elastic energy could be a driver for the elastic recovery of the fibre network during unloading.


Introduction
Paperboard is a widely used packaging material for many goods including food and beverages. In this context, mechanically stable packages are of high importance since they can ensure the quality of the contents and minimize waste. When transformed into packages, the mechanical properties of paperboard materials are altered by operations such as

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10570-023-05168-x. stretching, compression, creasing and folding or 3D forming. Numerical modelling is often employed to predict the mechanical response of the material structure when subjected to these operations based on continuum descriptions of the paper. Examples of continuum-scale material models for paperboard include that used by Wallmeier et al. (2015) to simulate deep drawing operations and the models by Xia et al. (2002) and Borgqvist et al. (2014) that have been the basis for simulations of industrial operations ranging from creasing (e.g. Nygårds et al. 2009) to large-scale package forming (e.g. Robertsson et al. 2018). Such models are usually calibrated on the macroscale response of the material from boundary measurements in standard mechanical tests and used to predict local changes within the thickness of the board. However, the accuracy of these descriptions and their ability to describe internal deformation and failure mechanisms is not well established. For more accurate understanding of the model validity and for further material model development greater understanding of the physical changes and mechanisms occurring within the material through in-depth experimental analyses are essential. In this context, fullfield experimental data with resolution at the level of the material microstructure are required to go beyond macroscale analyses.
Out-of-plane compression of paperboard is an example of a key process present in many industrial processing steps, such as calendering, creasing, deepdrawing and 3D-forming (e.g. Rodal 1989;Östlund 2017). As such, uniaxial out-of-plane compression testing is a standard approach for calibration of material models for paperboard. However, to date, physical interpretations of microscale compression mechanisms occurring in such tests are primarily based on macroscopic mechanical testing (stress-strain responses), curve fitting procedures or numerical modelling at different scales. In the very first stage of loading, surface roughness effects are considered to be significant (e.g. Stenberg et al. 2001;Kananen et al. 2018;Chen et al. 2020), causing surface deformation and a slow increase in stress with strain (Stenberg et al. 2001). Thereafter, pores and material voids within the paper start to collapse. Macroscopic stress-strain curves typically exhibit an exponential appearance (e.g. Nygårds 2008), which is attributed to the gradual decrease of void volume leading to a resistance to further deformation (Rättö and Rigdahl 1998). Rodal (1989) described two further stages of deformation after this initial elastic reduction of void volume; elastic-plastic buckling of the fiber network and, ultimately, fiber crushing at higher stresses (Rodal 1989). Other works propose bending of fibres around contact points and compression of lumens as possible additional compression mechanisms (e.g. Lundquist et al. 2004;Toll and Manson 1995). Pawlak and Keller (2005) suggested that fiber slippage, shear deformation of fibres and bond breakage also occur. However, experimental validation of these proposed mechanisms or other possible phenomena is lacking. From a modelling perspective, questions of interest include more detailed and experimentally verified understanding of the mechanisms occurring during the first onset of material compression as well as factors affecting the different slopes along the macroscopic load-unload curves, such as the influence and emergence of heterogeneity, the partitioning of strain between fibres and pores and the origins of permanent deformation.
High-resolution X-ray tomography in combination with in-situ mechanical testing provides the opportunity to image physical changes within samples in 4D (3D + time) while they occur. Analysis of such data have the potential to reveal active microscale mechanisms during in-situ compressive loading of paperboard, both qualitatively and quantitatively with image analysis. In a recent study by Wallmeier et al. (2021), X-ray tomography was used to visualize and analyse in-plane compression mechanisms in paperboard, such as delamination and wrinkle formation. For out-of-plane compression of pressboard samples, Joffre et al. (2015) used X-ray tomography and Digital Volume Correlation (DVC) to study material changes between the non-loaded and the maximum compressed state. They observed a uniform densification of the sample, as well as a correlation between the directions of the largest density gradients and the first principal strain component. The latter was interpreted as a flow of material from high density to low density zones. For other paper-based materials than pressboards, experimental full-field methods have not yet been applied to study 3D strain fields during uniaxial out-of-plane compression and their development over time have, to our knowledge, never been published for any paper-based material.
In the present study, we investigate the process of uniaxial out-of-plane paperboard compression in detail through in-situ testing with 4D synchrotron X-ray tomography. The objective was to image and quantify the microscale 3D structural evolution within the samples at high spatial and temporal resolution during load-unload cycles. The main aim of the paper is to provide experimentally-based interpretations of microscale mechanisms during compression of paperboard through quantified full-field structural and 3D strain evolution measurements that highlight when and where within the sample deformation occurs. A second objective of the experiments was to quantify the actual boundary conditions of the test at the microscale, i.e., the evolution of the contact area between the loading equipment and the sample to, thus, investigate the impact of non-perfect boundary conditions on macroscopic stress-strain curve. Experimentally supported interpretations of different deformation mechanisms can provide increased knowledge of the compressive behaviour of paperboard materials and contribute as a basis for advancing numerical modelling of paperboards at various scales ranging from continuum scale to fibre network modelling.

Experiments
Two types of commercial paperboards were tested with in-situ uniaxial compression during synchrotron X-ray tomography at high spatial and temporal resolution at the TOMCAT beamline, Swiss Light Source (Stampanoni et al. 2006). Circular samples (7 mm in diameter) were prepared from a singleply (1-ply) paperboard from Westrock (WR) and a three-ply (3-ply) paperboard from Billerud-Korsnäs (BK), respectively. Both paperboards have a surface layer of clay/latex coating and have been previously investigated under in-situ uniaxial tension with X-ray tomography and DVC analysis (Johansson et al. 2021(Johansson et al. , 2022. The 1-ply WR paperboard consist of bleached softwood and hardwood fibre kraft pulp and is produced at the Evadale Mill, USA. The sheet structure in the 3-ply BK paperboard, produced from kraft pulp at the Gävle Mill, Sweden, consists of bleached softwood and hardwood fibres in the top ply, unbleached softwood and hardwood fibres in the mid ply and unbleached softwood fibres in the bottom ply. Basic nominal properties of the materials can be found in Table 1. Two samples of each material were prepared and tested under load-unload uniaxial compression cycles with continuous deformation and X-ray tomography acquisition. The loading was performed using an inhouse developed loading device designed for in-situ use during X-ray tomography imaging. Figure 1a shows the device mounted on the rotation stage at the synchrotron beamline. The samples were placed on a polymer (polyether ether ketone, PEEK) bottom plate attached to a motor-driven metal piston (Fig. 1b). A transparent plastic (Polymetyl methacrylate PMMA) cylinder with a polymer (PEEK) top pin (3 mm in diameter) was then attached to complete the experimental setup (Fig. 1a). Before each test, the top pin was vertically adjusted with a screw until it was in contact with the sample surface (Fig. 1a). The device is controlled via a Labview program, which controls and records the piston displacement and the global force from a linear potentiometer and a 500N force transducer (U9C from HBM with accuracy class 0.2) mounted in-line with the piston below the sample.
During the continuous compression tests, the piston was displaced upwards with an approximate rate of 52 µm/min. The low displacement rate was selected to enable X-ray tomography imaging without significant motion artifacts. When a global force of -71 N was reached, the motion of the piston was shifted to the opposite (downward) direction until a global force of 0 N was recorded. During unloading, the displacement rate varied between 32 -45 µm/ min. Based on the cross-sectional area of the top pin (7.1mm 2 ), the samples were compressed to a maximum pressure of approximately 10 MPa.
The X-ray tomography imaging was performed with the GigaFRoST system (Mokso et al. 2017) at the TOMCAT beamline, designed for fast data acquisition with polychromatic radiation. The basic components of the system include a visible-light CMOS camera detector, a scintillator and an optical lens system with 6.8 × optical magnification resulting in a voxel size of 1.61 μm for a Field of View (FoV) of 3.25 × 3.25 mm 2 , which was centered on the sample (Fig. 1c). After acquiring a first 3D image of the sample, the compression test was started and synchronized with the imaging via a trigger signal from the camera to the Labview program. Each scan took 1 s to perform and consisted of 1000 projections collected during 180° rotation. The time interval between scans was 16 s and the experiment was completed once the sample had been unloaded to 0 N. This resulted in between 29-41 data sets per compression test. From the acquired projections, absorption contrast tomography images were reconstructed with the fast Fourier transform based gridrec algorithm (Marone and Stampanoni 2012;Marone et al. 2017). The reconstructed images were cropped to a region of interest (ROI) marked in Fig. 1c (3.25 × 0.88 × 3.25 mm 3 ) and postprocessed with 3D Gaussian filtering to remove some of the noise.

Segmentation
The reconstructed and postprocessed 3D images were segmented to enable extraction of the evolution of different material properties and boundary conditions during loading. Equipment, fibres, pores and background were segmented from each other during the process. Figure 2a shows an example slice through one of the 3D images at two different time points (before loading and at maximum compression). Because of similar intensities of fibres and the compression equipment (Fig. 2b), X-ray phase effects around the plates with similar intensities as pores, and variable contact between the equipment and the material during compression, the segmentation procedure consisted of several steps. In the first step, a segmentation mask representing the loading equipment was extracted. The c Schematic illustration of the experimental setup with the X-ray FoV and ROI of the reconstructed images extracted mask was used to obtain a correct interface between the sample surfaces and the equipment plates. Two separate image stacks, representing masks for the top pin and bottom plate respectively, were extracted. This segmentation was performed using slice-by-slice 2D region growing on a data set where the sample was not in contact with the top pin/bottom plate (load step 1 in Fig. 2a for the top pin). The segmented regions grew from defined seed points using the X-ray phase effects between the background and the pin/bottom as a stopping criterion. In a few 2D slices, local problems were encountered with leakage into the background caused by a locally blurred interface between the equipment and the background. Such regions were handled with morphological post-processing of the masks. Some final manual editing of the masks was performed in Fiji ImageJ, where some image slices with severe problems of region growing leakage were replaced with the neighboring image slice. From these segmented images, the top pin and bottom plate surfaces could be seen to be not perfectly flat and aligned at the microscale (see Fig. 2, which, together with the surface roughness of the sample, prohibited the use of simpler methods to extract the sample volume, e.g. a rectangular ROI).
The second segmentation step consisted of positional adjustments of the equipment masks to each sample. Sample changes between each experiment resulted in slightly different translation of the bottom plate as well as both rotation and translation of the top pin. Therefore, the two equipment masks were translated and rotated, then exported to match the first 3D image in each experiment.
In a third step in the segmentation, time vectors describing how the loading equipment translated vertically between each successive 3D image during each experiment were determined. While the displacement of the bottom plate was measured with the loading device, the top pin was not fixed rigidly and, therefore, moved vertically at certain time instances.
The equipment masks and time vectors for each experiment were used in the fourth segmentation step, where a mask representing the sample volume was extracted for each time step. A manually defined threshold (higher than in Fig. 2b) was first used to extract the paperboard fibres with the highest intensities since these, on average, had higher intensity Fig. 2 Example images illustrating the results of the binarization procedure. a Tomography images at before loading and at maximum compression. b Binary images from direct threshold-based segmentation of the fibres. c Final results after the stepwise segmentation procedure starting with extraction of equipment masks and sample volumes followed by fine-tuned fibre-and pore binarization compared to both the background and the polymer equipment. The resulting binary images did, however, also contain X-ray phase effects at the interfaces between the paperboard sample and the measurement equipment. These effects were removed for each time step by translation of the equipment mask according to the time vectors. Thus, the actual interface between the paperboard surface and the equipment was preserved at each time instance. The resulting fibre segmentation was, however, not accurate enough for quantification of material properties since only the brightest fibres were segmented. Furthermore, the pore space could not be separated from the background in this step. 3D morphological operations including erosion and dilation in several steps were, therefore, performed to fill the sample volume, slightly smoothing the outer sample surface and removing remaining random noise distributed across the binarized image stacks.
The fifth and last segmentation step consisted of a fine-tuned binarization of fibres and pores within the sample volume. In this step, the raw images were binarized again using an optimised manual threshold to separate fibres from pores (the same threshold as in Fig. 2b) whereafter the background and equipment were removed with the previously extracted sample volume mask. Figure 2c shows examples of the final results from the stepwise segmentation procedure where the equipment, fibres and pores are uniquely labelled and separated from the background.

Quantification of material properties and boundary conditions
The segmented images were used to calculate the evolution of several bulk sample properties; total sample volume, pore volume, fibre volume and porosity. Spatial 2D distributions of solid fraction in the sample out-of-plane direction were also quantified for each time step. For each pixel in the out-of-plane 2D direction, the proportions of voxels labelled as fibres with respect to the total number of voxels along the in-plane direction were calculated.
The segmented images were also used to quantify the contact area between the sample and the loading equipment at each time instance. The approach was to dilate the sample and the equipment masks whereafter the overlapping region was skeletonized. The total contact area was defined conservatively as the number of voxels where the sample was in contact with both the top pin and the bottom plate.

Digital volume correlation and strain field analysis
For full-field 3D quantifications of the evolution of the strain tensor field during the experiments, DVC analyses and strain calculations were performed using the open-source Python code Spam 0.5.2.1 (Stamati et al. 2020). The DVC analysis code compares 3D images from two different time steps and finds the relative 3D displacement of nodes evenly distributed across the image. Because of the inherent material texture and the size of the correlation window, we could downsample the reconstructed and postprocessed input image stacks to 25% to save computational time and RAM usage without affecting the resolution of the results. We used a node spacing of 64 × 64 × 32 µm 3 , which corresponds to 10 × 10 × 5 voxels in the downsampled images. The size of the correlation window is controlled by the inherent texture in the images and a correlation window width of 71 µm (11 voxels) was found to be optimal for our data. Larger correlation windows resulted in pixelated results while smaller windows led to noisier results.
The 3D displacement fields from the DVC analysis were used to define the full-field 3D strain fields in Spam. 3D deformation tensors, as well as the volumetric and deviatoric strain invariants, were computed within the framework of large strain theory in Spam. In this work, we computed the Hencky strain tensor components from the right stretch tensor U exported from the software: An initial analysis of the results showed that the 3D strain tensor is strongly dominated by the outof-plane normal component, zz , (which is the same direction as the applied load). In the following, the results are mainly visualized using the volumetric strain invariant. The evolution of strain fields was analysed over the entire experiment, both between each consecutive time step (incremental strain analysis) and in relation to the first (unloaded) time instance (total strain analysis).

Results
Initial average sample thickness and porosity were calculated from the image data (first time step) and are reported in Table 2. The thickness was larger in 035-BK compared to the other samples. The initial porosities were approximately the same in the four samples, just slightly lower in the WR samples. There are some uncertainties in the exact magnitudes of the porosity values since they are dependent on the threshold value used for the image processing as well as the resolution of the tomography images. However, for the following data analysis, only the relative changes over the loading will be considered.

Macroscopic stress-strain responses
Macroscopic stress-strain curves, shown in Fig. 3, were calculated from the measured force-displacement data. The cross-sectional area of the top pin was used to calculate engineering stress, while the initial average gauge lengths in each data set were obtained from the tomography data and used to calculate true (logarithmic/Hencky) strain. True stress was also calculated using the actual contact area at each load step quantified from the images (results of the contact area evolution are reported below). The numbered labels represent the positions on the loading curves where X-ray tomography data were acquired, hereafter referred to as time steps. In 032-WR and 035-BK, the loading of the sample began at time steps 4 and 5, respectively.
The stress-strain responses in Fig. 3 show that the WR-samples are stiffer than the BK-samples and there are small differences between the two samples of each type. For three samples in Fig. 3, the initial parts of the loading curves follow the typical exponential/J-shape commonly observed for paperboards (Nygårds 2008;Chen et al. 2020), whereafter the stress-strain response becomes linear up to peak load. However, the material behaviour for 038-WR appears to be linear during the entire loading curve, which might partly be caused by a too large preload on the sample with the initial, manual approach of the top pin. Calculated permanent strain values remaining in the samples after the load-unloading cycle are reported in Table 3 together with the tangent modulus calculated for different parts of the curves. Due to the more similar macroscopic responses of 032-WR and 033-BK, in terms of strain at maximum compression, these will be in focus in the following presentation of results.

Microscale boundary conditions
Detailed information on microscale boundary conditions of uniaxial compression experiments are usually not available and, therefore, not considered for the interpretation of macroscopic mechanical responses. Such boundary conditions have, however, been quantified from the X-ray tomography data in this study with the aim of evaluating their significance for the results. Figure 4a shows topological maps of the loading equipment, calculated from the segmented masks of the upper pin and bottom plate. It is clear that the equipment surfaces are not perfectly flat at the microscale. Most noticeable is the local cavity in the top pin and the locally elevated centre of the bottom plate. The relative elevation of the top pin decreases from the centre and radially outwards. These microscale variations mean that the sample is heterogeneously loaded, which will affect the outer shape change of the sample. Figure 4 b shows, as an example, in-plane variations in sample thickness at different time steps calculated from the segmented images of 032-WR. The sample is originally slightly thicker in the northeastern part (time step 1). At maximum compression (time step 15), the sample thickness is affected by the topology of the equipment. The minimum thickness is found at the location where the highest relative elevations of the top pin and the bottom plate meet. The sample is also generally thinner in the bottom part, compared to the top part, due to the non-perfect horizontal alignment of the bottom plate. After unloading to 0 N (time step 40), it is clear that there is a remaining imprint caused by the loading equipment. Figure 4 c shows corresponding data for 033-BK. The compressed zone is larger in the BK paperboard than the WR paperboard at maximum compression (time step 14) and after unloading (time step 35). Microscale heterogeneities of the loading equipment surfaces will also cause a changing equipmentsample contact area during the load-unload cycle, which affects the macroscopic stress values. The effects of considering the actual contact area for the macroscopic stress-strain curve during uniaxial compression of copy paper was previously analysed by Chen et al. (2020). They used carbon paper to estimate the actual contact area at different load steps. The contact area increased with loading, which was attributed to sample surface roughness effects, and the use of the measured contact area values for the macroscopic stress calculations had a major effect on the magnitude and shape of the stress-strain curve (Chen et al. 2020).
The quantified evolutions of the contact areas between the loading equipment and the samples from this study are plotted in Fig. 5a. As expected (given the shape of the equipment in Fig. 4 a), the contact area for all samples gradually increased with loading. The samples did not, however, not lose contact with the equipment at the same rate during unloading. The images in Fig. 5b show, as an example, 2D visualizations of sample-equipment contacts for the labelled load steps in Fig. 5a. At load step 1, the 032-WR sample rests on the elevated part of the bottom plate while not in contact with the top pin. The net (conservatively calculated) contact area is therefore 0.
During loading, the top pin gradually comes in contact with a larger, central part of the sample (load step 8) and at maximum compression (load step 15), the top pin is in (more or less) full contact with the sample. The data in Fig. 5a were used to calculate true macroscopic stress and the results are shown with 3D strain field evolution An analysis of the 3D strain tensor components showed that the deformation is strongly dominated by the out-of-plane normal component e zz (although some minor shearing is present in the e zx and e zy components, see Supplementary Information online). In the following, 3D strain fields are shown in terms of the volumetric strain invariant. Orthogonal 2D slices through the 3D fields of the total volumetric strain invariant at three different load increments are visualized in Fig. 6 for 032-WR. Increment 3-7 corresponds to the flatter increase in stress with strain in the macroscopic curve (Fig. 3). The first in-plane strain localization occurs where the local distance between the loading plates is at a minimum (compare load increment 3-7 in Fig. 6 with Fig. 4). In the outof-plane direction, this localization is located in the middle/bottom half of the sample. At maximum compression (load increment 3-15), the in-plane strain distribution essentially reflects the shape of the top pin (Fig. 4); overall, the amount of strain decreases radially outwards from the in-plane centre of the sample, and the point of minimum distance between the equipment plates, as well as the local cavity of the pin, are visible in the strain fields. The strains are well-distributed in the out-of-plane (thickness) directions from the top pin and through the centre of the sample. However, a thin layer close to the bottom plate appears to be relatively undeformed. After the load-unload cycle (load increment 3-40), permanent strain remains in the in-plane centre and out-of-plane upper part of the sample. Figure 7 shows corresponding data for 033-BK. The general patterns of strain field distributions are similar to the WR sample (Fig. 6) at the first strain localization (increment 1-3), at maximum compression (increment 1-14) and after unloading (increment 1-35). The strain fields are, however, consistently stronger and more widespread across the BK samples. The results in Figs. 6 and 7 confirm that local heterogeneities of the loading equipment have a major impact on the net development of 3D strain field distributions during uniaxial compression. The magnitude and in-and out-of-plane extent of the strain fields are, however, also strongly dependent on the out-of-plane structure of the sample material.
An incremental analysis of volumetric strain fields represents the relative change between consecutive load increments and, therefore, contain more detailed information on the deformation behaviour at specific time instances. Figure 8 shows volumetric strain field distributions for 032-WR across five selected increments marked on the macroscopic stress-strain curve. The strain fields for the three different increments along the loading curve show that, in the outof-plane direction at load increment 5-6, the volume strain localized along the top-and bottom-surfaces of the sample. At later increments (increment 9-10), the main localization is centred in the out-of-plane direction, but extends to the upper sample surface around the edges of the top pin. At load increment 14-15, the strain is fully localized in the out-of-plane centre of the sample. This evolving strain accumulation in the out-of-plane centre can also be seen in the in-plane slice through the strain field with a localized in-plane region first in the sample centre that  More details on the temporal development of the out-of-plane strain distribution are shown in Fig. 9 where data from the 12 first load increments are shown. These data show that the sample was initially compressed close to both sample surfaces when the applied stress was low. When the bottom plate was displaced further and the stress on the material increased, the main strain localization gradually shifted to the central part of the sample. This gradual shift can be observed in e.g. load increment 9-10, where the material in the in-plane centre is compressed with a higher stress compared to the material in the outer regions of the sample (due to varying contact area with the top pin). In the in-plane central parts of the sample, the strain fields are localized in the out-of-plane centre. Outwards in the in-plane radial directions, the strain distributions follow the shape of the top pin towards the upper sample surface. Around the edges of the top pin, small amounts of out-of-plane shearing also occur (see Supplementary Information online). Figure 8 shows that the outer in-plane regions of the sample expanded first during unloading, i.e., regions that had been compressed for the shortest amount of time/with the lowest local stress. The expansion at load increment 19-20 was diffuse, but localized to the outer in-plane regions and out-ofplane centre of sample. The last major change in the sample during unloading occurred at load increment 36-37. At this increment, the compressive strain decreased in the same (extended) in-plane region where the strain first localized in increment 5-6. Some material close to the bottom plate also expanded strongly over the same increment. Figure 10 shows a corresponding incremental analysis of the BK sample. As already visible in the macroscopic stress-strain curves, the amount of strain in the BK sample was higher compared to WR over a given stress increment (compare with Fig. 8). A relatively large in-plane localization can be seen to have occurred in the out-of-plane centre of the sample during an early load increment (increment 2-3). Like in the WR sample, strain also localized close to the sample surfaces in the out-of-plane direction. During load increment 7-8, the strain was distributed in a relatively wide region in the out-of-plane centre of the sample. At a later load increment (13-14), the strain distribution was more diffuse and the strongest strain values are seen in a thin layer in the out-of-plane centre. During unloading, the BK sample behaved in a similar way to the WR sample; the expansion started from the outer regions of the sample and the last major change in the strain distribution is related to the expansion of an extended region around the location where the first strain localization appeared, as well as close to the bottom plate in the out-of-plane direction. Figure 11 shows, as for WR in Fig. 9, temporal details on the out-of-plane development of the strain fields during the first 12 load steps. The shift of the main strain localisation to the out-of-plane centre occured quicker in BK compared to WR. After establishing in the out-of-plane centre, the main strain localization became gradually thinner with increasing load increments.
Evolution of sample properties. Figure 12a-c shows the evolution of sample-, fibreand pore volumes relative to the initial values, as functions of load step for all four samples. Fibre lumens are included in the total pore space volume whereas the fibre volume only represents the fibre wall material. As expected, the total sample and the pore space volumes decreased during compression (Fig. 12a, c). The data in Fig. 12b show that fibres were also compressed, with a magnitude that is around half of the pore space compaction. Figure 12d shows the porosity evolution in each sample (calculated from the pore-and sample volume data). These data show that the relative decrease in porosity was larger in BK samples than in WR samples at maximum compression (~ 15% and ~ 8%, respectively). After unloading, the BK samples recovered to ~ 90-92% of their initial porosities while the corresponding figure for the WR samples is ~ 96%.
It is of interest to examine the dependence of the sample properties in Fig. 12 on the displacement of the loading equipment. Figure 13a-d shows the sample properties plotted against displacement calculated from the images. Mostly, the data in Fig. 13a-c show almost linear variations with applied displacement at larger displacement values, which means that the evolution of the properties are closely related to the movement of the equipment plates. The curves do, however, deviate from this linear behaviour at low displacement values where the rate of volume change decreases significantly for the pore volume and the total sample volume.
A noticeable feature in Fig. 13a is the slower rate of decrease of sample volume in the beginning of the loading of 032-WR. The most likely explanation is that this sample was oriented with the clay coated surface (which is more resistant to deformation than the opposite side) downwards on the bottom plate. The more easily deformed surface of the sample was in contact with the bottom plate for the other three samples, which could explain why the pore space reduced at a higher rate; the sample-bottom plate contact area was generally larger than the sampletop pin contact area for most time steps. The small contact area between the top pin and the deformable surface of the paperboard during the early phase of the loading is a likely cause of the slower decrease in sample volume for 032-WR. A similar pattern can be observed for the same sample in the evolution of the pore volume with displacement (Fig. 13c), but, interestingly, not as clearly in the evolution of the fibre volume (Fig. 13b).
Linear curves were fitted to the data in Fig. 13b, c to estimate the rate of decrease of the pore and fibre volume in the samples during loading. The rate was evaluated over the first 4-5 load steps and over the range from load step 4 or 5 (depending on the sample) up to maximum compression (see Supplementary Information online, figures A7 and A8). The results show that the fibre and pore volume decrease rates were systematically lower for 032-WR than the other samples, likely because of the different sample orientation during the tests. For 038-WR, 033-BK and 035-BK, fibre and pore volumes decreased at approximately the same rate (1.3-1.4 mm 3 /mm and 1.3-1.6 mm 3 /mm, respectively) over the first 4-5 load steps. After this, the pore volumes decreased at a higher rate than the fibre volumes (2.9-3.5 mm 3 /mm and 2.3 mm 3 /mm, respectively). The fibre volume decrease rate was, in fact, the same for 038-WR both BK samples (2.3mm 3 /mm) after load step 4/5, while the pore volume decrease rate differed with 2.9mm 3 /mm in 038-WR and 3.3-3.5 mm 3 /mm in BK samples.
The results presented above suggest that uniform compression of the sample material, as a whole, occurred in the very first stages of loading when the applied force was relatively low. In this stage, fibres and pores were compressed at approximately the same rate. During a second stage, the pore volume started to collapse at a higher rate than the fibre volume. The pore volume decrease rate is higher in BK samples than WR samples, while the fibres compress at the same rate in both types of paperboards.  (Fig. 14a-b). At higher stresses, 032-WR and 035-BK show a tendency of approaching a shallower gradient in the rate of sample-and fibre-volume decrease (Fig. 14a, b, ~ 8 MPa to ~ 12 MPa). In the same samples during the first part of the stress release (~ 12 MPa to ~ 10 MPa), the volumes continue to decrease slightly. The reason for this decrease is that the displacement is still increasing over increments 15-19 in 032-WR and 19-24 in 035-BK (Fig. 13). This means that stress is being released in this last stage of sample compression, which suggests damage or microstructural-failure mechanisms are activated in the fibre networks.
The volumetric recovery of the sample, fibres and pores during unloading is a slow process (Fig. 14a-c).
In all samples except 038-WR, exponential recovery occurs between ~ 2 MPa to 0 MPa. The exponential behaviour is more pronounced in the BK samples. The behaviour of 038-WR clearly deviates from the other samples in Fig. 14a, b. With the same amount of stress, the sample-and fibre volumes of this sample decrease to a much lesser degree, which indicates that the fibre network or the fibres were stiffer in this sample. The load and unload behaviours of this sample are more linearly correlated with the applied stress compared to the other samples. When unloading began, 038-WR immediately started to expand, but at a slower rate compared to the compression rate.
The behaviour of the pore volume development (Fig. 14c) shows similar trends in both WR samples (unlike the fibre volume development). The evolution of the pore volume with true stress (Fig. 14c) is smooth and systematic for all samples. The loading curves are clearly steeper for BK samples, for which the pore volume reduced more with the same applied stress. In contrast to the fibre volume development, Fig. 12 Relative evolution of sample-and fibre-and pore volumes as well as porosity during the loadunload cycles there are no clear signs of shallower gradients of the curves during loading at high stresses. During unloading, the pore volume increased very slowly in all samples until an exponential recovery from ~ 2 MPa until 0 MPa. Towards the end of the load-unload cycle, the BK curves appear more exponential shaped than the WR curves.
The porosity evolution in Fig. 14d shows very similar behaviours for the two BK samples. The decrease in porosity occurred relatively fast during loading, especially between ~ 0 to ~ 4 MPa, and is close to linear. During unloading, the recovery is exponential from ~ 8 to 0 MPa. In the WR samples (Fig. 14d), the porosity decrease rate is slower during loading (compared to the BK samples) and the recovery during unloading is linear between ~ 10 to ~ 1 MPa.
In Fig. 15, the pore volume is plotted against the fibre volume for all samples. The relationship between these properties can be approximated as linear (possibly apart from 032-WR), which indicates that material properties independent on actual loading conditions (force and displacement) can be extracted from the data. Linear curves were fitted to the data points belonging to the first and second stage of loading, that is from the beginning of the loading until load step 4/5 and from load step 4/5 until maximum compression (for 032-WR, the last two datapoints during loading were excluded). The slope and R 2 -values of the linear functions are reported in the figure. During the first stage of loading, the slope is close to 1 in 038-WR and 035-BK while it is lower in 033-BK and insignificant in 032-WR (R 2 = 0.66). During the second stage, the slope is greater than 1 and larger in the BK samples than in the WR samples. During unloading, the BK samples follow the same linear curve while the WR samples unload with a similar slope (but shifted on the horizontal axis). The results have so far shown where within the sample deformation occurs in the different load stages, as well as how sample scale properties evolve during the process. The questions about how the micro structures evolve during loading remain; answering such questions could lead to a better understanding of the mechanisms of material deformation. For this purpose, the out-of-plane variations in solid fraction at three different time steps are shown in Fig. 16 for 032-WR and Fig. 17 for 033-BK. The differences in out-of-plane structures are evident before loading (time steps 4 and 1, respectively). The WR paperboard ( Fig. 16) consists of thin, alternating layers of relatively low and high solid fraction throughout the sample thickness. In a thin layer closest to the claycoated sample surface at the bottom of sample, solid fraction values can be seen to be slightly elevated. In the BK paperboard (Fig. 17), the 3-ply structure is clearly visible; relatively thick layers of higher solid fraction can be observed close to both sample surfaces and in the sample centre. These layers are separated by two relatively thick layers of low solid fraction.
At maximum compression in 032-WR (Fig. 16, time step 15), the thin layers of lower solid fraction were completely closed in the (in-plane) central part of the paperboard. The layered structure is still visible in the outer regions where lower local stress was experienced due to the shape of the top pin (Fig. 4). After unloading (time step 40), the layered structure was partly recovered, especially in the outer regions. Figure 17 shows a similar behaviour for 033-BK. At maximum compression (time step 14), all layers of low solid fraction were closed in the central parts of the sample and the structure appears more homogeneous. After unloading (time step 35), the layered structure was partly recovered.

Discussion
The results of this study have revealed detailed data on uniaxial paperboard compression tests that are unachievable with macroscopic testing, such as the development of the actual contact area between the loading equipment and the sample, strain localisation and the microstructural evolution. Concerning the contact area, the conversion from engineering to true stress with the actual contact area evolution, derived from the images, did, however, only change the magnitude of the stress values and had no dramatic effect on the shape of the stress-strain curves, in contradiction to previously reported results for copy paper (Chen et al. 2020).
For all of the samples, the 3D strain fields derived from the tomography images by DVC generally show the highest levels of strain in the outof-plane central part of the sample during compression, except during the very first stage of the loading when the strain primarily accumulated close to both sample surfaces. The observed in-plane distributions of volumetric strain in the out-of-plane centre of the samples are closely dependent on local imperfections of the loading plates. The first localization in the out-of-plane centre was seen to occur in each case where the distance between the top and bottom plates was minimum, whereafter the strain localization grew radially outward during the compression as the contact area grew.

Compression mechanisms
It is well known that uniaxial compression of paper materials leads to both reversible and irreversible thickness changes after unloading. With pulse loading experiments, Rättö and Rigdahl (1998) quantified creep and recovery behaviour of paper and observed a significant amount of permanent strain remaining in the samples after the stress pulses. Kananen et al. (2018) loaded paper handsheets in several compression cycles to remove irreversible effects, the main part of which occurred after the first cycle. After this, a curve fitting procedure suggested that the magnitude of the reversible strain is primarily determined by the initial porosity (Kananen et al. 2018). In our study, both the macroscopic load-unload curves and the total volumetric strain field evolution showed that a large part of the accumulated compressive strain is recovered when the applied force is removed. However, significant permanent strain remains after the load-unload cycle. A central question of interest, in this context, concerns which mechanisms lie behind the reversible versus irreversible compression of paperboard. The results from this study suggest that elastic energy stored during volumetric reduction of fibre walls might drive the reversible expansion of the fibre network during unloading. Our results show that both fibre and pore volumes are reduced during uniaxial compression (where the former quantity represents the volume of the fibre walls with lumens excluded). Compression of fibre walls has, to the best of our knowledge, not been suggested as a compression mechanism previously. With compression of fibre walls, elastic energy can be stored in the fibre network and the data show that at least part of the reduction in fibre volume is reversible. It is likely that also the reversible recovery of the pore space is driven by the expansion of compressed fibre walls, together with other possible, but experimentally unverified, mechanisms, such as, e.g., fibre bending.
For stresses up to ~ 8 MPa, the decrease rate of the fibre volume with true stress was equal in three out of the four samples. This shows that the fibres respond in the same way to applied stress, despite the different out-of-plane structures of WR and BK paperboards. In two of the samples (032-WR and 035-BK), the fibre volume decrease appeared to approach a plateau, whereafter the stress dropped despite further compression. This indicates that permanent deformation or damage of the fibre network occurred at higher stresses. A similar tendency of shallower gradients at high stresses was not observed for the pore volume as a function of true stress. It is not possible to determine from our data what kind of permanent changes to the fibre network occurred, but irreversible fibre deformation or fibre bond damage are likely hypotheses. According to Rodal (1989), a shift from linear elastic to elastoplastic behaviour in paper materials typically occurs at strain magnitudes of one to a few percent during compression. The elastoplastic behaviour is hypothesised to be a result of fibre network buckling and precedes a third regime where fibres start to collapse and stress increases very rapidly with strain (Rodal 1989). Such a third regime was not observed in the macroscopic stress-strain curves in this current work.
In one of the samples, 038-WR, the fibres were much stiffer compared to the other samples (Fig. 3, tangent modulus in Table 3). The volumetric change of the fibres as a function of true stress (Fig. 14) was much smaller, for this sample, while the pore volume behaved in a similar way as compared to the other three samples. Other differences were that the contact area with the loading plates decreased faster during unloading (Fig. 5a) and the amount of permanent (macroscopic) strain after unloading was significantly lower in 038-WR (Table 3). It is possible that the stiffness of the fibres partly prevented permanent fibre network damage to occur in this sample.
Although the fibre volume decreases significantly during compression (10-20% of their initial volumes at maximum compression), the largest net changes in volume occur in the pore space, which reduces twice as much as the fibre volume.
The solid fraction evolution (Figs. 16 and 17) showed that the main pore space reduction process appears to be closure of layers between fibres. The process is most clearly evident in the 3-ply paperboard, where two layers of lower solid fraction are clearly visible before the load-unload cycle. At maximum compression, these layers are barely visible in the central parts of the sample, but they recover partly during unloading. The initial out-of-plane structure of the 1-ply paperboard consist of thin and densely spaced layers of low solid fraction. These layers are interpreted as voids between layers of fibres, which are reduced during loading and are barely visible at maximum compression. Similar to the 3-ply structure, these layers partly recover after the unloading of the sample.
In general, we conclude that the out-of-plane structure of the studied paperboards has a significant influence on the compressive behaviour, as the WR and BK samples in our study had similar initial porosities at the sample scale. The initial presence of two layers with lower solid fraction within the 3-ply paperboard led, compared to the 1-ply paperboard, to greater and more distributed volume reduction in the out-ofplane direction and a larger global porosity reduction during compressive loading. Permanent strain also remained over a larger in-plane region compared to the 1-ply paperboard.

Evolution of deformation mechanisms during compression
Surface roughness effects leading to surface deformation have previously been suggested as a mechanism responsible for the flat increase in macroscopic stress with applied strain in the early phase of uniaxial compression (Stenberg et al. 2001). In this study, a shift in the out-of-plane localization of volumetric strain fields was confirmed in the early part of the loading process. Over the first four load steps in Figs. 9 and 11, the strain was seen to accumulate mainly close to the sample surfaces. These accumulations have large in-plane extensions despite low sample-equipment contact areas (Fig. 5a). In each case, the entire sample reacts to the low applied force during the beginning of the loading, not only the parts of the sample that are in direct contact with the equipment. This suggests that the surface deformation is not exclusively determined by sample surface roughness, at least not in single-or multiply commercial paperboards with coating. After the first load steps, the main deformation region shifts to the out-of-plane centre of the sample. That is, material close to the surfaces is deformed first whereafter the deformation shifts to the out-of-plane centre.
The rate of pore and fibre volume decrease as a function of displacement indicated a transition between (at least) two deformation phases: an initial phase (first 4-5 load steps) where the fibre and pore volumes decreased at the same (slower) rate, followed by a second phase where the decrease in pore volume occurred at a faster rate than the decrease of fibre volume. This suggests that when the strain localization shifts to the out-of-plane centre of the samples, the pore volume starts to decrease faster than the fibre volume. In future work, this phenomenon can be studied in detail though assessment of local poreand fibre-volume changes within these shifting strain localization zones.
In summary, the combined results of this work have resulted in new insights into changing material behaviour along the stress-strain curve. The DVC results showed that with low applied global stress, sample material close to both equipment-sample contact surfaces deforms. This stage also corresponds to a phase with a uniform rate of fibre-and pore volume compression. When the stress increases and the sample becomes more strained, the main region of deformation shifts to the out-of-plane centre where the pore volume is reduced at a higher rate than the fibre volume. The interplay of these two mechanisms and the gradual transition from one to the other appears to cause the changing slope of the macroscopic stress-strain curve during uniaxial compression.

Conclusions
This study has exploited synchrotron microtomography with high spatial and temporal resolutions to study the details of the microscale variations in boundary conditions and the microscale deformation mechanisms occurring during uniaxial compression of paperboard and to relate these to the macroscopic stress-strain response. Two types of commercial paperboards with different out-of-plane structures (1-ply and 3-ply) were tested in load-unload cycles. The combined results of the work provide new insights into the significance of varying microscale boundary conditions on uniaxial compression as well as experimental evidence of material changes occurring during compression. Furthermore, the results show that a gradual transition between different deformation mechanisms occur along the stress-strain curve.
3D volumetric strain field analysis shows material close to the sample surfaces is deformed during the initial stages of the loading, whereafter the deformation zone shifts to the out-of-plane sample centre. In a first loading phase, fibre-and pore volumes compress uniformly at the same rate while the pore space is compressed at a higher rate in a second phase. The fibre volume decrease rate is equal in the two paperboards studied (WR and BK), while the pore space decrease rate is greater in the 3-ply BK paperboard compared to 1-ply WR paperboard. The differences in the out-of-plane structure of the two boards explain why the net amount of compressive strain and reduction in global porosity was larger in 3-ply BK samples than 1-ply WR samples.
A key result from the work is the experimental verification of fibre wall volume reduction during uniaxial compression of paperboard, which has, to the best knowledge of the authors, not been discussed in previous research. This mechanism indicates that elastic energy can be stored in the fibres during compression that could drive the elastic recovery during unloading. The microscale 3D imaging reveals that the uniaxial compression of the paperboard samples causes reduction of both fibreand pore volumes, which is, to some extent, recovered during unloading, but neither retain their initial values immediately after the applied force has been removed (longer term recovery effects have not been studied). The mechanism of pore space reduction has been observed to be primarily by reduction in the spacing between fibre layers and this effect was observed in both the 3-ply and 1-ply paperboards studied. The net reduction of pore volume was observed to be around twice as large in magnitude compared to the net fibre volume reduction and the approximately linear relations found between these two parameters can a useful addition to models describing the deformation of paperboard under compression.
Funding Open access funding provided by Lund University. This work was funded primarily by Tetra Pak AB and by VIN-NOVA, Formas and Energimyndigheten through their common strategic innovation programme BioInnovation (VINNOVA Project Number 2017-05399). Additional funding was from the 'FORMAX-portal-access to advanced x-ray methods for forest industry', VR Project (Number VR 2018-06469).

Data availability
The original and processed datasets analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
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 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/.