Saturated Hydraulic Conductivity in Northern Peats Inferred From Other Measurements

In northern peatlands, near‐saturated surface conditions promote valuable ecosystem services such as carbon storage and drinking water provision. Peat saturated hydraulic conductivity (Ksat) plays an important role in maintaining wet surface conditions by moderating drainage and evapotranspiration. Peat Ksat can exhibit intense spatial variability in three dimensions and can change rapidly in response to disturbance. The development of skillful predictive equations for peat Ksat and other hydraulic properties, akin to mineral soil pedotransfer functions, remains a subject of ongoing research. We report a meta‐analysis of 2,507 northern peat samples, from which we developed linear models that predict peat Ksat from other variables, including depth, dry bulk density, von Post score (degree of humification), and categorical information such as surface microform type and peatland trophic type (e.g., bog and fen). Peat Ksat decreases strongly with increasing depth, dry bulk density, and humification; and increases along the trophic gradient from bog to fen peat. Dry bulk density and humification are particularly important predictors and increase model skill greatly; our best model, which includes these variables, has a cross‐validated r2 of 0.75 and little bias. A second model that includes humification but omits dry bulk density, intended for rapid field estimations of Ksat, also performs well (cross‐validated r2 = 0.64). Two additional models that omit several predictors perform less well (cross‐validated r2 ∼ 0.5), and exhibit greater bias, but allow Ksat to be estimated from less comprehensive data. Our models allow improved estimation of peat Ksat from simpler, cheaper measurements.

One of the most important peat hydraulic properties, and one of the most commonly measured, is saturated hydraulic conductivity, K sat (dimensions of L T −1 ). Low K sat is thought to play a central role in maintaining shallow water tables in peatlands by inhibiting lateral drainage (Ingram, 1982). As in any porous medium, peat K sat is also relevant to unsaturated hydrological processes because it defines the upper limit for unsaturated hydraulic conductivity (e.g., van Genuchten, 1980). Peatland hydrological processes are complex, and water retention is governed by a number of interacting feedbacks between hydrophysical, ecological, and biogeochemical processes across a range of timescales (Siegel & Glaser, 1987). Many of these feedbacks are mediated through changes in peat hydrophysical properties, including K sat .
Peat K sat generally declines with increasing depth (at least in the upper few decimeters), dry bulk density, and measures of the degree of peat humification such as the von Post hand-squeeze test (Päivänen, 1973). Fresh plant litter is commonly highly permeable but as it decomposes to form peat, pore spaces close and collapse and become increasingly disconnected. Strong vertical gradients are commonly observed in peat K sat , which can decline by several orders of magnitude across the uppermost few decimeters of a peat profile (e.g., Bourgault et al., 2018;Clymo, 1992;Holden & Burt, 2003;McCarter et al., 2020;Menberu et al., 2021;Ronkanen & Kløve, 2005). Peat hydraulic properties can be affected by both natural and anthropogenic disturbances. Direct human alteration of peatlands, particularly drainage for peat extraction, or conversion to agriculture or forestry, causes rapid aeration and decomposition and increased compression, leading in turn to reduced K sat (Prevost et al., 1997;Price, 2003). On the other hand, peatlands in which water tables have been artificially raised through deliberate or inadvertent blockage of downflow drainage structures can exhibit high K sat that reflects open pore structures in partially buoyant peat . Wildfire has been shown to reduce near-surface K sat (Holden et al., 2014), possibly due to large pores becoming clogged by ash (Mallik et al., 1984). Burning can also remove highly permeable, near-surface peat through combustion (Sherwood et al., 2013), particularly in sites that have undergone compound disturbances, in which wildfire follows another disturbance such as drainage, and which can enhance depth of burn (McCarter et al., 2021;Sherwood et al., 2013). Some studies have found peat K sat to vary depending on surface microform type (e.g., hummocks, lawns, and hollows) at spatial scales of a few meters (e.g., Baird et al., 2016;Branham & Strack, 2014;Morris et al., 2015), a phenomenon that likely reflects the influence of peat floristic composition upon pore structure (McCarter et al., 2021). It is unclear whether K sat also differs between peat trophic types, such as the distinction between wet, nutrient-rich, biodiverse fens and less-wet, nutrient-poor bogs. It has been postulated that climate exerts an independent control on K sat , possibly through its influence on peat floristic composition (Branham & Strack, 2014) or the enlargement of pores in continental sites through physical mechanisms such as desiccation cracking in hot, dry summers and deep freezing in winter.
Pedotransfer functions are statistical models used to infer hydraulic properties of sediments and soils from multiple independent variables, particularly grain-size distribution, organic matter content, and dry bulk density. The use of such functions has proliferated in recent decades, in large part because they allow representations of porous media in hydrological models to be parameterized using simple, cheap measurements, without the need for direct measurements of hydraulic properties (e.g., Wösten et al., 2001). Pedotransfer functions developed for mineral soils are inapplicable to peat because grain-size distribution can rarely be defined, and organic matter content is usually close to 100%. However, equivalent schemes are beginning to be developed for peats, often based on the reanalysis of large secondary databases (e.g., Wang et al., 2021). These models are commonly simple in structure and employ few predictor variables, such as dry bulk density and/or macroporosity. Menberu et al. (2021) demonstrated that a simple model exhibited high explanatory power (r 2 = 0.74) for specific yield. However, statistical models developed from studies of large databases are considerably less skillful for the estimation of peat K sat , with r 2 typically in the range of 0.3-0.5 for most peat types (e.g., Liu et al., 2020;Wang et al., 2021). Studies of K sat from individual peatlands indicate that much greater explanatory power can be achieved by using a combination of several independent variables, including continuous covariates such as depth, dry bulk density, and measures of peat humification (e.g., Päivänen, 1973); and categorical factors such as distinctions between hummock and hollow 3 of 18 microforms (Branham & Strack, 2014;Morris et al., 2015Morris et al., , 2019. A wealth of published peat K sat measurements exists and is accompanied by varying amounts of other data that could serve as predictor variables, presenting an opportunity to develop more skillful predictive equations that can be generalized beyond individual sites.

Aim and Objectives
We sought to use a meta-analysis of existing K sat measurements from northern peats to develop skillful predictive equations. We addressed the following specific objectives: 1. Compile a large database of peat K sat measurements and accompanying measurements that can be used as independent variables in a statistical analysis. 2. Develop models that can be applied at a range of spatial scales to predict unknown values of peat K sat , from simpler, cheaper measurements. 3. Evaluate the ability of the predictive models to simulate data points not used in their construction, as an estimate of their reliability in predicting unknown values of K sat .

Methodological Overview
We compiled a database of both published and unpublished K sat measurements from northern peatlands (i.e., those north of 45°N). For each K sat measurement, we also collected as many accompanying data as were available, both continuous measurements and categorical observations, which might be used as independent predictors. We fitted four linear models to describe log 10 -transformed K sat based on a selection of independent variables chosen objectively and judiciously in an effort to maximize model performance without overfitting. The four models each exclude certain combinations of candidate independent variables to mimic situations in which some field observations are unavailable or where they are inapplicable at the spatial scale under consideration. We assessed the performance of each model using a cross-validated variant of r 2 designed to test the ability of a model to predict unseen data points, and which is more conservative than standard r 2 .

Dependent Variable: Saturated Hydraulic Conductivity
We constructed a database of K sat measurements from 2,507 distinct peat samples. Each record in the database represents a single peat sample for which K sat has been determined, and also includes data fields to represent additional variables (described below) that accompany that K sat measurement. Our data set represents 47 distinct study sites across boreal and temperate-humid zones of Europe and North America ( Figure 1). We excluded tropical, subtropical, and southern-hemisphere peatlands, where measurements of hydraulic properties are comparatively sparse. We also omitted Russian peatlands, where data are less readily available. Our data set also omits permafrost peatlands of the Arctic and Subarctic, in which hydrological processes are complicated by the presence of permanent ice lenses and a seasonally variable frost table (Gamayunov et al., 1990;Nagare et al., 2012). Some of the 47 study locations are pristine peatlands, while others have undergone some level of disturbance. We omitted data from locations where the depth distributions of K sat and other peat properties were likely to have been severely disrupted by direct human alterations, such as the construction of reclaimed peatlands in former resource-extraction sites, and areas of peatlands that had undergone historical peat extraction.
Of the 2,507 K sat measurements, 2,390 are from our own published or unpublished data and were supplied directly by the authors based on original measurements. We digitized an additional 95 K sat measurements from published graphs, using the WebPlotDigitizer online tool (Rohatgi, 2020) and transcribed a further 22 data points from published tables of raw data. We converted the units of all K sat measurements to m s −1 .
The K sat determinations were made using a variety of measurement techniques, including in situ determinations from piezometer slug tests and tension-disc infiltrometers; and a variety of laboratory methods, including the modified cube method (Beckwith et al., 2003) and permeameters. In situ measurements from tension-disc infiltrometers and piezometer slug tests give an estimate of K sat that represents 3-dimensional flow into or out of the instrument in response to a head gradient; whereas laboratory methods typically provide directional measurements of horizontal (K h ) or vertical (K v ) hydraulic conductivity. For some samples, only either K h or K v were 4 of 18 available; in these instances, we took this single, directional measurement to be representative of K sat . For samples which had both K h measured in a single horizontal axis and K v , we took the geometric mean of K h and K v to be representative of overall K sat . In a small number of samples, for which K h had been measured in two horizontal axes in addition to a single measurement of K v (i.e., Cunliffe et al., 2013), we took the geometric mean of the three values. Replicate determinations of K sat (or K h and/or K v ) were available for some samples; in these cases, we took the median of the replicates to be representative of that sample and discarded the individual replicate measurements. Where replicate K h and K v determinations were available, we first calculated the between-replicate medians of K h and K v separately, before calculating the geometric mean between the two directions.

Continuous Covariates
Where available, we recorded three continuous covariates for each K sat record: sample midpoint depth relative to the ground surface, von Post score as an estimate of degree of humification, and peat dry bulk density. For in situ piezometer slug-test estimates of K sat , we took the midpoint depth to be the midpoint depth of the piezometer intake below the ground surface. Von Post scores occupy an ordinal scale between H1 (fresh, undecomposed plant litter) and H10 (black, amorphous, highly decomposed organic matter) and are therefore strictly not continuous. However, following the rationale of Morris et al. (2019), we treated von Post as a continuous predictor for the purposes of our statistical modeling. We converted all depths into units of m and all dry bulk density measurements into g cm −3 . Midpoint depth was available for all 2,507 samples in our database. However, the two other continuous covariates are unavailable for a large proportion of records. Missing values for dry bulk density or von Post score prevent those records from being used to fit models that include those variables as predictors, placing important limits on sample sizes for models that employ one or both of these covariates (see Section 2.7, below).

Categorical Factors and Harmonization
For all records in the database, we used original field notes and/or published descriptions to assign values for five categorical variables that describe the immediate vicinity of the sampling location or the broader characteristics of the site. These factors describe the surface microform type (e.g., hummock and hollow) from which the K sat determination was made; the trophic type of the wider study site (e.g., blanket peat, raised bog, and poor fen); the nature of any disturbance such as human alterations to water tables or recent burning; whether the sampling location was forested or open; and the measurement method used to determine K sat . Where insufficient information was available to identify a sample's status confidently with respect to one or more of these categorical factors, we assigned a value of "unspecified" rather than a blank. Doing so meant that, unlike for missing values of the continuous covariates described above, missing data for one or more categorical factors do not prevent records from being used to fit our models.
The categorical data were originally recorded by a variety of operators from different research groups and projects, according to different protocols and with different levels of detail. We harmonized the levels of each category according to a common set of descriptors to allow direct comparison between records. In many cases this harmonization involved removing detail from the original data and grouping records into simpler classifications. For example, where the original records indicated that samples had been taken from a "low hummock" microform type, we grouped these into a single microform category along with "hummocks." Similarly, "low lawns" and "wet lawns" were combined into a single category along with "lawns." A large number of records had no microform specified. In most cases this reflects an absence of clear hummock-hollow topography in the field. In total, the harmonized microform factor has four levels: hummock, lawn, hollow, and unspecified. It is unclear whether surface microforms persist for extended periods or whether they change during the course of peatland development (cf. Koutaniemi, 1999;Swanson & Grigal, 1988). Surface microform type may or may not, therefore, be a good descriptor of the floristic composition of subsurface peat. Ideally, we would have used a descriptor of peat type that relates directly to each K sat measurement, rather than surface microform, but information about peat type was unavailable for many samples. However, the large majority of our samples are from shallow peat layers where surface microform type and hydraulic properties seem likely to be most closely linked to subsurface peat type (cf. Baird et al., 2016): 87% of our data are from within 1 m of the surface, while 60% are no deeper than 50 cm. We proceeded with the surface microform factor and allowed our statistical modeling (described below) to indicate whether or not it improves estimates of K sat .
We harmonized descriptions of peatland trophic type into six levels: blanket peat, raised bog, poor fen, moderate fen, rich fen, and unspecified. We used the scheme described by Rydin and Jeglum (2013) to distinguish between poor, moderate, and rich fens. The only records for which trophic type is unspecified are the ones we digitized from the literature and where the original site descriptions do not specify trophic type, vegetation community, nor sufficient precision in spatial coordinates for us to establish the trophic type of the sampling location ourselves.
The original site descriptions that we collated included a wide variety of (sometimes subtly) different descriptions of human and natural disturbance. We harmonized these descriptions into seven categories, each of which might reasonably be expected to exert a distinct effect upon peat K sat ; and an eighth, unspecified category. Where no obvious human or natural disturbances seemed likely to have affected K sat , we described the record as representing 6 of 18 pristine peat. We considered sites to be burnt where obvious surface burn scars could be identified or where historical records indicate substantial fires in recent decades. A large number of samples come from undrained, largely intact locations in sites which had experienced drainage or other human-induced drying at the margin of the peat or at other distant locations, which may have affected the characteristics of the peat at the sampling location through long-range water table drawdown and wasting; we refer to these samples as having experienced marginal drainage. Where water tables had been lowered in the immediate vicinity of the sampling location, we refer to these samples as representing peat that has been dried; similarly, a small number of samples come from locations where water tables had been artificially raised, and we refer to these samples as wetted. Some Finnish peatlands represented in our database have received treated municipal wastewater or runoff from land-use practices such as forestry and agriculture, as part of efforts to retain nutrients and suspended solids; we refer to these sites as having been subject to water treatment. Finally, some samples represent peat that has experienced more than one of these disturbances, which we refer to as compound disturbance.
In all cases we recorded the measurement method used to estimate K sat . We grouped measurement methods into six harmonized categories: piezometer; modified cube method; tension-disc infiltrometer; multistep outflow; permeameter, including a variety of designs and laboratory setups; and the KSAT instrument (METER Group, Inc., Pullman, Washington, USA). We did not distinguish between rising and falling head tests in piezometers or permeameters because in some cases this information was unavailable.
Finally, we recorded whether the peatland in the immediate vicinity of the sampling point was forested or open.
In a small number of cases where the presence or absence of tree cover could not be determined unambiguously from reports of secondary data sets, we recorded a value of unspecified.

Climate Data: The Kerner Oceanity Index
For each of the 47 study locations represented in the database, we extracted gridded mean monthly temperature data from the CRU TS 4.04 database (Harris et al., 2020), averaged for the period 1961-1990. These gridded climate data have a spatial resolution of 0.5° latitude × 0.5° longitude. From the monthly temperature data, we derived a bioclimate index that can be used to represent continentality at each site, the Kerner Oceanity Index (KOI): where T oct and T apr are monthly average temperatures (°C) during October and April, respectively; and T max and T min are the average temperatures (°C) of the hottest and coldest months, respectively. Originally proposed by Kerner (1905), the index is predicated upon the greater thermal inertia of oceanic locations leading to cool springs and warm autumns; and the opposite in dry, continental locations where thermal inertia is lower. In mid to high latitudes of the northern hemisphere, high values of KOI indicate oceanic climates, while low values indicate continental climates. Although KOI varies on a continuous scale between −100 and 100, all K sat measurements from a single site will possess identical values of KOI, while almost every site possesses a unique value of KOI. This means that KOI is quasi-continuous in our analysis and is aliased with the identity of each study location. To avoid problems relating to unequal variance, we grouped KOI values into ordinal bands, thereby converting this variable into a categorical predictor. Following Stonevicius et al. (2018), we categorized KOI into three ordinal levels: continental and subcontinental (KOI < 10), oceanic (10 ≤ KOI < 20), and hyperoceanic (KOI ≥ 20). Our climatic characterization of each study site represents modern conditions, yet the deepest peat will have formed several millennia ago, during which time climate will inevitably have changed. However, Holocene climate change in northern latitudes is thought to have been modest (e.g., Marcott et al., 2013) compared to the much more dramatic changes during the preceding deglaciation of the northern hemisphere (e.g., Shakun et al., 2012). As such, we assume that our modern climate data are sufficient in most cases to provide an approximation of sites that have experienced continental, oceanic, and hyperoceanic climates.

Variable Transformations
Preliminary analysis of the database indicated that the K sat data vary across several orders of magnitude, are highly positively skewed, and exhibit highly nonlinear, heteroscedastic relationships to the three continuous 7 of 18 predictors. We log 10 -transformed K sat as well as midpoint depth and dry bulk density, which remedied these problems (Figures 2a and 2b). We did not transform von Post scores, which already exhibit a linear relationship to log 10 (K sat ) (Figure 2c). A debate exists as to whether it is preferable to fit nonlinear models to untransformed data in such situations (cf. Xiao et al., 2011;Tian et al., 2013). However, preliminary experimentation with nonlinear models indicated that they were numerically intractable for most combinations of our predictor variables, mainly because parameter estimates failed to converge. We discarded the possibility of using nonlinear models and proceeded to use linear models fitted to log 10 -transformed data.

Variable Selection and Model Calibration
We fitted four different linear models to predict log 10 (K sat ). We began by specifying a candidate pool of independent variables for each model that were to be considered for inclusion as fixed effects. For each model we omitted some independent variables from the candidate pool to reflect that model's intended application. Seeking an appropriate balance between predictive skill and model parsimony, we then used an iterative process, described below, to identify which variables from the candidate pool to include in the final model specification and which could be omitted.
Model 1 is intended to provide the most skillful predictions possible, for use in situations when all three continuous covariates-log 10 (depth), log 10 (dry bulk density), and von Post score-are available. All three are included in the candidate pool, but doing so comes at the cost of replication. All 2,507 records in our database possess values of midpoint depth but only 453 of those also possess both dry bulk density and von Post score, meaning that Model 1 utilizes just 18% of the records in the database (see Table 1 for further details). We made all categorical factors available to the candidate pool for Model 1, apart from the measurement method used to determine K sat , which we included only as a random effect. We intended for our models to predict unknown values of K sat , a situation in which the measurement method cannot be defined. At the same time, we sought to remove from our models any artifacts associated with measurement method, so that the fixed effects parameters in our models are free of methodological bias. Setting measurement method as a random-effect intercept means that a portion of the variance in log 10 (K sat ) is attributed to measurement method before the fixed effects of interest are estimated. Because measurement method is included only as a random effect, and therefore has no parameter estimates, we do not consider it further.
Model 2 is intended to allow rapid estimation of K sat from the simplest predictor variables that can be measured entirely in the field, without the need for laboratory work. Model 2 therefore omits log 10 (dry bulk density) from the candidate pool because it requires laboratory drying; but includes both log 10 (depth) and von Post score, and all categorical factors apart from measurement method, which is again included as a random effect. The requirement for von Post score means that Model 2 utilizes 762 records (see Table 1 for further details).
Model 3 takes a contrasting approach to that of Model 1, seeking to maximize replication by omitting log 10 (dry bulk density) and von Post score, leaving log 10 (depth) as the only continuous covariate. Doing so allows Model 3 to utilize all 2,507 records. As above, the candidate variable pool for Model 3 includes all categorical factors, apart from measurement method, which is again included as a random intercept. In Models 1, 2, and 3 we experimented with the two-way fixed-effect interactions depth × disturbance and depth × microform. These interactions are intended to capture any differences in the depth profile of K sat that occur in disturbed sites or beneath different microform types.
Model 4 is intended for prediction of K sat at larger spatial scales, such as in land surface scheme model tiles, where small-scale variables such as surface microform type or peatland trophic type cannot be defined. Despite their coarse horizontal resolution, land surface schemes commonly represent detailed depth-distributions of soil hydrophysical properties, so we included log 10 (depth) and log 10 (dry bulk density) in the candidate variable pool; we omitted von Post score, which is an operational variable not readily simulated by process-based land surface models. The requirement for dry bulk density means that Model 4 utilizes 883 records (see Table 1 for further details). We limited the categorical variables in the candidate pool to the two that apply at larger scales: KOI and the distinction between treed and open peatlands. We experimented with other categorical variables (measurement method, disturbance, trophic type, and microform) as potential random effects, again to remove any confounding influence from estimates of parameter values for fixed effects. All four models include a fixed-effect intercept term.
Having established the initial pool of candidate variables for each model, we used an iterative procedure to identify those variables which could be omitted from the final model without causing a significant decrease in predictive skill. For each candidate variable in turn, we calculated the change in the corrected Akaike Information Criterion (AICc) (Hurvich & Tsai, 1989) between a linear mixed-effects model that includes the variable under consideration and an equivalent model that omits it. We also removed variables that caused other problems such as multicollinearity. We used the lmer function from the R (R Core Team, 2014) package lme4 (Bates et al., 2015) to fit all models. We treated the change in AICc as a Chi-squared statistic, with degrees of freedom equal to the change in the number of model parameters. If the removal of a variable caused a nonsignificant increase in AICc 9 of 18 (significance threshold α = 0.05) or a decrease, then we omitted it from the final model specification (lower AICc indicates better model performance). At each iteration we also calculated the generalized variance-inflation factor (GVIF) for each variable, using the vif function from the R package car (Fox & Weisberg, 2019), corrected for the number of degrees of freedom (df) in the case of categorical variables: GVIF 1 2 (Fox & Monette, 1992 we compared parameter estimates between each iteration of model fitting, and between the four models, to check for any other indications of potential problems, such as instability of parameter values and/or directions. Given the size of our data set and its geographic extent, it might have been possible to develop separate, regional models of K sat for different locations. However, we sought to establish generalizable models of peat K sat that are independent of geographic location and which can be applied across the study domain (Europe and North America). For this reason, we not do present a regional analysis.

Model Validation and Performance Metrics
For each model we calculated the standard coefficient of determination, r 2 , by regressing predicted values of log 10 (K sat ) on the observed values. In order to assess the likely skill of our models in predicting unknown K sat values, we also calculated a cross-validated coefficient of determination, 2 cv , to assess performance when predicting data points not used in model specification. We split the records in the database randomly into five subsets of approximately equal size. We combined four of the subsets into a training set which we used to fit a model. We then used that model to predict log 10 (K sat ) for the fifth (unseen) subset only, which acts as a validation set, and calculated standard r 2 by regressing predicted on observed values. We repeated this five times, each time omitting in turn one of the five subsets from the training set and using it instead as a validation set. We took the arithmetic mean of the five r 2 values to give 2 cv . A large discrepancy between r 2 and 2 cv can indicate that parameter estimates are sensitive to a small number of data points, meaning that predictions made using the model may be unreliable (Shmueli, 2010). Final parameter estimates and p-values for all models are calculated from linear models fitted to the full data sets, not subsets.
Coefficients of determination as described above do not give a full picture of model performance. High r 2 indicates only that a high proportion of variance in the predicted values is explained by observed values; it gives no indication of whether predicted values agree with observed values. Models in which the linear relationship between predicted and observed values differs from the 1:1 line of perfect agreement are said to exhibit bias. Using the CCC function from the R package DescTools (Signorell et al., 2021), we calculated Lin's concordance correlation coefficient (ρ c ) to assess agreement between the linear relationship between predicted and observed values, and the 1:1 line, as a measure of model bias. Although ρ c can also be used as a measure of model fit, we use it primarily to indicate bias. Values of ρ c range between −1 and 1; values close to 1 indicate perfect concordance (no bias); values close to 0 indicate no concordance; and values close to −1 indicate perfect discordance. See Lin (1989) for further details.

Database Characteristics
Our measurements of K sat vary across nine orders of magnitude, with the highest value of 1.309 × 10 −1 m s −1 from a poor fen in the rock barrens of the Boreal Shield, in Ontario, Canada, and the lowest, 1.310 × 10 −10 m s −1 , from the Moor House blanket peat complex in the North Pennines, UK. A large majority of samples are from the uppermost meter of peat, with conspicuous clusters of measurements at depths of approximately 0.3, 0.5, and 1.0 m (Figure 2a). Only 20 measurements are from 3 m or more below the peat surface; the three deepest samples in the database, 6 m below the surface, are all piezometer slug tests from raised bogs in the United Kingdom. All 10 categories of the von Post scale are represented in the database, most of them approximately evenly so. The exceptions to this are the least humified category, H1, which only appears in 16 records; and the most humified, H10, which appears only five times (Figure 2c). Dry bulk density ranges between 0.01 and 0.30 g cm −3 , with a median of 0.09 g cm −3 (Figure 2b). More than 40% of the database, or 1,071 records, have an unspecified value for the categorical microform factor. The large majority of these represent peatlands that lack clear surface microtopography, meaning that hummocks, lawns, and hollows cannot be identified, rather than an absence of field observations. The unspecified categories make up relatively small proportions of the categorical factors that represent tree cover, trophic type, and disturbance.
With multiple independent variables, bivariate plots cannot convey independence of any apparent effect from other predictors. Nonetheless, the scatterplots and boxplots in Figure 2 suggest some potentially important trends in the database, which we are able to test within our formal modeling framework and which we report in 11 of 18 subsequent sections. In particular, K sat appears to be negatively related to log 10 (depth), log 10 (dry bulk density), von Post score, and increasing oceanity. Figure 2 also suggests the possibility of important differences in log 10 (K sat ) according to peatland trophic type, KOI, and levels of disturbance.

Variable Selection
The effects of the continuous and categorical predictors are qualitatively similar in all four models, with log 10 (K sat ) decreasing strongly and significantly with increasing log 10 (depth) and-in models that include them-von Post and log 10 (dry bulk density) (p < 0.001) (Tables S1-S4 in Supporting Information S1). Depth, dry bulk density, and von Post score are highly influential upon model performance; their omission causes large, significant increases in AICc (p < 0.001 in all cases) and large reductions in 2 cv .
Models 1, 2, and 3, which include categorical factors to represent local-scale site characteristics, all predict that log 10 (K sat ) increases almost monotonically along the trophic gradient from ombrotrophic to minerotrophic peat.
Fitted parameters indicate that blanket peat consistently exhibits the lowest log 10 (K sat ), followed by raised bog peat and then poor fen peat. Depending on the model, either moderate fens or rich fens have the highest log 10 (K sat ). Sample sizes in some models are small for some of these categories. In particular, blanket peat (n = 3) and moderate fens (n = 11) are underrepresented in Model 1 and those parameter estimates should be regarded with caution. However, Model 3, in which all trophic types are well represented by large samples, predicts a monotonic increase in log 10 (K sat ) along the trophic gradient, with the largest parameter estimate belonging to rich fens. In Model 3, the difference in parameter estimates between blanket peat and rich fens is 2.042, equivalent to more than two orders of magnitude difference in K sat on its original, untransformed scale.
All three models that include local-scale variables (Models 1, 2, and 3) also show a common effect of microform, with log 10 (K sat ) consistently and significantly lower in lawn peat than in hummocks or hollows (Tables S1-S4 in Supporting Information S1). Parameter estimates for lawns, which occupy an intermediate position with respect to water-table depth, are typically close to zero; only in Model 2 is the parameter for lawns significantly different from the reference category (unspecified microform). In contrast, the parameter estimates for peat beneath both hummocks and hollows are consistently and significantly greater than lawns (Tables S1-S4 in Supporting Information S1). The size of this difference in K sat on its original, untransformed scale is approximately an order of magnitude in Model 1 and about half that in Models 2 and 3. All levels of the microform factor are well represented in all models; even the smallest sample, hollows, is represented by 71 records in Model 1, giving us confidence in this finding of low K sat in lawn peat.
All four models include the climate index KOI, which provides a highly significant improvement in model performance according to changes in both AICc (p < 0.001) and 2 cv . Parameter estimates for KOI indicate that log 10 (K sat ) is significantly lower in hyperoceanic sites than in continental and subcontinental sites (see Tables S1-S4 in Supporting Information S1). The value of the parameter for hyperoceanic sites is large and negative in Models 1 (−1.451) and 4 (−1.498) (p < 0.001 in both cases); the absolute values are smaller in Models 2 (−0.666) and 3 (−0.342), albeit still highly significant (p < 0.005). Parameter estimates for oceanic sites are not significantly different from the reference category, continental, and subcontinental sites, in any of the four models.
All four models include a random intercept for measurement method, which led to significant improvements in AICc. As well as measurement method, Model 3 also includes a random intercept for disturbance, which led to a significant improvement in AICc. Model 4 includes random intercepts for measurement method and the two categorical variables that could not be included as fixed effects because of their small spatial scale: microform and trophic type. During the iterative model-specification process, it became apparent that the categorical factors representing disturbance and tree cover exhibited confounding effects upon other variables, so we omitted them as fixed effects from all models. For some models, inclusion of one or both of these factors led to a significant improvement in model performance according to change in AICc. However, in all models, most levels of the disturbance factor had no significant independent effect upon log 10 (K sat ) compared to pristine sites. The significance of each level of disturbance-or the lack thereof-and even the direction of their effects, were highly inconsistent between models. The effect of tree cover upon log 10 (K sat ) was also inconsistent between models, with some models indicating that tree cover increases log 10 (K sat ) and others indicating a decrease. We took these conflicting, inconsistent results to indicate that disturbance and tree cover are likely to represent over-fitted variables in our models.

of 18
Because our primary aim is to develop models that can be used to predict unknown values of K sat , we deemed it prudent to omit these variables in the interests of protecting model generality. The omission of disturbance and tree cover typically caused only small decreases in 2 cv . The fixed-effect interactions depth × disturbance and depth × microform both exhibited high GVIF values in Models 1, 2, and 3, indicating problematic levels of correlation with other predictors, so none of the final models includes either of these interactions. After omitting these problematic variables, the remaining fixed effects are all highly consistent between models and their values and directions are stable.

Fitted Models
Model 1, which utilizes all covariates, is the most skillful of the three models, with r 2 = 0.764 and 2 cv = 0.747. Model 1 has low bias, as indicated by a high value of Lin's concordance correlation coefficient, ρ c = 0.865 (see also Table 2). The final specification of Model 1 includes fixed effects for depth, dry bulk density, von Post score, trophic type, KOI, and microform and a random intercept for measurement method. The final form of Model 1 is as follows: where depth is midpoint depth below the peatland surface (m); dbd is peat dry bulk density (g cm −3 ); von Post is von Post score expressed as an integer between 1 and 10 (dimensionless); blanket, raised bog, poor fen, moderate fen, and rich fen are dummy variables that can take binary values of either 0 or 1 to represent a site's trophic type (the reference category is unspecified trophic type); oceanic and hyperoceanic are dummy variables, used to represent local climatic conditions according to the KOI (see Section 2.5 for thresholds; the reference category is subcontinental climate); and hollow, lawn, and hummock are dummy variables used to represent local surface microform type (the reference category is unspecified microform type). Note. Random effects are accounted for when estimating fixed-effect parameters but do not have parameter estimates themselves; dbd is peat dry bulk density; n is the number of records used in model fitting; k is the number of fitted model parameters, including the fixed-effect intercept term; AICc is the corrected Akaike information criterion; r 2 is the coefficient of determination when regressing predicted on observed values using ordinary least squares regression; 2 cv is a cross-validated coefficient of determination, used to assess a model's ability to predict unknown data points; and ρ c is Lin's concordance correlation coefficient, used to assess model bias. See Section 2 for full details. Model 2, which omits dry bulk density in order to allow rapid field estimation of K sat , is less skillful than Model 1, with r 2 = 0.656 and 2 cv = 0.640. Like Model 1, Model 2 also has low bias, with ρ c = 0.804 (see also Table 2). The final form of Model 2 is as follows: where variables are as specified for Equation 2.
Model 3, which utilizes all records in the database, is substantially less skillful than Models 1 and 2, with r 2 = 0.529 and 2 cv = 0.522 and exhibits greater bias, with ρ c = 0.667 (see also Table 2). The final form of Model 3 omits dry bulk density and von Post score but is otherwise similar to Model 1: where variables are as specified for Equation 2. cv is a coefficient of determination, calculated by regressing modeled values on measured values, using a five-fold cross-validation to indicate predictive ability for unseen data points; and ρ c is Lin's concordance correlation coefficient, used to assess model bias. See Section 2 for full details.

10.1029/2022WR033181
14 of 18 Model 4 is marginally less skillful than Model 3, with r 2 = 0.478 and 2 cv = 0.482 and exhibits a similar level of bias, with ρ c = 0.665 (see also Table 2). The final form of Model 4 is considerably simpler than the other three models, due to it including many fewer categorical levels, but its remaining parameter values are nonetheless similar: log 10 ( sat ) = −7.210 − 1.087 log 10 (depth) − 2.252 log 10 (dbd) − 0.275(oceanic) − 1.498(hyperoceanic), where variables are as specified for Equation 2.
Tables S1-S4 in Supporting Information S1 contain further details of all four models.
All four models exhibit some degree of directional bias, which causes them to underestimate the highest values of log 10 (K sat ) and overestimate the lowest values. This can be seen in Figure 3 as a difference in slope between the line of perfect concordance and the best-fit line between modeled and measured values of log 10 (K sat ). The magnitude of this bias is low in Models 1 and 2, as indicated by the high ρ c values, but is more pronounced in Models 3 and 4. Bias in all models passes through zero for intermediate values of K sat , between approximately 10 −5 and 10 −4 m s −1 .

Discussion
Some of our findings, particularly those relating to the roles of the three continuous covariates, are largely unsurprising, showing that peat K sat declines strongly and nonlinearly with increasing depth, degree of decomposition, and dry bulk density. These findings are in broad agreement with numerous previous observational studies from individual sites, the data from many of which are included in our database (e.g., Branham & Strack, 2014;Morris et al., 2019;Päivänen, 1973). Increasing dry bulk density is related to the constriction and closure of pores due to mechanical compaction and collapse (Clymo, 1992;Schlotzhauer & Price, 1999); while increasing levels of decomposition are associated with increasingly disconnected pores due to the biological breakdown of the peat matrix, particularly if also accompanied by compression (McCarter et al., 2020). Both dry bulk density and degree of decomposition typically increase with depth below the surface (McCarter et al., 2020;Morris et al., 2019), meaning that these three covariates are all positively correlated with one another (Figure 2, Table 3). Nonetheless, consistently low GVIF values for all three covariates indicate that each exerts an important independent control upon log 10 (K sat ), and that any multicollinearity is not so great as to reduce confidence in parameter estimates.
During the specification of all models, removal of log 10 (depth) led to a strong and significant decrease in model performance as measured by change in AICc, and led to important reductions in 2 cv . For the pragmatic purposes of predicting K sat , it seems that the inclusion of log 10 (depth) leads to an increase in accuracy, even in Model 1 which also includes both log 10 (dry bulk density) and von Post score.
The effect of KOI, which is manifest as an increase in K sat in more continental sites, may be attributable to a number of possible mechanisms. One plausible explanation is the development of macropores in continental peat due to greater depths of freezing during winter, although the pore-scale controls on peat hydraulic properties in cold environments are poorly understood (Hayashi, 2013;McCarter et al., 2020). Branham and Strack (2014) suggested that characteristic differences in peat floristic composition between continental and maritime sites may affect pore structures and hydraulic properties. Liu and Lennartz (2019) showed that woody peat has a higher K sat than Sphagnum peat of comparable dry bulk density, while Crockett et al. (2016) found that woody fen peat was more permeable than sedge fen peat. Although we excluded tree cover from our models due to numerical problems, KOI may contain some of this effect because continental sites are more commonly treed than maritime sites.
The strong and consistent changes in log 10 (K sat ) between microform types and along the ordinal gradient of peatland trophic types-increasing from blanket peat, through raised-bog and poor-fen peat, to moderate-and rich-fen peatare also noteworthy. Pore-scale studies of peat hydraulics paint a complex picture of the relationship between peat floristic composition and hydrophysical properties, mediated primarily through the susceptibility of different peat types to pore closure through compression and decomposition (McCarter et al., 2020). The remains of sedges and reeds that characterize fen peat are thought to decompose more readily than the Sphagnum peat found in bogs (  a larger proportion of small pores near the surface of a sedge-peat profile. Model 3 does not include von Post as a fixed effect, and the post hoc parameter estimates for microform and trophic type in Models 1 and 2 are strictly not independent of the parameter for the continuous variable von Post. As such, it may be tempting to attribute the strong effects of microform and trophic type to differences in peat compaction and/or humification according to floristic composition. However, the effect of minerotrophic peat types in Models 1, 2, and 3 is to increase log 10 (K sat ), rather than reduce it. Furthermore, low GVIF scores for log 10 (dry bulk density), von Post, microform, and trophic type in Models 1, 2, and 3 indicate that these variables are likely to represent separate sources of variance. It therefore appears that the effects of microform and trophic type upon K sat represent something other than different levels of peat humification and/or compaction. When poorly decomposed, the Sphagnum remains that characterize ombrotrophic peat have an abundance of large pores, reflecting the stem and branch structures of live Sphagnum (McCarter et al., 2020;Weber et al., 2017). However, hyaline cells in Sphagnum also lead to a large fraction of immobile pore space, the hydraulic effects of which may be retained even at high levels of decomposition (McCarter et al., 2020) and which may be partly responsible for the lower permeability of our blanket peat and raised bog samples. Additionally, although sedge peat may decompose more readily than Sphagnum peat, sedge peat has been shown to be less susceptible to the exponential loss of large pores during decomposition (McCarter et al., 2020;Rycroft et al., 1975). The orientation of sedge and reed remains also gives rise to the possibility of horizontal layering that may provide preferential flow paths (Baird & Gaffney, 2000;McCarter et al., 2020). Finally, at the hillslope scale, blanket peats form on slopes that may be as steep as 10°, meaning that peat may require a low permeability to retain enough water to prevent wasting. The development of low-K sat peat in blanket bogs may therefore be a result of autogenic feedbacks between peat accumulation, hydraulic properties, and water budget (cf. Holden, 2005Holden, , 2009Lewis et al., 2012;Waddington et al., 2015).
It has been shown that in some peatlands, preferential flow through pipes and macropores can be important components of the water budget (Holden & Burt, 2002;Rezanezhad et al., 2016). The small sample volumes of the individual K sat measurements in our data set mean that they cannot capture the site-scale influence of pipes, which are at least 1 cm in diameter. As such, we suggest that while our models are well suited to estimating K sat in a peat matrix, they may underestimate bulk, landscape-scale hydraulic conductivity in sites with well-developed preferential flow paths (cf. Glaser et al., 2021). Additionally, it seems possible that our models still contain some artifact of the measurement method used to determine K sat , despite its inclusion as a random effect, due to characteristic depth distributions of the various methods. For example, K sat values determined using piezometers are generally smaller than those collected using other methods (see Figure 2i), but they also come from the deepest peat where K sat is naturally low (see Figure S1 in Supporting Information S1), because slug tests can only be performed when the piezometer intake is fully below the water table and because piezometers can be readily installed to depths of several meters if desired. Laboratory methods, on the other hand, are mainly limited to more permeable, near-surface peats because recovery of intact samples from several meters' depth is usually impractical. Similarly, most of the determinations made using tension-disc infiltrometers are from blanket peats, in which K sat is low compared to other peatland trophic types. As such, it is difficult to judge whether some of the apparent artifacts of measurement method are, in fact, attributable in part or in whole to depth. Nonetheless, most measurement methods are represented at intermediate depths of ∼1 m ( Figure S1 in Supporting Information S1).
The omission of local-scale fixed-effect predictors from Model 4 allows it to be used for predictions at large scales, such as in land-surface schemes within Earth system models. Large-scale representations of peatlands within Earth system models have, until recently, relied on frameworks developed for mineral soils, with static soil properties, leading to a limited ability to simulate important feedbacks in peat development (Frolking et al., 2009). Efforts are ongoing to incorporate dynamic peat hydrophysical properties, including K sat , into land surface schemes (e.g., Chadburn et al., 2022), which our Model 4 may help to inform. Model 4 requires only depth, peat dry bulk density, and the climate index KOI to make predictions of peat K sat , with a level of predictive accuracy that may still be considered attractive ( 2 cv = 0.482).
The high predictive skill of our models, particularly Model 1, allows K sat to be estimated from simpler, cheaper variables with a degree of confidence that was previously not possible. Our models are more skillful than those previously developed from meta-analyses of peat K sat measurements across multiple sites (e.g., Wang et al., 2021), largely because they utilize several continuous and categorical predictors that previous studies did not consider. Consequently, our models contain a large number of parameters and may at first appear to be data-hungry and impractical. Although the categorical predictors each have several possible parameters, all but one are redundant in any given application. For instance, Model 1 has 14 fitted parameters 16 of 18 yet requires only six pieces of information to estimate K sat (depth, dry bulk density, von Post score, microform, trophic type, and KOI), while Model 4 requires only three (depth, dry bulk density, and KOI). Moreover, the categorical predictors used in our models (microform and trophic type) are simple and inexpensive to record in the field. Most researchers already collect this basic descriptive information during field sampling, and we suggest that doing so should be considered part of common good practice. The KOI can be estimated from freely available climate data, such as those used here. Dry bulk density and von Post score are more labor-intensive to collect, particularly the former, which requires more careful sample extraction and laboratory drying. The availability of both dry bulk density and von Post score make a large difference to our models' predictive abilities, and we encourage the collection of these data in situations where K sat may need to be estimated rather than measured directly. A full set of coefficients for each model is included in the Tables S1-S4 in Supporting Information S1, allowing readers to apply our models to their own data.

Conclusions
Our meta-analysis of a large catalog of published and unpublished peat K sat data yielded highly skillful predictive equations; our best predictive model explains three quarters of the variance in unseen values of log 10 (K sat ). Our models may provide valuable tools for researchers seeking to estimate peat K sat from simpler, cheaper measurements, in a manner similar to the application of pedotransfer functions in mineral soils. Our equations may also find utility in improving the representation of peatlands in land surface modeling. Dry bulk density and von Post scores are particularly valuable to the estimation of peat K sat and ideally should be measured if K sat is to be modeled. Even without these variables, acceptable estimates of peat K sat can still be made. Even our least skillful model, which is designed for use in large-scale land-surface schemes and necessarily omits some important variables, explains nearly half of all variance in log 10 (K sat ). Simple, commonly recorded categorical measurements such as surface microform and peatland trophic type also contribute to our models' high predictive performance compared to those from previous meta-analyses. It might be possible to extend our approach to peat hydrophysical properties other than K sat , such as those that describe unsaturated moisture-retention characteristics and, as datasets grow, to important locations not represented in our database, such as Russia and the tropics.

Data Availability Statement
All data reported in this article can be found in Data Set S1 and Data Set S2.