Peat macropore networks – new insights into episodic and hotspot methane emission

. Peatlands are important natural sources of atmospheric methane ( CH 4 ) emissions. The emissions are strongly inﬂuenced by the diffusion of oxygen into the soil and of CH 4 from the soil to the atmosphere. This diffusion, in turn, is controlled by the structure of macropore networks. The characterization of peat pore structure and connectivity through complex network theory approaches can give conceptual insight into how the relationship between the microscale pore space properties and CH 4 emissions on a macroscopic scale is shaped. The evolution of the pore space that is connected to the atmosphere can also be conceptualized through a pore network modeling approach. Pore regions isolated from the atmosphere may


Introduction
Peatlands are globally important modulators of hydrological and biogeochemical cycles (Gorham, 1991;Limpens et al., 2008). 25 Peatlands store vast quantities of water, and they may affect flooding patterns in watersheds (Holden, 2005). A large proportion of dissolved organic carbon in freshwater ecosystems originates from watersheds with peatlands (Kang et al., 2018;Asmala et al., 2019). Globally, approximately 400 Pg carbon is stored as peat, which corresponds to half of the carbon currently present in the atmosphere and a quarter or more of the total soil carbon stock (Rydin and Jeglum, 2013). Peatlands are major sources of carbon dioxide (CO 2 ) and methane (CH 4 ), thereby contributing to global warming (Frolking et al., 2011;Abdalla et al., 2016; 30 Leifeld et al., 2019). Management practices such as drainage and restoration affect the water table (WT) dynamics (Menberu et al., 2016;Evans et al., 2021) and therefore the air-filled porosity and the availability of oxygen (O 2 ) in peat (Waddington et al., 2015;McCarter et al., 2020). This, in turn, influences CO 2 and CH 4 emissions, as the oxidation of organic matter to CO 2 is stimulated under oxic conditions, whereas its reduction to CH 4 requires anoxic conditions. Such anoxic conditions prevail below the WT but also in microniches -anaerobic pockets -in unsaturated soil above the 35 WT (Silins and Rothwell, 1999;Deppe et al., 2010). These pockets form when the consumption of O 2 exceeds the transport, mainly diffusion of O 2 . If readily degradable organic carbon is also available, the anaerobic pockets may become microscale hotspots of CH 4 production (Wachinger et al., 2000;Hagedorn and Bellamy, 2011). As the peat CH 4 concentration increases above the atmospheric concentration, CH 4 may diffuse through the aerobic peat layer into the atmosphere if it is not oxidized by methanotrophic bacteria before reaching the peat surface (Whalen, 2005). Because the diffusion coefficient of CH 4 in air 40 is 4 orders of magnitude higher than in water (Ball and Smith, 2001), the air-filled porosity of the unsaturated zone largely regulates the diffusional CH 4 transfer in peat and therefore determines how much CH 4 is oxidized before reaching the peat surface.
Water retention characteristic is a fundamental soil property that links soil structure to water and aeration dynamics, redox conditions, and many accompanying biogeochemical processes (Bachmann and van der Ploeg, 2002;Lepilin et al., 2019). soil functions (Reddy and DeLaune, 2008;McCarter et al., 2020). The total volume of macropores generally decreases with depth because of a higher degree of decomposition in the deeper peat layers (Päivänen, 1973). Pore size distribution or the water retention characteristic alone does not provide information about the arrangement, connections, or topology of the pores.
Macropores form a complex network, where the individual macropores can be open and connected, dead-ended, or isolated 55 (Rezanezhad et al., 2016). The topology and structure of the network therefore regulates water, solute, and gas transport and, ultimately, biogeochemical processes in peat (see McCarter et al., 2020). Dead-ended and disconnected pores can block the gas transfer between peat and the atmosphere and promote the formation of anaerobic pockets above the WT (Knorr et al., 2009;. Thus far, the role of anaerobic pockets in the CH 4 processes has been neglected in simulation models (e.g., Fan et al., 2014), mainly because the pore network is difficult to characterize experimentally. 60 Pore characteristics have been earlier described by pore size distribution derived from the water retention characteristic (Laine-Kaulio, 2011;Weber et al., 2017;Lepilin et al., 2019) or by tortuosity indices derived from air permeability measurements (Laurén, 1997), both methods assuming homogeneous and isotropic soil structure (Beckwith et al., 2003). However, the CH 4 production and transport cannot be fully understood without considering the three-dimensional (3D) pore network structures in peat. X-ray imagery and complex network theory provide a very promising yet unstudied approach for describing these 65 structures and their effect on soil gas transport properties and mechanisms. X-ray micro-computed tomography (µCT) allows an explicit description of pore structure with resolution extending to micrometer scale (Perret et al., 1999;Blunt et al., 2013;Rezanezhad et al., 2016). Total porosity and pore size distributions can be determined directly from 3D images (Larsbo et al., 2014), and when the extracted pore space is represented as a 3D network of pores and pore throats, more detailed information about pore connections and topology can be obtained (Gostick, 2017). Methods of complex network theory are widely used 70 for quantifying multi-scale connectivity and transport processes in real-world networks (Newman, 2003). Network concepts such as clustering, centrality, tortuosity, isolation, and path lengths can be used to characterize the macropore network from the anaerobic pocket formation viewpoint. This may be a key in explaining the observed hotspots and episodic spikes of CH 4 flux, which are particularly difficult to explain in the current CH 4 models (e.g., Xu et al., 2016).
The aims of this study were (1) to introduce the µCT and complex network theory methods to analyze the characteristics of are 4.6 • C and 627 mm, respectively (Pirinen et al., 2012). The soil type is histosol dominated by Carex peat. The site was originally a mesotrophic fen classified as an herb-rich tall sedge birch-pine fen (Laine and Vasander, 1996). The dominating tree species, with a mean height of 20 m, are Scots pine (Pinus sylvestris L.) and Downy birch (Betula pubenscens Ehrh.) with an understory composed of Norway spruce (Picea abies Karst.). The stand volume is 230 m 3 ha −1 with a density of 2200 stems ha −1 . The ground vegetation is composed of dwarf shrubs with a coverage of 4 % (Vaccinium myrtillus L., V. vitis-idaea L.) and herbs (coverage 10.6 %) such as Dryopteris carthusiana (Vill.) H.P. Fuchs and Trientalis europaea L. The 90 moss layer is patchy and dominated by Pleurozium schreberi (Brid.) Mitt., Dicranum majus Turner, and D. polysetum Sw. In addition, Sphagnum girgensohnii Russow, S. russowii Warnst., and S. angustifolium (C.E.O. Jensen ex Russow) C.E.O. Jensen are present in moist patches. A detailed site description is available in Koskinen et al. (2014) and Bhuiyan et al. (2017).
Undisturbed peat samples were extracted into acrylic cylinders (diameter 50 mm, height 50 mm) using a sharp knife and scissors. The samples were collected from seven randomly located pits and three different depths (0-5 cm, 20-25 cm and 40-95 45 cm, hereinafter referred to as top, middle, and bottom layer, respectively). First, a 50-60 cm deep pit, with an undisturbed vertical face, was dug, and then the profile depth was measured with a ruler. Vertically oriented peat samples were extracted along the pit face paying attention to maintaining the undisturbed peat structure. The peat samples were located into plastic bags and transported to the laboratory, where the water retention measurement was started immediately.

100
The water retention characteristics of the cylindrical peat samples were measured in the laboratory using a pressure plate apparatus (Hillel, 1998). The samples were saturated with water, weighed, and placed into the pressure plate apparatus. The volumetric water content (θ, m 3 m −3 ) was then allowed to stabilize in the pressure plate apparatus under external pressures of 1, 3, 6, and 10 kPa for 1 week at each pressure level. These pressure levels are equivalent to soil matric potentials (Ψ) of −1, −3, −6, and −10 kPa. The sample mass (M Ψ , kg) was determined after each pressure level. At the end of the experiment, the 105 height and diameter of the sample were measured to determine the shrinkage, and then the samples were sent to µCT imaging (see Sect. 2.3). After the imaging, the peat samples were dried at 105 • C for 72 h to obtain the dry mass (M s ) of the sample.
The bulk density (ρ b , kg m −3 ) of the sample was determined from the dry mass and the volume of the sample in the saturated state (V sat , m 3 ). The volumetric water content θ was calculated at each matric potential in relation to the saturated volume of the sample as where ρ w (kg m −3 ) is water density. Total porosity f (m 3 m −3 ) was calculated as where the particle density ρ s was assumed to have the value of 1500 kg m −3 (Redding and Devito, 2006). Air-filled porosity f a was determined at each matric potential as the difference of total porosity and the respective volumetric water content.

115
Water retention properties were also characterized using the van Genuchten model (van Genuchten, 1980) where θ r is the residual water content, θ s is the saturated water content, and α and n are empirical fitting parameters. The fitting was performed for the average values of θ of the seven samples from each depth by applying the Levenberg-Marquardt algorithm in SciPy (Virtanen et al., 2020). Because the pressure range used in the water retention measurements was rather 120 narrow and therefore inadequate for a proper fit, the fitting procedure was simplified by setting the saturated water content to equal the total porosity and setting the residual water content to zero as suggested by Weiss et al. (1998).

Three-dimensional µCT imaging
In short, X-ray micro-computed tomography involves taking two-dimensional X-ray photographs of an object from multiple angles and then using a filtered back-projection algorithm to reconstruct the 3D volume of the sample. With the method it is 125 possible to get information about the internal structure of the sample noninvasively and to visualize the sample and assess its structural characteristics quantitatively. For example in soil studies, it is possible to determine the pore characteristics (porosity and pore size distribution), grain size distribution, and moisture distribution inside the sample (Taina et al., 2008;Helliwell et al., 2013).
Soil samples (see Sect. 2.1) were scanned in the micro-CT laboratory in the University of Helsinki with the GE Phoenix

130
Nanotom system. The final voxel (cubic 3D image element) size after reconstruction was 50 µm, and the data were stored in an unsigned 16-bit integer representation. The size of the resulting 3D images was 1142 by 1142 by 1152 voxels. Some darkening was observed in many of the images near the top and bottom of the cylindrical samples, which may have resulted from defects in µCT image reconstruction.

135
In the µCT image preprocessing stage, the 16-bit 3D grayscale images were converted to 3D binary images that represent the void and solid volumes of the samples. In this context, solid volume stands for the volume occupied by water or organic matter.
The conversion was done using the Python image processing packages SciPy ndimage (Virtanen et al., 2020) and scikit-image (van der Walt et al., 2014) and the image analysis toolkit PoreSpy (Gostick et al., 2019).
First, the 3D grayscale images were straightened through rotation and cropped to a size of 1000 by 1000 by 1000 voxels. A 140 cylindrical peat volume (height 1000 voxels, diameter 1000 voxels) excluding the acrylic cylinder was separated using PoreSpy.
Before noise filtering and binary segmentation, the 16-bit images were linearly mapped to an unsigned 8-bit representation.
The 16-bit to 8-bit mapping interval was selected for each image by visual inspection of the grayscale histogram so that the long, shallow tails of the intensity distribution, which were mainly generated by noise, were removed. The resulting 8-bit grayscale images were then filtered for noise reduction with a 3D median filter with a radius of 2 voxels (Fig. 1a). The intensity 145 contrast between air-filled regions and the regions containing water or organic matter in the images was often rather low ( Fig.   1b). Furthermore, not all the grayscale histograms were bimodal so that the intensity values for void and solid regions could not be readily distinguished from the intensity distributions. The segmentation of images into void and solid volumes was performed using the widely utilized Otsu's global thresholding method (Otsu, 1979) (Fig. 1c). Finally, isolated solid volumes were removed from the binary images using a method for finding disconnected voxels in PoreSpy.

Pore network extraction
Pore networks were extracted from the binarized 1000 3 -voxel images using a marked-based watershed segmentation method (Gostick, 2017) available in PoreSpy (Figs. 1d and 2). The extraction method has been designed to have good performance also for materials with a high porosity. It generates the topology of the pore network by dividing the void space into individual pore regions and determining the locations of pore throats, that is, the two-dimensional interfaces between adjacent pores, 155 and the connections between the pores. The method facilitates the subsequent determination of pore network geometry, which includes, for example, pore volumes, pore-to-pore distances, and throat diameters. Feature resolution is generally about twice the voxel size in µCT imaging (Stock, 2008), implying that the size of the smallest distinguishable feature was 100 µm.

Pore geometry
The size of a pore was characterized by its volume, which was determined by counting the number of voxels in an individual 160 pore region. The diameter of a pore was defined as the diameter of the largest sphere that fits inside the pore region. Similarly, the throat diameter was defined as the diameter of the largest circle that fits inside the throat region. Further, the equivalent pore diameter, defined as the diameter of a sphere with the same volume as the pore, was used in pore size classification in the study. Because the centroids of adjacent pores and the centroid of the throat between them were usually not collinear, the distance d between adjacent pores was determined as the sum of the distances between the centroids of each pore (p 1 and p 2 , The porosities of the binary images were calculated using the cylindrical sections of the 1000 3 -voxel binary images. In the calculations, it was assumed that the sample surfaces had been in level with the ends of the acrylic cylinder in the initial, 170 saturated state. The image porosities were determined as the ratios of the number of void voxels to the number of total voxels in the cropped image section covering the inner space of the cylinder. It was thus assumed that the shrinkage of the peat sample resulted in the displacement of air-filled void space within the sample into the space between the sample and the walls of the cylinder at −10 kPa matric potential, which ensured that the effect of shrinkage was included in the calculation. The extracted pore system can be partitioned into clusters of interconnected pores and a set of isolated individual pores.

175
The largest of these clusters, which is assumed to be the only one that extends through the applied sample domain, is defined as the pore network. The total pore system, including also smaller pore clusters and isolated pores, is thereafter referred to as total pore space. Network porosity is defined as the ratio of the sum of the volumes of individual pores in the network to the total volume of the applied sample domain. The volume of the total pore space is slightly smaller than the volume of the void space of the corresponding section of the binary image because the network extraction algorithm tends to discard some of the 180 smallest isolated pores especially adjacent to the borders of the image domain.

Water retention simulation
Water retention simulation was performed employing the algorithm for drainage percolation in the open-source pore network modeling package OpenPNM (Gostick et al., 2016). In the algorithm, a pore network is initially filled with a defending fluid (water). An invading fluid (air) enters the network through inlets located in a specified boundary region of the network and 185 gradually replaces the defending fluid under increasing external pressure as in a porosimetry experiment. Fluid invasion is access-limited, which means that the invading fluid can only enter the throats that are directly connected to the inlet. The pressure P needed to force the invading fluid to penetrate a throat and enter the adjoining pore is determined by the Washburn equation as where γ is the interfacial tension between the invading and defending phases (0.72 mN m −1 for an air-water interface), β is the contact angle of the invading fluid, and D is the diameter of the throat. The contact angle was assumed to be 180 • in the simulations, which gave the maximum limit of the throat entry pressure. The entry pressures corresponding to different contact angles can be readily estimated because the Washburn equation is linear with respect to cos(β). The air-filled porosity at each external pressure can be calculated as the volume fraction of air-filled pores. Water imbibition was simulated with the 195 percolation algorithm by using site percolation, in which fluid invasion pressure is controlled by pore diameters instead of throat diameters.
The network domain size used in the water retention simulations had to be as close to the total sample size as possible so that comparison with the measured retention curves would be reasonable. Thus, only 100 voxels, representing 5 mm slices, were excluded from the top and bottom of the sample images in order to exclude the roughness of the sample surfaces and the four subnetworks were performed using the same pressure steps in each simulation. The combined air-filled porosities at each pressure step were then calculated to represent the water retention characteristic in the total cylindrical network domain with a height of 800 voxels and a diameter of 1000 voxels. The maximum external pressure applied in the simulations, 2.88 kPa, was determined by the minimum throat diameter, which was 100 µm.

Network metrics 210
In addition to the effects of the vertical shrinkage, horizontal shrinkage created continuous void space between the sample and the cylinder walls in most of the sample images. In order to exclude those excess void regions, a centered cubic subregion with a side length of 600 voxels (30 mm), hereafter referred to as a subsample, was selected from each sample image for the analysis of pore size distribution and connectivity. Hence, the results characterized better the actual inner structure of a sample. The largest connected pore cluster, or the pore network, extracted from the cubic subregion was used in the network 215 connectivity analyses. The selection of only the largest cluster was justifiable because it contained on average more than 80 % of the total pore system volume and was the only cluster that extended through the network domain from top to bottom in each pore system. In order to determine the air-filled volume fraction of a pore network at different external pressures, we also performed drainage and imbibition percolation simulations for the cubic-domain networks.
We used the network analysis package NetworkX (Hagberg et al., 2008) for the estimation of network connectivity metrics.

220
The pore coordination number (the degree of a node in graph theory) gives the number of connections to an individual pore, or in other words, the number of throats emanating from a pore. The local clustering coefficient of a pore A is defined as the probability that two pores that are connected to A are also connected to each other. The network average clustering coefficient, the average of all local clustering coefficients (Watts and Strogatz, 1998), can be considered as the probability that two adjacent pores of a random pore are connected to each other. The closeness centrality of a pore is defined as the reciprocal of the average 225 of the shortest path lengths from the pore to every other pore in the network (Freeman, 1978). The calculated pore-to-pore distances were used as the edge weights in the calculation of the path lengths. The closeness centrality of the pore network was calculated as the average value of the closeness centralities of all the pores. A high network closeness centrality indicates that the overall global connectivity of the pore network is fairly high (van der Linden et al., 2019).
Geometrical tortuosity and betweenness centrality are, respectively, network measures related to the transport properties of 230 a spatial network in a certain direction and as a whole. To characterize properties related to the connectivity through the peat pore network between the opposite surfaces of the network domain, artificial boundary pores were added to the pore network to represent the interfaces at the surfaces of the samples using PoreSpy.
The geometrical tortuosity or path length tortuosity of a pore network is defined as the average value of the ratio of the lengths of the shortest paths between each pair of pores located at the opposite boundaries of the network domain to the 235 straight-line distance between the opposite boundary planes (Lindquist et al., 1996;Clennell, 1997). In the calculation of geometrical tortuosity, the shortest paths between boundary pores were determined using Dijkstra's algorithm. Geometrical tortuosity was determined separately for the vertical and for the horizontal direction. Both perpendicular horizontal directions were included in the calculation of horizontal geometrical tortuosity.
The betweenness centrality of a pore A is generally defined as the ratio of the number of shortest paths between all the 240 pairs of pores in the network that traverse A to the total number of the pairs of pores that do not include A (Freeman, 1977).
In this study, we also defined and quantified the top-bottom betweenness centrality by including only the shortest paths that connected the top boundary pores with the bottom boundary pores. As with closeness centrality, the betweenness centrality of a pore network was determined as the average value over the pores. A high network betweenness centrality suggests that the shortest paths through a network tend to be governed by a relatively small number of different routes (van der Linden et al., . In order to save calculation time, every tenth pore was used in the estimation of the network average closeness centrality, betweenness centrality, and geometrical tortuosity in networks with more than 8000 pores.

Statistics
We applied a one-way analysis of variance (ANOVA) followed by Tukey's pairwise multiple comparison test to determine the possible influence of depth on the water retention characteristics, pore sizes, and calculated network metrics. If residual 250 normality or variance homogeneity could not be assumed, a nonparametric Kruskal-Wallis test followed by Dunn's pairwise multiple comparison test or Welch's ANOVA followed by Games-Howell pairwise multiple comparison test, respectively, were applied instead. A paired sample t-test was applied to analyze the difference between vertical and horizontal geometrical tortuosity. The statistical analyses were conducted with the statistical function module in SciPy and the Python packages statsmodels (Seabold and Perktold, 2010), scikit_posthocs (Terpilowski, 2019), and hypothetical (Schlegel, 2020).

Air-filled porosity
The void fractions of the µCT images were generally somewhat smaller than the measured air-filled porosities of the samples (Fig. 3). The discrepancy was highest in the middle layer samples. The measured vertical shrinkage of the samples at −10 kPa matric potential decreased with depth, being on average 6.3 % in the top layer samples, 3.7 % in the middle layer samples, and network volumes, on the other hand, were 0.2-3.0 %, 6.3-35.1 %, and 6.6-27.8 % smaller than the total pore space volumes, as some of the void space consisting of small, isolated void regions was omitted during network extraction. The pore networks were the only connected pore clusters that extended vertically through the network domain in each pore system.

Water retention and air invasion dynamics
The total porosity of the peat samples differed significantly between sampling depths (Welch's ANOVA, F (2, 10.56) = 7.83, 270 p = 0.008), but the differences were generally rather small (Table 1). The between-sample variation in porosity was highest in the top layer. Also, the water retention characteristics differed significantly between depths (ANOVA or Kruskal-Wallis test, Table 1). Water content in all studied matric potentials was clearly lowest in the top layer and highest in the deepest studied peat layer.
The variation of air-filled porosity between samples was largest in the top layer and smallest in the bottom layer (Fig. 4). The 275 difference in the between-sample variation between depths was most pronounced under −1 kPa and −3 kPa matric potential conditions. The air-filled porosities derived from measurements and pore network simulations were rather coherent at different matric potentials for some of the samples, but considerable variation existed in the difference between measurement-based and simulated values in many samples at all depths. In the pore network simulations, the external pressure range extended only to about 3 kPa, which corresponds to the minimum throat diameter of 100 µm detected in the µCT imaging. In the samples with 280 no notable shrinkage in any direction, such as top layer samples 3, 6, and 7 and the middle layer sample 7, the percolation simulation matched the measured air-filled porosity values well at both −1 kPa and −3 kPa matric potentials. The hysteresis $LUILOOHGSRURVLW\GHWHUPLQHGIURP ZDWHUUHWHQWLRQPHDVXUHPHQW  effect is clearly seen in the drainage-imbibition simulations (Fig. 5). The air-filled porosity at a certain matric potential was higher during water imbibition than during drainage.

Pore and throat size distributions 285
The pore sizes of the connected pore networks obtained from the 600 3 -voxel subregions of the peat sample images ranged from 4 × 10 −4 mm 3 to 75 mm 3 , which corresponds to an equivalent diameter range of 0.09 mm to 5.2 mm (Fig. 6a). The median pore size of the top layer networks (median equivalent diameter 0.66 ± 0.06 mm) was significantly larger than that of the middle (0.58 ± 0.01 mm) and bottom (0.59 ± 0.02 mm) layer networks (Kruskal-Wallis test, H(2) = 9.67, p = 0.008). Also, the between-sample variation in pore size distribution was substantially higher in the top layer (the medians of the equivalent 290 diameters 0.579-0.765 mm) than in the deeper layers (medians 0.573-0.589 mm and 0.565-0.617 mm in the middle and bottom layers, respectively). The maximum pore sizes were generally highest in the top layer networks, but large pores also existed in some of the bottom layer networks (Fig. 6c-e). A larger total volume of a pore network was related to a larger average pore size in the top layer networks but not in the deeper layer networks.   The throats with the smallest detectable diameter (100 µm) were the most abundant at all depths, and the fraction of wider 295 throats decreased with depth (Fig 6b). However, the difference between the mean throat diameters at different depths was marginally nonsignificant (ANOVA, F (2, 18) = 2.82, p = 0.09).

Network porosity and connectivity metrics
The porosity of the pore network and all the pore network metrics differed significantly among depths (ANOVA or Kruskal-Wallis test, p < 0.05, Fig. 7). Network porosity and connectivity were clearly highest in the top layer. The connected pore 300 networks did not extend over the whole cubic network domain in most of the middle and bottom layer subsamples and in one of the top layer subsamples, which decreased the obtained network porosity. All the network metrics of the top layer differed significantly from those of the deeper layers (p < 0.05), whereas no significant difference was observed between the middle and bottom layer subsamples.
The vertical geometrical tortuosity was significantly higher than the horizontal geometrical tortuosity in the top layer net-305 works (mean difference 0.10; two-tailed paired sample t-test, t(6) = 3.67, p = 0.01), but no significant difference was found at the deeper layers. Furthermore, the variation of geometrical tortuosity in the vertical direction was slightly higher than in the horizontal direction (Fig. 7h,i). In the top layer, vertical geometrical tortuosity decreased with increasing network average coordination number, clustering coefficient, and closeness centrality and increased with increasing network average betweenness centrality, but the behavior was very different and inconsistent in deeper layers (Table 2).   The total volume of the connected pore space increased as a function of the number of pores in the network, but the variation was higher in the top layer subsamples than in the middle and bottom layer subsamples (Fig. 8a). Pore sizes tended to increase with increasing porosity in the top layer subsamples, whereas this connection was not found for the deeper layers. The network average coordination number, characterizing local network connectivity, increased with increasing porosity in the top layer subsamples, in which the connected pore space was rather evenly distributed in the network domain (Fig. 8b) Figure 8. Linear relationships among selected pore network properties and average pore network connectivity metrics.
The network average closeness centrality increased with increasing coordination number in the top layer except for the subsample with smaller network dimensions and largely also in deeper layers (Fig 8d). By contrast, high local network connectivity indicated a low betweenness centrality especially in the top and bottom layers (Fig. 8e). The connected pore cluster extended were rather uniform (Fig. 8g). In the deeper layers, horizontal and vertical tortuosity were unrelated. Also, vertical geometrical tortuosity and top-bottom betweenness centrality were rather well correlated in the top layer but not in deeper layers (Fig. 8h). Figure 9 describes the volume fraction of the connected, air-filled pore network of the total pore space at different external pressures. Values less than 1 indicate that a part of the pore space has been isolated from the surrounding volume and O 2 330 supply has ceased. Eventually, the isolation can lead to the formation of an anaerobic pocket. The volume fraction of the connected network was calculated for both imbibition (wetting) and drainage (drying). In imbibition, half of the network remained connected at 0.5 kPa in all layers, whereas in the drainage simulation this required 1 kPa in the top layer and 3 kPa in the middle and bottom layers. The layers behaved rather similarly in imbibition, but in drainage there were clear differences between the layers.

335
The fraction of the total pore network volume of the total pore space volume, which corresponds to the volume fraction at an external pressure of 3 kPa (Fig. 9h), was largest in the top layer subsamples, extending from 93.6 % to 99.9 %. In the deeper subsamples, 7.4 to 51.9 % of the total pore space was disconnected from the largest pore cluster. Thus, a significant fraction of macropore space was inactive in the drainage and imbibition simulations in the deeper layer subsamples, which indicates an even larger pore volume available for anaerobic pocket formation.

Evaluation of image and network analysis in macropore characterization
The combination of X-ray tomography, image analysis, and network analysis provides detailed information on pore structures, connections, and topology that cannot be obtained through traditional laboratory methods. These properties determine gas exchange in peat and are therefore essential in regulating biological activity. One of the aims of our study was to evaluate the 345 applicability of image analysis and network extraction methods to peat structure characterization. Overall, the image-derived porosities and the simulated water retention characteristics corresponded rather well with the laboratory measurements given the limitations imposed by the applied imaging resolution. The smallest pore diameter detectable in the µCT imaging was 100 µm, corresponding to a matric suction of approximately 3 kPa. This resolution is sufficient for accurately describing diffusional gas transport in peat soils, where matric suction typically remains low and pores smaller than 100 µm are generally 350 water-filled. By comparison, macropores have been found to dominate water and solute transport in peat (McCarter et al., 2020), and a similar resolution is also considered adequate for simulating hydrological transport processes in water-saturated peat (Gharedaghloo et al., 2018).
Good performance of image segmentation is a crucial prerequisite for a successful application of µCT image analysis and subsequent quantitative analysis tools (Iassonov et al., 2009). The binary segmentation stage succeeded fairly well in our study.

355
The sample void fractions obtained through µCT image analysis were in good agreement with the respective air-filled porosities derived from laboratory measurements (Fig. 3). The pore network extraction stage introduced another level of complication in macropore volume characterization. The porosities of the networks used in the water retention simulations were further diminished with respect to the corresponding image porosities, especially in deeper layers. The main reason for this was that a considerable fraction of pore space was dis-370 connected from the active pore network, especially in some of the middle and bottom layer samples. This was because narrower void connections may have not been detectable in the images, which may have resulted in some pore clusters or individual pores becoming isolated from the main connected cluster. In addition, some of the reduction of pore space connectivity and of the total, combined volume of connected pore space resulted from the division of the total network domain into four discrete regions.

375
The shrinkage behavior of peat further obscured the determination of network porosity. In the networks with a centered cubic domain, the average porosity was slightly lower in the middle layer subsamples than in the bottom layer subsamples (Fig. 7a), while the opposite was the case in the larger-domain cylindrical networks (Fig. 4). In addition to the spatial heterogeneity of air-filled porosity within the samples, the difference can be explained by the variation of horizontal shrinkage between depths.
Higher shrinkage in the middle layer samples decreased the air-filled porosity of the samples, which was reflected in the air-380 filled porosity of the cubic network domain, which did not include the void region between the shrunk sample and the cylinder walls.
Peat shrinkage also affected the results of the water retention simulations. Ideally, the simulated air-filled porosity at the maximum external pressure should have been close to the corresponding measured air-filled porosity at −3 kPa matric potential, which corresponds to the minimum detectable pore throat dimension in the images (Fig. 4). However, because of the 385 limited image resolution, the constructed pore geometry may not totally represent the actual void space geometry under the conditions of −3 kPa or higher matric potentials because the shrinking of the samples may have resulted in a decrease in void space. Thus, a fraction of the pore throats that were air-filled a −3 kPa matric potential may have shrunk so that they were not detectable in the µCT images constructed at −10 kPa matric potential. This may have generated disconnected pore space and decreased the total volume of the extracted pore network. Also, shrinkage may have decreased the dimensions of the pore 390 space so that a higher external pressure was needed for air invasion in the simulations. Conversely, the horizontal shrinkage of some of the samples created continuous void space near the cylinder wall at −10 kPa matric potential, and thus the extracted pore network contained pore space that had presumably not been present at higher matric potential conditions.

Network connectivity metrics related to peat structure and gas transport
Van der Linden et al. (2016) highlights the need for multiscale measures of pore space topology and connectivity and the 395 relationship of these descriptors to macroscopic transport processes in a porous medium, such as peat. Geometrical tortuosity is a structural characteristic of a porous medium (Clennell, 1997). It gives an estimate of the average path length through a pore network in a specified direction, which is a direct proxy for network transport efficiency. Therefore, it can be used as a benchmark measure of the applicability of different network metrics to characterize the transport properties of a network and to estimate the efficiency of macroscopic transfer processes in a network. We compared several kinds of network metrics with 400 the gas transfer capacity estimated by geometrical tortuosity. The pore coordination number and the clustering coefficient were used to characterize the local connectivity of a pore network, whereas closeness centrality and betweenness centrality were used to describe connectivity in the network scale. Furthermore, top-bottom betweenness centrality characterizes the positioning and shape of gas transport routes through the porous medium in a certain direction, thus also describing the efficiency of gas transfer through the medium. It is essential to determine how the local and global network connectivity measures are related 405 to the network gas transfer capability estimated by vertical geometrical tortuosity. According to our results, high average local pore connectivity is not always reflected as high global connectivity or gas transfer capability of the pore space. The distribution and spatial coverage of the connected pore space within the porous medium was found to regulate the applicability and comparability of the local and global measures and also their relationship to porosity and tortuosity.
Generally, a higher network porosity implied higher local pore connectivity characterized by the coordination number and 410 the clustering coefficient (Fig. 8b,c). However, if the connected pore space was concentrated in a smaller region in the network domain, which was the case in most of the deeper layer subsamples, local connectivity was rather high even though the porosity of the network was relatively low. Likewise, the relationship between local and global connectivity measures was largely dictated by the shape and topology of the network. In general, higher local network connectivity was reflected as a higher average closeness centrality (Fig. 8d), which means that the path lengths between pores shortened when more alternative 415 paths were available. The average betweenness centrality mainly decreased with increasing average coordination number at all depths, which indicates that the shortest routes between pairs of pores were spread out more widely when the number of pores and their connectivity was higher (Fig. 8e). However, some variation existed especially at deeper layers, where the porosity was lower and the spatial extent of the network was typically smaller.
The location of a pore largely determines its centrality in a spatial network because the probability of a pore being part of 420 the shortest route between two other pores is highest near the centroid of the network (Barthélemy, 2011). High abundance of dead-end pores or pores located near the extremities of the network decreases the average betweenness centrality because these pores do not belong to the shortest paths between other regions. In addition, the closeness centrality of a pore (the reciprocal of the average of the shortest path lengths from a pore to every other pore) is related to the dimensions of a network within its domain. If a network covers only a fraction of its domain or if it contains dense local clusters with narrow pore channels 425 between them, the average closeness centrality may become relatively high even if the average coordination number is rather low.
The ratio between the two local connectivity measures may reflect the structure of a pore network in regard to gas transfer efficiency. A high ratio of the network average clustering coefficient to the network average coordination number indicates that the number of three-pore loops is relatively high and may also imply that even larger clusters of interconnected nearby 430 pores are abundant in the network. This may be reflected as a high gas transfer capacity and good resilience to disturbances.
By contrast, a low ratio may result from a high abundance of long nonbranching pore conduits in the network, which may diminish the amount of alternative transport routes within the network and lead to a rapid suppression of gas transport if pores become clogged.
The geometrical tortuosity of a pore network in a certain direction is strongly dependent on the geometry of the connected  Top-bottom betweenness centrality is a network measure defined in this study to characterize the distribution of transport paths through a network. A high network average top-bottom betweenness centrality may imply that the optimal flow routes between the top and bottom of the network domain are governed by a small number of different paths and that transport through the network may therefore be sensitive to individual disturbances within the flow routes. Our results showed a correspondence between geometrical tortuosity and top-bottom betweenness centrality in the top layer networks, most of which had a spatially 445 rather even pore distribution (Fig. 8h). If the average shortest path length between the top and bottom boundary pores is low (the vertical geometrical tortuosity is low) and the boundary pores are located evenly on the top and bottom surfaces of the domain, the shortest paths between different pairs of boundary pores are also located more evenly throughout the network (the average top-bottom betweenness centrality is low). Because of a smaller number and more localized spatial distribution of boundary pores and more variable network topology and geometrical structure, there was evidently no correspondence between vertical 450 geometrical tortuosity and top-bottom betweenness centrality in the deeper layers.
A low geometrical tortuosity of a pore network suggests that gas transfer through the network is efficient in a specified direction. In the top layer subsamples with a higher network porosity and a stronger correlation between vertical tortuosity and network porosity (Fig. 8f), a lower vertical geometrical tortuosity was strongly related to network metrics values that indicated higher network connectivity (Table 2). By contrast, such a correlation was not found for any local or global connectivity metrics 455 in other layers. Thus, the applicability of network connectivity metrics to the description of network transport properties seems to be highly sensitive to the spatial distribution of the network within its domain. The average closeness centrality of a pore network has been found to be a proxy for The structural anisotropy of peat, characterized by the difference between vertical and horizontal geometrical tortuosity, can be used to estimate the diffusion capability of a gas in different directions through peat. In this study, the pore structure was slightly anisotropic in the top layer, but no anisotropy was found in deeper layers. Thus, structural anisotropy decreased with 465 lower network porosity and pore connectivity. However, the geometrical tortuosity was highly variable between samples in deeper layers because of the very heterogeneous spatial distribution of the pores in the middle and bottom layer networks. The variation of anisotropy between depths suggests that chains of pores existed between horizontally orientated, less degraded plant residues, giving rise to a low horizontal geometrical tortuosity. In deeper levels with more degraded and compacted peat, the horizontally oriented pore chains had fractured and collapsed, and the structural anisotropy had diminished. These findings 470 are in line with Kruse et al. (2008) and Liu et al. (2016) who found that the anisotropy of the hydraulic conductivity of fen peat decreased with increasing degradation. According to our results, the orientation of diffusion paths in pore networks is such that it does not promote or restrain gas transfer towards the atmosphere in deeper, more degraded peat layers. In less degraded peat near the surface, the impact of the orientation of plant litter on the primary diffusion direction may be more pronounced. However, higher pore connectivity and porosity in the less degraded peat may outweigh the hindering effect of

Conceptual implications for anaerobic pocket formation and methane dynamics
Emissions of CH 4 from northern peatlands are characterized by a large spatiotemporal variation (Abdalla et al., 2016;Rinne et al., 2018). The CH 4 efflux usually occurs episodically and in hotspots (Lai, 2009), and the supply of O 2 and gas diffusion conditions in the soil profile are the main abiotic factors affecting CH 4 emissions (Xu et al., 2016). We argue that the structure, 480 topology, and behavior of peat macropore networks above the WT can be a good candidate for explaining the spatiotemporal variation of peatland CH 4 emissions. The following conceptualization describes the pore network, anaerobic pocket formation, and CH 4 production and transport: When the pore network is internally connected and open to the atmosphere, the supply of O 2 to the soil is adequate and facilitates aerobic decomposition. When some of the pores become blocked by water, the O 2 supply to the soil gradually decreases, as the number of air-conducting pores decreases and the air flow paths become longer 485 and more tortuous. With increasing water content, part of the pore network becomes isolated and the O 2 supply is prevented.
Next, microbial activity consumes the trapped O 2 , and the microbial metabolism changes to alternative electron acceptors. In the final stage, the production of CH 4 onsets. An anaerobic pocket has formed. When soil water content decreases again, the network becomes connected, and the CH 4 trapped in the pocket is released and can be detected as a burst of CH 4 emission at the soil surface.

490
The existence of isolated anaerobic macropore clusters above the WT may also enable CH 4 production in newly aerated layers during a drying period (Knorr et al., 2008;. Large macropores are drained and exposed to O 2 as the WT declines, but pore clusters with narrower throats adjacent to the air-filled regions remain anoxic if diffusional O 2 supply to these water-filled pores is inadequate. If the isolated macropores are formed close to the peat surface, the time span that a substance needs to diffuse to the atmosphere is relatively short and, consequently, there is a limited 495 exposure to oxidation in the aerated layer. Thus, methanogenesis in isolated near-surface macropores may significantly increase the atmospheric emission of CH 4 (Estop- Aragonés et al., 2013). In addition, anaerobic pore clusters may be the regions where anaerobic microbes may survive above the WT and facilitate the onset of CH 4 production when the WT rises again (Kettunen et al., 1999).
Our drainage and imbibition simulations indicated that hysteresis may affect the evolution of the fraction of peat pore 500 space that is isolated from the atmosphere. The volume fraction of the air-filled pore space connected to the atmosphere at a specific matric potential was considerably smaller during drying than during wetting. This may further provide support for the conception that hysteresis may also affect the formation and destruction of anaerobic pockets and the temporal CH 4 dynamics in the unsaturated layer. Figure 9 describes conditions under which the anaerobic pockets are likely to occur. When the peat is drying, the macropore network remains largely water-filled until a matric potential of −2 to −3 kPa, which under hydrostatic 505 equilibrium corresponds to a distance of 20-30 cm from the WT. This promotes a high abundance of anaerobic pockets in the unsaturated layer close to the WT. The thickness of this layer is greater at deeper depths, where the average pore throat dimensions are smaller. In a hectare scale, this suggests that 2000-3000 m 3 of unsaturated peat is potentially active in CH 4 production. Under wetting conditions, the macropore network remains largely air-filled until the matric potential is higher than −1 kPa or the distance to the WT is less than 10 cm. This means that the thickness of the layer where anaerobic pockets can 510 be formed when the WT is rising is less than 10 cm. The large difference between the thickness of the pocketing layer under drying and wetting conditions is likely to cause a short-term mismatch between peatland CH 4 emissions and WT observations.
For example, the hysteresis effect may contribute to the lagged response of CH 4 flux to a rising WT that has been observed in peatland ecosystems (Kettunen et al., 1996;Goodrich et al., 2015).
The average pore volumes and pore throat diameters were smaller at deeper depths than in the near-surface peat (Fig. 6). The 515 higher volume fraction of large pores in the top layer compared to deeper layers is also seen in the air-invasion curves as a higher relative air-filled porosity at −1 kPa matric potential (Fig. 4). The reduction of pore space volume and pore dimensions with depth, which is due to increasing degree of decomposition of peat and higher compression by overlying matter (Rezanezhad et al., 2016), is a typical feature of peat soil (Rezanezhad et al., 2010;Weber et al., 2017;Gharedaghloo et al., 2018). Thus, the macropore network is drained at lower matric potentials in deeper peat than near the soil surface. In addition, the variation 520 of porosity and the variation of average pore volume between the samples were largest in the top layer (Fig. 8a). Increasing degree of decomposition deeper in the peat profile also results in the homogenization of peat and a smaller spatial variability of the pore size distribution (Weber et al., 2017). In our samples, virtually all macropores were connected at −3 kPa matric potential in the top layer, whereas roughly one-quarter of the macropore volume was still isolated in deeper layers (Fig. 9). As a result, the width of the layer favorable for anaerobic pocket formation is greater in deeper, more decomposed peat. However, 525 the lower macroporosity and smaller average volume of pores may reduce the total volume of anaerobic pockets and limit the capacity for methanogenesis in the pocketing layer at deeper depths.
In our study, the peat samples originated from a drained forested peatland site. The total porosity range of the studied peat layers was well in line with the values and the vertical variation reported for the near-surface layers of drained peatlands in the literature (Päivänen, 1973;Minkkinen and Laine, 1998). However, the peat soil of a drained peatland differs in physical 530 characteristics and pore size distribution from that of an undisturbed peatland (Liu et al., 2020). Pore deformation due to peat shrinkage induced by drainage is partly irreversible (Price, 2003), and the pressure of a growing tree stand may further compact the peat layer and increase its bulk density in deeper layers (Minkkinen and Laine, 1998). Also, the saturated peat layer below the WT may be further compressed due to the weight of overlying peat that has become drier and lost the support by buoyancy (Hooijer et al., 2012;Sloan et al., 2019). These processes decrease the macroporosity of the peat profile in drained peatlands in 535 comparison to undisturbed peatlands (Liu et al., 2020). Therefore, conditions for methanogenesis may be even more favorable in undisturbed peatlands because of a potentially larger volume of macropore space available for anaerobic pocket formation.

Conclusions
The network analysis of the peat pore system enabled by µCT imaging and a network representation of peat macropores demonstrated fundamental differences in peat pore structure and macropore characteristics between topsoil and deeper soil 540 layers. Pore space was more connected and routes through the peat matrix were less tortuous in the top layer than in deeper layers. Decreasing pore connectivity with depth was accompanied by a lower number of macropores, smaller macropore volumes, and narrower pore throats. This may indicate that the rate of gas diffusion in the air-filled pore space is reduced in deeper peat layers. The results also suggest that local and global network connectivity metrics might be used to estimate the efficiency of diffusional gas transfer in the air-filled pore space of peat. However, we highlight that connectivity metrics should 545 be evaluated with caution because not only pore connectivity but also the extent and spatial distribution of connected air-filled pore space regulate the gas transfer capabilities.
Pore network analysis may also provide new insights into the impact of pore structure and pore space connectivity on conditions regulating CH 4 production and transfer in peat. As the WT is generally close to the surface and low suction conditions prevail in peatlands, small pores remain continuously water-filled and diffusional gas transfer occurs in the air-filled macropore 550 network above the WT. We argue that the complex pore structure and the vertical variation in the pore characteristics of peat may promote the formation of anaerobic pockets above the WT during fluctuations of soil water content. When the WT finally declines and soil water content decreases, CH 4 produced in these pockets can be released rapidly via air-filled macropores.
In addition, hysteresis was found to regulate the thickness of the zone favorable for anaerobic pocket formation. Under the same WT, the pocketing can be distinctively different depending on whether the peat is wetting or drying. This may provide 555 an explanation for the observed hotspots and episodic spikes of CH 4 emissions in peatlands. Most importantly, a pore network representation of peat macropore structure enables the application of pore network modeling, which is a useful method for the pore-scale description of CH 4 production and transfer processes in peat and for the investigation of relations between peat pore structure and CH 4 dynamics.
Code and data availability. The data that support the findings of this study are available in Github [https://github.com/pjkiuru/macropore_ Author contributions. AL, TG, and IU developed the idea and designed the study. AL and MP collected the samples and performed the water retention measurements. TG and MR organized the µCT imaging and 3D reconstruction. PK processed the images, designed and 565 conducted the computations and simulations, and analyzed the data. PK and LK performed the statistical analysis. PK, AL, and MP wrote the manuscript with significant contributions from TG, VG, and IU. All other authors provided edits and comments on the paper. AL and MR are responsible for the funding acquisition.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This research has been supported by the Academy of Finland (grant nos. 325168 and 325169). LK holds a Marie the funding from the Academy of Finland to strengthen university research profiles in Finland for the years 2017-2021 (funding decision 311925). MR acknowledges SRC at the Academy of Finland (SOMPA, no. 312932) and EU-H2020 (VERIFY,no. 776810). This work used services of the Helsinki University X-Ray Micro-CT Laboratory, funded also by Helsinki Institute of Life Science (HiLIFE) under the HAIP platform.