Permafrost degradation at two monitored palsa mires in north-west Finland

,

Both peat plateaus and palsas are characterized by a thick, upheaved layer of peat and perennially frozen core (Zoltai and Tarnocai, 1975;Meier, 2015). The main differences are in their topographical profiles. Peat plateaus are usually extensive, flat-topped and less than three metres high, while peat mounds with a smaller areal extent but a more dome-shaped profile are deemed palsas. In the Fennoscandian context, both are often called palsas due to the similarity of their formation processes 25 and common occurrence within the same peatland areas (Seppälä, 2011;Meier, 2015). Additionally, smaller mounds identified as separate palsas can also be disintegrated parts of larger peat plateaus (Borge et al., 2017).
The rapid lateral degradation of palsas and peat plateaus during the past several decades has been reported from Fennoscandia (Zuidhoff and Kolstrup, 2000;Meier, 2015;Borge et al., 2017;Olvmo et al., 2020) and North America (Payette et al., 2004;Jones et al., 2016;Mamet et al., 2017). The occurrence of widespread permafrost degradation also before the time covered by 35 aerial photography time series is indicated by the presence of thermokarst ponds in the earliest available aerial photographs (Luoto and Seppälä, 2003;Payette et al., 2004;Borge et al., 2017).
The periods of most rapid lateral degradation vary depending on the studied region, methodology and temporal coverage of aerial photography time series used to map the extent of palsas and peat plateaus. A comparison of palsa degradation rates (area loss in % a -1 ) during different periods and respective climatic conditions has linked rapid lateral degradation to increases 40 in precipitation and temperature (Payette et al., 2004;Olvmo et al., 2020). On the contrary, Borge et al. (2017) did not find a correlation between the degradation rates and climatic parameters.
A delineation of the palsa extent from orthophotos does not provide information about the height and volume of palsas unless stereophotographs are available. High-resolution Light Detection and Ranging (LiDAR)-based Digital Elevation Models (DEMs) are therefore a great asset for (regional) palsa inventories (Wramner et al., 2015;Ruuhijärvi et al., 2022). The 45 increasing availability of Unmanned Aerial Systems (UAS) allows even more frequent and detailed mapping and change detection of both the extent and topography of palsas and peat plateaus (Martin et al., 2021;Fraser et al., 2022).
Active layer thickness (annually thawing and freezing ground layer underlain by permafrost) is one of the main parameters commonly used to monitor permafrost conditions (Brown et al., 2000). In the Northern Hemisphere, the ALT varies a from few decimetres in the High Arctic to several metres in the mountainous regions at mid-latitudes (Luo et al., 2016). An analysis 50 of the long-term trends showed an increase in the ALT, especially in the Qinghai-Tibetan Plateau (Luo et al., 2016), Greenland https://doi.org/10.5194/egusphere-2022-1173 Preprint. Discussion started: 14 December 2022 c Author(s) 2022. CC BY 4.0 License. and Scandinavia (Strand et al., 2021), West Siberia (Smith et al., 2021) and the interior of Alaska . The lack of a significant change in the ALT on the Alaskan North Slope and the Mackenzie River delta has been (partially) attributed to the effect of thaw consolidation following the melting of ice in the transient layer at the boundary between the active layer and the permafrost table (Shiklomanov et al., 2013;Streletsky et al., 2017;O'Neill et al., 2019). Some active layer 55 data are also available from palsas and/or peat plateaus (Åkerman and Johansson, 2008;Sannel et al., 2016;Mamet et al., 2017). More than decade-long ALT records from palsas are scarce, however.
In this research, we present the results of long-term monitoring of two palsa mires in the Kilpisjärvi area, north-west Finland inspected annually in late summer. The aim of our study is to investigate permafrost degradation dynamics and their relation to climatic parameters. In particular, we aim to answer the following questions: (1) How did the active layer, palsa height and 60 lateral extent change during the investigation periods? (2) Can the detected changes be related to climatic drivers and if yes, which?
The study consists of three overlapping periods of decreasing length and increasing data availability. Aerial orthophotos with various intervals and continuous meteorological observations are available from 1960 to the present. The annual active layer monitoring grids including palsa surface measurements with high accuracy Real-Time Kinematic (RTK) GNSS were 65 established on two palsas in 2007. For the last six years (2016-2021), the UAS data have enabled capturing changes in palsas at ultra-high resolution. This unique data set that includes multiple spatial and temporal scales allows comprehensive analysis of both lateral permafrost degradation and active layer dynamics of palsas in Finland.

Study sites
The study sites are located in the Kilpisjärvi area, the northwestern part of Finnish Lapland (Fig. 1a). The mean annual air 70 temperature measured at the Kilpisjärvi Weather Station is -1.38 °C (1990-2020FMI, 2022). The average air temperatures of the coldest (January) and warmest (July) months are -12.27 °C and 11.66 °C, respectively. The mean annual precipitation is about 514 mm, of which around 240 mm falls as snow (FMI, 2022).
Kilpisjärvi is located on the pre-alpine belt of the Scandinavian Mountains. The general landscape is characterized by flat summits typical for old peneplain surfaces and valleys with mires that have usually less than a 2 metre-thick peat layer (King 75 and Seppälä, 1987). This topography supports the formation of cold air ponds in the valleys and, in combination with low snow depths, favours the development of palsa mounds (King and Seppälä, 1987). South-southeast is the prevailing wind direction during the wintertime (November-April) as measured at the Kilpisjärvi Weather Station (FMI, 2022).

Peera
The Peera palsa mire is located south of Peerajärvi Lake, at 68°52'41" N, 21°04'41" E and around 460 m a.s.l. The studied mire area is ca. 5.2 ha and has four palsa mounds, of which the northernmost (hereafter 'Peera palsa') has been monitored in 90 this research (Fig. 1c). In 2021, the highest point of the Peera palsa was 2 metres above the surrounding mire, and the diameter of the palsa was around 50 metres. The vegetation cover on the palsa (Fig. 1d) varies from almost bare peat surface to dwarf birch shrubs of up to 50 cm in height. The average peat depth around the palsa is 1.3 ± 0.3 m (± SD).

RTK-GNSS and active layer measurements
The ALT monitoring sites were established at the end of August 2007. Palsa surface and edges were surveyed using the RTK-GNSS (Topcon Hiper Pro, II, and V) with 1 cm horizontal and vertical accuracies. Permanent control points established by the Finnish Geodetic Institute were used for the local RTK-GNSS base station installation. A 1 m long steel rod was used to 105 measure the depth of the active layer. In total, 515 points were measured at Peera (Fig. 2a), of which ca. 180 were placed every 2 m along 11 transects covering half of the palsa, indicated by the yellow rectangle in Fig. 1c. In Laassaniemi, 374 points were measured (Fig. 2c), of which ca. 240 were placed every 1 m along 18 transects. In the following years (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), the number of measured points and length of transects varied because many of the points were abolished after years of no permafrost (ALT > 1 m) observations and/or due to the collapse of the palsas. The measurements were done during the last week of August, 110 except in 2013, when the sites were visited on the 9 th and 10 th of September. In 2018, the active layer was not measured at exactly the same points as in other years due to an RTK-GNSS malfunction. Therefore, we interpolated the 2018 ALT values using the Natural Neighbour method and extracted the values from the interpolated raster layers at the locations of points in 2019.
The Natural Neighbour interpolation of the elevation information was used to create DEMs from each year's RTK-GNSS data. 115 We determined the edges of the Laassaniemi palsa and the western half of the Peera palsa for each year based on these DEMs and the ALT values. We then subtracted the base altitude (459.5 and 474 m a.s.l. for Peera and Laassaniemi, respectively), representing the water level next to the palsa from the DEMs to analyse the changes in palsa height and volume. All geospatial data were processed in ArcGIS Pro (ESRI, Inc.).
We delineated the areas not visibly affected by lateral degradation of permafrost to analyse the effects of climatic drivers on 120 the active layer without the effect of thermokarst at the edges of the palsas (top-of-palsa, Fig. 2). The ALT refers to this topof-palsa ALT in the rest of this article unless indicated otherwise. The trends in the ALT were analysed by fitting linear regression models using the least-squares method.
A correction for the thaw subsidence has been suggested for the ALT observations in ice-rich permafrost areas where the subsidence is slow and spatially uniform (Shiklomanov et al., 2013;Streletskiy et al., 2017). Despite the permafrost thaw in 125 palsa mires being relatively fast compared to the areas described by Shiklomanov et al. (2013)  we also tested whether using subsidence-corrected ALT values would improve the model fits. The correction for subsidence was done by calculating the sum of the annual mean ALT and total surface subsidence compared to the 2007 level.

Aerial and UAS time series
We used a combination of aerial orthophoto and UAS orthomosaic time series to assess the changes in the palsa extent between To quantify the lateral degradation of permafrost, palsa edges were delineated by visual interpretation of the orthophotos and orthomosaics. Distinguishing palsas from old, panchromatic orthophotos is challenging, especially in shadowed areas.
Therefore, the extent of palsas at the beginning of the time series is more ambiguous than in more recent years. In uncertain cases during delineation, the edge resulting in the lesser palsa area was prioritized to avoid overestimation. The time series where Astart and Aend are the total area of palsas at the start year (Ystart) of the respective period and at the end of it (Yend).
UAS DEMs, resampled to a 10 cm cell size, were used to analyse height and volume changes (2016-2021) of the two main 145 palsas, from which also the ALT and RTK-GNSS data were collected. Variations in the UA systems, settings used in the data collection and the devices used to collect the coordinates of the GCPs resulted in discrepancies in the DEMs of different years. Therefore, we used the palsa polygons as delineated from the 2016 orthomosaic to extract only the areas of the main palsas from the DEMs. We then used the minimum value within that area as the base altitude for the respective year.

Climate data and statistical analysis 150
The meteorological data in Kilpisjärvi were obtained from the official weather station of the Finnish Meteorological Institute (FMI, Fig. 1b). The data set includes daily values of precipitation (mm), snow depth (cm) and mean air temperature (°C) from 1960 to 2021 as well as wind speed (m s -1 ) and wind direction (°) from 1998 to 2021. We calculated several climatic parameters known to affect the active layer and/or permafrost degradation (Meier, 2015). These parameters are (square root of) thawing and freezing degree days (TDD and FDD), mean annual, summer and winter air 155 temperatures (MAAT, MSAT and MWAT), and precipitation during different seasons. We also examined the effects of the number of days with particularly low (< -10 °C) and warm (> 15 °C) air temperatures and the number of days with over 10 cm snow depth. We used the hydrological year starting from the beginning of September and ending at the end of August for all annual average values. TDD was calculated as the sum of positive daily temperatures from the start of the calendar year until the day of measurement of the active layer (for the comparison with the ALT values) or until the end of August (for the 160 comparison with annual area loss rates).
We also extracted the timing (day of year) and duration of seasonal snow cover, maximum and end of March snow depth and the fraction of precipitation as snow. As the snow accumulation on wind-exposed palsas is likely lower than measured at the station, we also estimated local snow depths at our sites using SnowModel (Liston and Elder, 2006) We used linear regression to assess temporal trends in the climatic parameters and to find which of the parameters have an 170 effect on the ALT, RTK-GNSS-based volume loss and annual area loss rate. For the analysis of climatic drivers of the lateral palsa degradation, we used average values for the same four time periods as in Sect. 3.2. All statistical analyses were done using R (R Core Team, 2021). The regressions of temporal trends and correlations were considered statistically significant if p ≤ 0.05.

Active layer thickness
The ALT fluctuates annually (Fig. 3). At both palsas, the number of values over 1 metre has increased, especially at the edges (Figs. 2b and 2d). During the study period of 2007-2021, the average (± SD) ALT at Peera was 60 ± 4.6 cm. A linear trend suggests a moderate deepening of the thawed layer of ca. 0.3 cm a -1 . The trend, however, is not statistically significant (R 2 = 0.07, p = 0.353), and the ALT has decreased after reaching ca. 70 cm in 2014. At Laassaniemi, the average ALT was 49 ± 6.5 180 cm, and the linear trend indicates a significant decrease in the thawed layer at the rate of ca. 1.2 cm a -1 (R 2 = 0.68, p < 0.001).

Climatic drivers of the active layer thickness
The R 2 and p-values of the regression analyses are presented in Table A1 relationship between winter air temperatures and ALT at Peera, and some correlation at Laassaniemi (Figs. 4d and 4e). The direction of the correlation at Laassaniemi suggests that the active layer is shallower after milder winters. The amount of precipitation affected the ALT only at Peera, and the correlation was significant only for ALT and the precipitation sum for September and October of the previous year (Fig. 4f). Timing of the snow cover onset had a significant negative correlation at both sites (Fig. 4g), while the thickness of the snow cover showed correlation with the ALT only at Peera (Figs. 4h and 4i). Surprisingly, the correlation at Peera was greater with the snow depth values measured at the weather 205 station than with the local snow depth estimated using SnowModel. The regression analysis of volume loss in relation to climatic parameters did not show any statistically significant correlations.
However, the results suggest some correlation with air temperatures (MAAT and MWAT) and spring rainfall at Peera, and some correlation with snow cover onset and duration at Laassaniemi. The subsidence-corrected ALT correlated very little with the climatic parameters. 210

RTK-GNSS-based surface elevation and volume change
Degradation of permafrost is noticeable at both palsas (Fig. 5), but has been relatively stronger at Laassaniemi than at the western half of the Peera palsa (Table 1). At Peera, most of the degradation has been taking place in the southwestern corner.
At Laassaniemi, the palsa has degraded from all sides, especially from its southern and northeastern edges. The top-of-palsa areas lowered by 30 cm (-2 cm a -1 ) and 74 cm (-6 cm a -1 ) at Peera and Laassaniemi, respectively, between 2007 and 2021. 215

Height and volume change based on UAS DEMs
The mean height of both palsas decreased by 20 cm between 2016 and 2021 corresponding to ca. -17 % change at Peera and 220 ca. -32 % change at Laassaniemi (Table 2). Overall, the lateral degradation was greater at Peera, while the relative volume change at both palsas was similar at ca. -33 % to -38 %. The 2016-2021 DEM difference shows that the northeastern part of the Peera palsa has been relatively stable compared to the southwestern half and the edges, where the thermokarst is the most active ( Figs. 6a and 6b). At Laassaniemi, the surface subsidence appears to be lowest at the southwestern side and gradually increases towards the northeastern edge of the palsa (Figs. 6c and 6d). At both sites, the apparent increase in the height at the 225 northern edges is due to vegetation growth and/or movement with the collapsed palsa surface.
The area of Laassaniemi palsa was in 2021 almost 40 m 2 , or 32.5 % larger when delineated using UAS data (122.2 m 2 ) than the area derived based on the RTK-GNSS and ALT measurements (82.5 m 2 ). This is due to the difference in how the extent of the palsa, and consequently the permafrost, is defined using these methods. Detection of the palsa edge from the UAS data relies on differences in surface topography and vegetation, whereas in the RTK-GNSS and ALT approach, the information 230 about the active layer affects the delineation of the palsa edge.

Palsa area change 1960-2021
The overall extent of palsas decreased from 2.18 and 0.39 ha in 1960 to 0.50 and 0.04 ha in 2021 at Peera and Laassaniemi, respectively (Fig. 7). The total palsa area loss between 1960 and 2021 has been -77.3 % at Peera and -89.7 % at Laassaniemi. 240 The average loss rates for the whole study period are -2.4 % a -1 at Peera and -3.6 % a -1 at Laassaniemi. Contrary to the accelerating degradation at Peera towards the end of the time series, the degradation rate at Laassaniemi has been the lowest during the 2012-2021 period (Fig. 8). No signs of anthropogenic disturbance, such as artificial draining or excavation, are visible in the aerial data.

Climate trends and lateral palsa degradation
The climatic averages for each of the four periods are summarized in Table 3. The mean annual air temperature increased significantly over the past 60 years by ca. 0.36 °C per decade. The average air temperatures have warmed especially after 1988 (Fig. 9a), and the increase has generally been stronger in the winter months than in the summer. The precipitation, both in the form of rain and snow, has also increased significantly by ca. 25 mm per decade. The interannual variability in precipitation 255 sum is most pronounced in 1985-2000 (Fig. 9b). The period 2000-2012 had the highest amount of precipitation as rain, while the highest average amount of snow was between 1985 and 2000. The snow cover duration was also on average the longest in 1985-2000 and has been decreasing since then.
At Peera, the annual area loss rates correlated the strongest with the parameters related to the air temperature and, to a lesser degree, with the precipitation sum and precipitation as snow (Table A2)   . Annual area loss rates (% a -1 ) of permafrost area within the Peera and Laassaniemi palsa mires based on digitized palsa edges from aerial photographs (1960,1985,2000,2012) 1961-1984 1985-1999 2000-2011 2012-2021 1961-2021   indicates snow cover duration). The four periods (I-IV) correspond to the periods used for the annual area loss rates.

Permafrost degradation 1960-2021
In Fennoscandia, the climatologically most optimal conditions for palsa existence are MAAT -3 --5 °C and annual precipitation sum of <450 mm (Luoto et al., 2004a). The MAAT was already -2.6 °C in Kilpisjärvi between 1960 and 1985 280 (Table 3). Over the past decades, both air temperatures and precipitation have increased significantly, especially during the wintertime. Following these climatic changes, palsas in the Kilpisjärvi area have experienced a rapid and continuous decrease in area over the past sixty years, and no new palsas have formed.
The degradation rates differ between our two sites and between the four time periods. At Peera, the trend of increasing area loss rates follows the increase in air temperatures and precipitation similar to the permafrost peatland near Hudson Bay in 285 Canada (Payette et al., 2004) and at the Vissátvuopmi palsa mire in Sweden (Olvmo et al., 2020). Based on the volume loss rates derived from the RTK-GNSS surveys, the permafrost degradation at Peera is particularly driven by the increase in winter air temperatures.
At Laassaniemi, the area loss rate peaked during 1985-2000, which is also the period with the longest snow cover duration and highest maximum snow depths. The RTK-GNSS-based volume loss appears to correlate mainly with the snow cover 290 duration as well. Opposite to Peera, the annual area loss rate was lowest during the last period (2012-2021) at Laassaniemi.
Palsas at Laassaniemi are generally much lower than at Peera and are surrounded by dwarf birch thickets and tall common cotton grass (Fig. 1f), which increases the possibility of errors in the delineation of palsa edges by visual analysis of the aerial data. Therefore, it is possible that the UAS-based palsa edge delineation overestimates the permafrost extent, which is also indicated by ca. 30 % larger palsa area compared to the RTK-GNSS and ALT -based extent (Fig. 6). 295 Variability in local conditions, such as microclimate, vegetation cover, hydrology, and palsa shape, may cause differences in degradation rates despite a similar regional climate (Borge et al., 2017). The increasing fragmentation and complexity of the palsas' shape in Peera can have a positive feedback effect on the area loss rate, and vice versa. The effect of local conditions is also demonstrated by an almost three times slower degradation rate at Vissátvuopmi (Olvmo et al., 2020) than at our sites, although the distance between the areas is only 10-20 km. A direct comparison with other area loss rates reported from North 300 America (Payette et al., 2004;Jones et al., 2016;Mamet et al., 2017) and Norway (Borge et al., 2017) is not possible due to the different ways of calculating the average annual change in per cent. However, when calculated using the same method as used here and in Olvmo et al. (2020), the long-term (~ 1950s to 2010s) degradation rates, for instance, at sites D6 and GF in Mamet et al. (2017) and in Goatheluoppal palsa site in Borge et al. (2017), are similar to ours at ca. 2-4 % a -1 (see supp. material in Olvmo et al., 2020).

Active layer and climate dynamics
The analysis of the active layer and climate dynamics over the past 14 years revealed two different cases. At Peera, the ALT varied annually following summer air temperatures, autumn precipitation, snow cover onset and snow depth. The correlation with the snow cover measured at the weather station is uncertain, however, as the snow accumulation is generally lower on the Peera palsa (Störmer, 2020) than at the station (FMI, 2022). There is no clear trend in the ALT between 2007 and2021, 310 although the height and area of the palsa were decreasing annually. At Laassaniemi, the annual variations in the ALT could be explained only by the differences in snow cover onset and TDDs. The active layer is becoming shallower, while the palsa surface is continuously subsiding and thawing is occurring at the edges.
The lack of correlation between the ALT and climatic parameters indicates that other factors, such as vegetation cover, hydrology, soil properties or palsa shape, are more important factors regulating the seasonal thaw at the Laassaniemi palsa. A 315 shallower active layer in lower palsas compared to higher ones in the same area has also been observed elsewhere in northern Finland (Rönkkö and Seppälä, 2003;Seppälä, 2011). Seppälä (2011) attributed the thicker active layer in higher palsas to a longer snow-free period, which allows more time to receive solar radiation and thaw. Therefore, the lack of any significant increase in the ALT at Peera and a significant decrease in the ALT at Laassaniemi could be (partially) attributed to the lowering of the palsas. However, more palsa and peat plateau sites including ALT and surface elevation monitoring as well as local 320 snow cover observations are needed to support these hypotheses.
A comparison of our results to the published active layer records at the geographically close locations in Sweden shows some similarities, but also contrasting trends. Similar to the observations in Tavvavuoma (Sannel et al., 2016) were also observed in Abisko (CALM, 2022). However, the deepening of the active layer in Abisko continued after 2014, 325 whereas at Peera, the ALT returned to similar values as before 2013. And at Laassaniemi, the active layer was even shallower in 2021 than after a particularly cold summer in 2012.
The strong negative correlation between the ALT and snow cover onset could be due to the longer time of peat saturation, which increases its thermal conductivity (Kujala et al., 2008) and promotes cooling of the ground in late autumn. The timing of snow accumulation has also been shown to have an impact on the ALT at a peat plateau in northern Sweden (Johansson et 330 al., 2013). Warm air temperature in November can delay the snow cover onset, which could explain the surprising negative correlation between the winter air temperatures, in which the November is also included, and the ALT at Laassaniemi. The changes in snow cover duration have been observed in various arctic regions, and an even stronger shift in snow-on and snowoff dates is expected with the projected climate change (AMAP, 2017). The effects of snow cover timing on the ground thermal regime and ALT have been so far investigated mainly in continuous permafrost environments (Ling and Zhang, 2003;O'Neill and Burn, 2017;Jan and Painter, 2020). Hence, the response of discontinuous and sporadic permafrost to changing snow cover duration warrants more research.
Correcting the ALT values for subsidence did not improve the linear model fits with thaw degree days or any other climatic parameters at our sites. This was not surprising as the rate of subsidence at our sites has been around two to six times greater than at the Alaska North Slope (Shiklomanov et al., 2013;Streletskiy et al., 2017) and even more compared to the Mackenzie 340 River delta in Canada (O'Neill et al., 2019). However, as noted by Martin et al. (2019), the permafrost thaw occurring from the bottom may also cause subsidence of palsas without an effect on the ALT. This suggests that the active layer may appear stable due to subsidence also in rapidly degrading sporadic permafrost environments.
Our results imply that using the ALT alone for assessing permafrost conditions can be insufficient in permafrost peatlands and, therefore, underline the importance of high-accuracy surface position measurements at the ALT monitoring sites. The 345 possibility that the active layer may remain stable or become shallower in the inner parts of palsas despite overall permafrost degradation should be considered when analysing ALT records from similar environments and in modelling efforts of biogeochemical, hydrological, and permafrost dynamics of palsa mires.

Conclusions
The analysis of the aerial photography time series indicates a rapid degradation of the permafrost area in north-west Finland 350 over the past 60 years. At the Peera and Laassaniemi palsa mires, the permafrost area in 2021 was less than 25 % of that in the 1960s. At the same time, air temperatures and precipitation have increased significantly. The accelerating lateral degradation at Peera coincides with increasing (winter) air temperatures, while the lateral degradation at Laassaniemi was greatest during the period of high winter precipitation in the 1990s.
Although we found no increasing trend for the ALT between 2007 and 2021, permafrost degradation is evident in the form of 355 lowering palsa surfaces and thermokarst at their edges. The active layer was deeper especially after winters with thick snow cover and warm summers at Peera, while at Laassaniemi, the relationships between ALT and climatic drivers are weak. Thinner active layer at both sites following a later snow cover onset in the preceding autumn indicates the complexity of feedbacks related to climate warming and consequent shifts in snow cover duration in arctic regions. Furthermore, the decreasing ALT at Laassaniemi despite the overall decrease in palsa area and height emphasises that the active layer thickness alone is not 360 enough to assess the (in)stability of permafrost.
Data availability. RTK-GNSS, ALT, and UAS data used in this study are available upon request from the authors.