Controls on Saturated Hydraulic Conductivity in a Degrading Permafrost Peatland Complex

Permafrost peatlands are vulnerable to rapid structural changes under climatic warming, including vertical collapse. Peatland water budgets, and therefore peat hydraulic properties, are important determinants of vegetation and carbon fluxes. Measurements of hydraulic properties exist for only a limited number of permafrost peatland locations, primarily concentrated in North America. The impacts of thaw‐induced collapse upon properties such as horizontal saturated hydraulic conductivity (Kh), and thus lateral drainage, remain poorly understood. We made laboratory determinations of Kh from 82 peat samples from a degrading Swedish palsa mire. We fitted a linear mixed‐effects model (LMM) to establish the controls on Kh, which declined strongly with increasing depth, humification and dry bulk density. Depth exerted the strongest control on Kh in our LMM, which demonstrated strong predictive performance (r2 = 0.605). Humification and dry bulk density were influential predictors, but the high collinearity of these two variables meant only one could be included reliably in our LMM. Surprisingly, peat Kh did not differ significantly between desiccating and collapsed palsas. We compared our site‐specific LMM to an existing, multi‐site model, fitted primarily to boreal and temperate peatlands. The multi‐site model made less skillful predictions (r2 = 0.528) than our site‐specific model, possibly due to latitudinal differences in peat compaction, floristic composition and climate. Nonetheless, low bias means the multi‐site model may still be useful for estimating peat Kh at high latitudes. Permafrost peatlands remain underrepresented in multi‐site models of peat hydraulic properties, and measurements such as ours could be used to improve future iterations.

here we focus on the hydraulic properties of the peat matrix.Peat K sat has been widely measured in boreal and temperate regions in samples with volumes on the order of 10-100s of cm 3 .These data have been used to develop multi-site pedotransfer functions, for which r 2 is typically between 0.3 and 0.5 (e.g., Lennartz & Liu, 2019;Liu & Lennartz, 2019;Wang et al., 2021), although a recent multi-site model by Morris et al. (2022) has demonstrated greater predictive power (cross-validated r 2 = 0.747).However, peat K sat has rarely been measured in permafrost peatlands outside of North America, and permafrost dynamics are not currently represented in multi-site pedotransfer functions of peat K sat (e.g., Morris et al., 2022), so it is unclear whether such models can be reliably applied in permafrost regions.
Existing K sat measurements from permafrost peatlands are concentrated in specific geographic locations, particularly around Scotty Creek Research Station, Northwest Territories, Canada (Ackley et al., 2021;Gharedaghloo et al., 2018;Nagare et al., 2013;Quinton & Baltzer, 2013;Quinton et al., 2008Quinton et al., , 2019) ) and the North Slope of the Brooks Mountain Range, Alaska (O'Connor et al., 2019(O'Connor et al., , 2020)).Scotty Creek is located in the discontinuous permafrost zone where peat-covered frost mounds are prevalent, termed palsas or peat plateaus depending on their areal extent (Zoltai & Tarnocai, 1975).Analyses of these peat plateaus have identified strong, vertical variations in near-surface K sat across several orders of magnitude: an upper zone (depth ≤0.1 m) of high K sat and a deeper zone (depth ≥0.2 m) of consistently low K sat (Nagare et al., 2013;Quinton et al., 2008Quinton et al., , 2019)).Vertical gradients in peat K sat are driven by changes in properties such as dry bulk density and peat humification (Morris et al., 2022;Päivänen, 1973).Deeper, older peats have typically been subjected to more decomposition and compaction, which weakens cellular structures and causes pores to collapse and become less connected (Pichan & O'Kelly, 2012;Rezanezhad et al., 2016).Other Canadian permafrost sites show similar depth relationships (Morison et al., 2017;Quinton et al., 2000Quinton et al., , 2004)), and dry bulk density has been identified as a universally strong predictor of K sat in organic soils (including peats) on the Alaskan North Slope (O'Connor et al., 2020).Palsas and peat plateaus are also susceptible to desiccation and burning, which can remove the most permeable, near-surface peats (Ackley et al., 2021).Furthermore, pore network structures in permafrost peatlands are likely to be influenced by peat floristic composition (McCarter et al., 2020) and annual freeze-thaw cycles (Liu et al., 2022).Thaw can reduce groundwater flows by lowering the water table into deep, low K sat peats, decreasing hydraulic gradients across permafrost landforms, and shifting the peat-forming vegetation toward more readily-decomposed sedges (Olefeldt et al., 2021;Quinton & Baltzer, 2013), although it remains unknown whether permafrost thaw has any direct impact on peat hydraulic properties.
Permafrost peatlands in Fennoscandia appear to be close to a climatic tipping point (Fewster et al., 2022).Palsa complexes in this region that are currently undergoing permafrost thaw exhibit a mosaic of desiccating, raised palsas with intact, underlying permafrost, and areas that have collapsed and become inundated (Swindles et al., 2015).Surface collapse would appear to have the potential to drive important changes in pore structure and peat permeability, and ongoing thaw in these sites provides an opportunity to study peat hydraulic properties at different stages of palsa degradation.To our knowledge, only Wetzel et al. (2003) have previously measured peat K sat in permafrost peatlands in Fennoscandia.They identified high dry bulk densities (up to 0.30 g cm −3 ) and declining K v with increasing depth in nine samples from Norwegian palsas.However, peat K h remains unstudied for palsas in this region, which limits our understanding for their hydrological functioning and lateral drainage following thaw.This is an important knowledge gap because peat K sat can be strongly anisotropic, varying by more than three orders of magnitude in the horizontal (K h ) and vertical (K v ) dimensions of some temperate peat samples (Beckwith et al., 2003).Furthermore, it remains unclear how generalizable existing, spatially-limited measurements from permafrost peatlands are to Fennoscandian palsas.
In this study, we aimed to characterize peat K h , and to establish its controls, in a degrading European palsa mire.Specifically, we addressed the following research questions: 1. Can combinations of easily-measured properties, such as depth, dry bulk density and the von Post hand-squeeze test for degree of humification, be used to predict peat K h in degrading palsas? 2. Can an existing model, trained primarily on boreal and temperate peats (Morris et al., 2022), reliably predict peat K h in a degrading permafrost peatland?3. Do desiccating and collapsed areas of the palsa mire exhibit characteristic differences in peat K h ?10.1029/2023WR035398 3 of 14

Study Site
Rensjön palsa mire, sometimes also known as "Railway Bog", is a degrading permafrost peatland in northern Sweden (68°05'12"N, 19°49'53"E; 488 m above sea level) (Figure 1), bordered by a regional railway line to the east, a watercourse to the south, and forested, mineral-soil uplands to the west.The site contains a mosaic of desiccating palsas with remaining permafrost, and collapsed areas and inundated fens that are permafrost-free, which characterizes the trajectory of permafrost peatland degradation in Fennoscandia (Swindles et al., 2015).Desiccating palsas have either bare peat surfaces, or surface vegetation primarily comprised of Cladonia sp., Empetrum nigrum, Andromeda polifolia, Betula nana, Salix sp. and Rubus chamaemorus.Collapsed palsas are dominated by brown mosses, feather mosses and Rubus chamaemorus, with some Sphagnum spp., Empetrum nigrum, Andromeda polifolia, and Betula nana.Thermokarst and open-fen areas contain Eriophorum spp., Carex spp., Sphagnum spp., brown mosses, and Lycopodium sp., while Salix sp. are prevalent along the site's permafrost-free, eastern margin.Palsas from Rensjön palsa mire have previously been sampled for palaeohydrological reconstructions, which have indicated peatland surface drying and reduced local C accumulation during the late-20th century (Sim et al., 2021).At the nearby Rensjön weather station (68°04'23"N, 19°50'06"E; 493 m above sea level), mean annual temperature is currently −0.7°C with mean monthly temperatures ranging from −13.3°C in January to 13.0°C in July, while annual unfrozen precipitation is 508.1 mm yr −1 (10-year averaging period: 2013-2022) (Swedish Meteorological and Hydrological Institute, 2023).

Field Sampling
We extracted shallow peat cores from 20 locations across Rensjön palsa mire in July 2022, following a similar method to Morris et al. (2019).We inserted a cylindrical section of polyvinyl chloride (PVC) downpipe of 0.5 m length × 0.1 m diameter carefully into the peat using the scissor-cut method of Green and Baird (2012) to minimize sample compression.Each core was then carefully excavated using a spade.Cores were retained within their PVC cylinders during transport to the laboratory at the University of Leeds, UK, where they were stored in an upright, refrigerated state at ∼4.0°C.Given the relatively small volumes of our peat core samples, our analyses best represent matrix-scale heterogeneity in peat K h and may somewhat underestimate the landscape-scale hydrological fluxes of the site, where the presence of voids, pipes, and preferential drainage channels may also become important.To investigate whether peat hydraulic properties vary between different stages of palsa degradation, we collected 10 cores from desiccating palsas with intact, underlying permafrost, and 10 cores from permafrost-free, collapsed palsas (Figure 1).Peat in inundated fen areas was too wet and loose to be recovered as intact samples.
At each coring location, we measured the active-layer thickness adjacent to the coring location using a ∼2 m-long steel probe.We also conducted a more extensive survey of active-layer thickness by taking 200 probe measurements across the site (Figure 1c).Active-layer thickness at coring locations on desiccating palsas ranged from 0.31 to 0.41 m (median = 0.36 m), which is similar to desiccating palsas across the site (0.22-0.46 m).No coring locations on collapsed palsas contained permafrost within the ∼2 m reach of the steel probe, which enabled deeper samples to be collected for K h determinations from collapsed areas than from desiccating palsa peats (Figure 3a).More broadly, some collapsed palsas across the site retained some permafrost at depths between 0.30 and 1.68 m, while 64 out of the 110 probed locations on collapsed palsas had no permafrost within ∼2 m of the peat surface.

Laboratory Analysis
We determined peat K h in the laboratory using a custom-built permeameter method.We froze each peat core prior to our analysis at ∼−18.0°C.We carefully cut open the PVC downpipes containing each frozen peat core using a thin-bladed angle grinder and checked the surface of each peat core for any cracks or voids.We used a pillar drill with a customized core drill bit to cut horizontal, cylindrical samples of 4.2 cm diameter from each frozen peat core.We then removed curved end faces from each subsample using a sharp knife.Samples were thawed overnight in cold storage at ∼4.0°C and were then sealed inside an acrylic split-cylinder (4.0 cm internal diameter × 12.0 cm length), which consisted of two half cylinders that allowed the sample to be inserted without damage.We fastened a thin, neoprene gasket to the edge of one of these half cylinders, to form a water-tight seal when the sample was enclosed.We lined the acrylic cylinders with petroleum jelly to prevent preferential flow along the acrylic-peat interface and used tensioned rubber O-rings to secure the two halves of the PVC split-cylinder together during K h determinations (Figure 2).
Each enclosed subsample was saturated with pool water from a blanket bog in northern England (Keighley Moor) with a similar acidity (pH range: 3.91-4.13)to our studied palsa peats (pH range: 3.95-4.42).
Once saturated, a Mariotte regulator was positioned above the sample to generate a hydraulic gradient, with an outflow tube positioned into a measuring cylinder with a "U" bend to control the outflow head (Figure 2).We took care at all times to ensure no air entered the saturated sample.We adjusted the height of the Mariotte regulator to ensure that the flow rate across the sample was similar for peats with different degrees of humification (e.g., shallow hydraulic gradients were required for highly permeable samples with open, connected pores).We ran pool water through the sample for an initial settling period of 30 min before running any tests.To measure peat K h , we recorded the volume of water released through the split-cylinder during time intervals of between 10 and 20 min and calculated K h from Darcy's Law, as described by Beckwith et al. (2003).Tests were repeated five times to evaluate any drift in measured K h over time (e.g., which may occur due to bubbles mobilizing and blocking flow at pore necks) and we deemed K h to have stabilized if at least three of the five K h measurements were within ±5% of the median; any measurements outside this range were considered outliers and were removed.Despite our careful approach, we cannot rule out the possibility of some pores clogging during the initial saturation of our samples, for example, by colloidal material (Beckwith et al., 2003).We temperature-corrected each K h measurement to 20°C (Thomas et al., 2016), to account for the temperature-dependence of fluidity.For each sample, we then calculated the between-replicate harmonic mean of K h from our standardised measurements.We used the harmonic mean to reduce the possible effects of any spuriously high values, which may have resulted from any preferential flow between the cylinder and the peat sample.In total, we measured peat K h in 82 samples, comprising 35 from desiccating palsas and 47 from collapsed palsas.
We measured peat dry bulk density (g cm −3 ) (DBD) by sampling known volumes of peat (between 1.0 and 2.5 cm 3 ) in peat-core offcuts immediately adjacent to each sample used for our K h determinations, and oven-dried these samples at 105°C overnight, following Chambers et al. (2011).We repeated this process three times for each sample and calculated the arithmetic mean of DBD.We also established the degree of peat humification for each sample using immediately adjacent peat-core offcuts, following the von Post classification method (Ekono, 1981).All laboratory measurements reported in this article are available in Fewster et al. (2023).

Model Development
Because we measured peat K h at multiple depths in each core, it is possible that our data set contains a hierarchical structure whereby measurements may exhibit greater variance between cores than within cores.Any such structured variance would violate the assumption of independent observations required by standard multiple regression.We therefore fitted a linear mixed-effects model (LMM) with core identity treated as the subject for a random intercept in an attempt to account for any between-core effects.
We developed a LMM to predict our measurements of peat K h using: three continuous fixed-effect predictors (depth, DBD and von Post score); one categorical fixed-effect predictor describing the stage of palsa degradation (two levels: desiccating and collapsed); and the random effect describing core identity.We considered von Post score to be a continuous predictor, following the rationale of Morris et al. (2019).Briefly, ordinal predictors with numerous discrete categories, such as our von Post measurements (range: H1 to H8), rarely have important deviations from linearity and multiple regression models are generally insensitive to any uneven spacing between ordinal categories (Pasta, 2009).Our K h measurements extended across four orders of magnitude, producing highly skewed, non-linear relationships with all continuous predictors.We therefore log 10 -transformed our measurements of peat K h , following Morris et al. (2022).We also log 10 -transformed depth and DBD, because residuals from preliminary models fitted with log 10 -transformed variables were more closely aligned with a Gaussian (normal) distribution in residual diagnostic plots.Our preliminary analyses found alternative variable transformations to be unsuitable.Our log 10 transformations resulted in approximately linear, homoscedastic relationships between our continuous predictors, including von Post score, and log 10 (K h ).Residuals for our final LMM, which includes log 10 -transformed depth as a predictor (see below for details), indicated less deviation from a Gaussian distribution (Shapiro-Wilk; p = 0.550; W = 0.987) than an equivalent model that includes depth as an untransformed predictor (Shapiro-Wilk; p = 0.108; W = 0.975).
We refined the specification of our LMM using a stepwise, forwards and backwards model building procedure, which balanced predictive performance with model parsimony, in R v.4.1.3(R Core Team, 2022) using the lmer function from the lme4 package (v.1.1-29)(Bates et al., 2014).We began by including all candidate fixed-effect predictors (log 10 (depth), log 10 (DBD), von Post score, and degradation stage), and core identity as a random effect.We then tested whether the addition of any two-way, fixed-effect interactions significantly improved model performance, before testing the removal of the remaining fixed effects.We considered interactions between palsa degradation stage (desiccating or collapsed) and each continuous, fixed-effect predictor: log 10 (depth) × degradation stage; log 10 (DBD) × degradation stage; and von Post score × degradation stage.Additionally, we tested two-way interactions between pairs of continuous, fixed-effect predictors: log 10 (depth) × log 10 (DBD); log 10 (depth) × von Post score; and log 10 (DBD) × von Post score.To evaluate the significance of each addition or removal, we followed the method of Morris et al. (2022) to calculate the difference in the corrected Akaike Information Criterion (AICc) (Hurvich & Tsai, 1989) for each model iteration.This approach considered AICc to be a Chi-squared statistic with degrees of freedom reflecting the difference in the number of predictors included in each model.We removed any predictor or interaction from our LMM where their omission led to a non-significant increase (p > 0.05) or decrease in AICc (lower AICc indicates improved model performance).Lastly, we assessed each pair of remaining predictors for any evidence of high collinearity, by calculating the Spearman's Rank correlation coefficient (r s ), and the generalized variance-inflation factor (GVIF) using the vif function from the R package car (v.3.1-0)(Fox & Weisberg, 2018).We deemed any predictor that exhibited a GVIF greater than √5 to be too highly collinear with at least one other predictor.In such cases, we excluded one of the affected predictors and recalculated GVIF scores.For our final LMM, we calculated standardised coefficients for each predictor by re-running our LMM with standardised predictors and response (Schielzeth, 2010), which we calculated using the std function in the sjmisc package (v.2.8.9) (Lüdecke, 2018).

Applicability of a Model Trained on Boreal and Temperate Peats
We evaluated the ability of an existing pedotransfer function for peat K sat by Morris et al. (2022), trained primarily on boreal and temperate peats, to predict our measurements of log 10 (K h ) at Rensjön palsa mire.We tested the best-performing model from Morris et al. (2022), which they termed Model 1, and which requires the continuous predictors log 10 (depth), log 10 (DBD), and von Post score (all measured in this study), and the categorical predictors trophic type (e.g., bog, fen), microform (e.g., hummock, hollow), and the Kerner Oceanity Index (KOI) (Morris et al., 2022).We initially classified the trophic type and microform of our samples as the unspecified category given by Morris et al. (2022), because palsa mires were not represented in the training data used to fit their model.However, given that our site exhibits ombrotrophic characteristics, including low peat pH (see Section 2.3, above) and a prevalence of lichens, shrubs, and mosses at the coring locations, we also tested the use of the raised bog classification for trophic type.Following Morris et al. (2022), we calculated the KOI of Rensjön palsa mire from gridded mean monthly temperature data from the CRU TS 4.04 climatology (spatial resolution 0.5° latitude × 0.5° longitude) (Harris et al., 2020), averaged across the period 1961-1990.This classified the local climate as oceanic (KOI = 10.8),albeit close to the KOI threshold for continental and subcontinental climates.

Evaluation of Model Performance
We evaluated the predictive performance of each model using the standard coefficient of determination (r 2 ), which we calculated from the linear relationship between measured and modeled values of log 10 (K h ).We also adjusted this value (   2  ) to account for the number of predictors in each model, using the Wherry formula-1 in Yin and Fan (2001).We assessed bias in model predictions using Lin's concordance correlation coefficient (ρ c ), which has a range of −1 to 1, and which quantifies the agreement of the linear relationship between measured and modeled values with a line of perfect concordance.Higher values of ρ c indicate lower bias (i.e., closer agreement), while 0 indicates no agreement and −1 indicates perfect disagreement (Lin, 1989).We calculated ρ c using the CCC function in the DescTools package (v0.99.47) (Signorell et al., 2022).
Our final LMM for Rensjön palsa mire indicates that log 10 (K h ) declines significantly with increasing log 10 (depth) (p < 0.001) and von Post score (p < 0.001) (Table 1).The removal of either log 10 (depth) or von Post score leads to a significant increase in AICc (both p < 0.001), so both predictors were retained.However, retaining both log 10 (DBD) and von Post score in our LMM resulted in GVIF scores greater than 3.4 for both predictors, indicating unacceptable levels of variance inflation in their parameter estimates.The predictors log 10 (DBD) and von Post score were strongly, positively correlated (r s = 0.818) (Figure 3f), while lower levels of collinearity were evident between log 10 (depth) and von Post score (r s = 0.387), and between log 10 (depth) and log 10 (DBD) (r s = 0.254) (Figures 3d and 3e).The omission of either von Post score or log 10 (DBD) from our LMM did not lead to a significant increase in AICc, but the omission of both variables did.For this reason, it was necessary to omit one of log 10 (DBD) or von Post score from our final LMM.Residuals for a LMM fitted with von Post score were marginally more closely aligned with a Gaussian (normal) distribution (Shapiro-Wilk; p = 0.550; W = 0.987) than an equivalent LMM fitted with log 10 (DBD) (Shapiro-Wilk; p = 0.467; W = 0.985), so we retained von Post score and omitted log 10 (DBD) for our final model.However, predictive performance is similar for these two LMMs, and either predictor could be used effectively for modeling peat K sat in degrading palsas.During the backwards removal of fixed effects, we found that the stage of palsa degradation (desiccating or collapsed) did not have any significant, independent effect on log 10 (K h ) and its removal caused AICc to decrease (see Section 3.3, below).
Our final LMM has normally distributed, unstructured residuals with low GVIF scores for log 10 (depth) (1.4) and von Post score (1.4), indicating that the issue of multicollinearity has been resolved.The standardised parameter coefficient for log 10 (depth) (−0.587) suggests that it exerts a stronger influence on log 10 (K h ) than von Post score does (−0.343).Model performance metrics indicate that our LMM explains a substantial proportion of the variance in measured log 10 (K h ) (r 2 = 0.605;   2  = 0.590) with low predictive bias (ρ c = 0.756) (Figure 4a).Our LMM generally underestimates near-surface, high values of log 10 (K h ) and overestimates deeper, low values (Figure 4a).
While testing the addition of two-way interactions in our model building process, we discounted the interaction between log 10 (DBD) and von Post score, because the inclusion of an interaction between these two highly correlated predictors (Figure 3f) risked introducing severe multicollinearity and parameter variance inflation to our LMM.No addition of any interactions improved model performance, so none were included in our final LMM.As expected, our data set exhibits some evidence of within-core dependency, indicated by a significant increase in model AICc (p < 0.05) if the random subject for core identity is removed.Our final LMM therefore includes core identity as the subject of a random intercept term.

Performance of an Existing K sat Model
The peat K sat model by Morris et al. (2022) accounted for a lower proportion of the variance in log 10 (K h ) at Rensjön palsa mire (r 2 = 0.528;   2  = 0.429) than our site-specific LMM (see Section 3.1, above).Setting the trophic type of the site to be unspecified resulted in model predictions that had substantial bias (ρ c = 0.604) and consistently underestimated measured log 10 (K h ), particularly for shallow peats with high log 10 (K h ) (Figure 4b).Changing the trophic type to raised bog removed some of this bias, increasing ρ c to 0.693 (Figure 4c).This level of bias is similar to that of our LMM (Figure 4a).For high measured values of log 10 (K h ), the model by Morris et al. (2022) more commonly overpredicted values for collapsed palsas and underpredicted values for desiccating palsas, although no discernible trends were evident in predictions of lower log 10 (K h ) values (Figures 4b and 4c).

The Effect of Palsa Desiccation and Collapse on Peat K h
The stage of palsa degradation did not exert a significant, independent control on measured log 10 (K h ) at Rensjön palsa mire.Measurements of K h were similar between desiccating palsas (median = 1.79 × 10 −5 m s −1 ) and collapsed palsas (median = 1.17 × 10 −5 m s −1 ) (Figure 3g).Contrary to our expectations, DBD was generally higher in peats from desiccating than collapsed palsas (medians of 0.17 and 0.11 g cm −3 , respectively) (Figure 3h), and the most highly humified peats were also sampled from desiccating palsas (Figure 3i).However, when all fixed effects were included during the stepwise model building process, there was almost no evidence of the categorical variable for degradation stage exerting an effect on log 10 (K h ) (p = 0.983) and the removal of this predictor caused model AICc to decline (an improvement in model performance).Furthermore, no interactions involving degradation stage (log 10 (depth) × degradation stage; log 10 (DBD) × degradation stage; or von Post score × degradation stage) gave any significant improvement in the performance of our LMM.

Horizontal Permeability of Degrading Palsas
Peat K h in degrading palsas near Rensjön evidenced similar negative relationships with depth, humification and DBD to previous studies of northern peatlands, with our highest K h values occurring in lightly decomposed, near-surface peats.While between-site comparisons of peat K h are somewhat complicated by different measurement techniques (Morris et al., 2022), our measured values for the upper 0.1 m of peat were at least one order of magnitude lower than those reported for peat plateaus at Scotty Creek, Northwest Territories, Canada (Gharedaghloo et al., 2018;Quinton et al., 2008), perhaps due to the comparatively lower DBD of unburnt peats at that site (Ackley et al., 2021).The higher DBD values for Rensjön palsa mire are comparable to previous DBD measurements from Fennoscandian palsas (Sim et al., 2021;Wetzel et al., 2003) and our maximum K h value (3.61 × 10 −4 m s −1 ) is similar to the highest K v value recorded for palsas near Kautokeino, Norway (2.94 × 10 −4 m s −1 ) (Wetzel et al., 2003).The minimum value of peat K h measured in our study (4.85 × 10 −7 m s −1 ) is similar to some previous measurements from anoxic peats on the Alaskan North Slope (O'Connor et al., 2020).
Log 10 (depth) exerted the strongest control on log 10 (K h ) in our LMM, as indicated by its high standardised coefficient, while the von Post measure of peat humification also exerted a highly significant control (Table 1).Peat depth accounts for some of the variation in both peat compaction and decay because the structural matrix of deeper, older peats is weakened by prolonged decomposition, which increases their susceptibility to compression under the weight of fresher, overlying peats (Branham & Strack, 2014;Ohlson, 1998).Alternatively, declining K h with depth may relate to differences in peat hydraulic properties between the oxic and anoxic components of the peat profile (O'Connor et al., 2020), increased pore tortuosity, and loss of pore connectivity (Gharedaghloo et al., 2018).Peat desiccation can also drive a vertical displacement of the peat profile through peat shrinkage and compression (Ackley et al., 2021;Hobbs, 1986) or by rendering raised palsa surfaces susceptible to aeolian erosion (Seppälä, 2003).Log 10 (depth) was positively correlated with von Post score (Figure 3e), although low GVIF scores suggest that any collinearity is not great enough to destabilize parameter estimates and that each predictor exerts a strong, independent effect on log 10 (K h ) in our LMM.
The high collinearity of log 10 (DBD) and von Post score (Figure 3f) was greater than reported for previous studies (Morris et al., 2019(Morris et al., , 2022) ) and prevented us from including both predictors in our final LMM due to unacceptable inflation of parameter variance.However, even after omitting log 10 (DBD), the strong explanatory power of our LMM (r 2 = 0.605) suggests that only one of these predictors is required to explain this portion of variance for log 10 (K h ) in degrading palsas.Indeed, our preliminary analyses identified a similar level of predictive performance from an equivalent LMM fitted with log 10 (DBD) instead of von Post score, and DBD has independently explained large portions of peat K sat variance in other permafrost sites (O'Connor et al., 2020).The retention of log 10 (depth) and von Post score, and removal of log 10 (DBD), means that our LMM, like Model 2 by Morris et al. (2022), can be used to make rapid field estimations of peat K sat at our site without the need for any laboratory work.
Estimated slope coefficients for log 10 (depth) and von Post score in our LMM (Table 1) closely align with those from the multi-site model by Morris et al. (2022), providing some support for the assertion by Quinton et al. (2008) that peats with a similar degree of decomposition and compaction evidence similar K sat irrespective of their geographic location.The generally high DBD of palsas at our study site means that our near-surface K h determinations are toward the lower end of values recorded for the upper 0.35 m of boreal and temperate peats (Morris et al., 2022).For any given site, a model fitted specifically to that site is always likely to demonstrate greater predictive performance than a multi-site model, trained on data from other sites and peatland types.It is not surprising, therefore, that the multi-site pedotransfer function by Morris et al. (2022) made slightly less skillful predictions of log 10 (K h ) at Rensjön palsa mire than our site-specific LMM (Figure 4), and had less explanatory power for this site (r 2 = 0.528) than for the boreal and temperate peatlands on which it was trained (cross-validated r 2 = 0.747).This reduced performance may, to some degree, relate to unrepresented processes specific to permafrost peatlands, such as their different floristic composition or freeze-thaw dynamics (discussed below).However, the multi-site model still explained approximately half of the total variance in unseen values of log 10 (K h ) in degrading palsas and exhibited similar levels of bias to our site-specific model, which may be considered adequate depending on the application.

Permafrost Thaw and Peat Hydraulic Properties
Contrary to our expectations, we did not find any significant difference in peat K h between desiccating and collapsed palsas (Figure 3g), nor any significant interaction involving the stage of palsa degradation.For the uppermost peat layers (<0.1 m depth), the range of peat K h was almost identical in desiccating and collapsed palsas (Figure 3a), despite samples from desiccating palsas often having higher DBD and slightly greater degrees of peat humification (Figures 3d,3e,3h,and 3i).Alternative, unmeasured factors, such as macroporosity, may have reduced any differences in the K h of these shallow peats (Holden et al., 2001), particularly because some deeper, lightly-humified samples from collapsed palsas did have higher K h than comparably-deep, desiccating palsa peats (Figures 3a and 3e).However, this comparison may be affected by the greater number of measurements from depths ≥0.2 m in collapsed palsas (n = 24) than desiccating palsas (n = 9).Raised, desiccating palsa surfaces were more commonly inhabited by ligneous shrubs than collapsed areas, where saturated, anoxic conditions restrict root growth (Limpens et al., 2021;Olefeldt et al., 2021).Ligneous roots can increase K sat , even in decayed and compacted peats, through the creation of macropores and relict root channels (Holden, 2009;McCarter et al., 2020;Surridge et al., 2005).Shallow, impenetrable frost tables in desiccating palsas force rooting zones to be near-surface and to extend laterally, meaning shrub roots may be less influential for the K h of deeper peats or K v .Alternatively, freeze-thaw cycles can increase the macroporosity of low K sat peats through the formation of ice lenses in micropores, which produce cracks upon expansion that subsequently increase peat K sat (Liu et al., 2022).Conversely, for peats with initially high K sat , freeze-thaw processes generally reduce peat K sat , because existing macropores can become blocked by the breakdown of large aggregates (Liu et al., 2022).
Given the similarity of peat K h between desiccating and collapsed areas of the palsa complex, thaw-driven changes to the hydraulic functioning of the site will more likely be driven by increases to active layer thickness and reduced topographical gradients than any direct alterations to peat hydraulic properties.Currently, lateral runoff from desiccating palsas into neighboring, collapsed areas and inundated fens is driven by the high K h of the uppermost, seasonally-thawed peats, steep hydraulic gradients, and the prevalence of near-surface, impermeable permafrost.As palsas thaw and active layers thicken, water tables may be initially lowered into deeper peats with lower K h , which could increase water retention as hydraulic gradients simultaneously decline (Quinton & Baltzer, 2013).A shift to persistently saturated, anaerobic conditions may explain why many collapsed palsa peats unexpectedly had lower humification and DBD, although this could also indicate that compaction following subsidence occurs deeper in the peat profile, below the active layer.Longer-term permafrost degradation may increase the hydrological connectivity of the site, because water tables in collapsed areas persist closer to the surface, where peat K h is highest (Figure 3a), and thaw removes impermeable, frozen barriers to groundwater flows (Connon et al., 2014;Haynes et al., 2018;Quinton & Baltzer, 2013).

Opportunities for Future Research
Our modeling demonstrates that simple, affordable measurements of peat properties can be used to predict K h in peatlands beyond the boreal and temperate regions to which many existing, pedotransfer functions are currently fitted (e.g., Liu & Lennartz, 2019;Morris et al., 2022).Generalizable models of peat K sat have the potential to provide skillful estimations for the hydrological functioning of remote Arctic peatlands without the need for extensive field or laboratory measurements, although we recommend that K sat measurements from permafrost peatlands should be better represented in the training data of multi-site models.During our analysis, we found that the Rensjön palsa mire was best classified as a raised bog (rather than an "unspecified" trophic type) in the model by Morris et al. (2022) (Figures 4b and 4c).Developing additional categorical descriptors to account for factors specific to permafrost peatlands, for example, the presence of permafrost landforms (e.g., palsas, peat plateaus, ice-wedge polygons), annual freeze-thaw cycles, and differences in surface vegetation, may further improve model predictions of peat K sat for Arctic regions, once suitable training data are available.
Further research is also required to explore spatial variability in peat hydraulic properties of permafrost peatlands in Fennoscandia and Siberia, because measurements there remain sparse, and other permafrost peatland types remain understudied, such as polygon mires (Helbig et al., 2013).We do not know when the palsas at Rensjön palsa mire began to thaw, but spatial chronosequences of palsa mires with different ages of permafrost degradation could be analyzed to assess whether the post-thaw stability of peat hydraulic properties changes with time.
Limited aerial imagery of our study site indicates that thermokarst drainage and collapse scar formation occurred at some point between 1960 and 2013 (Sim et al., 2021).For more intensively-studied sites, monthly observations of snow-cover, active-layer thickness, and disturbance (e.g., burning) may account for some unexplained variance in peat K h (Ackley et al., 2021).Lastly, future studies could attempt to measure peat K sat in newly terrestrializing areas of inundated fens, which represent the final stage of palsa degradation after collapse (Swindles et al., 2015) and are important for the longer-term hydrological functioning of the site following permafrost thaw.

Conclusions
We found that much of the variance in log 10 (K h ) of degrading Arctic palsas can be explained using simple measurements of peat properties, in a similar manner to previous models of boreal and temperate peatlands.Our LMM showed that log 10 (depth) exerted the strongest, independent control on log 10 (K h ), while the highly influential effects of von Post score and DBD were found to be highly collinear, meaning the latter was omitted from our final model.An existing pedotransfer function of peat K sat , fitted primarily to boreal and temperate peatlands, made slightly less skillful predictions of log 10 (K h ) in degrading palsas than our site-specific LMM, but demonstrated similarly low bias.Unexpectedly, near-surface peat K h did not significantly differ between desiccating and collapsed palsas, indicating that any post-thaw changes to their hydrological functioning will be more likely driven by reduced hydraulic gradients than altered peat properties.Permafrost peatlands should be more fully represented in multi-site models of peat K sat before such models are widely extended to Arctic regions.

Figure 1 .
Figure 1.The Rensjön palsa mire study site: (a) the location of the site within Sweden (inset) and the Kiruna Municipality; and aerial imagery (August 2022) of the site overlain with (b) the 20 coring locations, colored by the stage of palsa degradation; and (c) active-layer thickness across the site, measured using a ∼2 m-long steel probe.Land cover data sourced from www.lantmateriet.se/.Aerial imagery sourced from ESRI (2023).

Figure 2 .
Figure2.The experimental setup for measuring horizontal saturated hydraulic conductivity of peat.The head difference across the sample (ΔH) is calculated from the end of the Mariotte bubble tube and the end of the outflow tube positioned in the measuring cylinder.We repositioned the Mariotte bubble tube to adjust ΔH for peats with different degrees of humification (e.g., a lower ΔH was required for highly permeable peats with open, interconnected pores).

Figure 3 .
Figure 3. Peat hydraulic properties in degrading palsas.Top row: the response of K h to (a) depth, (b) dry bulk density, and (c) von Post score.Middle row: bivariate correlation plots showing collinearity between continuous predictors (d-f); r s is the Spearman's Rank correlation coefficient.Data points are colored by the stage of palsa degradation.Bottom row: variation in (g) K h , (h) dry bulk density, and (i) von Post score between stages of palsa degradation.K h , depth and dry bulk density are presented on a log 10 scale.In panels (g-i): centerlines indicate median values; boxes indicate the interquartile range (IQR), and whiskers extend up to 1.5 times the IQR beyond the upper and lower quartiles.

Figure 4 .
Figure 4. Performance plots showing measured values against modeled values for: (a) our linear mixed-effects model (LMM); and for an existing, multi-site model by Morris et al. (2022) with trophic type specified as (b) unspecified and (c) raised bog.Dashed line indicates perfect concordance between measured and modeled values; solid line indicates a linear regression between measured and modeled values; r 2 is the coefficient of determination; and ρ c is Lin's concordance correlation coefficient, used to estimate model bias.

Table 1
Summary of Fixed Effects for Our Linear Mixed-Effects Model to Characterize log 10 (K h ) in Rensjön Palsa Mire