Thermal erosion patterns of permafrost peat plateaus in northern Norway

Subarctic peatlands underlain by permafrost contain significant amounts of organic carbon and our 2 ability to quantify the evolution of such permafrost landscapes in numerical models is critical to provide 3 robust predictions of the environmental and climatic changes to come. Yet, the accuracy of large-scale 4 predictions is so far hampered by small-scale physical processes that create a high spatial variability of 5 surface ground thermal regime and thus of permafrost degradation patterns. In this regard, a better 6 understanding of the small-scale interplay between microtopography and lateral fluxes of heat, water and 7 snow can be achieved by field monitoring and process-based numerical modeling. 8 Here, we quantify the topographic changes of the Šuoššjávri peat plateau (Northern Norway) over a 9 three-years period using repeated drone-based high-resolution photogrammetry. Our results show that edge 10 degradation is the main process through which thermal erosion occurs and represents about 80 % of measured 11 subsidence, while most of the inner plateau surface exhibits no detectable subsidence. Based on detailed 12 investigation of eight zones of the plateau edge, we show that this edge degradation corresponds to a 13 volumetric loss of 0.13 ± 0.07 m 3 yr -1 m -1 (cubic meter per year and per meter of plateau circumference). 14 Using the CryoGrid land surface model, we show that these degradation patterns can be reproduced 15 in a modeling framework that implements lateral redistribution of snow, subsurface water and heat, as well 16 as ground subsidence due to melting of excess ice. We reproduce prolonged climate-driven edge degradation 17 that is consistent with field observations and present a sensitivity test of the plateau degradation on snow 18 depth over the plateau. Small snow depth variations (from 0 to 30 cm) result in highly different degradation 19 behavior, from stability to fast degradation. These results represent a new step in the modeling of climate-driven landscape development and 21 permafrost degradation in highly heterogeneous landscapes such as peat plateaus. Our approach provides a 22 physically based quantification of permafrost thaw with a new level of realism, notably regarding feedback 23 mechanisms between the dynamical topography and the lateral fluxes through which a small modification of the snow depth results in a major modification of the permafrost degradation intensity. In this regard, these 25 results also highlight the major control of snow pack characteristics on the ground thermal regime and the 26 benefit of improving snow representation in numerical models for permafrost degradation projections. assuming rotational rather than translational symmetry. This suggests that our simulations are indeed most realistic during the CED phase, whereas changes of the peat plateau a maximum of 5-10 cm on the plateau only triggers a 4-5 meters edge retreat within 100 years, a 10-20 cm cover fully degrades the plateau in 100 years and a 20-30 cm cover degrades it in 40 years. Our simulations reproduce the observed lateral edge degradation with a stable plateau top, but the final phase 559 of plateau degradation corresponds to complete plateau collapse with subsidence occurring throughout 560 the entire plateau. These results highlight the fast and high spatial variability of permafrost landscape evolution in response to climate change. They also show that the related microtopographic, thermal and hydrological 563 modifications can be represented in numerical models, opening the way for substantial improvements in simulating permafrost landscape evolution and its impact on greenhouse gas emissions from thawing permafrost.


29
Observations show that permafrost is warming at a global scale (Biskaborn et al., 2019). Its 30 thawing has major consequences on Arctic and Boreal ecosystems and landscapes (Beck et al., 2015;31 Farquharson et al., 2019;Liljedahl et al., 2016) and potentially represents an important climate feedback 32 through the degradation of thawed organic matter Schuur et al., 2009Schuur et al., , 2015. Carbon 33 emissions from permafrost regions towards the atmosphere are already observed (Natali et al., 2019); 34 field measurements show that these emissions are influenced by the timing of the active layer deepening 35 (Morgalev et al., 2017) and by the state of degradation of the permafrost terrains (Langer et al., 2015;36 Nwaishi et al., 2020;Serikova et al., 2018). Particularly, abrupt thaw of ice-rich permafrost is expected 37 to become a dominant process regarding carbon emissions, offsetting the potential ecosystemic negative 38 feedback that is expected for gradual thaw (McGuire et al., 2018;Turetsky et al., 2020). 39 As such, our ability to quantify and represent the physical evolution of permafrost landscapes is 40 critical to provide robust predictions of the environmental and climatic changes to come (Aas et al.,41 2019; Andresen et al., 2020;Teufel and Sushama, 2019). While permafrost affects about 14 millions 42 square kilometers throughout the Northern Hemisphere (Obu et al., 2019), the ground thermal response 43 to climatic signal and morphological changes of permafrost are governed by processes occurring within 44 a spatial scale of few meters (Gisnås et al., 2014;Jones et al., 2016;Martin et al., 2019;Way et al., 45 2018). Indeed, at the small-scale, Arctic and Boreal permafrost landscapes in low lands (such as peat 46 plateaus and polygonal tundra) are characterized by low amplitude (0-3 m vertically) and high frequency 47 artificial ground control points were acquired with a differential GPS (dGPS) to support the topographic 124 reconstruction. Images were processed, including the ground control points, using the Agisoft Photoscan 125 software (version 1.2.6). The final digital elevation models have a grid resolution of 0.1 m. To guarantee 126 a meaningful subsidence signal, we only considered subsidence values exceeding 5 cm in this study. 127 To frame our modeling approach with field constraints, we measured snow depth on the 128 Šuoššjávri peat plateau. Snow measurements for winter 2016 were presented in Martin et al. (2019). 129 Here we extend these measurements to winters 2017 and 2018 for the plateau top part only. Snow depth 130 was measured in the end of March with an avalanche probe at precise locations covering the top of the 131 plateau, as defined in Martin et al. (2019). The exact location of each measurement point under the snow 132 was localized using a dGPS system to ensure a 5-10 cm horizontal accuracy. 133

Quantification of thermal erosion patterns 134
Our drone-based photogrammetric approach (Sect. 3.1) enabled us to derive repeat DEMs of the 135 peat plateau for September 2015 and 2018. We computed the elevation difference (elevation from 2018 136 minus elevation from 2015) to quantify the spatial pattern of elevation changes of the plateau. From this 137 elevation difference map, we selected eight degrading zones (later referred as edge transect areas) 138 presenting a 10-30 m long and roughly straight lateral extent for comparison with modeling results (Sect. and per year (m 3 m -1 yr -1 ). Because elevation changes occurred in the mire due to water level variations 142 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

8/33
between the two dates, we relied on an estimation of the elevation of the plateau edge inflection point 143 (around 309.7-309.8 m asl, yellow color on Fig. 2B) to delimitate the plateau from the mire and thus 144 identify elevation changes associated to the plateau only. 145 Permafrost degradation creates different subsidence patterns depending on the size of the 146 permafrost landform. Small structures like palsas tend to sink entirely from the edge to the top, while 147 peat plateaus show stability of their top part and pronounced lateral retreat. To distinguish between these 148 two types of thermal erosion patterns, we introduce a so-called Horizontal vs Vertical (HvsV) shape 149 index that we can apply to both field observations and model results. The basic concept of the HvsV 150 index is illustrated in Fig. 3. To compute the HvsV index, the plateau edge elevation is averaged over 151 five points, from its very base (z1) to the first meters where its flat top is reached (z5). The elevation 152 difference between 2015 and 2018 (Δz) for z1, z2, z4 and z5 is then used as follow: 153 For the field observations, the HvsV index was obtained by first laterally averaging the slope of 155 each edge transect area using five to fourteen parallel elevation profiles across the zone, for the 2015 156 and 2018 DEMs. These two synthetic elevation profiles were then averaged into the five required points 157 (from the edge to the top of the plateau). The HvsV index was then calculated based on these synthetic 158 elevation differences (between 2018 and 2015). 159 For the simulation results (Sect. 4.2), a 10 m long window was used to capture the topography 160 from the base of the plateau to its flat top. These 10 m profiles where then averaged into the five required 161 points. Then, elevation differences were observed over three-year-long time periods to compute the 162 HvsV index.   169 We simulate the ground thermal regime and related topographic evolution of the Šuoššjávri peat 170 plateau using the CryoGrid model . CryoGrid is a land surface model designed 171 for permafrost modeling. It consists of a physically based description of 1D heat transfer in the soil 172 column, including freeze-thaw process of soil water/ice. The model features a simple snowpack module, 173 which includes heat conduction, dynamic buildup, melt, sublimation, water infiltration, and refreezing. 174 CryoGrid uses a surface energy balance module to calculate the ground surface temperature. The 175 turbulent fluxes of sensible and latent heat are calculated using a Monin -Obukhov approach (Monin 176 and Obukhov, 1954). Evapotranspiration is adjusted to soil moisture and the water uptake is distributed 177 vertically so that it decreases exponentially with depth. The soil moisture computation along the soil 178 column relies on a bucket scheme . 179 CryoGrid represents ground subsidence resulting from the melt of the excess ice in the ground 180 (thermokarst process, see Westermann et al., 2016). Subsidence calculation is based on soil stratigraphy 181 (ice content and natural porosity) and modifies the 1D vertical soil mesh when excess ice (soil grids for 182 which the ice content exceeds the natural soil porosity) melts. This functionality was first implemented 183 in Westermann et al. (2016) and later used in  to represent the transient 184 evolution of polygonal tundra landscapes caused by permafrost degradation for various hydrological 185 scenarios. In the present study, this scheme is used to account for the thermal erosion of the peat plateaus 186 and the microtopography changes induced by permafrost degradation within the plateau. 187 Following Nitzbon et al. (2019Nitzbon et al. ( , 2020, CryoGrid includes a parallel framework to 188 simultaneously compute several 1D tiles that can exchange numerical information at defined time steps 189 (quantities of water, snow and heat). This approach, denoted as laterally coupled tiling, allows us to 190 couple 1D tiles with different stratigraphies/topographies. With this method, the spatial variability tied 191 to Arctic landscape such as the topographic gradient between the center and the rim of polygonal tundra 192  or the stratigraphy differences between Yedoma and Holocene deposits (Nitzbon 193 et al., 2020) has been simulated. It also permits to calculate the lateral fluxes of snow, subsurface water 194 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

10/33
and heat that are small-scale key drivers of the ground thermal regime . In the 195 present study, laterally coupled tiling is used to simulate the Šuoššjávri peat plateau and its interface 196 with the surrounding wet mire and represent the subsidence and lateral fluxes along this transition (Sect. 197 3.2.4). 198

Modelling framework and sensitivity analysis to plateau snow depth
199 Snow depth is a major parameter influencing the ground thermal regime (cf. Introduction) and 200 assessing its influence on the degradation dynamics of peat plateau is one of the objectives of this study. 201 Because of microtopography-dependent wind-driven snow redistribution, snow depths over the peat 202 plateaus are usually shallow (10-40 cm, Martin et al. 2019). In the Cryogrid model, the efficiency of this 203 redistribution is adjusted with a parameter that sets, for each tile, the maximum snow depth which cannot 204 be transported towards the mire. Snow accumulated beyond this depth is systematically transported by 205 the model towards the lower-lying tiles, here representing the mire. 206 The present-day topography of the peat plateau is clear evidence that the long-term evolution is 207 not only governed by thermal erosion of the plateau edges, but also by a complex interplay of 208 thermokarst process in the interior of the plateau, with pond formation and drainage, as well as drainage 209 gully development and deepening. As a full three-dimensional simulation of these phenomena is beyond 210 the capability of present-day models, we focus modeling on the simplified situation of a laterally 211 symmetric peat plateau edge. Furthermore, we conduct simulations for steady-state climate forcing from 212 the period of the field observations (Sect. 3.3.3). This makes it possible to compare the magnitudes of 213 modeled volumetric plateau degradation over three-year periods with field observations for sufficiently 214 straight sections of the plateau edge (Sect. 3.2, Fig. 2). As field observations of snow depth reveal a 215 considerable spread of snow depths values on the plateau (that cannot be reproduced by modeling), we 216 investigate the sensitivity of modeled topography changes towards snow depths on the plateau (which 217 can be controlled by the above mentioned model parameter), using a realistic range of values between 218 which snow depths vary during the simulations (Sect. 3.1). We therefore designed four simulations 219 called: 0 cm snow, 5-10 cm snow, 10-20 cm snow and 20-30 cm snow. 220 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

Model setup 221
To simulate the edge thermal erosion of the Šuoššjávri peat plateau, we used the laterally 222 coupled tiling framework of CryoGrid  with 40 coupled tiles in a linear 223 configuration so that they can exchange snow, subsurface water and heat (Sect. 3.2.1). In the third 224 dimension, translational symmetry is assumed in order to represent the evolution of the 10 to 30 meters 225 long and roughly straight geometry of these specific zones of the plateau. These edge transect areas are 226 detailed in Fig. 2. 227 Fig. 4 provides details on the model setup. The wet mire is divided into three tiles with surface 228 elevation at 300 m asl. They are composed of a 3 m thick layer of unfrozen saturated peat above a 7 m 229 thick silty mineral layer that also extends below the plateau. Their respective widths are 50, 2 and 0.5 230 m. These values are significantly higher than those of the tiles initially hosting permafrost to ensure 231 stable unfrozen lateral boundary conditions in the mire side, as is observed in the plateau mire complex. 232 The largest and most external mire tile is linked to a hydrologic reservoir  to ensure 233 a steady water level in the mire at its surface (at 300 m asl), keeping it always saturated and preventing 234 ponding of surface water above 300 m asl. The peat plateau tiles are 0.3 m wide, giving a total width of 235 the plateau of 11.1 m. It stands at 302 m asl, making it 2 m high, in line with observations (Table 1). For 236 the initial steady state simulation (Sect. 3.2.5), the topographic transition between the mire and the 237 plateau was initially set as an arbitrary regular slope between 300 and 302 m asl as well as the active 238 layer depths, which were linearly interpolated between 0.9 m (at the contact with the mire) and 0.7 m at 239 the other hand of the plateau. The excess ice is distributed homogeneously along depth and the soil grid 240 ice content is set so that the plateau surface reaches 300 m asl (mire surface elevation) when the excess 241 ice is fully melted.  normal period. Additionally, the period was wetter than average, with 2015 and 2016 being, respectively, 255 33 and 50 % wetter than the normal period. As such, the simulated period gives a good opportunity to 256 study the response of ground surface temperature to an anomalously warm and wet 12 month-long 257 climatic conditions (Grinde et al., 2018;Heiberg et al., 2017). This is of particular interest for this study 258 because both of these two factors enhance permafrost degradation (Sect. 1). 259 The forcing data are looped to generate 100-year time series with steady-state climate forcing. 260 To achieve a realistic initial temperature profile also in deeper layers, a 100-year spin-up is performed 261 for all simulations using the 0 cm snow scenario, for which the peat plateau is stable (Sect. 4.2). Note 262 that the other snow scenarios cannot be used for model spin-up, as the plateau edge starts to retreat 263 instantly, so that a true steady-state cannot be reached. layer (3 % porosity, as in Westermann et al., 2013). 271 Over the Šuoššjávri plateau, the soil stratigraphy presents a significant spatial variability. While or more subsidence (i.e. the sensitivity threshold of the measurements), which consequently implies that 290 81 % of the plateau is stable during the observation period. The mean subsidence value (considering 291 values larger than 5 cm) is 17±15 cm (1σ) and the median 12 cm, with 1.2 % the total plateau surface 292 subsiding by more than 40 cm. The maximum observed subsidence is a one square meter patch in Zone 293 6, exhibiting between 1.5 and 1.7 m of subsidence over the 3 years. 294 We extracted the plateau edge over its outermost 2 m (i.e. the band of the 20 outermost pixels 295 delimiting the plateau). We find that this surface corresponds to one third of the total plateau surface, 296 but represents 77 % of the total subsidence (including the rims of the depressions within the plateau). 297 The distribution of subsidence values for the whole plateau and its edge only are presented in Fig. 6. 298 Due to the differences in peat plateau stratigraphies detailed in Sect. 3.3.3., the west side of the 299 plateau features higher subsidence values than the east side. On the eastern edge, ground subsidence is 300 lower due to the limited thickness of the peat layer, with mineral soil at a depth of less than one meter 301 below the surface. A description of the eight edge transect areas and their subsidence between 2015 and 302 2018 is presented in Table 1     year to year. We denote this phase as the "Initial Slope Adjustment" (ISA, Fig. 8 and 9). This phase lasts 344 for 40 years in the 5-10 cm snow simulation and 30 years for the 10-20 cm snow simulation. For the 345 same two simulations, the thermal erosion later affects the slope in a more uniform way and the plateau 346 edge retreats with a constant slope at a constant rate; during 60 years for the 5-10 cm snow simulation 347 and 20 years for the 10-20 cm snow one. We denote this phase the "Constant Edge Degradation" (CED, 348 Fig. 8 and 9). During the second half of the 10-20 cm snow simulation, both the edge and the top of the 349 plateau subside. We denote this phase as the "Plateau Collapse" (PC, Fig. 8 and 9). Contrary to the 10-350 20 cm snow simulation, the 20-30 cm snow simulation does not exhibit the phases of initial slope 351 adjustment and constant edge degradation but only the evolution corresponding to a plateau collapse.

358
The top panel of Fig. 9 presents the volumetric rate loss of the plateau for the three simulations. 359 For the 5-10 cm snow simulation, this rate is constant around 0.06 to 0.08 m 3 lost per meter of lateral 360 extension of the plateau (i.e. per meter of plateau circumference) and per year (m 3 m -1 yr -1 ) for the whole 361 simulation. For the 10-20 cm snow simulation, the volumetric loss during the initial slope adjustment 362 and the constant edge degradation phases shows a steady increase from 0.08 to 0.28 m 3 m -1 yr -1 . During 363 the plateau collapse phase, this rate steadily decreases to 0.12 m 3 m -1 yr -1 at the end of the simulation. 364 For the 20-30 cm snow simulation, the rate reaches 0.28 m 3 m -1 yr -1 over the first decades and stabilizes 365 at this value for 10-20 years before undergoing another rapid increase to 0.35 m 3 m -1 yr -1 , at which it 366 stabilizes until the end of the simulation. Further quantification of the thermal erosion is given in Sect.

396
The red square reminds that no subsidence was observed for the "0 cm snow" simulation.

Field measurements 399
At the western edge of the Šuoššjávri plateau, subsidence is highly variable, ranging from 0 to more 400 than 1 meter within 3 years. This pattern highlights the chaotic behavior of permafrost landscapes facing 401 degradation due to a positive feedback between subsidence, snow accumulation and water drainage. 402 When an initial perturbation, such as intense rainfall or extraordinary snow accumulation triggers 403 subsidence (Seppälä, 1988(Seppälä, , 2011, both the snow redistribution and the subsurface drainage towards the 404 mire are affected, which creates warmer surface conditions and, in return, triggers more subsidence. 405 Considering the complex geometry of the Šuoššjávri plateau edges, meter-scale variability of the snow 406 and hydrological conditions likely contribute to observed variability of ground subsidence. Besides, due 407 to the dependence of heat transfer to the surface area of the interface between the mire and the plateau, 408 the geometry of edges affects degradation speed. As such, zones 1, 2, 4 and 6 belong to convex features 409 of the plateau edges and show particularly high subsidence rates. 410 Our results confirm that edge degradation is a major degradation pathway of peat plateaus with 77 % 411 of the total subsidence occurring within the outermost 2 meters of the Šuoššjávri plateau. This result 412 shows consistency with Jones et al. (2016) who reported that 85 % of the degradation of forested 413 permafrost plateaus was due to lateral degradation along the margins. Estimating the elevation of the 414 inflection point of the plateau edge from the DEMs allows to quantify the plateau surface area at a given 415 time, even though in practice the precise positioning of such limit can be discussed. Between 2015 and 416 2018, we find that the Šuoššjávri plateau lost 3.2 % of its surface area. If we take the percentage of 417 annual loss rate of plateau area as 100 * (1 -Si+1/Si ) (where the fraction corresponds to the ratio between 418 the area at year i and the area at year i+1), then over an observation period of n years, the average ratio 419 can be expressed as: 420 Where r is the annual loss rate in % yr -1 , n is the number of year between the 2 observations, Sn the 422 plateau surface at the end of the observation period and S0 the plateau surface at the beginning of the 423 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

22/33
period. Applying Eq. 2, the aerial change in this study corresponds to a surface loss rate of 1.1 % yr -1 . 424 Reconstructing the Šuoššjávri peat plateau extent from 1956 to 2011 with aerial imagery, Borge et al. 425 (2017) observed annual loss rates that they compared to the 1956 extension. Using Eq. 2, we compute 426 the annual loss rate from their data to be 0. Plateau Collapse), the ISA and PC phases are less relevant than the CED phase for model to field 440 comparisons: the ISA phase is strongly affected by the initial response of the model to the applied 441 climate forcing and snow depths on the plateau, whereas the CED phase corresponds to the prolonged 442 edge retreat observed in the field, while the bulk of the peat plateau remains stable. The PC phase 443 corresponds to the sustained collapse of the plateau with ground subsidence in all parts, which is not 444 observed for the Šuoššjávri peat plateau, but regularly occurs for smaller circular palsas in the vicinity. 445 For such palsas, the assumption of translational symmetry inherent in the model setup is not valid and 446 simulations should be performed assuming rotational rather than translational symmetry. This suggests 447 that our simulations are indeed most realistic during the CED phase, whereas changes of the peat plateau 448 shape need to be taken into account to model the final stages of degradation. 449 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

The role of snow 450
Our simulations confirm the crucial role of snow on the ground thermal regime and peat plateau 451 degradation. The sensitivity of modeled ground temperatures towards snow depth is in good agreement 452 with the field measurements for the Šuoššjávri and the Iškoras peat plateau (40 km east of Šuoššjávri) 453 presented by Martin et al. (2019). For example, they showed that the measurement points on the plateau 454 exhibiting a mid-march snow depth smaller than 10 cm were associated with coldest mean annual 455 ground temperatures, and that 70 % of the plateau locations with stable permafrost had a March snow 456 depth smaller than 30 cm. 457 However, our idealized model approach assumes snow depths to be constant on the peat plateau, 458 which does not capture the significant spatial and interannual variability of snow depths on the plateau 459 observed in measurements (Fig. 7). In particular, the complex geometry of snow drift patterns (snow 460 accumulation) along the plateau edges, with snow drifts forming on lee sides, is not captured by the 461 simple snow redistribution model implemented in CryoGrid. Field observations show that snow drifts 462 along the plateau edges feature considerably higher snow depths than the surrounding wet mire, thus 463 introducing additional winter warming in the zone of maximum change. Additionally, persistent wind 464 patterns can strongly influence the distribution of snowdrifts. Unfortunately we do not possess field 465 measurements to discuss this parameter at the Šuoššjávri site. In CryoGrid, on the other hand, snow 466 removed from the plateau is evenly distributed over the entire mire, not taking edge effects into account. 467 Furthermore, our model assumes a fixed value for the snow density and thus snow thermal 468 properties, while this parameter in reality varies with e.g. snow depth and time, responding strongly to 469 synoptic conditions and imposing metamorphosis of the snowpack. Measurements of snow density in 470 Šuoššjávri showed that the snow on top of palsas is slightly less dense than in the mire. This could be 471 due to a thinner snowpack leading both to greater kinematic metamorphism and the formation of depth 472 hoar, notorious for high porosity and a real difficulty to measure density and high effective thermal 473 conductivity (Domine et al., 2016). A thinner snowpack also implies a lower overburden load and 474 therefore less compaction. Such limitations could be moderated by using more sophisticated snow 475 models taking snow microphysics and the transient evolution of snow density into account, such as 476 CROCUS (Vionnet et al., 2012) or SNOWPACK (Bartelt and Lehning, 2002). Yet, even these models 477 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

24/33
show limitations to reproduce the thermal characteristics of snow deposited in Arctic regions because 478 they do not account for the vapor fluxes in the snow pack, which significantly affect the snow thermal 479 conductivity profile (Domine et al., 2016). In this study, we present idealized numerical simulations of peat plateau thermal erosion that 483 reproduce the general patterns of edge retreat as observed in-situ through repeat digital elevation models. 484 We demonstrate that the snow depth on the plateau is a strong control for subsidence patterns and 485 dynamics. This result indicates that peat plateau systems will react sensitively to changes in the applied 486 climate forcing, not only regarding temperature but also regarding (snow) precipitation and windspeed 487 variations which all affects in turn snow pack building. In this regard, predictions regarding snow fall 488 for the coming decades are complex. An overall decrease of mean snowfall is expected at the global 489 scale (O'Gorman, 2014), consistently with observed trends over the past decades (Liston and Hiemstra, 490 2011) but this decrease will be accompanied by strong regional trend (Brown and Mote, 2009;Lader et 491 al., 2020) and the frequency and distribution of extreme snowfall remain unclear (O'Gorman, 2014). In 492 any case, it appears that for future modelling works, the accuracy of both the snowpack buildup and its 493 thermal properties at the small-scale (10-100 m) should be considered of major importance to robustly 494 simulate permafrost degradation. 495 Additionally, our implementation illustrates the idea of small-scale feedback mechanism on 496 permafrost degradation. The feedback between the dynamical microtopography and the lateral fluxes 497 shows how a limited increase in snow cover (when comparing the 10-20 cm snow and the 20-30 cm 498 snow simulations) results in a dramatically faster degradation rate. Such a sensitivity to minor 499 perturbation resulting into major modifications of permafrost degradation finds consistency with the 500 observed permafrost destabilization when punctually augmenting the snow depth with a fence (Hinkel 501 and Hurd, 2006), when implanting linear road infrastructures (Deimling et al., 2020) or due to the traffic 502 of heavy vehicle in Alaskan lowlands (Raynolds et al., 2020). 503 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

Spatiotemporal stability and degradation conditions 504
Our simple approach is clearly not able to capture the complex patterns of different subsidence rates 505 that are observed around the edges of the plateau (Fig. 5). On top of small-scale variations of ground 506 stratigraphy, excess ice contents and plateau heights, we suggest that the irregular plateau outline with 507 both concave and convex shapes also affects the lateral fluxes of heat, water and snow, which in turn 508 exert a control on the edge dynamics. While computationally demanding, our multi-tile approach could 509 be embedded in an ensemble framework to represent a range of edge geometries and other critical 510 parameters, yielding a range of different degradation scenarios and therefore capture the high spatial 511 variability of subsidence at the plateau scale. degradation (e.g. Borge et al., 2017). 517

518
Most Land Surface Models (LSMs) used to simulate the future response of permafrost to climate 519 changes still rely on simplified one-dimensional implementations of permafrost thaw dynamics which 520 ignores subsidence and only reflects gradual top-down thawing of the frozen ground (Andresen et al., 521 2020). Excess ice melt and the resulting microtopography changes exert a major control on the evolution 522 of hydrologic conditions, which in turn greatly influence the timing of permafrost degradation, as 523 demonstrated for polygonal tundra . Aas et al., (2019) presented a similar 524 approach for peat plateaus in Northern Norway. It is based on two tiles (one for the wet mire, one for 525 the plateau) and reproduces both climate-induced stability and degradation. However, in this approach, 526 the plateau subsides as a whole when a climate-related threshold is exceeded and excess ice begins to 527 melt. This contrasts with our field observations, which show ongoing edge retreat on decadal timescales, 528 while the plateau interior is largely stable. Our approach uses a larger number of tiles to explicitly 529 represent the temperature and soil moisture gradients across the plateau edges, which causes excess ice 530 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

26/33
melt to only occur in a narrow zone at the plateau edge, in agreement with our observations. Over longer 531 timescales, on the other hand, this process leads to the reshaping and finally the complete collapse of 532 the entire peat plateau. In ESMs, implementing a multi-tile approach is challenging due to its complexity 533 and computational demands. Yet, parameterized approaches could likely be developed, based on 534 sensitivity tests using our framework. In particular, future studies should investigate to what extent the 535 two-tile approach demonstrated by Aas et al. (2019) can emulate the results of the multi-tile model, 536 especially when averaging the results over the range of climatic conditions under which peat plateaus 537 occur in the sub-Arctic domain. However, a multi-tile setup would probably reveal different 538 hydrological regimes (zones of well-drained plateaus would remain until the end), which would in turn 539 affect carbon decomposition. 540

541
We present and compare field measurements and numerical modeling of thermal erosion patterns 542 of the Šuoššjávri peat plateau in Northern Norway. We use high resolution digital elevation models 543 derived from drone-based photogrammetry to quantify changes of surface elevations of the plateau 544 between September 2015 and 2018. Thermal erosion of the plateau edges is the main process through 545 which thermal erosion occurs and accounts for 80 % of the total measured subsidence (in terms of 546 subsided soil volume), while most of the total plateau surface exhibits no detectable subsidence. We 547 show that this retreat corresponds to a volumetric loss of 0.13 ± 0.07 m 3 m -1 yr -1 for the edge transect 548 zones we studied. 549 Using the CryoGrid land surface model we show that these degradation patterns can be reproduced 550 numerically in a framework that implements lateral redistribution of snow, subsurface water and heat, 551 as well as excess-ice-melt-triggered subsidence. Overall, the modeled volumetric change rates are in the 552 same order of magnitude as the measurements. Based on a steady-state climate forcing, our simulations 553 demonstrate the importance of lateral snow transport and resulting snow depths on the plateau. The 554 modelled peat plateau is fully stable when all snow on its top is transported towards the mire (0 cm snow 555 depth on the plateau), whereas its edges degrade at increasing rates with increasing snow depths. While 556 https://doi.org/10.5194/tc-2020-338 Preprint. Discussion started: 22 December 2020 c Author(s) 2020. CC BY 4.0 License.

27/33
a maximum of 5-10 cm on the plateau only triggers a 4-5 meters edge retreat within 100 years, a 10-20 557 cm cover fully degrades the plateau in 100 years and a 20-30 cm cover degrades it in 40 years. Our 558 simulations reproduce the observed lateral edge degradation with a stable plateau top, but the final phase 559 of plateau degradation corresponds to complete plateau collapse with subsidence occurring throughout 560 the entire plateau. 561 These results highlight the fast and high spatial variability of permafrost landscape evolution in 562 response to climate change. They also show that the related microtopographic, thermal and hydrological 563 modifications can be represented in numerical models, opening the way for substantial improvements 564 in simulating permafrost landscape evolution and its impact on greenhouse gas emissions from thawing 565 permafrost. Code availability. The model code and settings used for the simulations is available from 571 github.com/CryoGrid/CryoGrid3/tree/xice_mpi_palsa_newsnow. It will be permanently deposited upon 572 acceptance of the manuscript. The code is published under the GNU General Public License v3.0. 573 Data availability. Field data will be permanently deposited on archive.sigma2.no upon acceptance of 574 the manuscript.  Martin et al. (2019) compared to the field measurements from Martin et al. (2019) for the same region.

592
Values indicated with the letter n correspond to the number of field observations in Martin et al. (2019).

593
While MAGSTs from the present study are in perfect agreement with the measured ones on the 594 peat plateau, the model seems to underestimates them slightly (0.5°C too cold) in the mire. This in 595 mainly due to a too low water inflow from the reservoir which does not manage to always fully saturate 596 the wet mire. As a consequence, the uppermost layer of the mire partially dries during summer, which 597 shifts the heat flux during autumn from latent to sensible and imposes colder temperatures in the mire 598 during winter. 599 Active layers are overestimated by 20 to 30 cm for snow depths smaller than 10 cm. This 600 probably arises from the difference between the real and simulated snow conditions. While real 601 conditions can be erratic and show important variations during winter (from snow free to snow covered,