3D Strain Field Evolution and Failure Mechanisms in Anisotropic Paperboard

Experimental analyses of the 3D strain field evolution during loading allows for better understanding of deformation and failure mechanisms at the meso- and microscale in different materials. In order to understand the auxetic behaviour and delamination process in paperboard materials during tensile deformation, it is essential to study the out-of-plane component of the strain tensor that is, in contrast to previous 2D studies, only achievable in 3D. The main objective of this study is to obtain a better understanding of the influence of different out-of-plane structures and in-plane material directions on the deformation and failure mechanisms at the meso- and microscale in paperboard samples. X-ray tomography imaging during in-situ uniaxial tensile testing and Digital Volume Correlation analysis was performed to investigate the 3D strain field evolution and microscale mechanical behaviour in two different types of commercial paperboards and in two material directions. The evolution of sample properties such as the spatial variation in sample thickness, solid fraction and fibre orientation distribution were also obtained from the images. A comprehensive analysis of the full strain tensor in paperboards is lacking in previous research, and the influence of material directions and out-of-plane structures on 3D strain field patterns as well as the spatial and temporal quantification of the auxetic behaviour in paperboard are novel contributions. The results show that volumetric and deviatoric strain, dominated by the out-of-plane normal strain component of the strain tensor, localize in the out-of-plane centre already in the initial linear stress-strain regime. In-plane strain field patterns differ between samples loaded in the Machine Direction (MD) and Cross Direction (CD); in MD, strain localizes in a more well-defined zone close to the notches and the failure occurs abruptly at peak load, resulting in angular fracture paths extending through the stiffer surface planes of the samples. In CD, strain localizes in more horizontal and continuous bands between the notches and at peak load, fractures are not clearly visible at the surfaces of CD-tested samples that appear to fail internally through more well-distributed delamination. In-plane strain localization preceded a local increase of sample thickness, i.e. the initiation of the delamination process, and at peak load, a dramatic increase in average sample thickening occurred. Different in-plane material directions affected the angles and continuity of the in-plane strain patterns as well as the sample and fracture properties at failure, while the out-of-plane structure affected how the strain fields distributed within the samples.


Introduction
Improved understanding of the mechanical behaviour of paperboard materials is important for the packaging industry. This knowledge is used in the industrial process to reduce both the amount of wasted material and the risk of defects in final packages. In the food package industry, it is crucial to avoid defects in packages, as these might lead to reduced shelf-life of the food. To analyse the mechanical behaviour of paperboards in different steps of the package manufacturing process, continuum modelling can be used. This approach offers the possibility to identify and adjust crucial steps in the process where different types of defects are likely to appear, from creasing and folding to sealing of filled packages. Increasingly detailed and complex models are being developed that require higher fidelity experimental characterization for model validation.
Experimental mechanics traditionally deals with the extraction of material-specific mechanical parameters that can be used to build and calibrate continuum models, and to identify general mechanical behaviour of the sample, for example the type of material failure. Paper and paperboards can be described as viscoelastic plastic materials [1]. In-plane tensile tests of paperboards (ISO 1924-3) show a typical macromechanical behaviour of an initial linear phase (elasticity) followed by non-linear hardening (plasticity) until failure occurs with a drop in stress to zero. A gradually descending postpeak behaviour, that can be described as cohesive failure, is present in short samples [2,3]. However, despite a relatively well-understood macro-mechanical behaviour of papers and paperboards, challenging aspects of these materials is their anisotropic nature, through thickness heterogeneity and variability at the microscale.
Paperboards generally consist of fibres oriented in a preferred in-plane direction, referred to as the Machine Direction (MD). Since their structure is built up by either one or several layers (plies) with or without a thin top layer of printenhancing coating, paperboards exhibit properties differing by several order of magnitude in the out-of-plane (thickness) direction (ZD) compared to the in-plane direction. The material therefore behaves differently when tested in different material directions. Furthermore, microstructural variations within the same type of paperboard and material direction leads to variability in the mechanical behaviour of the same sample type. Due to these challenges, more knowledge is needed about the micro-and mesoscale mechanical mechanisms leading up to material failure. This requires more advanced, nontraditional experimental approaches that reveal the heterogeneity of the material response and the deformation mechanisms at a pertinent scale.
3D visualisation and quantification of material structures at a microscale is achievable through X-ray tomography. Furthermore, an effective tool used to quantify the development of strain fields in 3D based on X-ray tomography data is Digital Volume Correlation (DVC). With DVC, displacement of regions within the sample during loading is quantified and the full 3D strain tensor field is obtained through a continuum mechanical approach. DVC can, thus, be described as a mesoscale experimental approach and its usefulness for a better understanding of mechanical deformation-and failure mechanisms have previously been demonstrated for a range of different materials. Previous work includes, for example, the application of DVC to bone samples [4], geological materials such as sand and shale [5,6] and wood [7]. The application of the DVC method to paperboard materials and comprehensive analyses of the 3D development of the full strain tensor during loading is, however, not common, although it is highly relevant since the results can be used to analyse, evaluate and improve the performance of paperboard continuum models. For example, 2D experimental (in-plane) strain field mapping was used to evaluate the performance of a continuum-based constitutive paperboard model in Borgqvist et al. [8].
2D strain field localization patterns in different kinds of paper and paperboard samples subjected to tensile loading have been in focus in several previous studies. For example, plastic strains in silicon treated linerboard and newsprint paper samples were studied optically by Korteoja et al. [9], while Borodulina et al. [10] and Hagman and Nygårds [11] used digital image correlation (DIC) to investigate strain fields in isotropic handsheet paper and 3/5-ply paperboard samples, respectively. In several of these studies, diagonal in-plane strain field patterns were observed during the plastic phase of the loading. In Korteoja et al. [9], plastic strain patterns could be observed optically in notched Cross Direction (CD) samples while they were not distinguishable in MD-tested samples. At a microscopic scale, the observable strain patters in the CD-tested samples appeared to develop as micro failures around horizontally aligned fibres. These micro failures collectively created diagonal patterns across the samples at a macroscopic scale. Hagman and Nygårds [11] observed partly diagonal macroscopic strain streak patterns across rectangular 3-ply/5-ply paperboard samples in both MD-and CD-tested samples. The main difference between samples tested in MD and CD was that the strain patterns were weaker and more diffuse in MD. Borodulina et al. [10] described the damage localization in isotropic handsheet paper samples as always occurring at an angle (around 45°) to the loading direction, reminiscent of semi-ductile shear-band formation. Prior to failure, diagonal patterns of in-plane strain localization were visible in corresponding zones.
In-plane strain field localization patterns are generally believed to be controlled by the initial structure of the paper or paperboard, e.g. local variations in inter-fibre bond density and fibre orientations [10]. Hagman and Nygårds [12] concluded that in-plane strain field localization is determined by variations in the initial formation of the paper sample, a relationship that has been discussed also in previous studies (e.g. [13]). Hagman and Nygårds [12] suggested that zones with elevated plastic strain and deformation occur in areas of locally lower material density. It is often assumed that strain field localization starts to develop in the plastic (hardening) regime [10][11][12] which could indicate that it represents irreversible damage, e.g., broken fibre bonds or micro failures. Hagman and Nygårds [12] showed that elevated strain (in the loading direction) is associated with similar patterns of locally increased temperature during plastic loading indicating that energy is dissipated locally due to plastic deformation or evolving damage (e.g. fibre movements or breakage of fibre bonds). Measurements and analysis of acoustic emission during uniaxial tensile tests also suggest that fibre bonds start to break when the loading curve becomes nonlinear [14]. Fibre network modelling [10] suggests that fibre bond failures localize in areas of high local strain and not the other way around. Their results also support the already established belief that the non-linear response of the fibre network is caused by fibre level plasticity rather than properties related to fibre bonds [1,10].
Strain field mapping in 3D of paper or paperboard materials are rare, but a recent application of synchrotron tomography and DVC analysis of notched handsheet paper samples was presented by Golkhosh et al. [15]. Due to the known auxetic behaviour of paper materials, e.g. [16,17], Golkhosh et al. [15] stressed the importance of considering the out-ofplane strain component that is not obtainable through 2D DIC analysis. One of the main findings of their study was that observed out-of-plane thickness increase appeared to be inversely proportional to the strength of inter-fibre bonds in the samples (controlled by three different handsheet production methods). Straightening and separation of fibres was suggested as the mechanism responsible for the auxetic behaviour. The evolution of strain field patterns during loading was, however, not discussed in this paper; only two components of the strain tensor (normal strain components in the loading direction and in the out-of-plane direction) were considered at the last load increment before failure. Also Krasnoshlyk et al. [18] presented results from 3D synchrotron imaging of low-density paper samples subjected to tensile loading. DVC was applied and a few examples of displacement results were reported. However, fracture mechanisms were in focus in this study and properties such as the Fracture Process Zone (FPZ) and appearance of the crack path were analysed, mainly through local variations in fibre floc distributions and local spatial changes in sample thickness and porosity during loading. Fibre rotations and deformation of inter-fibre bonds were observed and suggested as the cause of the global plastic behaviour of the material [18].
It is clear that a need for comprehensive analyses of the 3D strain field evolution in different paper and paperboard materials exist, focusing on, in particular, out-of-plane deformation and failure mechanisms. In this study, the 4D (spatial 3D and time) development of meso-scale strain fields before failure plus deformation and failure mechanisms at the scale of fibres and fibre networks are analysed for different types of commercial paperboards tested in different material directions. Xray tomography and DVC analysis have been applied to study paperboard samples in different loading steps during in-situ tensile tests. In addition, image analysis was used to characterize the evolution of sample properties such as spatial variations in sample thickness and solid fraction, porosity and fibre orientation distributions. The results contribute to an improved understanding of the deformation and failure mechanisms in paperboards that can be used to improve the performance of paperboard continuum models and, in turn, can lead to reduced package defects as well as material and food waste.

Sample Preparation and in-Plane Tensile Tests
In this work, in-plane tensile tests were performed on samples from two different types of commercial paperboard; a singleply paperboard (WestRock, WR) and a 3-ply paperboard (Billerud-Korsnäs, BK). The WR material is a bleached single-ply paperboard with clay and latex coating, while the structure of the BK paperboard consists of two unbleached bottom and middle plies covered by a third top ply of bleached fibres with clay coating. Samples were prepared in both MD and CD directions from both types of paperboards. The samples were notched to encourage the failure zone to develop in a well-defined area that was the in focus of the X-ray imaging. Notches were cut symmetrically on both sides of the samples resulting in an effective width of 3 mm between the notches, see sketch in Fig. 1(a). The samples were mounted in a custom-built tensile test device using specially designed grips as described in Tryding et al. [3], see Fig. 1(c). By clamping the sample between a cylinder (3 mm diameter) and a flat surface, the grips ensure that the sample is held a long a well-defined contact line. Each sample was mounted in the tensile tester with a 14 mm macroscopic gauge length, defined by the distance between the contact lines of the grips. The tensile test device used is an in-house built apparatus designed for in-situ X-ray tomograph measurements. Figure 1(b) shows the boundary conditions for the tensile loading. Photographs and a close-up of a mounted sample is shown in Fig. 1(c).
The loading device was mounted on the rotationtranslation stage of the tomograph before the sample was fixed between the clamps. The temperature was fixed to 28°C inside the X-ray tomograph while the ambient room temperature was 20-21°C during the experiments. The relative humidity in the room was around 31 ± 4% with minor variations during a day. After being mounted in the X-ray tomograph, the tensile test device and the sample was left to settle for several hours before the experiments were started. The X-ray tomography field-of-view is shown in the sketch over the boundary conditions in Fig. 1(b).
The tensile tests were carried out with a constant displacement rate of 0.28 mm/min, measured on the moving grip of the device. The force and displacement values were set to zero at the start of the test and the force was increased in increments of 15 N for the MD-tested samples and 8 N for the CD-tested samples. After each load step, the tensile tests were interrupted, and X-ray tomography data were collected. This resulted in pauses with durations of around 1 h and 15 min. During these pauses, force relaxation occurred fast (after a few minutes), before the tomography was started. The defined force increments resulted in four to five load steps along the increasing slope of the force-displacement curves (load steps 00-04, where load step 00 represent the initial state of the sample). The following load step (04 or 05) was located at the peak in the force-displacement curve and often occurred after a smaller force increase compared to the previous constant load increments. Finally, X-ray tomography data were also collected for a post-peak load step (05 or 06) at the descending part of the force-displacement curves. The total time for the experimental procedure was around 8 h for each sample.
The experimental procedure was repeated at least twice for each sample type to study the repeatability of the results. Figure 1(d) shows a sketch of the 3D volumes obtained from the X-ray tomography. The blue and red planes, respectively, marks the in-plane and out-of-plane central images in the 3D volume that are used to visualize 3D-volumetric strain field data below.

X-Ray Tomography
The X-ray tomography was performed with a Zeiss Xradia Versa XRM520 at the 4D Imaging Lab, Lund University. After each load step, 1601 radiographs were collected over 360°rotation with an X-ray tube voltage of 140 kV, a power of 10 W and an exposure time of 1 s. The field of view was focused around the in-plane centre of the samples ( Fig. 1(b)) and the source-sample-detector distances together with the 0.4X optical magnification resulted in a voxel size of 4 μm in the reconstructed tomography images. The total scanning time of each data set was around 1 h 10 min and automatic tomographic reconstruction results from the Zeiss Xradia Reconstructor software were used.
The reconstructed image stacks were cropped to a region of interest (ROI). The chosen ROI was minimized in size whilst being large enough to contain the full sample throughout the loading sequence (e.g., accounting for slight initial realignments and/or the geometry changes that occurred throughout the test). Therefore, the exact dimensions of the ROI varied slightly between the different samples. The dimensions of the image stacks after post-processing were on average 4.04 mm × 3.85 ± 0.04 mm (1010 × 961 ± 10 pixels) in the in-plane direction of the sample and 0.60-0.64 mm (151-161 pixels) in the out-of-plane direction.

DVC and Strain Fields
DVC was applied to the image data to analyse of the full 3D strain tensor evolution in each sample during tensile loading. The post-processed X-ray tomography volumes were analysed in an incremental manner, so that the strain evolution across each consecutive load step in the experimental series was obtained. For this analysis, the open source Python code Spam 0.4.3 [19] was used. Spam follows a regular-grid, local DVC approach and uses the local gradient of the displacement field to calculate the full 3D strain tensor. Isotropic and deviatoric strain invariants, representing the volume and shear distortion of the sample, are calculated within the framework of either large or small strain theory depending on the user input. The results from these frameworks were nearly identical when applied to our data and the small strain approximation was used here. Hence, the volumetric and deviatoric strain invariants can be defined as: The volumetric strain invariant (ε vol ) quantifies the average (isotropic) volumetric change between neighbouring DVC analysis points, while the deviatoric strain invariant (ε dev ) represents the total magnitude of the strain that is not related to the mean volumetric change (including shear deformations). A node spacing of 20 voxels (~80 μm) and a correlation window size of 21 voxels (~84 μm) was used for the DVC analyses. For all data sets presented in this study, most voxels converged to a solution with the stopping criterion of less than 0.001 change of the deformation function. This resulted in low subpixel correlation errors and smooth displacement fields.
After ensuring the quality of the obtained displacement fields, the components of the strain tensor were obtained together with the volumetric and deviatoric strain invariants.
The accuracy and precision of the displacement values obtained with Spam was evaluated experimentally. Two X-ray tomography datasets of a paperboard sample were collected, differing only by a known horizontal rigid body movement of 100 μm (controlled through the translation stage of the X-ray tomograph). The tomography settings were the same for the two scans and similar (but not identical to) the settings used for the loading experiments; the optical magnification and the voxel size of the resulting image stacks were the same (4 μm), but the X-ray tube voltage was 80 kV and the exposure time was 3 s (in the main experiments these were adjusted to 140 kV and 1 s, respectively, to reduce the total data acquisition time for the loading experiments). The two data sets were correlated in Spam using the same DVC settings as reported above. The resulting X-, Y-and Z-displacement values in all voxels of the DVC volume were inspected. The mean Xdisplacement value was −99.6 ± 0.37 μm, which is very close to the target, rigid-body movement of 100 μm. The mean Yand Z-displacement values from the DVC results were − 0.14 ± 0.03 μm and − 0.96 ± 0.13 μm, respectively.

Quantification of Sample Structures
The initial structures of all samples at load step 00 were quantified from the tomography data in Matlab 2019a®. In a first step, a Gaussian filter was applied to the image stack to improve the signal-to-noise ratio. Thereafter, a global image threshold value was found with Otsus's method and the intensity values below this threshold were removed from the image stacks. This approach removed most of the random noise that was primarily visible in voxels representing air around the paperboard samples, as well as within the samples between the fibres. However, local areas of noise with higher intensity than the global threshold remained. Therefore, a morphological filtering was applied involving connected pixels below a manually chosen threshold of 50 pixels.
The filtered image stacks were binarized (all remaining intensity values were set to 1) and used to calculate the initial solid fraction of the samples in three principal directions. For each direction, a 2D image was obtained where each pixel value represents a solid fraction value. The solid fraction values were calculated as the total number of pixels containing fibres along the axes perpendicular to the 2D image, divided by the total pixel length of the sample along the same axes.
The presence of the notches was accounted for in the solid fraction calculation so that the presented results are independent of the sample geometry. By calculating the solid fraction values in 2D projections, the high resolution of the original images could be retained (the resulting solid fraction images represent projections through the 3D volume in three principal directions as in Fig. 2(a)). Solid fraction estimations performed in 3D would require calculation over small volumes and reduced resolution of the results when projected along the three principal directions of the sample.
The inverse of solid fraction is the porosity, that was also quantified from the binarized image stacks in a similar manner. The main difference is that the porosity values are presented as average values below, as opposed to as 2D images. To avoid border effects that could lead to erroneous or misleading average porosity values, the binarized images were first cropped to dimensions of 910 × 601 pixels in the inplane direction and 31 pixels in the out-of-plane direction.
The out-of-plane thickness of the samples was obtained in a similar manner as the in-plane solid fraction images, resulting in 2D images where each pixel represents the total number of pixels between the top and bottom of the sample. This approach allowed a quantification of spatial variations in sample thickness. From the thickness maps, mean thickness values (at five vertical locations) were extracted over 200 horizontal pixels centred on the sample to also enable quantitative plots of the average thickness change with time.
The incremental changes of out-of-plane thickness between different load steps were also calculated. To reduce the effects of small-scale variations between the samples in consecutive load steps (e.g. individual fibre displacements as well as global changes in sample geometry and displacements within the field-of-view), the level of detail was reduced in the thickness maps before the difference between load steps were calculated. This was accomplished by morphological filling of the binary images, resulting in image stacks representing the outer geometry of the samples with slightly reduced level of detail from which the differences in thickness between consecutive load steps were calculated.

Quantification of Fibre Orientation Distributions
It is known that fibre orientations affect the mechanical properties of paper materials. The fibre orientation distribution is, for example, correlated to the stiffness distribution in paper materials [20] and the elastic stiffness is greater in the fibre orientation direction [21]. Fibre reorientation during loading might also influence affect the microscale deformation mechanisms. It is therefore of interest in this study to quantify and compare fibre orientation distributions in the samples during loading.
The fibre orientation distribution was quantified from the filtered image stacks with an open source Matlab tool from the QIM project. The QIM Orientation Analysis tool [22] utilizes computation of the structure tensor to produce two parameter values for each voxel in the 3D image: the azimuth angle and the elevation angle. For the paperboard samples in this study, the azimuth angle is the in-plane orientation while the elevation angle is out-of-plane orientation, see illustration in Fig. 2 All fibre orientation data were post-processed and visualized in Matlab 2019a®. As shown in Fig. 2(b), the azimuth orientations were represented by orientations ranging from 0°( oriented vertically/along the loading direction) to 90°(oriented horizontally/cross the loading direction). The elevation angles are most often close to 0°in the paperboard samples due to their planar structure. From these orientation data, the total in-plane and out-of-plane fibre orientation distributions were calculated by the use of histograms with 1°bin width. This results in volumetric proportions of voxels oriented in each 1°i nterval from 0 to 90°in the in-plane direction (and approximately 0-15°in the out-of-plane direction). Because of the anisotropic structure of the paperboard samples, the sample shape affects the global orientation distributions differently in MD and CD. Therefore, the fibre orientation distributions were obtained for subsets of the 3D volumes cropped to a rectangular shape centred in the in-plane direction (601 × 601 pixels). In addition to the calculation and comparison of initial sample properties, the change in fibre orientation distribution during loading was also analysed. Table 1 summarize some of the average initial properties of the samples used in this study. The average initial sample thickness at the in-plane centre was around 370 ± 15 μm, with no major differences between WR and BK samples. The average out-of-plane porosity (57 ± 3%) is higher than the average in-plane porosity (51 ± 3%) for all samples.

Initial Sample Properties
In Fig. 3, the initial fibre orientation distributions in the inplane and out-of-plane directions (azimuth and elevation angle, respectively) are shown. As expected, samples loaded in MD generally have a larger proportion of vertically aligned fibres (in-plane direction) while the opposite is true for samples loaded in CD (Fig. 3). An exception is BK-CD-1 that has a relatively even distribution of orientations. In some cases, there are significant variations between different samples even though they were extracted from the same paperboard. In WR-CD-2, for example, the distribution of fibre orientations is more diverse compared to WR-CD-1. The out-of-plane orientation distributions ( Fig. 3) have larger data uncertainties compared to the in-plane orientation data since these measurements might be affected by the non-perfect vertical alignment of the entire sample in the 3D image stack. Nevertheless, the initial distribution of out-of-plane orientations in Fig. 3 show that more fibres in the WR samples are inclined in the out-of-plane direction compared to the BK samples. During loading, the change in fibre orientation distributions were generally small and no general trends that could be related to either the loading direction or the out-of-plane paperboard structure were indicated by the data. The results are presented in Appendix 3.
The data in Table 1 and Fig. 3 do not indicate any major differences in average initial properties between WR and BK samples. However, spatial variations in initial physical properties exist and differ between both the individual samples and the paperboard types. Figure 4 shows the calculated spatial variations in the initial volumetric solid fraction (in three principal directions, see Fig. 2(a)) as well as initial sample thickness for the WR samples WR-MD-3 and WR-CD-2 and the BK samples BK-MD-4 and BK-CD-1. Small-scale in-plane variations in solid fraction can be observed for all samples, where areas of higher solid fraction (e.g. the~0.5-1 mm long structures marked by arrows A in Fig. 4) represent structures, possibly flocs, with higher fibre content than the more porous areas of low solid fraction. The in-plane solid fraction is generally higher in WR-CD-2 compared the other samples (as also reflected by the average porosity in Table 1). In the out-of-plane directions in Fig. 4, the clay/latex coated top layers of both types of paperboards are clearly visible as sheets with solid fraction values close to 100% (e.g. arrows B in Fig. 4). In both WR samples, the solid fraction is clearly lower in the out-of-plane centre of the sample (e.g. arrows C in Fig. 4), while the three-ply structure is visible in the BK samples as sheets with higher solid fraction in the out-ofplane bottom, centre and top (arrows D in Fig. 4). Figure 4 also shows significant spatial variations in initial sample thickness. In WR-MD-3, minimum thickness can be observed in the central parts of the sample (rectangle E) while the thickness generally is larger closer to the notches (the steep increase in sample thickness along the right notch is, however, an effect of the sample preparation). In WR-CD-2, maximum thickness is found in the right-and bottom parts of the sample while a zone of minimum thickness is observable close to the left notch (rectangle E). The BK samples are generally thinner compared to the WR samples. In BK-MD-4, a diagonal band of lower thickness is observable approximately from the upper right corner of the sample to the lower part of the left notch (rectangle E). In BK-CD-1, elevated sample thickness is found in the upper right part as well as in a band-like structure running from the bottom of the left notch towards the centre of the sample (between rectangles E).
For brevity, in the following we focus on the results from the four samples shown in Fig. 4, while only selected data for the replica samples in Table 1 will be shown and commented.
Global Mechanical Response Figure 5 shows the force-displacement curves measured on all the samples in this study together with markers identifying the load steps along the curves. Also shown on the right axes in Fig. 5 are the stress-displacement curves, where the width and the initial mean sample thicknesses between the notches were used for the stress calculations.
When tested in MD, all samples exhibit higher tensile strength and lower strain at failure compared to the samples tested in CD. The WR MD samples follows a linear increase in stress up to around (at least) 45-50 MPa and 35-40 MPa, respectively, and fail at 63 MPa and 59 MPa, respectively. The tensile strengths of the BK MD samples differ more significantly from each other, with values of 51 MPa and 73 MPa, respectively. The non-linear part of the stress-

Strain Field Evolution
The Spam code returns two measures that can be used to assess the quality of the DVC result in each voxel; a return status value indicating whether the voxel converged or not as well as a subpixel correlation error. The latter value represents the absolute difference in grey levels between the two correlated images. Table 2 summarizes the average subpixel correlation error of all voxels in the volumes representing each sample and load increment (the grey level errors have been normalized with the grey level range of the 16-bit input images). The results show that the subpixel correlation errors are small, generally below 1% of the total image intensity range, and confirms successful correlation is achieved. Volumetric and deviatoric strain fields from the DVC analysis of WR-MD-4 are visualized in Fig. 6 (left and middle column), where each vertical subfigure represents a load increment between load steps 00-04. The DVC results are represented by slices through the centre of the 3D volumes in the three principal directions, see Fig. 1(d). Also shown in the right column in Fig. 6 are maps of the change in sample thickness between the corresponding load steps for each pixel in the plane. These data are biased by in-plane changes in the sample geometry between the concurrent load steps that is especially evident around the edges of the sample. However, general trends in the spatial evolution of sample thickness are captured by the images and can be compared to the strain field evolution, indicating that the observed strain concentrations are, to some extent, spatially correlated to out-of-plane thickening.
In Fig. 6 (left column), patterns start to develop in the volumetric strain field at early load increments. The first clear signs of strain localization in the in-plane direction develop close to both notches in load step 01-02 (arrows A and B). In the following load increment (02-03), it becomes evident that the strain localization A propagates in up-and downward oriented diagonal patterns from the left notch towards the centre (dashed lines). This strain localization intensifies significantly in magnitude in load increment 03-04, while the smaller localization B close to the right notch increases to a lesser degree in load increment 03-04 compared to the previous load increment. In the out-of-plane direction, the strain initially localizes in the central part of the paperboard (arrow C in load increments 01-02 and 02-03). In load increment 03-04, the strain localization at the left notch appears to affect the full sample thickness in the out-of-plane direction (arrow C). It is likely that irreversible deformation, associated with the initiation of the non-linear shape of the force-displacement curve (Fig. 5) is initiated before or during this load increment.
The volumetric strain field evolution is associated with an increase in sample thickness in the same in-plane regions of the sample, most clearly visible in load step 03-04 (rectangle D in Fig. 6, right column). A slight increase in sample thickness around the left notch can be observed already in load increment 02-03, however not as widespread across the sample as the volumetric strain field at the same instance. It appears as if the in-plane volumetric strain localizations precede the local out-of-plane sample thickening with more pronounced thickening and coalescence with the volumetric strain field in load increment 03-04.
The individual strain components of the 3D strain tensor (data from load increment 03-04 are shown in Appendix 1, Fig. 13) confirm that the volumetric strain field is dominated by the zz component, i.e. a thickness increase in the out-ofplane direction. Another strain component significant for the volumetric strain is the yy component that represents the elongation of the sample in the loading direction. The yy component is clearly localized in a limited area close to the left notch (around the strain localization marked by arrow A in Figs. 6 and 13) and, to a lesser degree, along the previously described diagonal strain patterns. The diagonal patterns are also associated with a slight negative xx-strain, i.e. a compression in the horizontal (in-plane) direction.
The patterns in the deviatoric strain field in Fig. 6 are similar to the volumetric strain field patterns but slightly lower in magnitude. The deviatoric strain field is clearly affected by the zz and yy components of the strain tensor (Fig. 13), but shear components are also of importance. The individual strain components show that xy-shearing occurs, where the upper left part of the sample is displaced towards the (in-plane horizontal) right and the lower left part towards the left (see   Figure 7 show three different planes, extracted at different depths across the sample thickness, of the 3D volumetric strain field at load step 02-03, in WR-MD-3 (top row) and the replica sample WR-MD-2 (bottom row). Each plane is approximately 90-100 μm thick and slice 1 represents the plane furthest away from the clay coated top layer of the paperboard. In the out-of-plane direction, the magnitude of the strain field is, in both samples, generally lower in slice 1. The data also show similar features of the in-plane strain distribution in both samples. Although the strain clearly localizes at both notches in WR-MD-2, the magnitude of the strain is, in both samples, clearly higher close to the notch/notches compared to the sample centre. Diagonal strain patterns from the notches towards the sample centre are distinguishable in both samples (marked with dashed lines in the central slices in Fig. 7).
The in-plane appearance of the 3D strain fields is clearly different when the same type of paperboard is tested in CD, see Fig. 8. Instead of following diagonal paths with clearly higher strains close to the notches, the strain localizes relatively evenly across the full sample width in load step 03-04. In most other aspects, the strain field evolution behaves in a similar manner as in MD; for example, strain field patterns Fig. 6 Volumetric (left column) and deviatoric (middle column) strain fields for different load increments in WR-MD-3 together with spatial variations in out-ofplane sample thickness in different load steps (right column). All strain fields are visualized with a colorscale ranging from −1 to 6% start to develop in the out-of-plane central part of the paperboard (arrows C) at early load increments and is associated with localized increase in sample thickness (dashed rectangles D) at later load increments. Furthermore, the volumetric strain is dominated by the yy and zz strain tensor components, and significant shearing is observable in the zy and zx components (Appendix 1, Fig. 14). In contrast to WR-MD-3, however, the yy strain component (Fig. 14) is not localized to limited zones close to the notches but clearly extends across the full sample between the notches. The xy strain component is slightly negative in a diagonal band running from the (in-plane) lower left to the upper right part of the sample, while it is slightly positive in the upper left and lower right corners, thereby creating an x-shaped shear band pattern. Figure 9 shows three different planes of the 3D volumetric strain field at load step 02-03, in WR-CD-2 (top row) and the replica sample WR-CD-1 (bottom row). In WR-CD-1, the volumetric strain localization is more inclined compared to WR-CD-2 since it propagated across the sample from the lower part of the left notch to the upper part of the right notch. Other than that, there are similarities between the replica samples. In contrast to the samples tested in MD (Fig. 7), the strain localizations in the CD-tested samples clearly extend across the full sample widths between the notches, while the uppermost and lower-most parts of the samples are less affected. The strain patterns propagate in relatively straight bands between the notches and no significant diagonal features are observable. Except the slightly lower average magnitude of the strain furthest away from the coated top layer (slice 1) in both WR-CD-2 and WR-CD-1, the volumetric strain fields appear relatively constant in the out-of-plane direction, which is a feature in common with the samples tested in MD (compare with Fig. 7).
A corresponding presentation of the strain field evolution results from the BK samples tested in MD and CD is given in Appendix A2 . In general, the results are similar to those described for the WR paperboard above; the strain field patterns are more diagonal in MD and straighter and more continuous between the notches in CD. General differences between the WR and BK results are that the latter strain field patterns are, in comparison, more patchy and evenly distributed across the BK samples and develop more gradually with increased loading. Furthermore, larger out-of-plane variations in the strain field patterns can be observed in the results from the BK samples, while the individual shear strain components, in particular the zy component, are less significant compared to WR samples.

Delamination
The results of the 3D strain field analyses consistently showed that the deformation of all samples is strongly dominated by expansion in the out-of-plane (thickness) direction. The data shows that the strain localizes in the out-of-plane central part of the paperboard at early load increments. This can be explained by the lower solid fraction in the central plane of the paperboard, indicating lower material strength. In the in-plane direction, there is, for all samples, a spatial agreement between the strain field localization and a local increase in sample thickness. The in-plane strain localizations do, however, appear to precede the thickening of the same areas.
The spatial thickness change in the MD-tested samples at load increment 03-04 is clearly localized to regions around initial strain localization points at the notches for both paperboard types (Figs. 6 and 15). Another behaviour is indicated by the data from the CD-tested samples (Figs. 8 and 18). Here, Fig. 7 Out-of-plane variation in the volumetric strain field in WD-MD-3 and the replica sample WR-MD-2 at load step 02-03 that represent an increment at the late stage of the linear part of the global force-displacement curves. Each plane is approximately 90-100 μm thick and slice 1 is located furthest away from the clay coated top layer of the paperboard. All strain fields are visualized with a colorscale ranging from −1 to 4% the increase in sample thickness does not appear to develop around well-defined initial strain localization points to the notches, but rather across the full sample width (in agreement with the strain field patterns).
Minor temporal differences in the delamination process between different types of paperboards are indicated by the data. In the WR paperboard samples (Figs. 6 and 8), notable changes in thickness start to appear at later load increments (in the non-linear regime), while the process appears to be more gradual in the BK samples where signs of sample thickening can be indicated already in the linear phase (Figs. 15 and 18). This pattern is also distinguishable in Fig. 10 that shows the average thickness evolution in five horizontal planes across the thickness maps as a function of load step (and, on the right axis, the applied force). Figure 10 also shows that the increase in sample thickness of the WR MD samples is small and gradual up to load step 04, whereas the WR samples tested in CD shows a slightly larger increase in sample thickness already between load steps 03 and 04. In all samples except BK-CD-1, the slight sample thickening during the linear and/ or plastic regime is followed by a steep increase in sample thickness when the samples fail.
Failure and Fracture Figure 11 shows, for both WR-MD-3 and WR-CD-2, the spatial distribution of solid fraction at load step 05, i.e. at the peak stress, as well as a visualisation of the fracture after failure (visualized as the average intensity through the image stack from load step 06). Figure 12 shows a corresponding image for the BK samples BK-MD-4 and BK-CD-1. Clear differences between the samples tested in MD and CD, respectively, can be  Figs. 11 and 12). In the CD-tested samples, the decrease in solid fraction and increase in out-of-plane thickness appears more evenly distributed at peak load compared to the MD-tested samples. In WR-CD-2, slightly lower values of solid fraction can be observed at the locations where the fracture develop as well as in the out-of-plane centre (arrows B in Fig. 11). Such subtle spatial differences are more difficult to identify in BK-CD-1 (Fig. 12).
Also the characteristics of the fractures differ between MD and CD. The fracture in WR-MD-3 (Fig. 11) develops in a Fig. 9 Out-of-plane variation in the volumetric strain field in WR-CD-2 and the replica sample WR-CD-1 at load step 02-03 that represent an increment at the late stage of the linear part of the global force-displacement curves. Each plane is approximately 90-100 μm thick and slice 1 is located furthest away from the clay coated top layer of the paperboard. All strain fields are visualized with a colorscale ranging from −1 to 4% Fig. 10 Average thickness at different sample location plotted as a function of load step for all samples in this study. The global applied force at each load step is also shown for comparison. The thickness labels refer to the pixel number in the z-direction, i.e. the in-plane vertical location. Since all samples have a total vertical height of 1010 pixels, the data represent five planes through the approximate vertical centre half of the sample zigzag pattern relatively straight across the sample, from the initial strain localization points A to B (Fig. 6). Also in BK-MD-4 ( Fig. 12), the fracture path has a well-defined and angular appearance. At the left notch, the fracture coincides with the initial strain localization point A (Fig. 15). A comparison with the initial properties in Fig. 4 shows that the first part of the fracture in BK-MD-4 is located in a zone with minimum sample thickness (Fig. 4, rectangle E). The fracture then sharply turns downwards at an angle of around 45°and the remaining part of the fracture runs towards the right notch in an area where the initial thickness was elevated and the initial solid fraction was slightly lower than average (Fig. 4).
In both CD-tested samples, the fractures appear straighter and less angular. In WR-CD-2, an initial in-plane strain localization is clearly visible in the lower part of the right notch in load step 01-02 (arrow A in Fig. 8) coinciding with the location of the fracture in Fig. 11. In this zone of the sample, the initial thickness was relatively large while the initial solid fraction was relatively low (Fig. 4). At the left notch, the strongest strain localizations appeared close to the bottom of the notch (arrow B in Fig. 8). However, based on the appearance of the mean intensity image the fracture in Fig. 11, the fracture does not seem to reach, nor propagate towards, strain localization B. Instead, Fig. 11 shows that another fracture developed in the upper part of the left notch, within a zone of the sample where the initial thickness of the sample was at a minimum ( Fig. 4, rectangle E). It therefore appears as if the fractures initiated from weak zones in the material structure along the notches, i.e. spatial locations that are either thinner or less dense than the average sample, and not necessarily from the points where the strain initially localized most strongly. In BK-CD-1, the exact location of the fracture in Fig. 12 is slightly uncertain due to some sample movement during the X-ray tomography at load step 06, leading to lower quality of the mean intensity image for this sample. It is, however, clear that the fracture seemed to have developed along or above the uppermost of the early strain localization bands (Fig. 18).
In summary, the results indicate that the failure and fracture process generally differs between samples tested in MD and CD, where the fractures in the former are sharply defined already at peak stress and the main damage (in terms of decreased solid fraction/increased porosity) seem to be clearly localized in a limited region around the fracture, while the upper and lower parts of the sample appears almost unaffected. In contrast, the decrease in solid fraction/increase in porosity is more evenly distributed in the CD-tested samples at the peak stress and it is difficult to define the fracture at this point. These data indicate that the CD-tested samples might fail internally (at peak load), before any fracture appear at the surfaces of the sample.

Discussion
In our DVC data, patterns in the volumetric and deviatoric strain fields clearly appear in an early phase of the deformation, i.e. during the linear part of the global mechanical response. This is in contrast to previous observations in 2D where the first streak patterns in strain fields appear after the beginning of the plastic regime [11]. According to Hagman & Nygårds [11,12], random strain variations were present in their data during the initial linear stress-strain phase and clear strain patterns did not emerge until the plastic regime. In the data from Korteoja et al. [9], it is unknown when macroscopic strain patterns started to develop since the optical method used relies on plastic deformation and micro failures to identify strain localization. Differences between the appearances of the strain patterns from our DVC analysis and previous DIC results are, however, not unexpected since the full 3D strain tensor and internal variations, not just surface strain fields, is considered in the present study. Nevertheless, early strain localization occurs in all strain tensor components in our study.
A consistent result among all of the samples studied here is that the out-of-plane (zz) strain component is clearly stronger in magnitude than the in-plane (xx and yy) strain components (as well as all the shear strain components). The out-of-plane component therefore dominates the patterns in the volumetric strain fields in our study, while this component is not detectable with 2D DIC methods.
In low-density papers, previous research has shown that the fracture propagation path appears to be controlled by variations in the paper formation; the crack appears to develop in non-straight patterns thorough low-density zones in between flocs [18]. In our data, the roughness of the crack appears to be also related to the material anisotropy; straight crack paths were observed for samples tested in CD while rough and irregular cracks formed in samples tested in MD. The strain fields and thickness maps both indicate more localized nucleation of fracture process zones in the MD-tested samples.

Strain Field Patterns -An Indication of Deformation Mechanisms
The results in this study have shown that there are general differences in the types of strain field patterns that develop in MD-and CD-tested samples, respectively. Diagonal strain streaks propagating from the notches towards the in-plane centre were present in all MD-tested samples. In the WR MD samples, these streaks appeared at angles of around 45°t o the loading direction and their localization was relatively independent on load increment (i.e., the main strain field patterns appeared already in load increment 01-02, where after the strain tended to intensify in magnitude in already established patterns). In the BK MD samples, the directions of the diagonal stain streaks ranged between approximately 25-45°, with small angle variations between different load increments. In contrast to previous DIC studies in 2D [9,11], diagonal strain streaks were not observed in the CDtested samples in our study. In the MD-tested samples, however, the appearance and angles of the observed diagonal strain streaks (especially in WR MD samples) are similar to what has previously been reported and described in 2D DIC studies [9][10][11].
In all the CD-tested samples in this study, the strain streaks were generally less inclined, often aligned perpendicular to the loading direction, and propagated more continuously across the samples from notch to notch. The volumetric and deviatoric strain fields were not lower in magnitude in the in-plane centre of the samples, as observed in all the MDtested samples. The individual strain tensor components show that significant elongation in the loading direction (yy strain) occurred across the full sample width in CD, in contrast to the MD-tested samples where this strain tensor component is of a significantly lower magnitude and spatially limited to distinct strain localization points on the notch edges. The angles of the strain streaks in WR-CD-2 ranged between 70 and 90°to the loading direction, while they propagated at an angle of 65°to the loading direction in WR-CD-1. In the latter sample, however, the inclined strain streaks are of a clearly different character compared to the diagonal strain streaks in the WR MD samples, since they in WR-CD-1 propagate across the sample from one notch to the other. In both BK CD samples, the strain streaks generally appeared across the samples at angles between 80 and 90°to the loading direction, although some strain streaks also propagated at slightly more inclined angles.
There are several possible mechanical explanations to the appearance of diagonal strain patterns in the MD-tested samples, one being the structural shape of the notched samples themselves. The rounded shape of the notches are, however, designed to reduce the emergence of purely sample shapedependent strain patterns. Furthermore, the fact that the diagonal strain streaks appeared only in the MD-tested samples and not in the CD-tested samples suggests that they are at least partly caused by material anisotropy. In general terms, the angles of the MD strain streaks (close to 45°) could indicate both ductile deformation behaviour and brittle shear failure. Tryding et al. [3] showed that the stable sample length parameter, that depends on the elastic modulus and the post-peak slope of the global stress-displacement response, can be used as an indicator of if a paperboard sample behaves in a more ductile-or a more brittle-like manner. The stable length differs between samples tested in MD and CD and the diagonal strain field patterns in the MD-tested samples in this study might, together with the global force-displacement response, therefore indicate a more brittle behaviour compared to the CDtested samples. The effects of more or less brittle/ductile behaviour on strain field patterns are, however, outside the scope of this study and require more systematic experiments with varying sample lengths. Below, we interpret the differing strain field patterns in MD and CD paperboard samples in terms of micro-mechanical considerations.
Plastic elongation of paper materials is considered to be caused by plastic elongation of the fibres in the paper network, while macroscopic failure is caused by breakage of inter-fibre bonds [1]. In an ideal paperboard sample tested in CD, the majority of the fibres are aligned perpendicular to the loading direction. When subjected to a loading force, these fibres separate from each other (as the fibres aligned in the loading direction elongate plastically), leading to stress/strain localization along the fibre bonds oriented cross the loading direction. This can explain the appearance of the strain fields that develop in the WR samples tested in CD, i.e. a sharply defined band of elevated strain across the sample width. When the paperboard is tested in MD, plastic elongation of the fibres oriented along the loading direction leads to increased vertical distance between inter-fibre bonds and stress/strain localization at these points. The strain field patterns developing in MD-tested samples are likely to become more variable compared to CD, since the vertical displacement of the inter-fibre bonds depends on the plastic elongation of individual fibres. When tested in CD, the occurrence of vertical/horizontal inter-fibre bonds in the loading direction is much denser compared to when tested in MD. Therefore, it is more likely in this case that the fewer fibres aligned in the loading direction behave "collectively" when the material is tested in CD, due to their strong bonding to fibres aligned cross the loading direction.
In addition to the differences between MD-and CD-tested samples, there are also clear differences in strain field pattern between the two paperboard types WR and BK. The significant out-of-plane variations in strain patterns observed for BK-CD-1 are likely caused by the 3-ply structure of the BK paperboard. Although some out-of-plane variation exists also in WR samples, the transition between the slices in the 3D volume appears to be less smooth in the BK samples. Furthermore, the individual strain tensor components indicated that out-of-plane shearing between planes appears stronger in magnitude and therefore more important in WR samples compared to BK samples. The in-plane strain field patterns in the BK samples were, generally, more disperse and distributed across the samples. This indicates that there is a larger variability in material properties within the BK paperboards that causes the strain to distribute across several locations.

Volumetric Expansion, Delamination and Fracture Mechanisms
The in-plane evolution of sample thickness in the MD-tested samples suggests that the delamination of the sample initiates in a FPZ that is well-defined in the in-plane direction, likely through plastic elongation of fibres followed by breakage of inter-fibre bonds and inter-fibre separation. Thereafter, the same delamination process appears to spread to neighbouring areas, limited by the in-plane strain field distribution. MDtested samples appear to fail more abruptly due to strongly localized strain, deformation and delamination in strain localization points close to the notches. At peak load, the solid fraction images show that the fracture appears to be severe and extend beyond the stiffer top-and bottom sheets (including the clay coated top layer), despite the fact that the delamination process occurs mainly in the out-of-plane central part of the paperboard.
In CD-tested samples, the increase in thickness before failure occur across the full sample width suggesting that many inter-fibre bonds break simultaneously across the sample and thus immediately initiates the delamination process in a larger region or FPZ. The results support the microscale interpretation of the CD strain fields above, i.e. that the continuous strain fields from notch to notch arise as a result of strained inter-fibre bonds along fibres aligned cross the loading direction. The CD-tested samples appear to fail more gradually through a wide-spread delamination process, active across the entire sample in between the notches. The solid fraction images show that at peak load, the fracture is barely visible at the surface of the sample and the top-and bottom sheets of the paperboard appears to be intact. This suggests that these samples break internally when enough inter-fibre bonds have failed, and this continuous process could also explain the more gradual global post-peak behaviour of the CD-tested samples. Furthermore, the relatively straight fractures appearing visually after peak load also suggest a gradual separation of fibres aligned cross the loading direction, rather than the abrupt rupture and crack propagation behaviour indicated by the data from the MD-tested samples.
The main difference in the delamination behaviour of the samples in this study is clearly related to the main fibre orientation, i.e. if the sample is tested in MD or CD. The different out-of-plane structures of the paperboards (WR and BK) appears to only have a minor influence on the behaviour. There are, however, slight indications of a more continuously developing delamination process in the BK samples. This could be related to the 3-ply structure that likely infer a larger variation in inter-fibre bond strengths within the samples.
The main temporal trend in the delamination process of all samples is a dramatic increase in sample thickness associated with peak stress and failure (Fig. 10). The acoustic emission data representing breakage of fiber bonds in Isaksson et al. [14] indicate a similar trend during uniaxial tensile loading, with an abrupt increase in the number of acoustic events occurring close to peak load and failure [14]. It is therefore likely that breakage of inter-fibre bonds is a direct cause of the delamination and increase in sample thickness that has been shown here to be localized in zones of high out-of-plane normal strain.
Previous research has indicated that fibre-scale movements, such as torsion [18] or straightening of fibers [15] leads to inter-fibre bond deformation and breakage, fibre detachments and a global increase in paper thickness. The mesoscale approach used here to quantify change in fibre orientation distribution via the structure tensor resulted in small variations in fiber orientation distributions with tensile loading. Although our results provide some indications of microscale mechanisms, a successful segmentation of the entire fibre network together with tracking of fibre-and fibre-bond movements and strain evolution is probably required in order to fully understand the micro-scale mechanisms responsible for the material deformation.

Conclusions
This work has shown that in paperboards subjected to tensile testing, the 3D strain fields are clearly dominated by the outof-plane normal strain component. It is previously known that delamination of paperboard is a key mechanism leading up to failure of the sample during tensile loading, and the experimental 4D approach used here has enabled this auxetic material behaviour to be quantified both spatially and temporally. In-plane strain localization patterns (in the out-of-plane central parts of the samples) precede a local increase of sample thickness, i.e. the initiation of the delamination process. At peak load, a dramatic increase in average sample thickening occurred.
The in-plane distribution of the strain fields differed significantly between samples tested in MD and CD, respectively. The patterns, that emerge already in the linear part of the global force-displacement curve, have a more diagonal appearance in MD and generally propagate from the first strain localisation points at the notches and towards the inplane centre of the samples. In CD, the in-plane strain field patterns are more continuously distributed across the samples between the notches. These differences in strain field patterns might be controlled by the distribution and displacement/ failure of inter-fibre bonds in the fibre network. However, more experiments with varying sample lengths are needed to understand and account for the effects of sample length and brittle/ductile material response on the strain field patterns emerging during tensile loading.
After failure, appearances of the fractures were dependent on the material direction; the fractures were angular in MDtested samples and relatively straight in CD-tested samples. The MD sample failed abruptly with fractures clearly visible on the sample surface at peak load, while the CD-tested samples appeared to fail more gradually. At peak load, the CDtested samples increased more evenly in thickness with no clearly visible fractures on the sample surfaces, suggesting that they might initially have failed internally through delamination.
Similar deformation and failure behaviour were observed for the single-ply WR and 3-ply BK paperboards. However, the strain localization patters were more dispersed across the BK samples, larger out-of-plane variations were observed, and the strain tensor shear components were less significant than in the WR samples. The delamination process also appeared more gradually during loading in the BK paperboard.
We conclude that X-ray tomography, image analysis and DVC analysis of strain field evolution in 3D are effective tools to gain further understanding of the micro-and meso-scale deformation and failure mechanisms in different types of paperboards and different material directions. In future work, a similar approach will be applied to in-situ tensile loading during continuous synchrotron tomography imaging.
Strain field evolution, BK samples Figure 15 shows, for the BK sample tested in MD, the evolution of the volumetric and deviatoric strain fields together with the change in thickness. Already in load increment 00-01, diagonal patterns are visible in the volumetric and deviatoric strain fields (marked with dashed lines in the volumetric strain field). These patterns originate from both notches and propagate in up-and downward oriented diagonal patterns towards the centre of the sample where they meet and create an approximately squarelike pattern. In load increment 01-02, the strain fields extend to neighbouring areas and the diagonal patterns become more diffuse with the maximum strain localizations close to both notches, marked by arrows A and B. At the same time, parts of the initial diagonal strain localizations from load step 00-01 start to dominate both the volumetric (marked by dashed lines) and deviatoric strain fields, and become more intense with each load increment up to load step 03-04. These localizations consist of diagonal patterns propagating from strain localizations A and B and downwards in a south-southeastern direction, at an approximate angle of 45-65°from horizontal. In the out-ofplane direction, the volumetric strain field localize clearly in the central part of the sample (marked by arrows C) and remain there up to load increment 03-04.
The change in thickness (Fig. 15, right column) indicates a gradual thickness increase, relatively evenly distributed over the sample. In load increment 03-04, maximum thickness increase occurs in the central/lower part of both notches (marked with dashed rectangles D). The individual components of the volumetric strain field (Fig. 16) is similar to the observations made for the WR sample tested in MD; the zz component (representing a thickness increase in the out-of-plane direction) is dominant while the yy component (elongation in the loading direction) localizes in small areas close to both notches (stain localizations A and B). Some weak indications of the squarelike pattern of the volumetric and deviatoric strain fields at load step 00-01 can be distinguished in the xx strain component.
The deviatoric strain fields in Fig. 15 (middle column) are similar to the volumetric strain fields, but smaller in magnitude. The shear components of the strain tensor for this sample (Fig. 16) indicate a weak x-shaped shear band pattern in the xy component. The zyand zx components, that indicate shearing between planes in both in-plane directions, are generally of a similar magnitude as the xy component (i.e. relatively weak compared to both WR samples in Figs. 13 and 14). The zy component strain is diffuse and the identification and directions of zy shearing planes are not clear. The zx strain component field is also relatively weak and diffusely distributed, except for clear in-plane localizations of negative strain at the left notch and positive strain at the right notch. In the out-of-plane direction, the negative strain appears to localize in a bottom sheet of the sample. These patterns thus indicate that an upper part of the sample (close to the coated top layer) is displaced from a bottom part in an in-plane rightward direction. Figure 17 shows three different slices through the 3D volumetric strain field at load step 02-03, for BK-MD-4 (top row) and the replica sample BK-MD-1 (bottom row). The out-of-plane variation of the strain field patterns appears to be more significant in the BK samples tested in MD compared to the WR samples, which is likely caused by the 3-ply structure of the BK paperboard. Figure 17 shows that the volumetric strain clearly localizes in the out-of-plane centre (slice 2) in BK-MD-4. On either side of this plane (slice 1 and slice 3), it is clear that the strain magnitude is generally lower in large parts of the sample. Furthermore, the locations of elevated strain, as well as their patterns and directions, differ between the different slices (marked with dashed lines). Similar observations can be made from the data in BK-MD-1, where the appearance of the volumetric strain field varies significantly in the out-of-plane direction.
The volumetric and deviatoric strain fields developing in BK-CD-1 during loading (Fig. 18) have differences and similarities to both the BK MD samples and the WR CD samples. In contrast to BK-MD-4, but like WR-CD-2, the strain localizations are mainly developing in horizontal directions (marked with dashed lines). However, the distribution of strain localisation is generally more diffuse and widespread across the sample in the in-plane direction and more variable in the out-of-plane direction than WR-CD-2. In this aspect, the results from BK-CD-1 are more similar to the BK samples tested in MD.
In Fig. 18, weak early strain localization is visible at load increment 00-01, where slightly elevated volumetric and deviatoric strain appears close to the left notch (arrow A) and continues across the central/lower part of the sample towards the right notch. In load increment 01-02, two clear horizonal strain bands appear across the sample, from the localization at the left notch (arrow A) towards two points on the right notch (arrows B). These patterns intensify with further loading (02--03 and 03-04) and areas of elevated strain spread to neighbouring areas. Load increment 03-04 corresponds to a stage on the non-linear part of the global force-displacement curve (Fig. 5).
The initial thickness of BK-CD-1 is smaller in the lowest part of the sample (see Fig. 4), and the first localization of strain bands in Fig. 18 appears to be related to this geometrical feature. The lowermost of the two main bands marked in load increment 01-02 appear close to the border towards a thicker zone in the sample (compare with Fig. 4). The visualisation of Fig. 15 Volumetric (left column) and deviatoric (middle column) strain fields for different load increments in BK-MD-4 together with spatial variations in out-ofplane sample thickness in different load steps (right column). All strain fields are visualized with a colorscale ranging from −1 to 6% the thickness change (Fig. 18, right column) shows a localized zone of thickness increase along this strain band (rectangle D). Furthermore, an area of significant compression (rectangle E) is present below the strain band. Moving though the slices of the strain field 3D volumes in the out-of-plane direction (data not shown), it becomes clear that the location of the compression zone (rectangle E) is associated with negative volumetric strain in the planes furthest away from the clay coated top layer. The same area is also characterized by elevated deviatoric strain that indicates shearing of the sample in the initial stages of the loading. This strain localization appears clearly already at load increment 00-01, remain in 01-02 but becomes less significant in load increments 02-03 and 03-04, in agreement with the thickness change in Fig. 18. In load increment 03-04, the volumetric strain is elevated in large patches spread across almost the entire sample as are the thickness increase, marked by rectangle D in load increment 03-04.
As with the other samples in this study, the individual components of the strain tensor (Fig. 19) confirm that the volumetric strain field is dominated by the zz component and to a lesser degree the yy component. The yy component is, like the zz component, elevated across the sample but with a significantly smaller magnitude.
The individual components of the deviatoric strain field indicate a weak x-shaped shear band feature in the xy component, as previously described for WR-CD-2 and BK-MD-4. Like the other samples, the zy and zx components indicate shearing between different planes in the sample, although the magnitude of all the off-diagonal components are small in this sample. In BK-CD-1, the zy component is slightly negative in an upper plane in the out-of-plane direction (marked with arrow F in Fig. 16) and slightly positive below; this indicates that a plane close to the coated top layer is displaced downwards (in the loading directon) in relation to a plane further away from the top layer. The strain localizations at the right notch (marked by arrows B) are elevated in the zy strain field. The zx component reminds of the behaviour in BK-MD-4, indicating shearing between planes in a horizontal direction with positive strain at the right notch and negative at the left notch. Figure 20 shows three different slices of the 3D volumetric strain field in BK-CD-1 at load increment 02-03 that represent a late stage of the linear part of the force-displacement curve (Fig. 5). Since the strength of the replica sample BK-CD-2 is significantly lower (Fig. 5), data from this sample at load increment 01-02 are shown for comparison with BK-CD-1 in Fig. 20. The noisy appearance in the upper right corner of both datasets in slice 1 is a result of a slight inclination of the samples in the 3D volume, causing the clay coated top layer of the paperboard to appear in this region (where the image correlation failed due to the lack of texture in the clay coated layer). Figure 20 shows that the general patterns of the volumetric strain fields (marked by dashed lines) are of a similar character in BK-CD-1 and BK-CD-2. The directions of the strain bands are mainly horizontal or inclined and extend across the sample. For both samples, the strain magnitude is generally lower closer to the coated top layers (slice 1). Although the variation between different slices in the out-ofplane direction is smaller than in the BK samples tested in MD (Fig. 18), they appear more variable in the out-of-plane direction compared to the WR samples (Figs. 7 and 9). Compared to the WR samples tested in CD (Fig. 9), the strain localization patterns in the BK CD samples at load increment 03-04 in BK-CD-1 are less clearly localized and smaller in magnitude.  Figure 21 shows changes in fibre orientation distribution in both WR MD samples during loading. Unexpectedly, the behaviour differs between these two replica samples. In WR-MD-3, a very slight gradual increase in the proportion of vertical fibres (close to 0°orientation) can be observed from load step 00 to 04. At peak stress in load step 05, the proportion of vertical fibres increase significantly at the expense of horizontal and aligned fibres (with orientations greater than 25°). The out-of-plane orientations (azimuth angle) do not change much with loading.
In WR-MD-2, there is a significant decrease in the proportion of vertical fibres between load step 00 and 01 that is counter-intuitive given the loading direction. With further loading (between load step 02 and 03), fibres instead start to realign back from horizontal to vertical, i.e. towards the loading direction, before another decrease in the proportion of vertical fibres occur at load step 04. This is followed by a large increase in the proportion of vertical fibres at load step 05, i.e. at peak load when the sample failed. In both WR-MD-3 and WR-MD-2, the proportion of vertical fibres appears to decrease slightly with further separation after failure (load step 06).
An initial behaviour similar to WR-MD-2 can be distinguished in WR-CD-2 (Fig. 22), where a slight decrease in the proportion of vertical fibres occurs between load step 00 and 01. This is followed by a gradual realignment of fibres from horizontal back to vertical between load steps 02 and 04 and a more significant increase in the proportion of vertical fibres in load step 05. With further displacement after failure (load step 06), the fibre orientation distribution becomes flatter, indicating that a relatively large part of originally aligned and horizontal fibres realign in the loading direction, possibly by a fibre pull out mechanism.
The behaviour of WR-CD-2 during the initial stage of loading is different from the replica sample WR-CD-1 (Fig. 22). In WR-CD-1, the proportion of vertical fibres increases constantly during loading, and the difference between each consecutive load step is slightly larger than in WR-CD-2. However, the fibre orientation distribution becomes flatter at load step 06 also in this sample, indicating similar behaviour after failure in both WR CD samples. Figure 23 shows the change in fibre orientation distribution in the BK MD samples during loading. In BK-MD-4, no significant difference with loading can be observed, except an abrupt out-of-plane reorientation of fibres at peak load (load step 05). Again, the replica sample behaves in a different way; in BK-MD-1, the fibre orientation distribution is constant in load steps 00, 01 and 02. After this, the proportion of vertical fibres decreases in load step 03, again counter-intuitive given the loading direction. This is followed by a large increase in proportion of vertical fibres in load step 04 (peak load) and then another decrease in load step 05 (after failure).
In BK-CD-1 (Fig. 24), a slight and gradual increase in the proportion of vertically aligned fibres occur during Fig. 17 Out-of-plane variation in the volumetric strain field in BK-MD-4 and the replica sample BK-MD-1 at load step 02-03 that represent an increment at the late stage of the linear part of the global force-displacement curves. All strain fields are visualized with a colorscale ranging from −1 to 4% loading. No abrupt changes occur when the sample failed at load step 05. In the replica sample BK-CD-2, there is a slight initial decrease in the proportion of vertical fibres between load steps 00 and 01, accompanied by an increase in the proportion of horizontal fibres, i.e. an initial behaviour similar to WR-CD-2 and WR-MD-2. Also like in WR-CD-2 and WR-MD-2, the proportion of vertical fibres then increase gradually with further loading in BK-CD-2.
In summary, an initial decrease in the proportion of vertical fibres during the first stage of loading (between load steps 00 and 01) was observed in four out of eight samples. These four samples were all of different types and there are no obvious similarities between the samples. Thus, several samples showed indications of a counter-intuitive initial redistribution of fibre orientations within the samples, effectively leading to less fibres being aligned in the loading direction after the first load step. These results suggest a fibre network mechanism taking place during the linear phase of the global mechanical behaviour. Straining and straightening of vertical fibres might also induce straightening of fibres aligned cross the loading direction, as the fibre-to-fibre bonds are displaced throughout the network. The drying process during paperboard production leads to anisotropic shrinking and internal stresses throughout the fibre network [23]. Therefore, there are likely more curved fibres present in a dry paperboard compared to initial state in the network during the production process, when the paper mass is still in a water saturated state. The initial change Fig. 18 Volumetric (left column) and deviatoric (middle column) strain fields for different load increments in BK-CD-1 together with spatial variations in out-ofplane sample thickness in different load steps (right column). All strain fields are visualized with a colorscale ranging from −1 to 6% in fibre orientation distribution observed during loading might suggest an initial straightening of the entire fibre network, possibly a mechanism reversing the shrinkage process during the material production. This mechanism could possibly explain why these data patterns were observed randomly in this study, with no apparent relation to neither initial fibre orientation distribution nor the initial structural properties of the samples.
With further loading beyond load step 01, a gradual increase in the proportion of fibres oriented in the loading direction was observed in most of samples in this study. This behaviour could be related to straightening or plastic straining of initially curved fibres. Furthermore, in most of the samples a significant increase in the proportion of fibres aligned in the loading direction occurred during and after the failure of the sample. This suggests that loose fibres are pulled out from the network.