Drivers of Holocene palsa distribution in North America

Palsas and peat plateaus are climatically sensitive landforms in permafrost peatlands. Climate envelope models have previously related palsa/peat plateau distributions in Europe to modern climate, but similar bioclimatic modelling has not been attempted for North America. Recent climate change has rendered many palsas/peat plateaus in this region, and their valuable carbon stores, vulnerable. We fitted a binary logistic regression model to predict palsa/peat plateau presence for North America by relating the distribution of 352 extant landforms to gridded modern climate data. Our model accurately classified 85.3% of grid cells that contain observed palsas/peat plateaus and 77.1% of grid cells without observed palsas/ peat plateaus. The model indicates that modern North American palsas/peat plateaus are supported by cold, dry climates with large seasonal temperature ranges and mild growing seasons. We used palaeoclimate simulations from a general circulation model to simulate Holocene distributions of palsas/peat plateaus at 500-year intervals. We constrained these outputs with timings of peat initiation, deglaciation, and postglacial drainage across the continent. Our palaeoclimate simulations indicate that this climate envelope remained stationary in western North America throughout the Holocene, but further east it migrated northwards during 11.5e6.0 ka BP. However, palsa extents in eastern North America were restricted from following this moving climate envelope by late deglaciation, drainage and peat initiation. We validated our Holocene simulations against available palaeoecological records and whilst they agree that permafrost peatlands aggraded earliest in western North America, our simulations contest previous suggestions that late permafrost aggradation in central Canada was climatically-driven. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Peatlands in the northern circumpolar permafrost zone have been estimated to cover more than 3.5 million km 2 , and to contain 277.3e338 Pg of soil organic carbon (Tarnocai et al., 2009;Hugelius et al., 2014). While frozen, this carbon store is rendered largely inert, although increasingly widespread observations of permafrost thaw across northern high latitudes have raised concerns as to the future of this vast carbon stock under warming climates (Camill, 2005;Schuur et al., 2008Schuur et al., , 2015. Rising temperatures can exacerbate decomposition of permafrost peatlands by increasing surficial drying and reactivating microbial processes (McCalley et al., 2014;Swindles et al., 2015;Singleton et al., 2018). Once underlying ice thaws, the development of saturated, sedge-dominated fens strongly promotes carbon release as methane (CH 4 ) (Frolking et al., 2011;Wu and Roulet, 2014), a potent greenhouse gas (Collins et al., 2013). However, the insulating properties of peat mean that permafrost peatlands may thaw more slowly than mineral-soil and sedimentary permafrost, enabling their persistence in warming climates as relict permafrost peatlands (Vitt et al., 1994;Halsey et al., 1995;Smith and Riseborough, 2002). It is possible that peat accumulation rates might also increase in response to warmer growing seasons and persistently saturated conditions (Estop-Aragon es et al., 2018;Taylor et al., 2019). On the other hand, losses of organic carbon from deep, previously frozen peat layers may prevent net accumulation for several decades or even centuries after thaw (O'Donnel et al., 2012;Jones et al., 2017).
Permafrost in peatlands can be identified by the presence of perennial frost mounds, known as palsas or peat plateaus depending on their spatial extent (Brown et al., 1970;Zoltai et al., 1993;Laberge and Payette, 1995;Sepp€ al€ a, 2011). These landforms are widely believed to have developed in response to mid-Holocene cooling and are thought to be highly climatically sensitive, being primarily distributed towards the furthest extents of the modern discontinuous permafrost zone (Luoto et al., 2004a;Kuhry and Turunen, 2006;Treat and Jones, 2018). Where palsas/peat plateaus have thawed in response to recent climatic changes, increased emissions of nitrous oxide (N 2 O) and CH 4 have been recorded (Voigt et al., 2017;Kirkwood et al., 2019). Field observations suggest that palsas/peat plateaus aggrade under cold, dry climatic regimes that enable deep winter frosts to persist through the thaw season (Sepp€ al€ a, 1986(Sepp€ al€ a, , 2011. In addition to a suitable climate, Sphagnum-dominated communities are thought to be necessary for palsa/peat plateau development, as are sufficient peat depths in order to insulate ground-ice lenses Beilman et al., 2001;Luoto and Sepp€ al€ a, 2002;Jorgenson et al., 2010). At a site level, the presence of trees and shrubs may also protect permafrost by providing shade against insolation and by intercepting precipitation Shur and Jorgenson, 2007), although snow trapping by spruce stands can sometimes exacerbate permafrost degradation (Jean and Payette, 2014). Recent observations of degrading palsas/peat plateaus across northern high latitudes indicate clearly that climate warming is rapidly changing these ecosystems (Swindles et al., 2015;Borge et al., 2017;Mamet et al., 2017). The impacts of climatic change on these landforms at broad spatial scales may be better understood by relating their present distributions to modern climate data. Previous use of climate envelope models to predict palsa/peat plateau distributions has been restricted to Fennoscandia in northern Europe, where these landforms have been found to exist within a narrow climate envelope (Luoto et al., 2004a;Fronzek et al., 2006;Parviainen and Luoto, 2007;Aalto and Luoto, 2014;Aalto et al., 2017). In particular, Luoto et al. (2004a) found that Fennoscandian palsas/peat plateaus exist in optimum climatic ranges for mean annual temperature (À4.99 C to À2.87 C), seasonal temperature range (22.7 Ce27.7 C) and frost number (0.58e0.63) (a normalised freeze/thaw ratio: Anisimov and Nelson, 1997). An optimum threshold for mean annual precipitation of no more than 445 mm yr À1 was also identified. Later studies have demonstrated the important influence of the seasonality of precipitation, annual temperature extremes and growing degree days on palsa distributions (Parviainen and Luoto, 2007;Aalto et al., 2017). Specifically, reductions to seasonality and vegetation cover often increase the degradation of palsas/peat plateaus (Sepp€ al€ a, 2011), although these relationships are complex. While the modern climate space of palsa/peat plateau distributions has been extensively researched in Fennoscandia, it is unclear whether similar relationships apply elsewhere. In particular, the vast permafrost regions of North America and Siberia display much larger environmental gradients, but are understudied in comparison to northern Europe.
Attempts to ascertain broad-scale bioclimate relationships for North American permafrost peatlands have focussed on linear, temperature-dependent relationships and have been spatially limited to Alberta, Saskatchewan and Manitoba in western Canada, where permafrost peatland distributions are well constrained (Vitt et al., 1994(Vitt et al., , 2000Halsey et al., 1995;Beilman et al., 2001). Within this region, peat plateaus are found in mean annual temperatures below À1 C (Vitt et al., 1994), although under certain temperature regimes permafrost and non-permafrost peatlands are both present (Halsey et al., 1995;Beilman et al., 2001). This may reflect the importance of additional climatic factors (e.g., precipitation), and the role of local hydrological setting, topography and vegetation cover. Although detailed climate envelope modelling has yet to be attempted for North American palsas/peat plateaus, sufficient records of the distributions of these landforms now exist to make such analyses feasible. Furthermore, the spatiotemporal dynamics of Holocene ice sheet retreat (Dyke et al., 2004;Tarasov et al., 2012), terrestrial drainage (Teller and Leverington, 2004) and peatland initiation (Jones and Yu, 2010;Packalen et al., 2014;Morris et al., 2018) in North America are all well constrained. In Siberia, palsas/ peat plateaus are dispersed more broadly than in other regions (Kirpotin et al., 2011), and their distribution is not well defined, particularly in remote Eastern Siberia. The dearth of available data from Siberian permafrost peatlands currently makes meta-analyses impractical, so here we focus on North America.
Attempts to reconstruct Holocene distributions of North American permafrost peatlands have focussed exclusively on core-based, palaeoecological evidence, with hydroseral transitions in plant macrofossil records primarily used to indicate permafrost developments. Zoltai et al. (1995) have suggested that the southern limit of permafrost peatlands in western continental Canada was located 500 km further north at c. 6.0 ka cal. BP than at present. This idea is well-cited but is based on only a single temporal snapshot. Treat and Jones (2018) conducted the first genuine time series analysis of permafrost aggradation/degradation at large spatial scales in North American peatlands. They found that permafrost developed earlier in western North America than in the east. However, Treat and Jones (2018) focused solely on 6.0 ka cal. BP onwards, meaning that broad-scale distributions of palsas/peat plateaus during the early Holocene remain poorly understood. Additionally, only a limited number of palaeoecological records exist in some parts of North America, which makes regional comparisons more difficult. It therefore remains unclear whether North American palsa/peat plateau distributions have shifted in response to Holocene climate in a manner that may be thought of as analogous to reconstructions of moving treelines (Sorenson, 1977;Rochefort et al., 1994;Payette et al., 2002).
In this study, we aim to simulate the Holocene spatiotemporal dynamics of palsas/peat plateaus in North America and constrain the earliest possible timing at which these landforms could have developed. To do so, we used palaeoclimate simulations from a general circulation model to drive a climate envelope model fitted to modern climate-palsa relationships. We validate our simulated Holocene distributions of palsas/peat plateaus against available palaeoecological evidence from the literature. Specifically, we address the following objectives: i) To establish the modern climate space that maintains North American palsas/peat plateaus. ii) To simulate the spatiotemporal development of this palsa/ peat plateau climate envelope across North America throughout the Holocene, based on palaeoclimate simulations. iii) To constrain these Holocene predictions using broad-scale, non-climatic factors, including the timing of North American deglaciation, terrestrial drainage and peat initiation.

Overview
We developed a catalogue of 352 North American palsas and peat plateaus and related their locations to present-day climate variables, using binary logistic regression to predict the presence or absence of these landforms in a spatial grid across the continent. In doing so, we identified an envelope of climatic conditions that are suitable for the maintenance of these landforms. We conducted all of our analyses at a spatial resolution of 0.5 latitude Â 0.5 longitude, to match the resolution of our modern and palaeoclimatic datasets. The logistic regression model (LRM) predicts the likelihood of palsa/peat plateau occurrence in each 0.5 Â 0.5 grid square based entirely on climatic driving variables. We then applied the LRM to palaeoclimate simulations from the HadCM3 coupled atmosphere-ocean general circulation model to estimate changes in the distribution of the palsa/peat plateau climate envelope at 500-year intervals throughout the Holocene. We further constrained our estimates of Holocene permafrost peatland distributions in North America using: i) an interpolated map of peat initiation dates; ii) dates of deglaciation from the GLAC-1D ice sheet reconstruction (Tarasov et al., 2012;Ivanovic et al., 2016); and iii) dates of drainage of locations that were once occupied by postglacial lakes and seas, also using GLAC-1D. Our final constrained simulations seek to identify the earliest possible timing of palsa/ peat plateau presence for each grid cell in northern North America.

Modern distribution of permafrost peatlands
We collated a database of modern palsa/peat plateau locations in North America by performing a structured literature search, using the search terms "palsa", "peat plateau", "permafrost plateau" and "permafrost peatland" in combination with "Canada", "USA", and the names of selected Canadian provinces and territories (e.g. "Quebec", "Ontario") in Google Scholar. Where available, we studied evidence such as site descriptions and photographs to corroborate the presence of palsas/peat plateaus at each location. Palsas that were described as "lithalsas", "mineral palsas" or "inorganic palsas" were rejected because they are unlikely to have a substantial peat cover, and because permafrost beneath mineral soils exhibits different responses to climate (Pissart, 2002;Wolfe et al., 2014). We recorded the original coordinates of each palsa/peat plateau, and where necessary converted these to decimal degrees. Where palsas/peat plateaus were reported without exact coordinates, we digitised their locations from site maps to the nearest 0.5 Â 0.5 grid cell. Our catalogue comprises a binary map of the presence and absence of palsas/peat plateaus for each grid cell within the study area (Supplementary Dataset S1). We limited our study area to 44e80 N and 52e168 W, encompassing all recorded palsas and peat plateaus in North America. We omitted Greenland from our analysis due to its small proportion of unglaciated land area throughout the Holocene. We omitted grid cells that contain less than 50% land coverage because this matched the land mask applied to the modern climate data.

Modern climate data
We extracted monthly North American climate data from the CRU TS 4.02 dataset and averaged values from 1961 to 1990 to represent modern climate (Supplementary Dataset S2). This climatology is interpolated from climate stations across North America to a 0.5 Â 0.5 spatial resolution (Harris et al., 2014). Some palsas/peat plateaus in the database may be relict features that formed under earlier climates, before the recent influence of human induced warming, and may now exist under climatic conditions in which new palsas/peat plateaus cannot form (Halsey et al., 1995). Our use of 20th century climate data, when the magnitude of human-induced climate change was less than at present, was intended to help reduce any disequilibrium between the observed palsas/peat plateaus and the modern climate data. This time period also matches the baseline period used to bias-correct the general circulation model output (see below). During the model construction process we experimented with fitting our model to a preindustrial HadCM3 climate simulation (see Morris et al., 2018 for details) rather than to modern climate data. However, we found that models fitted using the modern climatology performed slightly better in predicting the current distribution of palsas/peat plateaus than did the simulated preindustrial climate, so we proceeded with the CRU TS 4.02 dataset.

Palaeoclimate simulation data
The palaeoclimate data (11.5e0.5 ka BP) were obtained from the 50-year climate means from a suite of equilibrium-type simulations performed with the UK Met Office Hadley Centre's HadCM3 coupled atmosphere-ocean general circulation model with dynamic vegetation (Gordon et al., 2000;Pope et al., 2000;Valdes et al., 2017) (Supplementary Dataset S3). The simulations were run at 500-year intervals from 11.5 ka to present; and are the same as those used by Morris et al. (2018). HadCM3's simulated climate was driven by changing boundary conditions e including solar insolation, atmospheric trace gases (carbon dioxide, methane, and nitrous oxide), and parameters relating to ice sheet evolution (ice sheet geometry, land-sea mask, global orography and bathymetry) e that were updated between each 500-year timestep, informed by the latest palaeo-evidence for each time interval (including Berger, 1978;Loulergue et al., 2008;Lüthi et al., 2008;Schilt et al., 2010;Veres et al., 2013;Peltier et al., 2015). Thus, these equilibrium-type snapshot experiments are kept broadly in line with the current Paleoclimate Model Intercomparison Project Phase 4 (PMIP4) protocol for transient simulations of the last deglaciation (Ivanovic et al., 2016), and the HadCM3 climate model predicts the climate response to these changing forcings. We downscaled the climate means from the model's native resolution (3.75 longitude Â 2.5 latitude) to a 0.5 Â 0.5 grid using bi-cubic spline interpolation; and performed regional bias correction against modern instrumental data spanning 1961e1990 CE (New et al., 1999). The methods for setting up and processing these simulations, including the downscaling and bias-correction, are more comprehensively described by Morris et al. (2018;Supporting Info). In common with Morris et al. (2018) only terrestrial climate data are used here.

Climatic predictors
We identified candidate climate variables that have previously been proposed to control the distribution of palsas and peat plateaus at coarse spatial scales in Fennoscandia (Luoto et al., 2004a;Fronzek et al., 2006;Parviainen and Luoto, 2007;Fronzek et al., 2011;Aalto and Luoto, 2014;Aalto et al., 2017). These predictors (Table 1) were all derived from the mean monthly temperatures and total monthly precipitation values we extracted from the CRU TS 4.02 climatology and HadCM3 palaeoclimate simulations.
We calculated mean annual temperature (MAT) from the monthly temperatures for each year. The seasonal temperature range (TRANGE) index is the difference between the warmest and coldest mean monthly temperatures and is identical to the continentality index presented by Luoto et al. (2004a). Degree days are the annual time integral of monthly temperatures relative to a specified threshold. Thawing degree days (TDD) and freezing degree days (FDD) are the time integrals of monthly temperatures above or below 0 C respectively, while Growing Degree Days (GDD 5 ) is the time integral of monthly temperatures above 5 C. We combined TDD and FDD to calculate the frost number (FROST), a normalised freeze/thaw ratio previously used by Anisimov and Nelson (1997) to predict permafrost status at coarse spatial scales: We calculated three precipitation metrics. We summed each month's total precipitation to give total annual precipitation (APREC). We calculated water precipitation (RAINFALL) and snow precipitation (SNOWFALL) as the total annual precipitation that occurred during months with MAT above and below 0 C, respectively.
We omitted a number of variables that have previously been linked to permafrost peatland distributions. For example, the mean January and July temperatures, previously used by Parviainen and Luoto (2007), may not necessarily always relate to the coldest and warmest months, respectively. Similarly, the seasons defined by Luoto et al. (2004a) for summer and winter precipitation used somewhat arbitrary choices of months (Carcaillet and Richard, 2000;Harrison et al., 2003).

Statistical modelling
We used binary logistic regression to develop a predictive model of palsa and peat plateau occurrence from the candidate pool of modern climatic predictors. Logistic regression models are commonly used to relate continuous predictors to binary observations within species distribution studies and bioclimatic research (Koutsias and Karteris, 2000;Pearce and Ferrier, 2000). We conducted all statistical analyses in IBM SPSS Statistics 23.
We began by selecting six predictors that have previously been linked to palsa/peat plateau distributions at coarse spatial scales: FROST, MAT, TRANGE, RAINFALL, SNOWFALL and GDD 5 . To limit multicollinearity amongst predictors we did not include FDD or TDD, because FROST combines both of these into a dimensionless freeze/thaw ratio. Additionally, we used SNOWFALL and RAINFALL rather than APREC, enabling the insulating effects of snow cover (Aalto et al., 2017) and dry surficial soils (Seppӓlӓ, 2011) to be represented separately.
A Spearman's Rank correlation matrix (Supplementary Dataset S5) indicates that all six climatic predictor variables are significantly correlated with one another, a situation known as multicollinearity. Multicollinearity is common in climatic datasets, but care must be taken when interpreting the results of models fitted using multicollinear predictors (Mac Nally, 2000). We followed a rigorous fitting procedure to construct a model that balances parsimony with predictive skill, meaning that the overall predictions of the model are likely to be robust when taken as a whole. However, multicollinearity means that the relative importance of individual predictors should be interpreted with caution.
To begin with we entered all six predictor variables simultaneously (block entry). In common with Luoto et al. (2004a), we also included the squared form of each variable (FROST 2 , MAT*, TRANGE 2 , RAINFALL 2 , SNOWFALL 2 and GDD 5 2 ) because the probability of palsa/ peat plateau presence is unlikely to increase monotonically with each climate variable, but rather exhibit an optimum window. To preserve the sign of negative temperatures in its quadratic term, we calculated MAT* as the product of MAT and its absolute value, |MAT|. We used a stepwise backwards-deletion procedure to remove non-significant predictors sequentially until all remaining variables were significant predictors of palsa/peat plateau presence (Student's t; p < 0.05 threshold), and which, if removed, significantly reduced the model's predictive performance according to deviance (-2LL) scores (p < 0.05 criteria). We chose a backwards deletion approach because forward additive approaches can be prone to suppressor effects, where the presence of existing predictors may boost the significance of new variables (Field, 2013). After fitting main effects, we tested a variety of first-order interaction terms. Once our final model had been constructed, we calculated standardised parameter coefficients (b s ) for each predictor variable using the method of Menard (2011).
We evaluated the final model's performance by comparing predicted and observed values for presence and absence using three metrics: accuracy, informedness, and the area under the curve (AUC) of a receiver operating characteristic (ROC) plot (following Powers, 2011). Accuracy evaluates the proportion (0e1) of correctly classified cases (both presence and absence) (Powers, 2011). Informedness is a deep measure of how consistently a model predicts a case by considering both presence and absence, and indicates how informed a prediction is compared to chance (Powers, 2011). Informedness values range from 1 (all cases classified correctly) through 0 (random predictions) to À1 (all cases classified incorrectly).
The predictions of a LRM are a set of probabilities that the condition represented by the response variable (palsa/peat plateau presence in our case) is true for any given combination of the predictor variables. These predicted probabilities are usually rounded to 0 or 1 based on the 0.5 threshold. However, doing so fails to account for a subject's rarity and can lead to erroneous predictions in situations such as ours with rare observations (only 211 of our 9545 grid squares, or 2.2%, contain palsas/peat plateaus) (Fielding and Bell, 1997;Pearce and Ferrier, 2000). Increasingly, ecological modelling studies have begun to choose classification thresholds other than 0.5 that optimise appropriate evaluation criteria (Nenz en and Araújo, 2011;Wong et al., 2018;Kanagaraj et al., 2019). Because of the comparative rarity of our observed palsas/peat plateaus, we used a lower predicted probability threshold (see section 3.2, below) to determine palsa/peat plateau presence, chosen so as to maximise informedness. Where informedness is unaffected by rarity, some other evaluation metrics, such as markedness (see Powers, 2011) would be dominated by the high proportion of grid cells without palsas/peat plateaus.
AUC is a measure of model performance that is unaffected by prevalence and which, unlike accuracy and informedness, is independent of classification thresholds (McPherson et al., 2004;Marino et al., 2011). AUC scores range from 0.5 (random classification) to 1 (perfect classification) (Pearce and Ferrier, 2000).

Non-climatic constraints of palsa/peat plateau distributions
Our LRM predicted the climate space of modern palsa/peat plateau distributions but did not account for non-climatic controls. To constrain our predictions of palsa/peat plateau distributions further, we filtered our results for each grid square using several binary plausibility criteria, including: i) the presence of peat; ii) the absence of continental ice; and iii) the necessity for a terrestrial (rather than marine or lacustrine) environment.

Modern peatland distribution
We used PEATMAP (Xu et al., 2018) to estimate modern peatland distribution in North America. We rasterised the catchment-scale PEATMAP shapefiles (available from: http://archive.researchdata. leeds.ac.uk/251/) and reclassified to produce a binary raster map of peatland presence/absence for each of our 0.5 Â 0.5 grid cells. To improve our estimate of modern peat distribution, we supplemented the PEATMAP dataset with lithology data from northern Alaska's most recent ecological land classification report (available from: http://alaska.portal.gina.alaska.edu/catalogs/9549landscape-level-ecological-mapping-of-northern) (Jorgenson and Grunblatt, 2013). Finally, we retroactively classified a small number of grid cells that were not identified as containing peat by either PEATMAP or the additional data for Alaska, but where our catalogues of modern palsas/peat plateaus (see above) and/or peat initiation dates (see below) indicated extant peatlands. Overall, we estimated that 4533 grid cells (47.5% of the total grid cells within our study area) presently contain peat (Fig. 1).
Palsas have been suggested to develop primarily on fine-grained sediments, which display lower thermal conductivities than coarser grained sediments and bedrock, and this encourages the development of segregated ice (O'Neill et al., 2019). However, many North American peatlands e not only permafrost peatlands e have developed on fine-grained sediments, due to their poor drainage (Tarnocai et al., 2011). Furthermore, palsas have been observed to overlie coarse-grained sands in Bugristoye, Russia (Blyakharchuk and Sulerzhitsky, 1999). As such, it is unclear whether the presence of fine-grained sediments exerts a genuinely independent control on palsa development. Our inclusion of a modern peat presence filter removed the majority of geologically unsuitable grid cells, so we did not include a separate geological filter.

Holocene changes in peatland distribution
We constructed a catalogue of radiocarbon dates of peatland initiation, by combining a number of recent meta-analyses (MacDonald et al., 2006;Gorham et al., 2007;Jones and Yu, 2010;Korhola et al., 2010;Packalen et al., 2014;Treat et al., 2016Treat et al., , 2019Morris et al., 2018;Treat and Jones, 2018) (Supplementary Dataset S4). We recorded uncalibrated radiocarbon ( 14 C) basal dates and the associated laboratory error terms for all available North American sites. All dates were then recalibrated using the OxCal software by recording the median intercept value from the IntCal 13 curve at 95.4% confidence (OxCal, 2018); except for the calibrated dates presented by Gorham et al. (2007). Gorham et al. (2007) did not provide uncalibrated dates, but we took their reported calibrations (performed using CALIB 3.0: Stuiver and Reimer, 1993) to be broadly comparable with our own recalibrations, given the 500year temporal grain of the HadCM3 mean climates. The final database contained basal ages for 1756 individual peatlands across North America.
We grouped the calibrated basal dates into our 0.5 Â 0.5 grid cells and discarded all but the oldest in each cell so as to constrain the earliest possible dates at which palsas/peat plateaus might have begun to form. Our catalogue has basal dates for 921 of the 4533 grid cells that contain peat. We used nearest-neighbour interpolation to assign calibrated initiation dates to peatland grid cells without their own initiation dates. Where multiple basal dates were equidistant from the target grid cell we assigned the oldest date.

Deglaciation and post-glacial drainage of North America
We used the GLAC-1D reconstruction of global ice sheet evolution to estimate the date of deglaciation for each of the 9545 0.5 Â 0.5 grid cells in our study area, which is derived from a Bayesian calibration of the 3D Glacial Systems Model (Tarasov et al., 2012). The original glacial reconstruction for North America is calculated at 0.5 latitude Â 1.0 longitude resolution at 100-year time intervals and utilises physically-based constraints including geological evidence of ice-sheet margins, relative sea level data, and modern rates of uplift. We determined the time since deglaciation as the most recent 100-year time interval in which the grid cell of interest was ice-free. Although our HadCM3 palaeoclimate   (Xu et al., 2018) and ecological land classification data from northern Alaska (NOAK) (Jorgenson and Grunblatt, 2013). Additional sites were identified from locations of modern palsas/peat plateaus and peatland basal dates (this study). Black dots represent grid cells with observed palsas/peat plateaus located within. simulations included the ICE-6G_c ice sheet reconstruction (Peltier et al., 2015), we preferred the use of GLAC-1D to estimate time since deglaciation because its incorporation of physical evidence, a dynamic ice sheet model, and climate forcings improved the spatiotemporal resolution of our calculations . We also used GLAC-1D to estimate the date of terrestrial drainage for grid cells that were previously submerged by isostatic depression and post-glacial lakes, such as Lake Ojibway and Lake Agassiz (Teller and Leverington, 2004). We used reconstructed surface elevations to date land emergence as the most recent 100-year time interval when each previously submerged grid cell first reached an elevation above sea level.

Modern distribution
The 352 palsas and peat plateaus in our catalogue occupied 211 distinct grid cells, equivalent to 2.2% of the entire study area (Fig. 1). Of these 352 sites, 181 have been documented previously in published databases of permafrost peatlands (e.g., Treat et al., 2016;Zoltai et al., 2001), while to our knowledge 171 have not previously been collected into a single database. Modern palsas/peat plateaus span almost the entire longitudinal breadth of North America, with 98.5% located in a relatively narrow latitudinal band between 51 and 70 N. Modern palsas/peat plateaus are present in grid cells with cold (median MAT ¼ À5.15 C), dry (median APREC ¼ 421 mm yr À1 ), continental (median TRANGE ¼ 39.82 C) climates (Table 2).

Statistical modelling
Our final climate envelope model, fitted using modern climate data from the CRU TS 4.02 climatology (Fig. 2), predicts the distribution of climate space suitable for palsas/peat plateaus in North America (Table 3). Our LRM included RAINFALL, MAT, and GDD 5 as significant predictors, alongside the quadratic forms of these variables (MAT*, RAINFALL 2 , and GDD 5 2 ). We also included the significant, first-order interaction term TRANGE Â RAINFALL, along with the significant main effect TRANGE. We chose a classification threshold for the final LRM of 0.028 which maximised model informedness (0.624) when predicting the modern palsa/peat plateau distribution from modern climate data (Fig. 3).
Considering all grid cells within our study area (n ¼ 9545), the overall accuracy of our model using the chosen classification threshold was 77.2%. Considering only those grid cells with observed palsas/peat plateaus, 85.3% were correctly classified by our model; of those grid cells without observed palsas/peat plateaus, 77.1% were correctly predicted. A high AUC score of 0.865 indicated the high probability that, independent of the chosen classification threshold, our model would assign a higher predicted probability to a randomly-chosen grid square that contained an observed palsa/peat plateau (presence) than to a randomly-chosen grid square without (absence). This means that our model performed well in distinguishing the presence and absence of North American palsas/peat plateaus, and that the model's high accuracy was not simply a consequence of the rarity of true positive cases.
Our model classified a total of 2321 North American grid cells (24% of the study area) as currently within the climatic envelope suitable for palsas/peat plateaus (Fig. 3b). The widespread presence of false positive predictions within our results was expected, at least in part, because of the likely existence of some as-yetunrecorded palsas/peat plateaus across remote, poorly mapped areas such as portions of the Canadian Arctic. The modern climate envelope broadly follows the band of observed palsas/peat plateaus, forming a continuous band across Canada covering much of northern Quebec, the Hudson Bay Lowlands (HBL), northern Saskatchewan and Manitoba, the MacKenzie River Basin and the interior of Alaska (Fig. 3).
Our LRM indicates that the probability of palsa/peat plateau presence declines strongly with increasing MAT, but rises with increasing RAINFALL and GDD 5 . Additionally, the presence of these landforms is significantly negatively associated with the interaction term TRANGE Â RAINFALL. The relationships for MAT, RAINFALL and GDD 5 in our predictive model are all curvilinear, indicating that North American palsas/peat plateaus primarily exist within an optimal climatic envelope ( Fig. 4).

Holocene distribution of climate space
Simulated changes in Holocene distributions of the climate space suitable for palsa/peat plateau presence during 11.5e0.5 ka BP are presented in Fig. 5. Our simulations showed contrasting spatiotemporal dynamics of the climatic envelope between eastern and western North America. West of 115 W across Alaska, British Columbia, Yukon, much of the Northwest Territories, and Alberta, our simulations indicated that the climate envelope has remained relatively stationary throughout the Holocene. East of 115 W the climate envelope was located much further south at the opening of the Holocene than its present position, forming an arc through the modern Laurentian Great Lakes. This envelope shifted steadily northwards across eastern Canada until 6.0 ka BP when it reached close to its present position, centred on approximately 55 N. Although the centre of this band has remained relatively stationary since the mid-Holocene, our simulations indicated that it has gradually contracted during the late Holocene, with the number of North American grid cells classified as within the climate envelope declining from 2959 at 6.0 ka BP (31% of the study area) to 2321 (24%) at present. Fig. 6a shows interpolated, calibrated peat initiation ages within the 4533 grid cells that currently contain peat. Early peatland establishment (before 10.0 ka BP) occurred in non-glaciated

Deglaciation and post-glacial drainage
Estimated dates of terrestrial deglaciation and drainage, calculated using the GLAC-1D ice sheet reconstruction, are presented in Fig. 6b and c respectively. GLAC-1D replicated the collapse of the Cordilleran and Laurentide ice sheets well, including the deglaciation of western North America mainly before 10.0 ka BP, whilst the Laurentide ice sheet persisted across much of eastern Canada until approximately 5.5ka BP (Dyke et al., 2004). Following ice sheet retreat, our GLAC-1D simulations captured the development of post-glacial Lakes Agassiz, McConnell and Ojibway, which initially submerged vast regions of North America before draining by approximately 8.0 ka BP (Smith, 1994;Teller and Leverington, 2004). The model also captured isostatic rebound in wake of the retreating Laurentide ice sheet, which has caused continued land emergence from 8.0 ka BP until present along the coastlines of Hudson Bay, James Bay and the Canadian Arctic Archipelago (Hillaire-Marcel and Fairbridge, 1978).

Holocene palsa/peat plateau distribution
Our final Holocene simulations suggested that palsa/peat plateau distributions were restricted from following the moving climate envelope in eastern North America due to the timings of peat initiation, deglaciation and postglacial drainage (Fig. 7). Across our North American study area the number of 0.5 Â 0.5 grid cells classified by our LRM as suitable for palsas/peat plateaus increased dramatically throughout the Holocene, from 124 (1.3% of the total study area) at 11.0 ka BP to a maximum of 1812 (18.9%) at 1.0 ka BP (Fig. 8a). Our results indicate that the number of grid cells suitable for the presence of palsas/peat plateaus expanded most rapidly from 11.0 to 10.0 ka BP and from 8.5 to 6.5 ka BP, and that the suitable climate envelope may have contracted slightly between 5.5 and 5.0 ka BP. We find that the spatial distribution of suitable grid cells for North American palsas/peat plateaus is likely to have remained relatively constant since 3.0 ka BP. The areal extent of grid cells suitable for palsas/peat plateaus has expanded throughout the Holocene at a rate proportional to the expansion of peatlands (Fig. 8b).  3. Predicted modern distribution of the suitable climate envelope for palsa/peat plateaus in North America, expressed as; (a) a continuous probability scale, with the transition from yellow to light blue indicating the threshold value (0.028); (b) a map of peatlands classified as those within the climate envelope for palsa/peat plateaus (predicted probability ! 0.028) and those outside of the envelope (predicted probability < 0.028). Black dots represent grid cells that contain observed palsas/peat plateaus. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Modern palsa-climate relationships
The modern climate envelope predicted by our model as suitable for North American palsas/peat plateaus is largely similar to those established in northern Europe, albeit with somewhat broader ranges of temperature. Some North American palsas persist even where modern MAT is slightly warmer than 0 C (Table 2; Fig. 4), which illustrates the insulating capabilities of peaty soils (Parviainen and Luoto, 2007;Halsey et al., 1995). These soil properties enable certain palsas to persist beyond the southern limit of discontinuous permafrost, defined for North America as the À1 C isoline of MAT (Brown, 1970;French and Slaymaker, 1993).
The negative relationship between MAT and palsa/peat plateau probability is as expected, but the lower MAT threshold of À11.32 C for predicted palsa/peat plateau presence in our model was much colder than the MAT range of À4.42 Ce0.42 C reported for Fennoscandia (Parviainen and Luoto, 2007). In Fennoscandia, it has been suggested that extremely cold environments prevent the formation of palsas/peat plateaus because slow accumulation rates restrict the development of sufficient peat depths (Luoto and Sepp€ al€ a, 2002;Luoto et al., 2004a). However, our results imply that peat depths can become adequate for palsa/peat plateau formation if a low MAT coincides with a large annual temperature range, such as in the Northwest Territories, leading to growing seasons conducive to peat accumulation. Perhaps somewhat counter-intuitively, extended growing seasons may also encourage permafrost aggradation within peatlands by encouraging tree and shrub growth. Canopy shading reduces incident radiation, thus thinning the active layer; while canopy interception results in shallower snow depths, thereby reducing insulation during winter months R€ onkk€ o and Sepp€ al€ a, 2003). Indeed, the prevalence of wooded palsas in North America (Cyr and Payette, 2010) likely explains the greater importance of GDD 5 within our climate envelope model because palsas in Fennoscandia are often treeless (Zoltai and Tarnocai, 1971;Luoto et al., 2004b).
In common with previous research, our results advocate dry environments as optimal for palsa/peat plateau persistence, although a minimum RAINFALL threshold of 112 mm yr À1 appears to indicate that palsas/peat plateaus cannot form in the driest climates of Arctic deserts. The significant, negative response of palsa/ peat plateau presence to the interaction term TRANGE Â RAINFALL further highlights the importance of aridity. SNOWFALL is not significantly related to the modern distribution of palsas/peat plateaus at the 0.5 Â 0.5 grid scale of our analyses. In individual sites, snow depth is known to be related to palsa/peat plateau development (Sepp€ al€ a, 1994), but our results suggest that snow depths are controlled by phenomena acting at smaller spatial scales than our 0.5 grids, including the direction of prevailing winds and canopy interception from tree stands (Zoltai and Tarnocai, 1971;Sepp€ al€ a, 2011). Moreover, our model performed less well in highlymaritime, coastal regions where palsa/peat plateau presence is heavily influenced by the wind-scouring of local snow cover (Way et al., 2018).
The broader climate envelope occupied by North American palsas/peat plateaus compared to those in Fennoscandia is perhaps indicative of the greater environmental gradients across the study area. Fennoscandian palsas/peat plateaus are distributed north of the Arctic Circle (i.e., north of 66.55 N), across an area of approximately 67,000 km 2 with a climate that is more continental than the rest of Fennoscandia (Luoto et al., 2004a), but which is less continental than much of northern North America. By contrast, the distribution of North American palsas/peat plateaus extends much further south to 45 N and spans a greater range of climates, from the maritime Newfoundland and Labrador to the highly continental Canadian Arctic, Prairies and Alaskan interior. These environmental gradients render climate envelope models fitted to European palsa/ peat plateau distributions unsuitable for the prediction of North American distributions.
Due to the coarse spatial scale of the climate data used in this study, we necessarily omitted a number of local-scale factors that have previously been shown to be important to palsa/peat plateau formation at the sub-grid scale, including topography, vegetation, and hydrology (Sepp€ al€ a, 1986; Aalto and Luoto, 2014). Consequently, some of our model's false negatives (i.e., failure to predict some observed palsas/peat plateaus) may reflect isolated palsas/ peat plateaus that exist largely in response to local factors, rather than grid-scale climate, such as in alpine valleys where strong, local temperature inversions enable so-called 'warm permafrost' to persist (Lewkowicz and Coultish, 2004;Throop et al., 2012). In certain cases, permafrost can also persist in unsuitable climates as isolated patches of "ecosystem protected" permafrost, where overlying trees, shrubs and mosses maintain cold ground temperatures (Shur and Jorgenson, 2007). Alternatively, some false negatives may represent palsas/peat plateaus that were already in disequilibrium with the surrounding climate during the 1961e1990 baseline period. The possible inclusion of relict landforms in our catalogue may mean that our model overestimates the temperature suitable for palsa/peat plateau presence, but the magnitude of any such overestimate is impossible to evaluate. We attempted to minimise this effect by using climate averages from 1961 to 1990, a period slightly before present. Additionally, the lag time between the development of suitable climates and the development of palsas/peat plateaus is poorly understood (Camill and Clark, 2000;Luoto and Sepp€ al€ a, 2003). Regardless, our predicted upper MAT limit of 0.19 C for palsa/peat plateau presence is lower than that observed in Fennoscandia (0.42 C) (Parviainen and Luoto, 2007), giving us confidence that our climate envelope model is not substantially affected by this issue. It is also possible that our model overestimates the modern distribution of palsas/peat plateaus where large numbers of false positives exist (i.e. grid cells without observed landforms that are predicted positive), for example in the Great Slave Lowlands and in northern Quebec. Despite the presence of peat and suitable climatic conditions, recent geological surveys suggest that these regions are predominantly underlain by tills and bedrock (Geological Survey of Canada, 2014), substrates that are significantly less frost susceptible than fine-grained sediments (O'Neill et al., 2019). However, local observations indicate that lithalsas in the Great Slave Lowlands (Wolfe et al., 2014) and fen palsas near Scheffervile, Northern Quebec (Doolittle et al., 1992) are underlain by highly frost-susceptible silts and clays. Where false positives are present in our modern simulation, palsa/peat plateau distributions may therefore be discontinuous, developing only where fine-grained sediments are present, although this is impossible to test given the resolution of existing geological maps.

Spatiotemporal dynamics of Holocene palsa/peat plateau distributions
Our Holocene simulations indicate locations across North America that were climatically suitable for palsa/peat plateau formation; that contained peat; that were free of continental ice; and that had drained of any postglacial water bodies. These locations could therefore have supported palsas/peat plateaus once other favourable local conditions had developed, such as the establishment of moss communities and the development of sufficiently deep peat . Our simulations, therefore, provide estimates for the earliest timing of palsa/peat plateau establishment.
We validated our continental-scale predictions of Holocene palsa/peat plateau distributions for North America against local, multiproxy palaeoecological evidence of permafrost aggradation. Broadly, our simulations agree with the available evidence that palsas/peat plateaus developed earlier in western North America than in the east (Tarnocai and Zoltai, 1988;Treat and Jones, 2018). Our simulations suggest that western North America could have supported palsas/peat plateaus since the early Holocene, where a stable climate envelope has persisted since 11.5 ka BP, and where early peatland expansion was possible in ice-free environments . In eastern North America our predicted climate envelope shifted northwards during 11.5e6.0 ka BP, but once nonclimatic constraints are taken into account, our simulations showed only isolated examples of suitable grid cells for palsa/peat plateau presence prior to 7.5 ka BP. Our results, therefore, suggest that the high mobility of the suitable climate envelope during the early Holocene was not sufficient alone to cause substantial redistribution of these landforms, and that North American palsas/peat plateaus have, for the most part, not occupied regions outside their current range since 11.5 ka BP. North American peatlands steadily aggraded permafrost during the Neoglacial (starting c. 4.3 ka BP) and Little Ice Age (LIA) (c. 700e150 years BP) (Vitt et al., 1994;Arlen-Pouliot and Bhiry, 2005;Wanner et al., 2011;Fillion et al., 2014;Treat and Jones, 2018), when global reductions in incident solar radiation caused cooling across the continent (Mann et al., 2009;McKay et al., 2018). Our simulations reflect this, with peak numbers of suitable grid cells consistently recorded after 3.5 ka BP (Fig. 8a).
Our results suggest that palsa/peat plateau development in western North America was primarily constrained by the timing of peat initiation. The widespread expansion of peatlands in this region during 11.5e6.0 ka BP has mainly occurred within our simulated climate envelope, suggesting that these peatlands have been climatically suitable for palsas/peat plateaus since initiation. Our results indicate that palsas/peat plateaus could have begun to establish amongst the early peatlands located within non-glaciated regions of Alaska and the MacKenzie River Basin, an assertion supported by palaeocological records that date peatland permafrost aggradation at c. 8.0 ka cal. BP (Ovenden, 1982;Tarnocai and Zoltai, 1988;Tarnocai, 2010). Our simulations indicate that by the mid-Holocene, the distribution of grid cells suitable for palsa/peat plateau presence within western North America approximately matched the present-day configuration, with minor differences being driven primarily by late peat initiation. Notably, our model suggests that after 1.0 ka BP northern Alaska became less climatically suitable for palsas/peat plateaus, possibly due to regional reductions in GDD 5 . However, the general stability of our simulated climate envelope is in contrast to the proposition by Zoltai et al. (1995) that palsas/peat plateaus migrated south after c. 6.0 ka cal. BP in response to climatic cooling, which they based on extensive core-based evidence of permafrost aggradation in peatlands from western and central Canada. The 500-year temporal resolution of our palaeoclimate simulations was too coarse to simulate changes in palsa/peat plateau distributions during the LIA, but distributions of bog landforms in Alberta, Saskatchewan and Manitoba indicate that permafrost peatlands reached their maximum extent at this  time (Vitt et al., 1994;Halsey et al., 1995;Beilman et al., 2001).
Regional palaeoecological records generally indicate that permafrost aggradation lagged the development of suitable conditions predicted by our modelling (e.g., Treat and Jones, 2018). This lag time likely accounts for the time required for sufficient peat to accumulate and for suitable moss communities to establish. Other important controls on permafrost aggradation that were not included in our simulations include vegetation cover, topography and site hydrology (Sepp€ al€ a, 1986;Zoltai et al., 1995;Aalto and Luoto, 2014). Comparisons with core-based reconstructions are further complicated by sparse sampling of remote regions and the well-cited absence of direct ecological indicators of permafrost aggradation (Kuhry and Turunen, 2006;Oksanen, 2006;Treat et al., 2016;Treat and Jones, 2018). Based on the few paleoecological records of North American palsas/peat plateaus that are available (see Treat and Jones, 2018), considerable variation is observed in the time that permafrost aggradation lags behind our simulations, ranging from a few centuries to several millennia. The timings presented by our simulations should therefore be interpreted with caution. Nevertheless, we take confidence from the fact that our predictions concur with the broad spatiotemporal patterns indicated by existing palaeoecological evidence.
Although the simulated climate envelope migrated northwards across eastern North America, we suggest that late deglaciation, drainage and peat initiation restricted the formation of palsas/peat plateaus there until at least the mid-Holocene. Our results indicate that palsas/peat plateaus may have formed amongst early peatlands located north of the Great Lakes between 10.5 and 8.5 ka BP, an area that is now too warm for these landforms, although no palaeoecological records exist that might be used to test this hypothesis. Our simulations indicate that northern Ontario has been climatically suitable for the presence of palsas/peat plateaus since approximately 7.5 ka BP, with regional pollen records suggesting that permafrost developed immediately in newly deglaciated/ drained areas (Allard and Seguin, 1987). Late drainage of the modern Hudson Bay Lowlands delayed peat initiation, causing the number of suitable grid cells south of Hudson Bay to increase only gradually, despite stable climatic conditions. Core-based analyses from this region are particularly scarce, although the initial aggradation of peat plateaus near Ennadai Lake, Nunavut at c. 4.5 ka cal. BP (Sannel and Kuhry, 2008) and near Churchill, Manitoba at c. 2.2 ka cal. BP (Kuhry, 2008) are consistent with our simulations. In southern Quebec, suitable climatic conditions for palsas/peat plateaus developed as early as 7.0 ka BP, before the climate envelope progressed northwards across newly deglaciated regions where peatlands have since expanded. Although sparsely distributed, local plant macrofossil records from northern Quebec provide no evidence of permafrost prior to 5.5 ka cal. BP (Payette, 1988;Treat and Jones, 2018), which is when our simulated climate envelope reached its present position and became stationary. The total number of suitable grid cells in northern Quebec remained relatively constant in our constrained simulations from 3.0 ka BP onwards, during which time palaeoecological records indicate steady permafrost aggradation within peatlands in response to continued Neoglacial cooling and the onset of the LIA (Arlen- Pouliot and Bhiry, 2005;Fillion et al., 2014;Treat and Jones, 2018).

Conclusions
Modern distributions of North American palsas/peat plateaus are well described by a logistic regression model driven by mean annual temperature, rainfall, growing degree days and the interaction between rainfall and seasonal temperature range. Our results demonstrate that North American palsas/peat plateaus occupy a climate envelope that is largely similar to that in Fennoscandia, but which exhibits a broader range of temperature regimes. We  used our logistic regression model to simulate Holocene distributions of palsas/peat plateaus across North America, by estimating the earliest 500-year timestep that each grid cell became suitable for palsa/peat plateau presence according to climatic and nonclimatic factors. Our results suggest that western North America became suitable for palsa/peat plateau establishment during the early Holocene, with expansion of these landforms primarily limited by the timing of peat initiation. The suitable climate envelope for palsas/peat plateaus in eastern North America has shifted dramatically northwards since the early Holocene, but the distributions of these landforms are likely to have been heavily constrained by the timing of peat initiation, deglaciation and drainage. In particular, late peat initiation likely prevented palsa extents from following the moving climatic envelope in central and eastern North America, meaning that palsas/peat plateaus developed later in eastern North America than in the west.

Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.