Rapid Vegetation Succession and Coupled Permafrost Dynamics in Arctic Thaw Ponds in the Siberian Lowland Tundra

Thermokarst features, such as thaw ponds, are hotspots for methane emissions in warming lowland tundra. Presently we lack quantitative knowledge on the formation rates of thaw ponds and subsequent vegetation succession, necessary to determine their net contribution to greenhouse gas emissions. This study sets out to identify development trajectories and formation rates of small‐scale (<100 m2), shallow arctic thaw ponds in north‐eastern Siberia. We selected 40 ponds of different age classes based on a time‐series of satellite images and measured vegetation composition, microtopography, water table, and thaw depth in the field and measured age of colonizing shrubs in thaw ponds using dendrochronology. We found that young ponds are characterized by dead shrubs, while older ponds show rapid terrestrialization through colonization by sedges and Sphagnum moss. While dead shrubs and open water are associated with permafrost degradation (lower surface elevation, larger thaw depth), sites with sedge and in particular Sphagnum display indications of permafrost recovery. Recruitment of Betula nana on Sphagnum carpets in ponds indicates a potential recovery toward shrub‐dominated vegetation, although it remains unclear if and on what timescale this occurs. Our results suggest that thaw ponds display potentially cyclic vegetation succession associated with permafrost degradation and recovery. Pond formation and initial colonization by sedges can occur on subdecadal timescales, suggesting rapid degradation and initial recovery of permafrost. The rates of formation and recovery of small‐scale, shallow thaw ponds have implications for the greening/browning dynamics and carbon balance of this ecosystem.


Introduction
The Arctic is warming twice as fast as the global average. It is expected that by the end of the 21st century, the global area of near-surface permafrost will decrease by 24 ± 16% (IPCC scenario RCP2.6) to 69 ± 20% (IPCC scenario RCP8.5) (Meredith et al., 2019). This thawing of permafrost exposes previously frozen soil organic carbon to microbial decomposition, resulting in the release of greenhouse gases (GHG). This constitutes a positive feedback to climate change (Schuur et al., 2015). The expected release of GHG from thawing permafrost is in the order of 10s to 100s gt carbon by the end of the 21st century (Meredith et al., 2019). There is, however, considerable uncertainty about the exact magnitude of this "permafrost carbon feedback" (Meredith et al., 2019;Schuur et al., 2015). Part of this uncertainty revolves around the relative importance of two thawing mechanisms: gradual versus abrupt thaw. Gradually increasing thaw depths and higher temperatures have resulted in increased vegetation cover and productivity (Tape et al., 2006), most notably in shrubs (Myers-Smith, Elmendorf, et al., 2015), which may partly offset GHG emissions from thawing permafrost (Euskirchen et al., 2009). Satellite monitoring of the Normalized Difference Vegetation Index (NDVI), a measure of vegetation greenness and photosynthetic capacity, confirms that the Arctic has shown a greening trend over most of the past decades (Epstein et al., 2018). In contrast, abrupt permafrost thaw occurs when excess ground-ice in the permafrost melts, causing a collapse of the soil surface (thermokarst) (Schuur et al., 2015). Such processes may counteract the widely observed Arctic greening trend (Abbott et al., 2016;Phoenix & Bjerke, 2016;. Since 2011, local browning (a decrease in NDVI) has diminished or even halted the overall Arctic greening trend (Epstein et al., 2018). Abrupt rather than gradual thaw may predominate in large parts of the Arctic (Olefeldt et al., 2016;Schuur et al., 2015) and has been found to occur increasingly throughout Arctic regions (Farquharson et al., 2019;Fedorov et al., 2014;Frost et al., 2018;Jorgenson et al., 2006;Jorgenson & Grosse, 2016;Liljedahl et al., 2016;Raynolds et al., 2014;Walter et al., 2006). Abrupt thaw can lead to formation and expansion of surface water in continuous permafrost (Jorgenson et al., 2006;Raynolds et al., 2014) and is expected to increase the area of small lakes by 50% by 2,100 under IPCC scenario RCP8.5 (Meredith et al., 2019). It thereby enhances GHG emission through increased wetting and disturbance of the existing vegetation (Abbott et al., 2016). However, most studies on permafrost carbon dynamics have focused on gradual thaw, and abrupt thaw processes are currently not considered in climate models (Meredith et al., 2019;Schuur et al., 2015;Van Huissteden & Dolman, 2012). Thaw ponds are an example of features resulting from abrupt thaw of which we have limited understanding. Thaw ponds are isolated depressions of less than 1 ha which form due to local melting of ground-ice. In these local depressions, snow and water accumulate, posing a positive feedback on thawing (Nauta et al., 2015). Thaw ponds are demonstrated sources of methane, as opposed to shrub-dominated tundra vegetation (Nauta et al., 2015;Van Huissteden & Dolman, 2012;Van Huissteden et al., 2005. Thaw ponds are abundant throughout ice-rich lowland tundra, but small waterbodies often go undetected in remote sensing analysis of thermokarst features (Campbell et al., 2018;Grosse et al., 2008;Raynolds et al., 2019). In contrast to the deep, ditch-shaped ponds resulting from ice-wedge degradation, melting of various other types of ground ice (e.g., ice lenses) may result in shallow, irregular shaped thaw ponds (Van Huissteden, 2020). Smaller, shallow ponds may be especially vulnerable to climate-induced changes in hydrological regimes (Campbell et al., 2018;Wolfe et al., 2011). Such thaw ponds have received little scientific attention compared to other abrupt thaw features (e.g., thaw lakes, ice-wedge degradation, retrogressive thaw slumps). This warrants research into trends in shallow thaw pond area, the processes that trigger formation of such thaw ponds and their long-term fate.
Several studies have focused on long-term development of thermokarst features related to ice-wedge degradation (Jorgenson et al., 2015;Kanevskiy et al., 2017;Liljedahl et al., 2016) or thermokarst bogs in subarctic regions (Myers-Smith et al., 2008;Payette et al., 2004;Robinson & Moore, 2000;Turetsky et al., 2007). These studies suggest the existence of a coupled vegetation succession and permafrost degradation-recovery mechanism, described by Kanevskiy et al. (2017) as "quasi-cyclic" succession. Such successions are characterized by soil subsidence due to melting of ground ice, followed by formation of open water features, which are in turn colonized by wetland vegetation (e.g., graminoids and peat mosses). Vegetation succession, especially colonization by peat mosses (Sphagnum), may improve insulation of the permafrost, which allows for permafrost recovery (Seppälä, 1988;Zoltai & Tarnocai, 1975).
We investigate whether shallow arctic thaw ponds in north-eastern Siberia, a large and understudied lowland tundra region, follow the quasi-successional succession model described above. As a novel aspect, we apply dendrochronological dating to assess the presence of Arctic dwarf shrubs across stadia of thaw pond development and to assess timespans of formation and recovery. Specifically, we explored (i) whether thaw ponds in lowland tundra display a dominant vegetation succession trajectory, (ii) if permafrost degradation and recovery are coupled to vegetation dynamics, and (iii) on which timescales vegetation succession and permafrost degradation and recovery in thaw ponds occur. We hypothesized that the following vegetation succession occurs: following soil subsidence and ponding of water in the new depressions, the initially dominant dwarf shrub vegetation drowns and is replaced by open water. Sedges are expected to be the first species to recolonize such open water areas, followed by formation of Sphagnum moss carpets. Additionally, we expect Sphagnum moss to be associated with permafrost recovery and to be an important microsite for the reestablishment of dwarf shrubs. Thus, we expect to find a (quasi-)cyclic succession mechanism similar to that described for ice-wedge degradation.

Study Site
Vegetation composition and associated abiotic conditions in shallow thaw ponds of various ages were studied in the "Kytalyk" Nature Reserve in north-eastern Siberia (70°49′N, 147°29′E) approximately 30 km north-west of the town of Chokurdakh. This area is situated in a region that displays frequent indications of browning rather than greening between 1982 and 2018 (Epstein et al., 2018). North-eastern Siberia is a region with high abundance of thermokarst wetlands and lakes (Olefeldt et al., 2016). Our study site is characterized by a shallow active layer overlying continuous permafrost of over 300 m thickness (Van Huissteden et al., 2005. Mean annual temperature is −13.4°C, with an average July temperature of 10.3°C . Mean annual precipitation is 196 mm, with 76 mm falling in June to August (1981-2010) (Trouet & Van Oldenborgh, 2013). Our study location is situated in a drained thaw lake basin or "alas," surrounded by remnants of Pleistocene yedoma deposits, floodplains, and thaw lakes. Our site represents an older episode of lake drainage, evident from the multiple levels of alases in the area ( Figure 1a). Paleoecological data suggest that peat accumulation has taken place in this lakebed for at least 4,000 years (Teltewskoi et al., 2016). Holocene alases are abundant throughout coastal eastern Siberia (Fedorov et al., 2018;Morgenstern et al., 2013) and cover approximately 111,000 km 2 in coastal Yakutia alone (Fedorov et al., 2018). The lakebed vegetation is characterized by slightly elevated shrub patches dominated by Betula nana, lichens, and mosses, interspersed with waterlogged depressions characterized by Eriophorum angustifolium, Carex spp., and Sphagnum spp. (Siewert et al., 2015). The Circumpolar Arctic Vegetation Map classifies this site as a mix of tussock-sedge tundra, erect dwarf shrub tundra, and sedge-moss wetland. These vegetation classes are typically found in low elevation areas with relatively high summer temperatures and are representative of a major part of coastal north-eastern Siberia (Raynolds et al., 2019).
Around Chokurdakh, ground ice content is high (75 vol-% on average in the top 1 to 2 m of soil) (Iwahana et al., 2014;Wang et al., 2019), and in recent decades a doubling in area of thaw ponds at the Kytalyk site has been observed (van Huissteden et al., 2017). Apart from these shallow, irregular shaped thaw ponds (Figure 1d), polygonal thermokarst features related to ice-wedge degradation ( Figure 1e) and diffuse drainage gullies (Figure 1c) also occur in the study area.

Thaw Pond Selection
For our research we selected isolated nonpolygonal thaw ponds (<100 m 2 ) surrounded by dwarf shrub vegetation (see Figure 1f). Thaw ponds were identified based on a series of commercial high-resolution satellite images (GeoEye I from 2010, WorldView II from 2015 and 2017, and WorldView III from 2018; specifications in supporting information Table S1) and inspected in the field. We observed both recently formed and older ponds. The latter frequently showed indications of terrestrialization. To represent temporal variability in thaw pond formation, we selected thaw ponds from three age classes based on their presence and extent in satellite images from previous years. Ponds that were not present yet in 2010 were labeled "Young ponds" (Y, n = 12). Ponds that were present in 2010 and showed an increase in the extent of open water from 2010 to 2018 were labeled "Old Increasing ponds" (OI, n = 13). Ponds that were present in 2010 and showed a decrease in the extent of open water from 2010 to 2018 were labeled "Old Decreasing ponds" (OD, n = 15). We assumed that these pond age classes represent successional stages from young ponds resulting from recent soil subsidence due to permafrost degradation to older ponds in which the open water area is decreasing due to terrestrialization. Table S22 contains an overview of the expansion trends of the 40 selected ponds. All thaw ponds were situated within a matrix of shrub vegetation with shallow thaw depths (±20 cm) to exclude external hydrological influences such as drainage. Satellite images were coregistered based on manually selected control points using the geo-referencing toolbox and pansharpened with the Gram-Schmidt algorithm in ArcMap 10.3.1 (ESRI, 2016). Pond boundaries were mapped manually as polygons based on visual distinction between open water and vegetated pixels in the pansharpened satellite images. While this manual delineation introduces a degree of subjectivity, we expect the effects of pond expansion and terrestrialization on the spectral data to exceed variation resulting from differences in interpretation (Table S22). Approximate lateral expansion rates of ponds were calculated based on changes in bounding box geometry of pond polygons.

Field Measurements 2.3.1. Vegetation Composition
To be able to identify dominant vegetation succession mechanisms, vegetation composition in ponds from the three age classes (Young, Old Increasing, and Old Decreasing) was studied. In all 40 selected thaw ponds, two perpendicular transects were laid out over the pond perimeter ( Figure 1f). Transect directions were selected to most accurately represent vegetation composition of the pond. The total length of the transect was measured as the distance from one pond margin with intact, Betula nana dominated vegetation to another (i.e., intact shrub vegetation outside the pond perimeter was not considered). Transitions in vegetation composition along the transect were marked with pickets ( Figure 1g). A vegetation transition is defined here as a transition from one dominant plant functional type (PFT) class to the other (Table 1). Boundaries were distinguished based on dominance of either of the PFT classes defined in Table 1 determined by visual cover estimates. Using this set-up, vegetation composition was quantified as total length between pickets for each class present within the thaw pond, expressed as a percentage of the total length of the transects. Measurements were conducted over the course of several days in late July 2018.

Abiotic Conditions
To assess the coupling of vegetation succession with the development of permafrost degradation and recovery, abiotic variables (microtopography, height of the water table above the permafrost, and thaw depth) were measured under various PFTs in thaw ponds of the three age classes. In a subset of 21 thaw ponds, 50 × 50 cm quadrats were placed on sites reflecting the PFT classes from Table 1 (Figure 1h), with the exception of "grass," "other moss," and "drowning shrubs," which were either only present in small and narrow areas (< 0.5 m width) or in few ponds (n < 5). The dominant PFT was determined by visual cover estimates. In the center of each of the 112 quadrats, thaw depth (referred to hereafter as TD), was measured using a pointed metal probe with cm scale indication. TD was defined as the distance between the moss surface and frozen soil as measured in early August, before maximum thaw at the end of the summer season. In our study site, temperatures typically fall below zero again in late August-early September (Parmentier, Van Der Molen, et al., 2011). As such, reported values should be interpreted as instantaneous thaw depths rather than an active layer thickness, which is generally defined as the top soil layer subject to annual freezing and thawing (Shur et al., 2011) or the depth of maximum seasonal penetration of the 0°C isotherm (Muller, 1947). The water table (hereafter: WT) was measured as the height of the water table above the frozen soil. In noninundated sites, a peat corer was used to determine whether there was a water table above the
Note. In order of hypothesized vegetation succession.

10.1029/2019JG005618
Journal of Geophysical Research: Biogeosciences permafrost. No WT was recorded in absence of a water table. Thickness of the organic layer (OLT) was measured as the distance between the moss or soil surface and mineral soil in a peat core. No OLT was recorded in case no mineral soil was visible (i.e., OLT > TD). Lastly, relative surface elevation (hereafter: RSE) of a quadrat was measured as millimeter elevation relative to a local reference level using an optical leveler (Kompensator-Nivellier NI 025, VEB Carl Zeiss, Jena, Germany). The reference level was calculated as the average elevation of four positions outside the pond perimeter on intact shrub vegetation, at 2 m distance of the pond margin. The permafrost position (PP) was calculated as RSE − TD, representing the depth at which frozen soil occurs (before maximum end-of-summer thaw) relative to the reference level. Measurements were conducted once per quadrat over the course of several days in early August 2018, except OLT, which was measured in early August 2019.

Shrub Distribution
To assess whether vegetation succession can lead to reestablishment of dominant dwarf shrub vegetation, an inventory was made of presence of living dwarf shrubs (yes/no) in all quadrats inside the pond perimeter.

Dendrochronological Analysis
To gain an insight in the timescales associated with degradation of permafrost, vegetation succession, and reestablishment of shrub vegetation, the ages of shrubs, their year of establishment, and death were determined using dendrochronology. In a subset of 15 thaw ponds (five from each age class), stem samples of Betula nana were taken for tree-ring analysis. Three types of Betula nana shrubs were sampled: (i) living shrubs from shrub vegetation surrounding the thaw pond, (ii) dead shrubs from within the thaw pond (including terrestrialized areas), and (iii) living shrubs (established from seeds or ramets) from within the thaw pond (including terrestrialized areas, e.g., on a moss carpet or other substrate). For each shrub type, two replicates per pond were sampled, if present, resulting in a total of 66 samples (30 living, outside ponds; 24 dead, inside ponds; and 12 living, inside ponds). From all sampled shrubs, three stem sections of 2 cm length were clipped at the field site and preserved in a 3:1 mixture of vodka and glycerol. Sections were cut at the root-shoot interface, at the moss surface, and above the moss surface just below the first branch (Li et al., 2016). The root-shoot interface was identified as a slight thickening in the stem diameter, often followed by a bend in the belowground stem (see also Figure S1).

Data Analysis 2.4.1. Vegetation Composition
Differences in PFT composition among thaw ponds of various age classes were assessed using multivariate analysis of variance. Cover percentages of both transects ( Figure 1c) were averaged per pond. Since the assumption of normality was violated and percentage data were zero-inflated, PERMANOVA using Bray-Curtis distance was used to test for differences in the covers of the PFTs among the three pond age classes (Anderson, 2001). Prior to PERMANOVA analysis, a betadisper test using the average distance of group members to the group spatial median (Anderson, 2006) followed by Tukey's Honest Significant Difference analysis was used to confirm multivariate homogeneity of dispersion among pond age classes.
To correct for an unbalanced design with small sample sizes, a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n n − 1 ð Þ ð Þ p correction (Stier et al., 2013) was used for the betadisper test. Pairwise PERMANOVAs for each pair of pond age classes were conducted as a post hoc comparison, using Bonferroni adjustment of p values to correct for multiple testing. Lastly, SIMPER analysis (Clarke, 1993) was used to assess contribution of specific PFTs to difference in overall PFT composition among pond age classes. SIMPER analysis was restricted to 1,000 permutations and used Bray-Curtis distance.

Abiotic Conditions
Differences in TD, WT, RSE, and PP among the PFTs in thaw ponds of the three age classes were analyzed using Linear Mixed Effects Models (LMMs) using Restricted Maximum Likelihood (REML). "Pond age class" and a quadrat's "PFT class" were used as categorical predictor variables. Since a water table was absent in 95% of the Intact Shrub quadrats (no water table above the permafrost), observations of WT in this PFT class were omitted from analysis. TD was log-transformed prior to analysis to improve residual diagnostics. All models were checked for significance of PFT class, pond age class, and their interaction term and for random slope and/or intercept for individual ponds. First, the optimal random structure was determined using Likelihood Ratio tests (LRTs) on a full (or as full as feasible) model. Then, the optimal fixed structure was determined using backward selection. P values of predictor variables were calculated on nested models, using F tests with Kenward-Roger approximation for fixed effects and likelihood ratio tests (LRTs) for random effects. Models were selected based on predictor p values, Akaike's Information Criterion (AIC), normality and homoscedasticity of residuals, and absence of patterns of residuals against grouping factors and fitted values. Normality of residuals for LMMs was determined visually. As a post hoc test, Tukey contrasts and Bonferroni adjusted p values were calculated. Lastly, the root-mean-square error (RMSE) of the candidate models derived from leave-one-out cross validation (LOOCV) was compared to check for overfitting.

Shrub Distribution
A binomial GLMM with Laplace Approximation of Maximum Likelihood was used to model presence (yes/no) of shrubs as a function of PFT class and pond age class, with Pond ID as a random factor. Observations in "Intact Shrub" and "Open Water" quadrats were omitted from statistical analysis due to complete presence or absence of living shrubs. The same procedure for model selection as described for Abiotic Data was adopted. Success rates for the prediction of presence of Betula nana were calculated from confusion matrices based on LOOCV.

Shrub Ages
Differences in age among shrub types (living & inside the thaw pond perimeter, living & outside the thaw pond perimeter, and dead & inside the thaw pond perimeter) were assessed using a Poisson GLMM with Laplace Approximation of Maximum Likelihood, using pond age class and shrub type as fixed effects and Pond ID as a random effect. To correct for age differences resulting from the time elapsed since drowning of dead shrubs relative to living shrubs, the GLMM was performed on year of establishment of the shrub, expressed as number of years since the earliest year of establishment found in the data set. The same procedure for model selection as described in paragraph 2.4.2 was adopted.
All analyses were carried out in R version 3.5.1 (R Development Core Team, 2019). The adonis2, betadisper, TukeyHSD, and simper functions from the vegan package (Oksanen et al., 2015) were used for vegetation composition analysis. The lme4 package (Bates et al., 2014) was used for (G)LMMs. For F tests with Kenward-Rogers approximation, the KRmodcomp() function from the pbkrtest package (Halekoh & Højsgaard, 2014) was used. LRTs were implemented using the base R anova() function. The emmeans() function from the emmeans package (Lenth, 2018) was used for pairwise contrasts. Residuals of GLMMs were transformed to standardized space prior to visual inspection and tested for overdispersion, uniformity, and zero-inflation using the DHARMa package (Hartig, 2017).

Selected Ponds
Based on the time series of satellite images, we selected ponds of contrasting age, size, and expansion dynamics. Twelve young ponds were selected, which were assumed to be at most 8 years  (Figure 2). Mean field transect lengths, including terrestrialized areas with sedge and Sphagnum dominated vegetation, were 3.61 ± 1.26 m for Young ponds, 6.41 m ± 2.57for Old Increasing ponds, and 7.79 ± 3.26 m for Old Decreasing ponds. Figure 2a suggests that ponds may reach highly variable open water extent before starting to terrestrialize.

Vegetation Composition
Vegetation composition as derived from transect distances differed among thaw ponds of the three age classes and showed a progression over time from dead shrubs toward sedge and Sphagnum dominated vegetation. Pond age class explained 31% of variance in vegetation composition compared to variance among individual ponds (pond age class: F = 8.18, df = 2, p = 0.001, PERMANOVA). Pairwise PERMANOVAS showed that all pairwise combinations of ponds differed significantly in vegetation composition (Young vs. Old Increasing: F = 5.26,p = 0.012,Young ponds vs. Old Decreasing ponds: F = 17.15,p = 0.003,Old Increasing vs. Old Decreasing: F = 4.43,p = 0.003). Young ponds were characterized by a relatively high cover of dead shrub vegetation, low cover of sedges, and in most cases absence of Sphagnum moss. Old Increasing ponds contained open water more frequently, increased cover of sedges and Sphagnum, and decreased cover of dead shrubs compared to Young ponds. Old Decreasing ponds contained the highest cover of Sphagnum, and low cover of drowning shrubs (Figures 3a and  3b). The difference between Young and Old Increasing ponds was mostly determined by dead shrubs Detailed visual cover estimates in quadrats indicated that dwarf shrub vegetation consisted mainly of Betula nana and occasionally Salix pulchra and Vaccinium vitis-idaea with a dense moss and lichen layer. Sedge vegetation consisted mainly of Eriophorum angustifolium. Sparse sedge cover was also present on Sphagnum dominated sites and especially Carex aquatilis was often found in association with Sphagnum (Table S21, Figure S3). Sphagnum species were almost exclusively S. obtusum and S. squarrosum, with incidental presence of S. subsecundum and S. balticum (Magnússon, pers. obs.). Sedges and Sphagnum were generally absent in intact shrub quadrats (Table S21).  Table S22 for images of individual ponds per year). Y = Young ponds; OI = Old Increasing ponds; and OD = Old Decreasing ponds. Shaded areas represent 95% confidence interval per pond age class. (b) Observed lateral expansion rates of open water area per pond age class. One decreasing pond displayed a positive lateral expansion while its net area decreased, due to development of a complex shape.

Shrub Distribution and Ages
Analysis of shrub distribution and ages indicated that young Betula nana individuals are found in thaw ponds, predominantly on Sphagnum carpets. Out of 66 shrubs sampled, 43 were cross-dated with sufficient statistical certainty (see Text S1) to determine their age and the year of establishment. Correlation of ring-width series among these 43 shrubs was strong (average correlation coefficient = 0.52 ± 0.10). Ages of living shrubs from outside the pond perimeter and dead shrubs inside the pond perimeter ranged from 30 to 84 years (mean = 46.23 ± 13.00, n = 30) and did not significantly differ (p = 0.058), suggesting that pond formation occurred in intact shrub vegetation. Living shrubs from inside the pond perimeter ranged from 2 to 18 years in age (mean = 7.83 ± 4.76, n = 13) and hence were significantly younger than dead shrubs or living shrubs from outside the pond perimeter (p < 0.001 for both contrasts), suggesting they established after pond formation ( Figure 4a). As expected, dead shrubs in Young Ponds generally died more recently than those in Older ponds (Figure 4a). See Tables S9-S10 for model specifications and pairwise contrasts.
Analysis of shrub distribution (recorded per quadrat) showed that living dwarf shrubs inside the pond perimeter were predominantly found in quadrats dominated by Sphagnum, and occasionally also in quadrats dominated by sedges and dead shrubs (Figure 4b). Pairwise contrasts indicated that they were more likely to be found on Sphagnum than on sedge (diff. = 0.46 ± 0.18, p = 0.035), other contrasts were not significant. See Tables S6-S8 for model specifications, RMSEs, and pairwise contrasts. Where present, dwarf shrubs had

Abiotic Conditions
We found that abiotic conditions in thaw ponds are significantly related to the dominant vegetation type (PFT class) in the pond. When we ordered dominant vegetation types according to their association with the three pond age classes (Figure 3), this resulted in a clear pattern of permafrost degradation and subsequent recovery ( Figure 5). Pond age class explained only a small additional amount of variance in TD, WT, and PP. Relative Surface Elevation (RSE) and PP decreased progressively from intact shrub to dead shrub to open water (Figures 5a and 5g), coinciding with an increase in WT and TD (Figures 5c and 5e). Under Sphagnum vegetation, and to some extent sedges, RSE, PP, and Organic Layer Thickness (OLT) increased again, whereas TD and WT decreased. No large differences were observed along the pond age class gradient (from Young to Old Increasing to Old Decreasing), however, in Old Increasing ponds PP was lower ( Figure 5h) and WT and TD were higher (Figures 5d and 5f)   Tables S11-S18 contain model specifications, RMSEs, and pairwise contrasts for all abiotic variables.

Discussion
The aim of this study was to determine whether thaw ponds display quasi-cyclic successional vegetation development linked to permafrost degradation and recovery and to provide a preliminary estimate of the

Journal of Geophysical Research: Biogeosciences
timescales associated with thaw pond formation and recovery. Overall, our results from three lines of evidence (vegetation composition in thaw ponds of various age classes, abiotic conditions related to various PFT classes, and dendrochronology on dwarf shrubs) suggest that formation and terrestrialization of thaw ponds can occur rapidly in shrub-dominated tundra. Colonization by sedges was already observed in young ponds with ages of up to 8 years. Though a vegetation succession trajectory toward Sphagnum dominated vegetation in older, decreasing ponds was observed and recruitment of shrubs on Sphagnum carpets was demonstrated, it remains uncertain whether the previously dominant shrub vegetation can fully reestablish in thaw ponds and on what timescale. In the paragraphs below, we will elaborate on this interpretation of our findings.

Permafrost Degradation
Our results suggest that thaw ponds may form and develop rapidly in the north-eastern Siberian lowland tundra. Young thaw ponds (up to 8 years in age, see Table S22) were found to be mostly composed of dead shrubs (Figure 3). Presence of open water (Figure 3) and associated deeper thaw depths ( Figure 5) were most common in older but still expanding ponds. Interestingly, we found that several of the young ponds already started declining again in open water extent between 2015 and 2018 ( Figure 2, Table S22). This suggests that complete degradation (the transformation from intact shrub vegetation to open water with dead shrubs) can occur within a time span of under 8 years but may continue toward more advanced degradation on longer timescales. Similarly, formation of ponds was observed by Li et al. (2017) and Nauta et al. (2015) within 5 years following removal of shrubs. However, some degree of interannual variability in the open water extent of ponds cannot be ruled out, as the water level in thaw ponds also depends on specific weather conditions (Campbell et al., 2018;Jones et al., 2009;Plug et al., 2008;Wolfe et al., 2011). We found elevation differences in thaw ponds of around 50 cm compared to intact tundra vegetation, which is relatively low compared to elevation differences reported for thermokarst troughs (i.e., ice wedge degradation) by Kanevskiy et al. (2017), ranging from 50 to 110 cm. Li et al. (2017) found an even smaller elevation difference of up to 30 cm within a time span of 8 years after thaw pond formation following removal of shrubs. Existing data on the time span of degradation in the case of ice-wedges indicate that for 50 cm soil subsidence to occur, about 10 years would be required (Jorgenson et al., 2015;Kanevskiy et al., 2017). In the case of ice-wedge degradation, such initial degradation may continue into advanced degradation on timescales of up to hundreds of years with soil subsidence of over 3 m (Kanevskiy et al., 2017). We did not find evidence of such magnitudes of subsidence in thaw ponds, indicating that rates and magnitudes of subsidence associated with thaw pond formation may be comparable to shallow or initial forms of ice wedge degradation. Lateral expansion rates of ponds suggest that initial expansion of Young thaw ponds (mean = 0.52 m/year) exceeds that of Old Increasing ponds (mean = 0.33 m/year). For a mean observed transect diameter of 7.79 m for Old Decreasing ponds (including already terrestrialized areas), a time span of pond formation of approximately 20 years would be required based on observed pond lateral expansion rates. However, ponds were observed to attain highly variable sizes before starting to terrestrialize (Figure 2). Observed variability in soil subsidence and degradation rate may be attributed to different initial ground ice content. The amount and configuration of ground ice may also differ between the observed shallow, irregular shaped thaw ponds (ice lenses) and thermokarst troughs (ice wedges) (Van Huissteden, 2020). Interestingly, dead shrubs that died recently still occur in older decreasing ponds (Figure 4a). This confirms our observation that thaw ponds may not grow and decline uniformly (Table S22). Overall, our results for thaw pond degradation in this lowland tundra ecosystem correspond well with reported rates of shallow ice-wedge degradation, indicating that they form on relatively short timescales compared to thermokarst troughs. However, this claim would need to be validated by longer term monitoring.
Dendrochronology on dead shrubs was used to more accurately identify the onset of thaw pond formation.
Evaluating the year of death of dead shrubs, it became evident that in many cases, timing of thaw pond formation (as derived from satellite imagery, e.g., before 2010 or before 2015) predated the year of drowning of shrubs. This indicates that dead shrubs identified in older thaw ponds in this study likely represent more recent expansions of older thaw ponds. It is likely that the oldest parts of an older thaw pond are currently overlain by moss carpets, overgrowing remains of dead shrubs. Therefore, the years of death of the sampled dead shrubs may not represent the earliest onset of permafrost degradation, unless woody material is preserved and can be recovered from the oldest parts of a thaw pond. In addition, Arctic dwarf shrubs,

10.1029/2019JG005618
Journal of Geophysical Research: Biogeosciences especially in older individuals, may show extremely narrow or completely absent rings due to reduced growth activity (Hallinger et al., 2010;Wilmking et al., 2012). Our study demonstrates that cross-dating of dead Betula nana shrubs is feasible. However, the detected year of death of dead shrubs should be considered a minimum estimate.

Aquatic Vegetation Succession in Thaw Ponds
Results from vegetation composition analysis confirm the hypothesized quasi-cyclic vegetation succession in our shallow thaw ponds. Older ponds are colonized by wetland vegetation, most notably E. angustifolium and Sphagnum mosses (Figure 3, Table S21). Sphagnum mosses, in particular, were clearly associated with shrinking ponds (Figure 3). Earlier studies confirm this association between encroaching vegetation and a decline in pond or lake area (Andresen & Lougheed, 2015;Frost et al., 2018;Jorgenson et al., 2015). We found mean lateral shrinkage rates of Old Decreasing ponds of 0.43 m/year, indicating that terrestrialization of ponds occurs at similar rates as their formation (an estimated 20 years based on observed transect diameters). This suggests a rapid transition of thaw ponds toward bog vegetation, despite highly variable shrinkage rates (Figure 2b). This observed succession and associated abiotic development correspond well to the initial degradation (dead shrubs), advanced degradation (open water), initial stabilization (sedges), and advanced stabilization (Sphagnum) phases reported by Jorgenson et al. (2015) and Kanevskiy et al. (2017) for ice-wedge degradation, and with vegetation successions derived from peat stratigraphy in boreal thermokarst ecosystems (Kuhry et al., 1993;Myers-Smith et al., 2008;Payette et al., 2004;Robinson & Moore, 2000;Turetsky et al., 2007).
Observed sedge vegetation consisted mainly of Eriophorum angustifolium (Table S21). E. angustifolium is known as a rapid colonizer of disturbed, waterlogged sites (Phillips, 1954). Occasional presence of tussock sedge E. vaginatum in ponds (Table S21) may be interpreted as a relic of pre-thaw vegetation. Seven out of 12 of our young thaw ponds contained varying amounts of E. angustifolium dominated sedge vegetation ( Figure 3, Table S21), which indicates that initial colonization by sedges can occur within a time span of 8 years or less (the time difference between the earliest satellite images and the field measurements). In the case of ice-wedge degradation, such initial colonization and associated increase in surface elevation in a thermokarst trough may occur within a decade after degradation (Jorgenson et al., 2015). Li et al. (2017) found establishment of sedges in thaw ponds resulting from the removal of B. nana shrubs within several years of thaw pond initiation. Additionally, narrow zones of grass vegetation, predominantly Arctagrostis latifolia, were found in Young and Old Increasing thaw ponds, often in association with dead shrubs and subsiding soil. This indicates that grasses may represent a short interval in the thaw pond vegetation succession. The early presence of grasses can be explained by their ability to utilize nutrients from recently thawed permafrost through deeper roots , while their lower tolerance to water explains the temporary character of the grass vegetation.
Sphagnum was most clearly associated with old decreasing ponds (Figure 3). Observed Sphagnum species in our study site are almost exclusively S. obtusum and S. squarrosum (Magnússon, pers. obs., see also Iturrate-Garcia et al., 2016). Both species are associated with relatively wet and mesotrophic to eutrophic sites (Daniels & Eddy, 1990). Occasional presence of S. subsecundum and S. balticum was observed mainly in old decreasing ponds (Magnússon, pers. obs.), which may indicate a shift toward mesotrophic to oligotrophic conditions (Daniels & Eddy, 1990;Kuhry et al., 1993). Earlier studies of peat stratigraphy in boreal thermokarst bogs indicate distinct succession of Sphagnum species associated with wetter conditions such as S. angustifolium and S. riparium upon degradation toward increasingly oligotrophic species associated with drier conditions such as S. fuscum (Camill et al., 2009;Kuhry et al., 1993;Robinson & Moore, 2000;Zoltai, 1993) and subsequent reestablishment of woody vegetation (Camill et al., 2009;Robinson & Moore, 2000;Routh et al., 2014;Treat et al., 2016;Zoltai, 1993). In our study region, Sphagnum species associated with oligotrophic, drier sites such S. lenense and S. rubellum (Daniels & Eddy, 1990) may be found in polygonal tundra (Figure 1e) but were not observed in the shallow, nonpolygonal thaw ponds analyzed in this study (Magnússon, pers. obs.). This suggests that the studied thaw ponds are in an early stage of Sphagnum bog succession. Our results correspond with an observed trend of replacement of woody vegetation by Sphagnum vegetation upon permafrost degradation in boreal ecosystems (Jones et al., 2017;Jorgenson et al., 2010;Myers-Smith et al., 2008). The relatively high summer temperatures (paragraph 2.1) of the study area due to a strongly continental climate and the relatively shallow nature of the studied 10.1029/2019JG005618

Journal of Geophysical Research: Biogeosciences
ponds could explain the rapid terrestrialization of ponds in this ecosystem and similarity to processes described in boreal ecosystems. High summer temperatures have been associated with enhanced vegetation encroachment in thaw lakes (Roach et al., 2011) and succession from sedge to Sphagnum dominated vegetation in thermokarst bogs (Myers-Smith et al., 2008). In addition, high summer temperatures may enhance evaporation in shallow, isolated ponds (Campbell et al., 2018;Wolfe et al., 2011), further promoting vegetation encroachment through a lowering of the water table (Payette et al., 2004;Roach et al., 2011).

Permafrost Recovery
Our results suggest that vegetation succession in thaw ponds is strongly related to permafrost dynamics. Sites characterized by colonizing wetland vegetation showed an increase in surface elevation and decrease in thaw depth and water table relative to sites with dead shrubs and open water, indicating initial recovery of permafrost. This effect was more pronounced under Sphagnum than under sedges. Differences in abiotic conditions among PFTs have also been described in polygonal tundra in the vicinity of the Kytalyk reserve by Minke et al. (2009) andDonner et al. (2007). Although the presence of PFTs can be interpreted as the result of abiotic gradients (Donner et al., 2007), changes in surface elevation and permafrost position can also be a result of vegetation development and peat formation (Minke et al., 2009;Shur et al., 2011). For instance, Sphagnum growing above the water remains dry during the summer, providing insulation to the permafrost below (Seppälä, 1988). Low thermal conductivity of dry Sphagnum mosses was demonstrated by Van Breemen (1995) and Soudzilovskaia et al. (2013). This improved insulation can lead to refreezing of thawed soil, as well as aggradation of ground ice resulting in frost heave (Shur et al., 2011). Frost heave further contributes to terrain rise above the water table and thereby to insulation (Seppälä, 1988). Our results confirm the tight coupling between permafrost status and vegetation cover and indicate that vegetation succession, most notably toward a Sphagnum dominated stage, may play a key role in recovery of permafrost and microtopography. However, manipulation experiments would be necessary to demonstrate causal effects. The proposed thaw pond succession mechanism is visualized in Figure 6. . Schematic representation of our proposed mechanism of thaw pond formation (degradation) followed by vegetation succession (recovery). Through melting of ice-rich permafrost, local depressions are formed in which the initially dominant dwarf shrub vegetation drowns due to waterlogging. This is associated with soil subsidence, increased thaw depths, and formation of open water. The open water can be colonized by water-tolerant species such as sedges, which is associated with slight decreases in water table and thaw depth and slight increases in surface elevation. Colonization by Sphagnum mosses is associated with larger decreases in water table and thaw depth, and an increase of the surface elevation, resulting from organic matter accumulation and permafrost recovery. Sphagnum carpets can provide a substrate for reestablishment of dwarf shrubs, although it is presently unclear whether this initial reestablishment also leads to complete reestablishment of the former shrub-dominated vegetation.

Journal of Geophysical Research: Biogeosciences
We found changes in permafrost position of +18 cm under sedge vegetation and +32 cm under Sphagnum vegetation, although the permafrost position remains significantly deeper than under intact shrub vegetation (Figure 5g). Quantification of ice content would be a valuable topic of future research to determine the relative influence of ground ice aggradation and frost heave as opposed to refreezing of previously thawed soil during successional development. Comparison of the increase in the relative surface elevation over the successional gradient from open water to Sphagnum (+23 cm) and the increase in organic layer thickness (+6 cm) would suggest that peat accumulation constitutes a relatively small component of surface heave, and that frost heave due to ice aggradation below the organic layer may play an important role in the recovery of the terrain. This is corroborated by Wang et al. (2019), who observed a positive relation between ground ice content in the top 20 cm of permafrost and moss layer thickness and relative surface elevation in the Kytalyk lakebed site. An observed decrease in thaw depth (−9 cm) over the successional gradient from open water to Sphagnum, despite the increase in organic layer thickness, could suggest the formation of an ice rich intermediate permafrost layer (Shur et al., 2011). Detailed soil stratigraphical data may provide clearer insights into the relative role of peat accumulation, refreezing and ice aggradation.
The recovery rates observed in our shallow ponds seem high compared to those reported for thermokarst ponds formed after ice-wedge degradation. Complete recovery of permafrost after ice-wedge degradation has been reported to take tens to hundreds of years (Jorgenson et al., 2015;Kanevskiy et al., 2017). The fast recovery rates observed in our study may be related to the limited subsidence compared to thermokarst troughs (see paragraph 4.1). Early colonization by sedges and associated changes in abiotic conditions indicate that initial recovery of permafrost can already be initiated within subdecadal timescales. Shrinking rates of older decreasing ponds suggest that complete terrestrialization of ponds can occur within 20 years (see paragraph 4.2). However, timescales associated with the full recovery of the permafrost underneath, if this occurs at all, remain uncertain. In addition, pond dynamics depend on pond characteristics such as size and depth and hydrological influences such as drainage, precipitation, and evaporation (Campbell et al., 2018;Jones et al., 2009;Plug et al., 2008;Wolfe et al., 2011). Small, shallow ponds in particular may be prone to evaporative loss of surface water in a warming climate (Campbell et al., 2018;Wolfe et al., 2011). Lastly, it is unclear to what extent degradation or recovery of permafrost in thaw ponds is entirely cyclic. Kanevskiy et al. (2017) propose the use of the term "quasi-cyclic" development for ice-wedges, since recovery of ice-wedges often leads to a complete rearrangement of ground ice patterns and accumulation of organic matter which in turn affects the ice-wedges' susceptibility to future thaw. It is unclear whether such rearrangements also occur in thaw ponds.

Reestablishment of Dwarf Shrub Vegetation?
Our study is the first to explicitly look at reestablishment of dwarf shrubs in thaw ponds. At our study site, the living shrubs (both established from seeds or from ramets) found within thaw ponds were significantly younger than living shrubs outside of the thaw pond in intact shrub vegetation. This indicates that recruitment of dwarf shrubs, both through sexual reproduction and clonally, occurs in thaw ponds. Similarly, Huebner and Bret-Harte (2019) found increased shrub recruitment on bare soil in retrogressive thaw slumps compared to undisturbed tundra vegetation, which was attributed to reduced competition by surrounding vegetation and creation of suitable conditions for germination. However, whereas Huebner and Bret-Harte (2019) found the highest recruitment densities on more recently disturbed sites (1-10 years since disturbance), we found the majority of recruitment on Sphagnum carpets in older, decreasing thaw ponds. This may be explained by the ability of Sphagnum species to grow above the water table, generating the moist but unsaturated substrate needed for recruitment of dwarf shrubs (Bell & Bliss, 1980;Huebner & Bret-Harte, 2019). Frost heave may contribute substantially to the generation of unsaturated conditions (Seppälä, 1988).
It remains unclear, however, to what extent this recruitment can lead to a full recovery of the original shrub-dominated vegetation. The year of establishment of recruiting Betula nana within the thaw pond perimeter provides insight into the timescales associated with this potential recovery of shrubs. Most recruiting shrubs were up to 9 years in age, with outliers of up to 18 years. This age range (2-18 years) of recruiting shrubs and their low cover compared to intact shrub quadrats (Table S21) indicate that a potential full recovery likely may require multidecadal timescales after terrestrialization, although no reliable estimate can be made based on available data. It may also indicate an underrepresentation or absence of 10.1029/2019JG005618 Journal of Geophysical Research: Biogeosciences sites in an advanced state of vegetation recovery. Additionally, recruitment rates of most Arctic species, including Betula nana, may strongly depend on climatic pulses, and mortality among recruiting shrubs is high (Büntgen et al., 2015). This indicates that presence of young recruitment does not necessarily mean successful recolonization of a thaw pond. Often, the year of establishment of recruiting Betula nana predates the death of dead shrubs from another site within the same pond, indicating again that both processes may occur simultaneously within the same thaw pond. Further tree-ring analysis of recruiting shrubs in thaw ponds may provide clearer insight into the timescales associated with the potential reestablishment of dwarf shrubs.

Implications for the Arctic Greening Trend and Carbon Balance
The balance between formation and recovery of thaw ponds on a landscape scale may have profound implications for the vegetation composition and GHG balance of lowland tundra ecosystems. Our results seem to suggest that thaw pond formation may occur very rapidly, whereas no sites in advanced state of recovery of the original dwarf shrub vegetation were observed. Such an imbalance could provide a potential explanation for local tundra browning observed around the study region by Epstein et al. (2018). Formation of open water features and subsequent death of shrub vegetation (Li et al., 2017) is associated with decreased NDVI . As such, thermokarst has been suggested as a driver of Arctic browning (Phoenix & Bjerke, 2016;. Conversely, colonization of a pond by aquatic vegetation may cause a local increase in NDVI (Zona et al., 2010) or other greenness indices such as Tasseled Cap Components . Thus, the balance between thaw pond formation and recovery over a larger spatial scale may affect Arctic greening/browning dynamics.
The balance between permafrost degradation and recovery also affects the carbon balance of the lowland tundra ecosystem. Earlier research at the Kytalyk site has indicated that thaw pond formation can shift the lowland tundra ecosystem from a methane sink to a methane source (Nauta et al., 2015). Methane emissions were found to be higher for Eriophorum dominated sites than for drier, Sphagnum dominated sites (Van Huissteden et al., 2005. Eriophorum spp. and Carex spp. enhance methane emissions through gas exchange via aerenchyma (Ström et al., 2005(Ström et al., , 2012. Conversely, Sphagnum may reduce methane emissions directly through methane oxidation facilitated by symbiosis with endophytic methane-consuming bacteria and presence of aerated surface peat Sundh et al., 1995). Presence of Sphagnum above the water table may further reduce GHG emissions by facilitating permafrost recovery through its insulation capacity (Seppälä, 1988;Soudzilovskaia et al., 2013). Lastly, Sphagnum may act as an ecosystem engineer during peat growth, enhancing soil carbon stocks (Myers-Smith et al., 2008;Robinson & Moore, 2000;Siewert et al., 2015;Van Breemen, 1995). Increasing organic layer thickness from open water to sedge to Sphagnum mosses ( Figure 5i) suggests that peat accumulation occurs during the observed succession. While there is considerable debate about the carbon source/sink dynamics of thermokarst features over the long term (Robinson & Moore, 2000;Turetsky et al., 2007;Van Huissteden & Dolman, 2012), studies from boreal peatlands suggest that a full recovery of pre-thaw carbon stocks through peat accumulation may require centuries (Jones et al., 2017;Turetsky et al., 2007). Lastly, pond dynamics and associated GHG dynamics may depend on weather patterns determining precipitation and evapotranspiration (Jones et al., 2009;Plug et al., 2008) and drainage processes (Gorham, 1991;Regmi et al., 2012;Zona et al., 2010). This may also explain the observed variability of expansion and terrestrialization rates of the studied ponds across different time intervals (Table S22). Overall, thaw pond formation results in the release of methane, whereas advanced vegetation succession, especially colonization by peat-building Sphagnum mosses, is associated with lower methane emissions and net carbon storage.
To estimate the net contribution of thaw ponds to GHG emissions over the course of their development it is necessary to (i) determine net amounts of carbon emission and uptake associated with permafrost degradation and recovery, respectively (Van Huissteden & Dolman, 2012); (ii) determine whether there is a net trend of thaw pond formation and associated shrub decline or recovery and associated shrub expansion; and lastly, (iii) how the rates of these two processes relate to climate. Additionally, efforts should be made to determine exactly how abundant small-scale thaw ponds are across the Arctic since small waterbodies tend to be overlooked in remote sensing inventories (Grosse et al., 2008;Raynolds et al., 2019).

Conclusion
Our study is the first to characterize vegetation successions in thaw ponds in the Siberian Arctic. In this lowland tundra ecosystem, thaw ponds can develop in shrub-dominated lowland tundra within subdecadal to decadal timescales, resulting in depressions with dead shrubs, standing water, and thawing permafrost. Early colonizers such as sedges can establish within several years of formation and complete terrestrialization may occur within 20 years. Especially colonization by Sphagnum moss is associated with permafrost recovery, and Sphagnum carpets are important microsites for the reestablishment of dwarf shrubs. The succession mechanisms identified in shallow, irregular shaped thaw ponds in the Siberian lowland tundra correspond well with earlier work on ice-wedge degradation and boreal thermokarst bogs, although degradation and initial colonization seem to occur on relatively short timescales and thaw ponds show relatively shallow subsidence compared to ice-wedge degradation. Main challenges for future research will be to assess whether the relative rates of thaw pond formation and recovery will result in a net trend of tundra wetting and to determine the net GHG emission of thaw ponds over the course of their formation and recovery. The balance of tundra greening and wetting will have profound implications for the contribution of Arctic lowland tundra areas to GHG emissions.