Modelling Tundra Ponds as Initiators of Peat Plateau Thaw, Northern Hudson Bay Lowland, Manitoba

. Frozen peat in permafrost regions poses a potential source of increased greenhouse gas production should these deposits thaw. Ponds on frozen peat plateaus in northern Manitoba are numerically modelled as heat sources to determine their potential to promote thaw. Modelling indicates that anticipated climate warming of approximately 2 °C between 2020 and 2050 will produce taliks up to few metres thick beneath ponds a few tens of metres across. However, active-layer thickness in the subaerial parts of peat plateaus will not increase beyond the peat thickness. These findings assume 1) a climate warming rate under a moderately effective intervention in greenhouse gas production, 2) pond freezing regimes that represent both rapid ice formation and ice formation delayed by rapid snow accumulation and 3) snow thermal conductivities that anticipate snow conductivity increase during the freeze interval. These conditions and properties may turn out to be less conducive to talik expansion than the values that will actually occur. Despite these uncertainties, peat plateau pond sizes and plateau margin positions can be monitored to ascertain the onset of accelerated thawing.


INTRODUCTION
As one of the world's largest accumulations of soil organic carbon, the Hudson Bay Lowland (HBL) holds a potential for increased greenhouse gas (GHG) generation under climate warming (Keller et al., 2014;McLaughlin and Webster, 2014). Of particular interest is the permafrost part of the HBL in northern Manitoba and Ontario (Fig. 1a), where the frozen condition protects this peat from decay. However, countless lakes and ponds perforate the frozen peat of this region, particularly in northern Manitoba, implying that water bodies on permafrost peatlands in the HBL coincide with permafrost thaw and release of carbon dioxide (CO 2 ) and methane (CH 4 ) from decomposition of the underlying peat Laurion et al., 2010). Although pond formation or expansion in northern Manitoba is not apparent during the past 70 years of overall gradual air temperature warming at Churchill and Gillam, ongoing erosion at the edges of the larger water bodies continues to expose peat to biologic transformation (Macrea et al., 2004;Dyke and Sladen, 2010;Sannel and Kuhry, 2011). Expansion of ponds in peatlands has been reported for the coastal region of eastern Hudson Bay (Payette et al., 2004) and the HBL in northern Ontario (Kirkwood et al., 2019) where mean annual air temperatures (MAAT) are approximately 2 to 4°C warmer than Churchill. With continued climate warming, new ponds are likely to form on peat plateau surfaces and very small ponds enlarge in the northern Manitoba part of the HBL.
The carbon budget of Arctic peatlands is increasingly a research focus, given the anticipation that climate warming will enhance the role of these peatlands as a CH 4 and CO 2 source. A synthesis of findings by the IPCC (Meredith et al., 2019) suggests a potential soil carbon release from the entire northern circumpolar permafrost region in the order of 100 Gt by 2100, approximately 20% of this carbon store. However, the most recent studies incorporated in this synthesis show a possible overall uptake of carbon from increased peatland productivity. For the HBL, approximately 30 Gt of carbon has accumulated since the beginning of the Holocene (Packalen et al., 2014). Packalen et al. interpret GHG production from the HBL as having increased over the last 3,300 years, coincident with a marked decrease in peatland initiation and a generally wetter climate. Wet conditions and enlargement of fens in the Manitoba part of the HBL since about 2,200 years BP is also interpreted based on palynology (Kuhry, 2008). Perhaps the myriad larger lakes with peat shorelines in this part of the HBL formed during this time interval.
Despite the complexity of peatland carbon budgets, it seems likely that ponds and fens will play an increasing role as GHG sources under climate warming. Studies of the HBL indicate that shallow ponds and fens in peat produce considerably more GHGs than bogs or forest (Hamilton et al., 1994;Rouse et al., 1995Rouse et al., , 2002Hanis et al., 2013). Should ponds on peat plateaus initiate or expand, GHG generation is likely to increase, with further generation on conversion of peat plateaus to fen (Kirkwood et al. 2021). In addition, degradation of peat plateaus and their transformation to fens has ecological implications due to the disappearance of polar bear denning habitat and elimination of caribou forage (Richardson et al., 2005;Piercey-Normore, 2010).
In a previous paper, Dyke and Sladen (2010) proposed snowbank insulation along the edges of frozen peat plateaus in Wapusk National Park as a locus of peat plateau thawing that would be promoted by climate warming. Artificial snowbank creation is known to degrade permafrost (Hinkel and Hurd, 2006). Thermal equilibration under this kind of disturbance and under steady climatic conditions can take decades (O'Neill and Burn, 2016), suggesting that climate warming will accelerate thaw. Ponds may also serve the same function by reducing or eliminating the penetration of freezing into the ground in proportion to the latent heat extraction required for freezing the water column (Mackay, 1963). Dredge and Nixon (1992) hypothesized a mechanism for developing thermokarst lakes on peat plateau surfaces, based on expansion of water-filled depressions above thawing ice wedges. Roy-Léveillée and Burn (2017) show talik formation is possible beneath water depths as little as 20 cm adjacent to snow-trapping shoreline banks. Thus, shallow lakes, even if presently maintaining permafrost below, may act as additional sites for the initiation of thaw under climate warming. This, and accelerated thaw along plateau margins may hasten shrinkage of peat plateau extent that has already been detected (Dyke and Sladen, 2010). The objectives of this study are to analyze ponds as heat sources and revisit the role of plateau-margin snow accumulations in accelerating the thaw of peat plateaus under climate warming. We apply numerical modelling by establishing thermal boundary conditions that take into account the freezing and thawing of ponds and insulating quality of snow cover.

Terrain and ground temperatures
Our analysis applies to the continuous and discontinuous permafrost region of the HBL in northeastern Manitoba, between Churchill and the Nelson River (Fig. 1b). The surficial geology of this region is well described by Dredge and Nixon (1992), who note the almost ubiquitous cover of peat over Tyrell Sea sediments in this region. Dredge and Mott (2003) note the transition from fen-dominated to peat plateau-dominated wetland with distance inland from Hudson Bay. The ground thermal regime in this region is associated with peatland type. In addition to the review of ground thermal regime in Dyke and Sladen (2010), continued subsurface and ground surface temperature records from fen and peat plateau sites established between 2007 and 2009 in Wapusk National Park document the thermal regime associated with these two contrasting environments (Sladen, in prep.).
Frozen peat plateaus in the Manitoba HBL cover a belt 50 to 100 km wide between the coastal fen zone and forested bog inland (Dredge and Nixon, 1992;Fig. 1b). The plateaus are typically polygonalized with ice wedge troughs and exhibit peat thicknesses increasing inland from Hudson Bay to as much as 3 m (Dredge and Mott, 2003). The permafrost peat is ice-rich with visually estimated volumetric ice contents commonly greater than 50 % and ice-rich mineral soil immediately below. The plateaus stand up to 2 m above adjacent fens, affording exposure to wind that, with the low stature of tundra-like vegetation typical of the plateau surfaces, restricts snow depths to about 10 cm (Kershaw and McCulloch, 2007). This condition promotes cooling to produce the lowest mean annual ground temperatures (MAGT) in the region of -3 to -5˚C. Fens typically bound the peat plateaus and are often elongated, paralleling drainages. Slightly protected and often with standing water, fens are warmer with MAGTs of about -2˚C. The frozen peat plateau belt is within the continuous permafrost area that Dredge and Nixon (1992) map as occurring north of the tree line that extends diagonally across northern Manitoba to the estuary of Nelson River, (Fig. 1b). Modelling by Zhang et al. (2012) indicates a similar distribution.
Within Dredge and Nixon's (1992) continuous zone, patches lacking permafrost exist. Stands of spruce often coincide with the edges of plateaus, adding to the snowtrapping efficiency of the banks at plateau edges. MAGTs average about 1˚C adjacent to a peat plateau embankment prone to snow accumulations of up to 2 m, and never drop below 0˚C to a depth of at least 5 m (Fletcher Lake site on Fig. 1b, Fig. 2). A transition to discontinuous permafrost in the southern half of this region is also indicated by fen MAGTs of about 3˚C near York Factory, just to the south of Wapusk National Park (Sladen et al., 2009).
Ponds and lakes are myriad in the northern Manitoba part of the HBL, probably numbering into the hundreds of thousands (Sannel and Kuhry, 2011) and ranging in size from ice-wedge trough filings to maximum lengths of about 2 km. Depths of larger lakes range up to about 2 m in the northern-most part of this region (Duguay and Lafleur, 2003). Dredge and Nixon (1979) distinguish between ponds and lakes. Ponds are metres to several tens of metres across, entirely within peat, less than 1 m deep and typically less than 0.5 m deep (Symons et al., 2014). Lakes are several hundred metres to greater than a kilometre in length, up to 2-3 m deep, with mineral sediment shores and peat banks up to 3 m high. The largest and deepest lakes undoubtedly have through-permafrost taliks, given that ice thicknesses are unlikely to exceed 2 m (Duguay and Lafleur, 2003). We have only one lake temperature record in the region, a 3-year record from below the depth of ice growth. The mean annual temperature for this record averages 4.0˚C, which supports the existence of through-permafrost taliks for lakes. Ponds, however, are shallow enough that taliks may exist but do not penetrate through permafrost. Lack of taliks beneath ponds is supported by winter electrical conductivity transects that show the lowest conductivities over ponds, interpreted as absence of any thawed ground (Dyke and Sladen, 2010).
The water bodies of interest in this study are the ponds that dot the surfaces of peat plateaus (Fig. 3). Ponds filling the depressions in the centers of ice-wedge polygons are often present in addition to those in ice-wedge troughs. The merging of these central ponds with ice-wedge troughs may begin Dredge and Nixon's (1979) process of pond amalgamation.

Climate
Air temperature records have been kept at Churchill, Manitoba since 1943 (Environment Canada, 2020). Environment Canada 30-year temperature normals are available starting with 1941-1970, ending with 1981-2010. For these five intervals they are, successively, -7.3, -7.2, -7.1, -6.9, and -6.5˚C, showing a gradual although slightly accelerating increase. A linear regression on the yearly MAAT ( Fig. 4a) shows an overall warming of 1.9˚C or 0.024˚C /yr, although there is considerable inter-annual variation (MAAT extremes are -9.8 to -3.2˚C from 1943 to 2020) and a general decrease in annual means until the mid-1960s. Following this decrease, overall warming averages 0.040˚C /yr from 1970 through 2020.
Year to year variability is particularly large since about 1995, with three years having MAATs between -3 and -4˚C.
If the Churchill climate data are replotted as yearly freezing and thawing degree-day indices (FDD and TDD), the variability since 1995 remains apparent, but only in the freezing indices (Fig. 4b). The coldest post-1995 freezing indices maintain much the same decrease as over the entire period of record but the post-1995 warmest freezing indices show a distinct jump to warmer values. Abrupt paleolimnological changes in the HBL and unprecedented fish die-offs since the early 1990s coincide with this onset of unusually warm freezing indices (Gunn and Snucins, 2010,  . Examples of small ponds on peat plateaus; a) ponds on the peat plateau hosting the Fletcher Lake peat plateau thermistor cable, b) ponds occupying the centers of ice-wedge polygons on a nearby plateau. Rühland et al., 2013). Rühland et al. (2013) suggest that diminishing Hudson Bay ice cover, which has historically retarded the warming otherwise prevalent in the Arctic, is allowing accelerated warming in bordering lands.

Peat thermal properties
Under climate warming, peat plateaus and raised bogs are the terrain features most likely to retain permafrost (Halsey et al., 1995). The wide range in volumetric moisture contents promoted by the fibrous, organic nature of peat enables a seasonal, permafrost-preserving switch of thermal conductivity. Drainage and drying of this elevated peat insulate permafrost in summer. Increasing wetness in fall results in a several-fold increase in thermal conductivity that promotes winter cooling (Brown, 1966). This change is responsible for peat plateaus retaining permafrost in the southern fringe of the discontinuous zone (Zoltai, 1972) and exhibiting the coolest mean annual ground temperatures in the study region.
Thermal conductivity measurements demonstrate the pronounced dependency of conductivity on water content (Farouki, 1981;Kujala et al., 2008). Because of the porous nature of peat, saturated volumetric water contents can approach 100% with thawed or frozen thermal conductivities approaching those of pure water or ice (0.605 and 2.23 W/ m˚C, respectively, Geoslope, 2014). The thermal conductivity of saturated unfrozen peat approaches about 0.5 W/m˚C but drops off steadily with drying to values of less than 0.1 W/ m˚C for volumetric water contents of less than about 20% (O'Donnell et al., 2009). Saturated frozen values range from 1.2-1.8 W/m˚C, diminishing to less than 0.2 W/m˚C for volumetric water contents of less than about 20%.
Volumetric water content below the water table for saturated peat in peat plateaus is 80-90% typically (Zoltai and Tarnocai, 1971;1975). However, most of the insulating value of peat is contained above the water table, where water contents can fall to around 10% or less. This seasonal variability can produce ratios of frozen to unfrozen thermal conductivity of 4 or more (Shur and Jorgensen, 2007), maximizing the phenomenon that gives peat its permafrostretaining character. extremes of 1526 and 1170 TDD, respectively. The peaks and troughs for each index are characterized as a sawtooth distribution, shown on the left half of Figure 5a, with the climate warming trend, discussed below, added.
To simulate climate warming, a linear warming trend predicted by global circulation models (GCMs) is added to the freeze and thaw index sawtooth. Although the averaged 1970-2020 warming trend could simply be extended to predict future air temperatures, GCM development over the last few decades offer a more deterministically grounded basis for anticipating future climate. The Scenarios Network for Alaska and Arctic Planning (SNAP) forecast employs 5 GCMs, including that developed by the Canadian Center for Climate Modelling and Analysis (SNAP, 2018). Forecasts for specific locations are produced by downscaling, a statistical interpolation of GCM predictions. The SNAP forecast for Churchill for the time interval 2040-2049 is a MAAT of -4.0˚C using their low range representative concentration pathway (RCP) of RCP 4.5. We selected this optimistic forecast expecting that it would provide a likely minimum climate warming impact.
The MAAT for the latest Environment Canada 30-year climate normal (1981-2010) is -6.5˚C. Therefore, from the middle of this climate normal interval (1995) to 2045 gives an average MAAT increase of 0.050˚C/yr. This rate of increase is 0.010˚C/yr higher than the 1970-2020 average increase and is used in our simulation. In terms of temperature indices, the freezing index decreases from 3500 FDDs in 1995 to 2918 FDDs in 2045. The thawing index over the same interval increases from 1177 TDDs to 1477 TDDs (Fig. 5a). The yearly increment of this trend is added to each yearly freeze or thaw index on the sawtooth extrapolation. The warming sawtooth trends, extended to 2050, are plotted on the right half of Figure 5a. Lastly, to provide a continuous air temperature boundary condition for the modelling, mean monthly temperatures from the 1981-2010 climate normal are extrapolated in proportion to the indices plotted in Figure 5a to give the air temperature projection used for all boundary conditions (Fig. 5b).
Boundary conditions must be specified for the peat plateau surface, the bottom of a pond, and under the snow accumulation in the fen adjacent to a plateau edge. For the subaerial parts of the simulation (plateau surface and adjacent fen), n-factors (Lunardini, 1981) are used. N-factors are an empirical means of accounting for environmental conditions that determine the temperature at the ground surface, the most significant of which is likely to be snow cover. Because thermal properties of moist ground change significantly between the thawed and frozen condition, n-factors are designated according to the moisture phase state. The thawing n-factor (n t ) is the ratio of the air thawing index to the ground surface thawing index. The freezing n-factor (n f ) corresponds to freezing conditions. We calculate n-factors using our ground surface temperature records for a peat plateau, a plateau-edge snow accumulation and a fen, together with an air temperature record in the vicinity of these installations (Sladen, in prep.).

Modelling approach
We use the commercially available temperaturemodelling software TEMP/W (Geoslope, 2014) to predict temperatures below the surface of a frozen peat plateau, as modified by a pond or a snow bank against the side of a plateau. The air temperature boundary condition is supplied by projecting the 1995 to 2020 accentuated air temperature variation, from 2020 to 2050. However, instead of applying MAATs, we apply successive freezing and thawing indices, a characterization that captures the observed post-1995 variation seen in Figure 4b. To define the variation, the freezing and thawing index extremes over the 1995-2020 time interval are selected and averaged to define the peaks and troughs of a cycle for both indices. The time intervals between the selected high and low values are also averaged to give the repeat interval of the cycle. Freezing index maxima or minima occur on average every 6 years, the high extreme being 3984 FDD and low extreme 2820 FDD. Thawing indices show less variability, with maxima or minima occurring every 4 years with high and low Our lake temperature record is in water too deep to represent the thermal regime at the bottom of a shallow plateau pond. The TEMP/W software treats a water layer as another solid, unrealistic considering that mixing is likely to be prominent in shallow water bodies. Instead, we synthesize the pond-bottom thermal regime. A pond depth of 30 cm is assumed, this shallow depth implying that a large portion of the pond-bottom thermal regime will be determined by cooling through relatively thin ice in contact with the pond bottom. Ice-free temperatures on the pond bottom are expected to closely follow air temperatures, given the active mixing likely to occur in such a shallow water body. Water temperatures measured for ponds and small lakes across Wapusk National Park (Symons et al., 2012(Symons et al., , 2014) support this contention. An average temperature of 17.1˚C was measured for 42 lakes and ponds, during which time the mean air temperature averaged 15.4˚C. The following year 21 ponds and lakes averaged 12.1˚C with a mean air temperature during the measurements of 12.4˚C. Our one lake temperature record Fig. 6. The four pond-bottom temperature boundary conditions used for the forward modelling, derived from the air temperatures in Figure 5b, depending on the pond ice snow regime and the snow thermal conductivity. Only the top graph includes the above-freezing part of the boundary condition as it is the same for all snow conditions. shows little of the hysteresis with air temperature normally observed in lakes (Piccolroaz et al., 2013), confirming that shallow pond water will match air temperature closely.
To specify the pond-bottom temperature when ice is in contact with the bottom, the temperature required at the pond bottom to extend freezing to the depth that is predicted is used. This temperature and the freezing time remaining after freezing reaches the pond bottom provides the index for freezing into the pond bottom. Two modes of pond freezing are analyzed, 1) freezing occurs with a constant snow depth established from the onset of pond freeze and 2) a constant snow depth is established once pond freezing reaches the pond bottom. These two modes reflect the character of snow accumulation on water bodies in the vicinity of Churchill (Brown and Duguay, 2011). Snow depths for our simulation are assumed to instantaneously reach 30 cm, on the expectation that following snowfall, windblown snow quickly fills in the recess of a pond surface below the plateau surface.
Yearly air temperature cycles as specified in Figure  5b are applied starting with the onset of freezing. For the plateau surface and the area covered by the snow bank, appropriate n-factors are used to determine the ground surface temperatures. For the pond bottom, once freezing has started, the temperature is assumed to remain at 0˚C until pond ice is calculated to reach the pond bottom. Then the freezing index remaining at the pond bottom is applied over the remaining freeze interval. Thawing the pond ice again maintains the pond bottom temperature at 0˚C until the pond ice thawing is complete. Thawing the 30 cm of pond ice requires 70 TDDs, determined using an empirically derived heat transfer coefficient (Ashton, 1983). The remaining air thawing index determines the pond water temperature until the commencement of the next freezing interval. The pond-bottom thermal regimes for the two freezing modes under two snow thermal conductivities (as selected below) are shown in Figure 6.
Two types of modelling are performed: 1) a hindcast to ascertain which of the two pond freezing modes and two snow thermal conductivities (selected as described below) result in the observed field condition of no sub-pond talik and 2) forward modelling that applies the SNAP-based climate warming. The hindcast is performed using the Churchill air temperature record for the interval 2008-2018. It also supplies the initial condition for the forward modelling, with the yearly freezing and thawing cycles, as specified in Figure 5b, then applied through 2050. Two kinds of cross section are investigated (Fig. 7): 1) a peat plateau with ponds 10 m, 20 m and 40 m in size, Figure 7a, and 2) a peat plateau-fen margin with adjacent snow bank, Figure 7b. Models with one and two ponds are analyzed, the two-pond version to examine if closely adjacent ponds can accelerate thawing of the intervening plateau surface. Peat thickness is set at 1.5 m, as measured at the Fletcher Lake site.

Material thermal properties
Thermal conductivities for plateau peat are estimated by using TEMP/W to find values that reproduce the 2008-2011 record from a 4.5 m-deep thermistor cable at the Fletcher Lake peat plateau site (location on Fig. 1b). We start with our n-factors determined for plateau sites and determine peat conductivities by finding values that reproduce peat plateau active-layer depths. N-factor values were determined from air-ground surface temperature records available for 3 to 5 years for peat plateau, fen, and snowbank locations (Table 1). Back-calculating with four years of ground surface temperature records gives consistent peat conductivity values for thawing to the active layer, averaging 0.26 W/m˚C. We obtain higher thawed conductivity values for thawing to depths less than the active layer, ranging between 0.50 and 0.87 W/m˚C, the highest values corresponding to the shallowest thaw depths.  Table 2 for thermal properties). Pond sizes of 10, 20, and 40 m are modelled. b) A fen adjoining a peat plateau. Fen peat is 30 cm thick. N-factor zones simulating decreasing snow depths away from the peat plateau edge are indicated. All material properties and conditions are listed in Table 2. Both models include a geothermal heat flux of 0.05 W/m 2 . c) Comparison of the Fletcher Lake site peat plateau ground temperature record with modelling. The range of annual extremes, 2008-2018 (shown in grey) and the mean for this time interval (white dashed line) are compared with the modelled mean annual extremes for different freezing n-factors. The thawed n-factor used is 0.90.  Figure 7c, along with the modeled mean minimum and maximum temperatures for the same time interval, using the selected peat thermal conductivities of k t = 0.30 W/m˚C and k f = 1.3 W/ m˚C. Modelling does not match both the minimum and maximum means using only a single freezing and thawing n-factor and conductivity. Rather, varying fits can be obtained with different freezing n-factors, as shown in Figure 7c for n f = 0.50, 0.67 and 0.90. We choose to have the best overall fit for the mean maximum and minimum temperatures. Therefore, for modelling, n f = 0.67 and n t = 0.90 produce satisfactory correlations with the mean minimum and maximum temperatures. Peat beneath ponds is assumed to be always saturated and is assigned thawed and frozen thermal conductivities of 0.50 and 1.50 W/m˚C, respectively (Table 2).
Snow thermal conductivities are based on the regression analysis of Sturm et al. (1997), calculated using snow densities of Kershaw and McCulloch (2007) and selected based on experience of the authors. We select 0.25 W/m˚C to represent wind-packed snow and 0.40 W/m˚C to represent snow that has undergone thaw. The other necessary parameters, i.e. material volumetric water contents, pond bottom peat and mineral soil thermal conductivities, are listed in Table 2.

Hindcast modelling
Only the pond models using the snow conductivity of 0.40 W/m˚C, with either of the snow cover conditions, show little or no talik development beneath the pond bottom. These selections are the only ones consistent with the evidence of no pond taliks from the geophysical surveys over ponds in early spring (Dyke and Sladen, 2010). Using the snow conductivity of 0.25 W/m˚C results in additional thaw of 1 m or 2.5 m for pond freezing without and with snow cover, respectively (Fig. 8). These results suggest that the higher snow thermal conductivity is the more appropriate value to use for extending the modelling forward in time.
The plateau margin with snow bank exhibits a talik under that segment of the adjacent fen assigned n f = 0.07. The talik is trough-shaped, thinning abruptly towards the next fen segment assigned n f = 0.18 (Fig. 9a). A temperature vs. depth plot through this talik compares well with the average of four years of observations from the Fletcher Lake snow bank thermistor cable (Fig. 9b). The hindcast model is initiated with the temperature distribution arrived at after 30 years of applying the 1981-2010 climate normal, at which time the talik extends to a maximum depth of 12.5 m. After the 11-year hindcast, the depth has increased only to 12.7 m. This slight increase suggests that the plateau   margin talik is slightly out of equilibrium and deepening with the present climate but that a minor increase in n f would result in talik shrinkage. Predicted peat plateau thaw depths range between 0.6 and 0.9 m during the 11-year modelling interval, consistent with active-layer depths recorded at the Fletcher Lake site.

Forward modelling
Plateau ponds: The most pronounced modeled change through 2050 for all pond sizes, freeze-back conditions and snow thermal conductivities is the increase in thaw depth beneath ponds. Figures 10a and b show the sub-pond talik development for the two onset of snow conditions. The 40 m pond size, early accumulation of snow and low snow thermal conductivity all result in the deepest talik development of 6.3 m below the pond bottom (Fig. 10b top). For the 40 m pond with complete pond freeze before snow accumulation and the high snow conductivity, thawing extends at most only 1.3 m below the pond bottom and no talik forms (Fig. 10a bottom).
Maximum thaw depths by 2050 beneath ponds for all pond sizes and both freeze-back conditions are shown in Figure 11a. Noteworthy is the constant and minimal thaw depth for all pond sizes for no snow accumulation during pond freeze and 0.40 W/m˚C snow conductivity. All other pond sizes and snow conditions result in thaw extending into mineral soil, implying talik development.
The progression of thaw depth for 40 m ponds is shown in Figure 11b. Smaller pond sizes are not shown as the thaw progression is essentially proportional to the 2050 depths shown in Figure 11a. For no snow accumulation during pond freeze and 0.40 W/m˚C, there is essentially no thaw progression beyond the depths reached during the hindcast interval. For snow on the pond ice as the ice forms, maximum thaw begins at the same depth but increases markedly through 2050. The same two snow accumulation conditions for the low snow conductivity also result in marked thaw depth increases through 2050. For the peat plateau surface, there is no complete thaw of the peat thickness but thaw depths do increase to a maximum of 1.1 m (not shown). If the thawed plateau peat conductivity is increased to 0.5 W/m˚C, as might be the case with an increasingly moist active layer, then thawing extends almost to the base of the peat layer.
Talik development beneath ponds implies that, given the probable ice-rich nature of the deeper peat and top-most mineral sediments beneath the peat, subsidence and pond deepening may accompany climate warming. An arbitrary pond deepening to 60 cm over the interval 2020-2050 is simulated by increasing the pond depth in increments proportional to the yearly thaw depth increase shown in Figure 11b for the case of pond freezing before the onset for ponds with the pond freezing during snow cover. Isotherms are labeled in degrees Celsius. Fig. 11. a) Maximum thaw depths for all pond sizes, snow thermal conductivities, and pond freezing regimes in 2050. Grey dots represent ponds freezing during snow cover, black dots, ponds freezing before snow cover. Power law regression curves have been fitted to each pond regime (dotted lines); note that the curves converge to the peat plateau thaw depth, as would be expected as ponds reach zero size. b) Progression of thaw for 40 m ponds, 2020-2050. Progression for a pond freezing during snow cover is shown in grey, for a pond freezing before snow cover in black of the higher conductivity snow cover. With thawing in the remaining frozen peat, the assumed 80% volumetric water content for pond-bottom peat would result in roughly 30 cm of thaw subsidence. The pond bottom temperature boundary condition and resulting thaw depth increase are shown in Figure 12. Despite the almost complete elimination of pond-bottom freezing by 2050, only a modest increase in maximum thaw depth occurs.
Peat plateau edge: Forward modelling results in the plateau edge talik deepening 1.0 m by 2050 (Fig. 13a). We do not know how far the snow bank extends from the edge of the peat plateau. If the associated model segment with n f = 0.07 was extended further from the edge of the peat plateau, the talik would be deeper. However, the modeled talik depth at the end of the hindcast time interval is consistent with the snow bank thermistor record (Fig. 9b) so, to this distance at least, the modelled talik extent is probably accurate. There is no talik initiation beneath adjacent segments of the adjoining fen with assigned n f = 0.18 and 0.28.
Terrain between two plateau ponds: Only the modelling result for the 10 m width of peat plateau between two 40 m ponds is presented (Fig. 13b). No permafrost degradation in the intervening width of peat plateau is produced, irrespective of snow conductivity and accumulation mode. Thaw depths are essentially the same as for peat plateau extending away from single ponds. Results for smaller ponds are not presented as they produce correspondingly less thermal interference.

DISCUSSION
Our modelling results suggest that, under the climate warming projected for Churchill, peat plateaus in northern Manitoba will remain intact through 2050. However, under plateau ponds, taliks initiate and deepen, depending on the timing of snow cover relative to pond freezing and on the snow thermal conductivity. Talik initiation and deepening is promoted by low snow thermal conductivities due to the reduced subsurface cooling associated with increased snow insulation. Deepening is also promoted by early onset of snow on pond ice. On the other hand, peat plateau surfaces show no talik development, even when thawed peat thermal conductivity is increased to simulate peat saturated throughout the summer.
The accuracy of these results depends on the accuracy of the SNAP climate change prediction. Our selection of RCP 4.5 implies an optimistic outlook for the reduction of greenhouse gas output whereby emissions peak around 2040, then decline to about half that value by 2100. The 30-year climate normal for our projected 2021-2050 air temperature regime is -4.5˚C. This value represents an accelerated increase from the last reported value of -6.5˚C for 1981-2010. If we note 1981-2010 climate normals plotted for northern Manitoba (Fig. 1c), a 2˚C increase in all locations shifts both Gilliam and Thompson to well above the -3.5˚C MAAT that Halsey at al. (1995) have determined to be the maximum temperature for the maintenance of frozen bog. Disappearance of permafrost at bog edges is presently occurring at Thompson (L. Dyke, personal observation, 2019), suggesting that this phenomenon would be occurring north of Gillam by 2050. In fact, signs of peat plateau retreat are documented for southern Wapusk Park (Dyke and Sladen, 2010).  Accuracy of the talik prediction also depends on the appropriate choice of snow thermal properties and sequence of pond freezing. We feel that the forward modelling based on complete pond freezing before snow cover and the higher snow thermal conductivity offers the most accurate climate warming prediction. This is because the initial observed condition of little or no talik is best reproduced by these conditions.
Although there is uncertainty regarding the response of frozen ground to climate warming in northern Manitoba, the results of this analysis offer guidance for continued observation. Peat plateaus are forecast to undergo increased summer thaw but not the initiation of taliks. However, this increased thaw may be accompanied by subsidence enough to induce shallow ponding or expand existing ponds on plateau surfaces. This expansion may lead to the lateral extension of taliks already forming under existing ponds. Regardless, pond expansion will be a sign that peat plateau degradation is underway.
To the extent that comparison is possible, our results are consistent with other climate warming projections for arctic ponds and peat regions. Langer et al. (2016) predict thaw depths beneath tundra ponds of similar size to ours for RCP 4.5 but for a site in northern Siberia with a present MAAT of -13˚C. For 2050 their thaw depth beneath a 0.3 m depth pond is slightly over 1 m and by 2100 is slightly under 2 m. For a tundra surface, thawing to approximately 0.5 m occurs until 2070, at which time thawing encounters ice-rich sediment and accelerates as ponding develops. Liljedahl et al. (2016) describe ice wedge degradation and the deepening of ice wedge troughs over the last few decades for several high arctic sites. Despite a warmer MAAT for our area, trough deepening does not appear to be happening, presumably because of the insulating quality of the peat. Zhang (2013) models peatland thaw specifically for Wapusk National Park through 2200 (but not beneath ponds) using an intermediate CO 2 emissions projection. Zhang's (2013) thicker peat site (0.24 m) is considerably thinner than our selection. By 2050 a roughly 1.5 m thaw depth has increased very slightly but markedly accelerates with talik formation at the end of the century. Although these studies are not directly comparable, both do suggest that continued warming beyond 2050 will result in accelerated pond thaw and subsidence followed closely by peat plateau thaw and subsidence.
Retreat of peat plateau margins under the applied climate warming is not forecast by the modelling. Although marginal ponds, as reported in Dyke and Sladen (2010), are not included in the modelling, their presence is unlikely to affect the modeled talik development, simply because the thick snow cover dominates the winter thermal regime. Other environmental factors, particularly tree growth along plateau margins, which extends deep snow onto the plateau surface, will ultimately contribute to plateau margin subsidence. In Northwest Territories, forested peat plateaus margins are retreating under a present MAAT of about -3˚C (Baltzer et al., 2014). In northern Norway frozen peat plateau margins are retreating under a MAAT of -2˚C for the interval 1967-2019 with yearly MAATs of as high as -0.1˚C (Martin et al., 2021). Under this highest MAAT, Martin et al. (2021) forecast plateau margins to retreat roughly 10 m over the next century with very thin snow cover but essentially collapse with a 20-30 cm snow cover. Palsas at four locations in the same general part of Norway have shrunk under climate warming that has raised MAATs from a range of -2.9 to 0.1˚C  to a range of -1.8 to 1. 1˚C (1991-2014) (Borge et al., 2017). These observations support the expectation that as northern Manitoba temperatures rise toward the same range, peat plateau retreat and collapse will ensue.

CONCLUSIONS
Frozen peat plateaus in northern Manitoba are maintaining their frozen state under the present climate. Climate warming is predicted to continue at an accelerated rate with the 30-year climate normal in 2050 rising 2˚C from the climate normal for 2010. Individual MAATs will have risen to as high as -2.3˚C if the interannual variability characterizing the last 20 years at Churchill continues. Although peat plateaus remain frozen under this climate trend through 2050, degradation will have commenced in the form of talik initiation beneath plateau ponds. Any warming additional to that prescribed in the modelling will induce additional degradation which should be readily apparent as expansion of plateau ponds.