Keywords: permafrost, paleolimnology, greenhouse gases, remote sensing

In Central Yakutia (Eastern Siberia, Russia), persistent freezing temperatures, flat topography, and a lack of glaciation during the Last Glacial Maximum (~20,000 years ago) have created optimal conditions for the proliferation of deep, ice-rich permafrost. Past and present climate change and other disturbances can have significant effects on permafrost landscapes and the global carbon budget. In areas of ice rich permafrost, degradation can result in the formation of thermokarst (thaw) lakes, which are hotspots for biological activity and greenhouse gas (GHG) emissions. The three axes of this thesis combine provide insights into 1) past lake dynamics, 2) spatial and temporal heterogeneity in GHG concentrations and diffusive fluxes from thermokarst lakes, and 3) large scale analysis of recent (since 1960s) lake dynamics in response to local climate trends and human activities.

Titre : Reconstitutions paléoenvironnementales et caractérisation des gaz à effet de serre dans les systèmes aquatiques pergélisolés de la Yakoutie centrale (Sibérie) Mots clés : pergélisol, paléolimnologie, gaz à effet de serre, télédétection Résumé : Yakoutie centrale se caractérise par un pergélisol profond et riche en glace. Les changements climatiques passés et présents et d'autres perturbations peuvent avoir des effets importants sur les paysages de pergélisol et le bilan carbone mondial. Dans les zones de pergélisol riche en glace, la dégradation peut résulter en la formation de lacs thermokarst, qui sont des 'hotspots' pour l'activité biologique et les émissions de gaz à effet de serre (GES). Les trois axes de cette thèse sont 1) la dynamique passée des lacs, 2) l'hétérogénéité spatiale et temporelle des concentrations de GES et des flux diffusifs des lacs thermokarst, et 3) l'analyse à grande échelle de la dynamique récente (depuis les années 1960) des lacs en réponse aux variations climatiques locales et aux activités humaines.

Pour comprendre le développement des lacs, l'accumulation de sédiments et de carbone organique, et les changements de la productivité primaire, dans le contexte des changements climatiques du Pléistocène final et de l'Holocène, une analyse paléolimnologique multi-proxy d'une séquence de carottes de sédiments du lac Malaya Chabyda en Yakoutie centrale a été réalisée (max âge ~14 cal kBP). Les propriétés sédimentologiques et biogéochimiques dans la section la plus profonde de la carotte suggèrent un environnement lacustre principalement influencé par la végétation terrestre, où l'accumulation de carbone organique a pu être relativement faible, beaucoup plus élevée que la moyenne mondiale moderne. La section médiane de la carotte a été caractérisée par une productivité primaire plus élevée dans le lac, une sédimentation beaucoup plus forte et une grande augmentation de la livraison de carbone organique (OC). Les conditions dans la section supérieure de la carotte suggèrent une productivité primaire élevée dans le lac et des taux élevés d'accumulation de CO, avec des conditions environnementales stables.

Les changements dans le caractère et la quantité de l'accumulation de CO peuvent avoir des implications importantes pour la libération future de GES. Ces résultats aident à comprendre l'histoire du développement des lacs dans les paysages de pergélisol dans le contexte des changements climatiques passés.

Pour analyser l'hétérogénéité spatiale et temporelle des concentrations de GES et du flux diffus, les concentrations de CO2 et de CH4 dissous dans les lacs thermokarst d'un site d'étude en Yakoutie centrale ont été mesurées sur quatre saisons. Les lacs formés au cours de l'Holocène (lacs alas) sont comparés aux lacs qui se sont développés au cours des dernières décennies. Les résultats montrent des différences dans les gaz à effet de serre dissous entre les types de lacs et les saisons. Les lacs peu profonds situés dans des dépressions d'alas hydrologiquement fermées ont agi comme des réservoirs de CO2 et de fortes sources de CH4 diffusif pendant certaines saisons. Les lacs thermokarstiques récents étaient des sources extrêmement élevées de CO2 et de CH4 diffusifs, avec une accumulation considérable de gaz à effet de serre sous la couverture de glace (hiver) ou dans de l'eau les plus profondes (été). Les flux diffusifs mesurés dans les lacs thermokarst de ce paysage typique de taïga alas de la Yakoutie centrale sont parmi les plus élevés présentés dans les régions arctiques et subarctiques.

Enfin, une analyse de télédétection a été menée pour étudier le développement récent (depuis les années 1960) des lacs et quantifier les différences dans l'histoire du développement entre les types de lacs. Cette étude a indiqué que certains types de lacs sont plus sensibles aux changements annuels des précipitations. Ces informations seront utilisées pour étendre les mesures in situ des gaz à effet de serre à une plus grande unité de paysage.

To understand lake development, sediment and organic carbon accumulation, and changes in primary productivity, within the context of Late Pleistocene and Holocene climate change, a multi-proxy paleolimnological analysis of a sediment core sequence from Lake Malaya Chabyda in Central Yakutia (Eastern Siberia, Russia) was conducted. Age-depth modeling with 14 C indicates that the maximum age of the sediment core is ~14 cal kBP. Three distinct sedimentary units were identified within the sediment core. Sedimentological and biogeochemical properties in the deepest section of the core suggests a lake environment mostly influenced by terrestrial vegetation, where organic carbon accumulation might have been relatively low, although high compared to the global modern average. The middle section of the core was characterized by higher primary productivity in the lake, much higher sedimentation, and a strong increase in organic carbon (OC) delivery. Conditions in the upper section of the core (< 376 cm; < 9.0 cal kBP) suggest high primary productivity in the lake and moderate OC accumulation rates, with stable environmental conditions.

Changes in the character and quantity of OC accumulation can have important implications for future GHG release. These results help to understand the developmental history of lakes in permafrost landscapes within the context of past climate change.

To analyze spatial and temporal heterogeneity in GHG concentrations and diffusive flux, dissolved CO 2 and CH 4 concentrations in thermokarst lakes in a study site in Central Yakutia were measured over four seasons. Lakes formed over the Holocene (alas lakes) are compared to lakes that developed over recent decades. The results show striking differences in dissolved greenhouse gases (up to two orders of magnitude) between lake types and seasons. Shallow lakes located in hydrologically closed alas depressions acted as CO2 sinks and strong sources of diffusive CH4 during some seasons. Recent thermokarst lakes were moderate to extremely high sources of diffusive CO2 and CH4, with considerable accumulation of greenhouse gas under the ice cover (winter) or in the deepest water layers (summer). The diffusive fluxes measured from thermokarst lakes of this typical taiga alas landscape of Central Yakutia are among the highest presented across Arctic and subarctic regions.

Lastly, a remote sensing analysis was conducted to investigate recent (since 1960s) lake development and quantify differences in developmental history between lake types. This study indicated that certain lake types (unconnected alas lakes) are more sensitive to annual changes in precipitation. This information will be used to upscale in situ greenhouse gas measurements to a larger landscape unit. This thesis combines paleolimnological analysis, in situ observations, and remote sensing analysis to provide a comprehensive understanding of lake development in permafrost landscapes and the contribution of these lakes to the global carbon budget.

5

Colombe de la paix. Pablo Picasso 1949

Résumé en français

Unconnected alas lake in May 2019 near Syrdakh Village (Sakha Republic, Russia) Cette thèse combine des analyses paléolimnologiques, des observations et mesures de GES in situ et des analyses de télédétection afin de fournir une compréhension complète du développement des lacs dans les paysages de pergélisol et de la contribution de ces lacs au budget global du carbone. Elle consiste en une introduction générale (premier chapitre), suivie de trois chapitres concernant les principaux axes méthodologiques de la thèse. Deux de ces chapitres sont publiés sous la forme d'articles de revue à comité de lecture, et un troisième sera soumis sous peu.

Le pergélisol occupe plus de 20 millions de kilomètres carrés et représente 24 % de la couverture terrestre au sein de l'hémisphère nord (Brown et al., 1997 ;Obu et al., 2019). Il est particulièrement abondant en Sibérie, en Alaska et au Canada, et son étendue spatiale, son épaisseur et sa teneur en glace dans le sol peuvent varier considérablement d'un paysage à l'autre (Grosse et al. 2013 ;Strauss et al. 2017). Selon les conditions environnementales environnantes passées et présentes, le pergélisol peut s'étendre sous la surface du sol de quelques mètres à 1 500 mètres. La plupart du pergélisol qui existe aujourd'hui s'est formé pendant les périodes glaciaires froides passées (depuis ~100 000 ans) (Vonk et al., 2013b). Le pergélisol continu et riche en glace formé à la fin du Pléistocène (~ 20,000 à 50,000 ans avant aujourd'hui) que l'on trouve en Sibérie orientale, dans certaines parties de l'Alaska et à l'ouest du Yukon (Canada) est appelé " Yedoma " (Schirrmeister et al., 2013 ;Strauss et al., 2017). Le pergélisol de type Yedoma est généralement plus ancien que le pergélisol des régions englacées au dernier maximum glaciaire (ex. Canada, Europe du Nord, Sibérie occidentale).

Les études sur la dynamique des paysages de pergélisol au cours de l'Holocène (~ les derniers 10,000 ans) et de ces dernières décennies ont montré que les zones dominées par un pergélisol riche en glace sont très sensibles aux changements de température et à d'autres perturbations locales (Grosse et al. 2013 ;Ulrich et al. 2019). L'augmentation des températures de l'air et du sol ainsi que les activités humaines telles que la construction de routes ou l'agriculture peuvent entraîner le dégel du pergélisol, ce qui peut avoir des répercussions étendues sur l'hydrologie locale et régionale (Biskaborn et al. 2019). Le dégel du pergélisol peut libérer des quantités substantielles de matières organiques et minérales dans les écosystèmes aquatiques, entraînant de profonds changements dans leur biogéochimie et leur rôle dans le cycle mondial du carbone (Vonk et al., 2015;Opfergelt, 2020). On estime que les sols des paysages de pergélisol nordique contiennent deux fois plus de carbone que l'atmosphère (Hugelius et al. 2014) et sa libération pourrait avoir un impact significatif sur le cycle mondial du carbone et le climate. Le taux et l'ampleur actuels de l'augmentation de la température dans l'Arctique sont disproportionnés par rapport aux moyennes mondiales, la température annuelle moyenne de l'air devant augmenter de 5,4 °C au cours du prochain siècle en l'absence d'un effort mondial significatif et dirigé pour réduire les émissions de gaz à effet de serre [START_REF] Ipcc ; H.-O. Pörtner | IPCC Special Report on the Ocean and Cryosphere in a Changing Climate[END_REF]. Cette évolution annoncera une période de changements dynamiques dans les paysages de pergélisol.

L'une des principales voies de dégradation du pergélisol dans les zones riches en glace est le déclenchement de processus thermokarstiques, qui aboutissent finalement à la formation de nombreux lacs dans les régions où la topographie est plate (Grosse et al. 2013). Les processus thermokarstiques commencent lorsque des perturbations telles que le réchauffement, destruction par l'abattage, ou par le feu des forêts provoquent un approfondissement de la couche active (la couche qui dégèle chaque été), ce qui entraîne la subsidence et l'érosion thermique de la surface du sol.

Une fois formés, les lacs thermokarstiques modifient profondément le régime thermique local du sol, augmentant parfois les températures des sédiments environnants jusqu'à 10°C au-dessus de la température annuelle moyenne de l'air (Brouchkov et al. 2004). Les processus thermokarstiques, en particulier la présence d'un lac ou d'un étang, augmentent considérablement le taux de dégel du pergélisol par rapport à ce qui serait prédit par les seules augmentations de la température de l'air (Brouchkov et al. 2004 ;Schuur et al. 2015). La poursuite de l'expansion latérale du lac peut provoquer l'affaissement des pentes (" effondrements de dégel rétrogressifs " ou 'retrogressive thaw slumps') et une subsidence et une érosion supplémentaires de la surface du sol (Séjourné et al. 2015 ; Bouchard et al. 2017). L'expansion et l'approfondissement se poursuivent jusqu'à ce que la profondeur du lac dépasse la profondeur du pergélisol riche en glace ou que la propagation du dégel sous le lac soit entravée par les propriétés isolantes des sédiments lacustres, ce qui met fin à l'approfondissement. Après cette phase, la taille et la profondeur du lac sont contrôlées principalement par le bilan précipitations-évaporation [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]French 2017 ;Grosse et al. 2013).

Plus précisément, les lacs thermokarstiques sont considérés comme des 'points chaud' biogéochimiques qui jouent un rôle important dans la transformation de la matière organique du pergélisol et donc sur le climat mondial (Walter Anthony et al. 2016). Les sédiments anoxiques au fond de ces lacs sont des sites de production de CH4 qui peut être libéré par ébullition et/ou diffusion. Une libération importante de CO2 à partir des lacs thermokarstiques des hautes latitudes a également été documentée [START_REF] Abnizova | Small ponds with major impact: The relevance of ponds and lakes in permafrost landscapes to carbon dioxide emissions[END_REF]Hughes-Allen et al., 2021). La première partie de la thèse concerne la compréhension du développement des lacs, de l'accumulation de sédiments et de carbone organique, et des changements de la productivité primaire, dans le contexte des changements climatiques du Pléistocène final et de l'Holocène. Les méthodes étaient centrées sur une analyse paléolimnologique multi-proxy d'une séquence de carottes de sédiments du lac Malaya Chabyda en Yakoutie centrale (âge maximum ~14 cal kBP). Trois unités sédimentaires distinctes ont été identifiées dans la carotte sédimentaire. Les propriétés sédimentologiques et biogéochimiques dans la section la plus profonde de la carotte suggèrent un environnement lacustre où l'accumulation de carbone organique (CO) a pu être relativement faible (moyenne ~ 100 g CO m -2 a -1 ), quoique plus élevée que la moyenne mondiale moderne. Cette CO est venu principalement par la végétation terrestre, plutôt que par la production primaire dans le lac même. La section médiane de la carotte (584-376 cm; 12.3-9.0 cal kBP) a été caractérisée par une productivité primaire plus élevée dans le lac, une sédimentation beaucoup plus forte et une grande augmentation des apports de carbone organique (moyenne ~300 g CO m -2 a -1 ).

Les conditions dans la section supérieure de la carotte (< 376 cm ; < 9,0 cal kBP) suggèrent une productivité primaire élevée dans le lac et des taux élevés d'accumulation de CO (moyenne ~ 200 g CO m -2 a -1 ), avec des conditions environnementales stables. Les changements dans le caractère et la quantité de l'accumulation de CO peuvent avoir des implications importantes pour la libération future de gaz à effet de serre (GES). Ces résultats aident à comprendre l'histoire du développement des lacs dans les paysages de pergélisol dans le contexte des changements climatiques passés. Cette étude constitue le second chapitre de cette thèse.

Pour analyser l'hétérogénéité spatiale et temporelle des concentrations de GES et du flux diffus, les concentrations de CO2 et de CH4 dissous dans les lacs thermokarst d'un site d'étude en Yakoutie centrale ont été mesurées sur quatre saisons. Les lacs formés au cours de la première partie de l'Holocène (lacs 'alas') sont comparés aux lacs qui se sont développés au cours des dernières décennies. Les résultats montrent des différences dans les concentrations et flux de GES dissous (jusqu'à deux ordres de grandeur) entre les types de lacs et les saisons. Les lacs peu profonds situés dans des dépressions d'alas hydrologiquement fermées ont agi comme des puits de CO2 et d'importantes sources de CH4 dissous pendant certaines saisons. Les lacs thermokarstiques récents étaient des sources extrêmement élevées de CO2 et de CH4 dissous, avec une accumulation considérable de GES sous la couverture de glace (hiver) ou dans l'eau profonde de l'hypolimnion (été). Les flux diffusifs mesurés dans les lacs de thermokarst de ce paysage typique de taïga de la Yakoutie centrale sont parmi les plus élevés detectés dans les régions arctiques et subarctiques.

Cette étude représente le troisième chapitre de la thèse.

Enfin, le quatrième chapitre de cette thèse présente une analyse de télédétection qui a été menée afin d'étudier le développement récent (depuis les années 1960) des lacs et quantifier les différences dans l'histoire de leur développement entre les types de lacs. Les lacs situés dans des dépressions d'alas hydrologiquement fermées sont les plus nombreux et constituent la plus grande partie des lacs de la zone d'étude en termes de superficie. Les lacs thermokarstiques récents se trouvent souvent dans des zones où l'activité humaine est évidente, comme des terres agricoles récemment défrichées ou à proximité de routes. Certains lacs thermokarstiques récents se sont formés en l'absence apparente d'activité humaine, probablement en raison de l'augmentation des températures. Par ailleurs, les résultats indiquent que certains types de lacs sont plus sensibles aux changements annuels des précipitations. Les lacs situés dans des dépressions d'alas hydrologiquement fermées ont été plus affectés par les changements de précipitations que les lacs thermokarstiques récents ou les lacs d'alas connectés, considérés comme davantage 'résilients'.

Les références pour cette section peuvent être trouver à page 31.

General Introduction 1 Permafrost definition and distribution

Permafrost, which is ground or rock that has remained at or below 0 ºC for at least two consecutive years, is differentiated based on its spatial extent (continuous (90-100%), discontinuous (50-90%), sporadic (10-50%), and isolated (0-10%)), as well as its thickness, and the percentage of ground ice present. Most permafrost areas are located in Eurasia (65%) with the rest occurring in North America and Greenland (Figure 1). Approximately 20-25 percent of land area in the Northern Hemisphere is underlain by permafrost (Brown et al., 1997;Obu et al., 2019). Most permafrost areas lie within the continuous permafrost zone (54%), the remaining 46% is split almost equally between discontinuous, sporadic, and isolated permafrost zones (Figure 1) (Tarnocai et al. 2009). Continuous permafrost makes up the majority of the permafrost area in Russia, which itself contains a large portion of global permafrost (Figure 2). The continuous and ice-rich permafrost of late Pleistocene age that can be found in Eastern Siberia, parts of Alaska, and the western part of the Yukon territory in Canada is called 'Yedoma' (Schirrmeister et al., 2013;Strauss et al., 2017). (From Stolbovoi V. and I. McCallum. 2002) Depending on past and present surrounding environmental conditions, permafrost can extend below the ground surface from a few meters up to 1,500 meters. Most of the permafrost which currently exists today was formed during past cold glacial periods (since ~100,000 years ago) (Vonk et al., 2013b). Yedoma is generally older (formed ~50-20 kBP) than non-Yedoma permafrost. Areas where deep, continuous permafrost is found (i.e., eastern Siberia) generally remained unglaciated during these cold periods, as opposed to vast areas of Canada, northern

Europe and western Siberia. Without the insolating effects of a glacier or ice sheet, extreme cold was able to penetrate deep into the Earth's surface and permafrost formed at substantial depths (Strauss et al., 2017).

Permafrost can form either when existing soil or rock already in place is exposed to persistent freezing temperatures ('epigenetic') or simultaneously with the deposition of sediment that is being frozen ('syngenetically') (Grosse et al., 2013;French, 2017). Permafrost developing in alluvial or aeolian sediments that are being deposited is an example of syngenetic permafrost (French, 2017). During permafrost formation, some of the ground frozen during the winter does not thaw completely in the summer and the freezing layer propagates downward each year, becoming thicker each winter.

Common features of permafrost landscapes

Permafrost is typically overlain by an 'active layer', which freezes and thaws seasonally (Figure 3; Figure 4). The depth of the active layer can vary between <1 m and 20 m and is controlled by many factors including exposure to solar radiation, changing temperature and precipitation patterns, and human disturbances such as removing forest for agriculture (Kelley et al., 2004;French, 2017). In taiga landscapes, the active layer is generally shallower in permafrost areas covered by forest, which insolates and protects the ground from solar radiation. When forest cover is removed, by wildfire or human activity, there is usually a rapid and dramatic deepening of the active layer (Kelley et al., 2004). A common feature in lowland areas (0-300 m above sea level) of ice-rich and continuous permafrost are ice wedges. Ice wedges are massive, wedge shaped (with the widest part of the wedge nearest the surface) aggregations of ice (Figure 4; Figure 5a) (Grosse et al., 2013).

Environments favorable to ice wedge formation are poorly draining lowlands that are abundant in unconsolidated sediments. Ice wedges can be quite prolific in suitable environments. In central and northern Siberia, ice wedges are considered to be the dominant type of subterranean massive ice, constituting between 20 and 50 percent of total ice volume in the upper 5-10 m of permafrost (Vonk et al., 2013b;Fritz et al., 2016;French, 2017). In central Siberia, ice wedges can form up to depths of 50 m (Vonk and Gustafsson, 2013;Fritz et al., 2016;Ulrich et al., 2017b). Ice wedges form when freezing temperatures during the winter cause the soil (including the active layer) to contract, forming cracks (Figure 5b). Water seeps into these cracks as spring brings warmer temperatures, melting of seasonal snow, and rain. Once in contact with the still frozen soil below the active layer, the water freezes and expands. The cracks continue to form in the same place every winter and grow each spring and more water flows into the cracks and expands. As this cycle continues, the lower density of the ice wedge compared to the surrounding soil and the accumulation of mass to the ice wedge can cause a dome or ridge to form above the ice wedge.

Termed 'ice-wedge polygons', landscapes rich in these features can have a distinct hummocky or blocky pattern (Figure 4; Figure 5c). These hummocks are related to the aggregation of ice at the top of the ice wedge displacing the surrounding ground (Fritz et al., 2016;French, 2017). If water begins to pool in the troughs between ice wedges, this can act as a catalyst for 'thermokarst processes' (Figure 6). Exposure of the top of an ice wedge polygon can induce melting of the wedge and accelerate thermokarst initiation. Melting increases the depth of the troughs between hummocks, exposing more ice and inducing more melting. Snow will accumulate in the deepened troughs during winter, insulating the ground from cold air temperatures. In spring and summer, rain will also puddle in the deepened troughs, creating small ponds and further accelerating melting of the ice wedges. Eventually, through ground collapse and the coalescence of the small ponds, a more substantial lake may appear (Figure 6). Lakes formed from these thaw events are called 'thermokarst' lakes (sometimes called 'thaw' lakes) (Katamura et al., 2006;Kokelj and Jorgenson, 2013). This lake will continue to deepen and expand laterally into the surrounding permafrost through thermal erosion, soil compaction, and thaw slumping (Grosse et al., 2013;Séjourné et al., 2015). Even small lakes can have a profound effect on the surrounding thermal regime and hydrological patterns. The heat storage capacity of water is two times higher than ice and four times higher than dry ground. This high heat capacity, the low albedo of lake surfaces compared to the surrounding landscape, and their efficient absorption of long-wave solar radiation, can cause lakes to increase surrounding sediment temperatures up to 10º C above the mean annual air temperature (MAAT) and contribute to permafrost thaw even during cold winter months when thaw would normally be negligible (Grosse et al., 2013). The active layer under these lakes can become quite deep as heat absorbed during the summer is transferred into the surrounding permafrost, causing substantial thawing. Eventually, an area of perpetually unfrozen ground can form beneath the lake, called a 'talik' (sometimes called a 'thaw bulb') (Figure 4; Figure 6). During the summer and winter months, heat is transferred from the lake into the surrounding permafrost and talik. During the autumn and early winter, however, the temperature of the talik can be higher than the temperature of the lake. During these periods, heat is transferred from the talik into the surrounding permafrost as well as the lake waters at the water-sediment interface. The persistent transfer of heat from the lake and talik into the adjacent permafrost causes a positive feedback cycle of lake expansion and continued permafrost thaw (Grosse et al., 2013;French et al., 2017).

Areas with high ground ice content and substantial ice-wedge formations, such as Yedoma deposits, are particularly susceptible to warming temperatures and other disturbances as ice melt can cause rapid and localized surface subsidence, thermal collapse, and thermokarst lake development (Grosse et al., 2013;Vonk et al., 2013b;Natali et al., 2021). In areas with low ground ice content and discontinuous to isolated permafrost distribution, permafrost thaw can progress steadily and gradually, as a top-down process (Natali et al., 2021).

Anthropogenic climate change in the Arctic and sub-Arctic

Past and present climate change related to natural perturbations in the Earth's temperature as well as anthropogenic induced climate change can have profound effects on permafrost landscapes and the distribution of permafrost globally (Boike et al., 2015;Grosse et al., 2016;Biskaborn et al., 2019b). Not only are these landscapes particularly vulnerable, but the current rate and magnitude of temperature rise is greater in northern latitudes (Arctic and sub-Arctic) compared to global averages (Pörtner et al., 2019;Czerniawska and Chlachula, 2020). Long term records indicate a 0.07º C annual temperature increase and an overall shortening of the winter season since 1900 (Czerniawska and Chlachula, 2020). The MAAT in Central Yakutia (Sakha Republic) has increased from -9.6º C in 1980 to -6.7º C in 2019. April has experienced a shocking 5º C (2020) temperature increase since 1980. A record breaking 38º C was recorded in Verkhoyansk, Yakutia in June 2020 (Kharuk et al., 2021). Spring snow cover has been disappearing 3.4 days earlier per decade over the pan-Arctic terrestrial region and climate models predict decreases in snow cover duration between 10-20% by 2050 (Callaghan et al., 2011). Changes in average snow depth, however, are highly spatially heterogeneous with some areas of Eurasia experiencing increasing snow depth totals (Callaghan et al., 2011). Increasing temperatures and aridity have contributed to a proliferation of large and destructive wildfires in recent years (Kharuk et al., 2021). Climate warming, the shortening of winter snow cover, increased frequency and magnitude of forest fires, and anthropogenic removal of forests (for pastoral practices and logging) have reduced the albedo of the landscape and increased its exposure and vulnerability to solar radiation (Badmaev et al., 2019;Czerniawska and Chlachula, 2020). All of these factors can contribute to permafrost degradation. The susceptibility and response of permafrost landscapes to these disturbances depends largely on permafrost conditions such as ground ice content, permafrost distribution, and ground temperature (Grosse et al., 2011). As described above, permafrost degradation in areas of high ground ice content and continuous permafrost often results in thermokarst lake formation. Not only do these lakes profoundly change local thermal and hydrological regimes, but they can also act as biological hotspots for the emission of greenhouse gases (GHG), notably carbon dioxide (CO2)

and methane (CH4) due to the release of stored ancient carbon (Walter Anthony et al., 2016;Turetsky et al., 2020;Hughes-Allen et al., 2021).

Permafrost and local and global carbon cycles

Permafrost carbon stores consist primarily of the remnants of terrestrial vegetation such as leaves and root detritus which have accumulated in the perennially frozen soil over thousands of years (Davidson and Janssens, 2006;Vonk and Gustafsson, 2013;Schuur et al., 2015). Lakes store not only terrestrial material but also organic matter that was produced in the aquatic ecosystem by primary producers including algae. Abundant organic matter is created during the productive, but short, growing season typical of northern latitudes. However, prevalent anaerobic conditions, poor soil drainage, persistent freezing temperatures, and cryoturbation result in generally low carbon decomposition rates in permafrost environments (Davidson and Janssens, 2006). Cryoturbation, which is unique to permafrost landscapes, transfers surface soil organic carbon (SOC) into deeper soil layers until it reaches the top of the permafrost table, where cold temperatures shelter it from mineralization (Umakant et al., 2022). However, incorporation of this organic matter into permafrost sediments is highly heterogenous and depends on many factors including: soil moisture levels, vegetation cover, soil formation, cryoturbation and sedimentation rates, subsurface microbial decomposition, freezing and thawing conditions in the active layer, and aggradation of permafrost (Grosse et al., 2011). For example, high levels of organic matter accumulation and preservation (~33 kg C m -3 ) were observed in areas of syngenetic permafrost development with ice wedge growth during late Pleistocene interstadial periods. Analysis of the organic matter in these deposits reveal generally high total organic carbon (TOC) and other proxies indicative of low decomposition rates. Fluvially dominated and wind-driven accumulation during the same time period resulted in comparatively low levels of organic matter preservation (~7 kg C m -3 ) (Schirrmeister et al., 2011). The characteristics of this organic matter (i.e., total organic carbon, total inorganic carbon, carbon/nitrogen ratio) depend largely on whether the deposits originated during temperate interglacial or interstadial periods, or during harsh glacial or stadial periods (Schirrmeister et al., 2011).

Considering the high spatial heterogeneity of SOC in permafrost as well as the large global extent of permafrost affected areas (~21 million km 2 ), it is not surprising that it is difficult to constrain the total amount of carbon stored in permafrost landscapes (Figure 7) (Tarnocai et al. 2009). These estimates are further complicated by the fact that most in situ samples are taken from relatively shallow depths (0-200 cm) (Tarnocai et al., 2009;Hugelius et al., 2014). An initial estimate of permafrost SOC by [START_REF] Post | Soil carbon pools and world life zones[END_REF] based on 30 samples to a depth of 100 cm concluded that the surface layer (<100 cm) of permafrost-affected soils globally contained approximately 200 Pg of organic carbon. These estimates have since been improved by additional in situ sampling and improved spatial modeling of SOC in permafrost areas. Tarnocai et al. (2009) estimate that the northern circumpolar permafrost region contains a total of 1672 Pg of organic carbon (~1030 Pg from surface layers (0-300 cm) and ~650 Pg from deep Yedoma layers) (Figure 7). These totals account for about half of the global subterranean organic carbon pool and approximately twice the amount of carbon which currently exists in the atmosphere (Tarnocai et al., 2009;Hugelius et al., 2014). There is still a lot of uncertainty associated with these estimates. The SOC in Yedoma deposits is particularly poorly constrained. For example, deep cores from Northern Siberia and Alaska yielded organic carbon pool estimates of approximately 10 + 7/-6 kg m -3 (Strauss et al., 2013). A 22 m-deep core in Central Yakutia (Yukechi) yielded a much lower value of organic carbon content of ~ 5 kg m -3 (Windirsch et al., 2020), while another Central Yakutian study (Spasskaya Pad/Neleger site) of a shallow core (2 m) showed a notably higher organic carbon content of 19 kg m -3 for the top two meters of larch forest-covered Yedoma deposits (Siewert et al., 2015). The mobilization of this stored carbon could significantly impact global climate.

The organic carbon in permafrost soils remains sequestered by freezing temperatures until warming and/or disturbances to permafrost landscapes initiates permafrost thawing [START_REF] Hugelius | Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps[END_REF]Strauss et al., 2017). Warming can be related to natural climate oscillations, such as the warm temperatures experienced during the early Holocene, or mostly anthropogenic in origin, such as the warming temperatures currently being recorded globally (Biskaborn et al., 2019b;Pörtner et al., 2019;Kharuk et al., 2021). The mobilization of ancient carbon has the potential to contribute significant quantities of CO 2 and CH 4 to the atmosphere and amplify positive feedback to global air temperature increase ('permafrost carbon feedback') (Walter Anthony et al., 2014;Bouchard et al., 2015;Schuur et al., 2015;Hughes-Allen et al., 2020). Once thawed, the biolability of permafrost carbon (the amount of carbon that is accessible to an organism for uptake) depends primarily on its chemical character and its physical environment (Ewing et al., 2015). The chemical character of carbon is related to the origin of the carbon (terrestrial plant (C3 or C4) vs aquatic vegetation or algae) and the extent to which the carbon was decomposed before and during burial and incorporation into the permafrost table (Kuhry et al., 2010;Ewing et al., 2015). It is possible that the permafrost carbon stored in areas of shallow, discontinuous to sporadic permafrost was originally accumulated in the absence of permafrost. In this case, the carbon likely experienced substantial decomposition and decay before permafrost aggradation. This is observed in several subarctic palsa (low permafrost mounds) and peat plateau deposits which initially developed during the Holocene Thermal Optimum as permafrost free marshlands (Kuhry and Turunen, 2006).

Cryoturbation and loess deposition in Yedoma, however, provide mechanisms for the rapid incorporation of relatively undecomposed organic material into permafrost (Zimov et al., 2006;Kuhry et al., 2010;Schirrmeister et al., 2013). Studies have shown that the carbon from Yedoma in particular is highly biolabile once released (Schirrmeister et al., 2011;Vonk et al., 2013a;Windirsch et al., 2020).

Once thawed, permafrost carbon can be readily decomposed by microbes to CO2 and CH4, especially when released into aquatic systems, such as thermokarst lakes (Elder et al., 2018). In oxic environments, more CO2 will be produced, whereas CH4 production will be dominate in anoxic environments (Walter Anthony et al., 2010;Turetsky et al., 2020;Hughes-Allen et al., 2021). CH4 and CO2 are both powerful GHGs. CH4 in particular is 23 times more effective in absorbing longwave radiation than CO2 on a 100-year scale (~60 times more effective on a 20-year scale) (Schimel et al., 2002). CH4 and CO2 can be emitted from aquatic systems to the atmosphere CH 4 , which is predominately produced in anoxic sediments, can be transferred from the sediment either by ebullition or by diffusion (Bastviken et al., 2004;Walter Anthony et al., 2010;Sanches et al., 2019). Ebullition in particular is a direct pathway for CH4 from the sediment to the atmosphere without significant oxidation of the methane in the water column. When released by diffusive flux, however, CH 4 enters the water column where it may encounter oxic sediment or water. Once in contact with these oxic substrates, a large portion of the CH4 will be oxidized to CO2 by methane-oxidizing bacteria. Once the CH4 reaches the upper mixed layers of the water column, it will be released by diffusive flux. Diffusive flux is controlled by the difference in methane concentration between the water and the atmosphere, and on the rate of exchange between the water and the air, which is largely regulated by local wind speeds. Lake stratification can result in particularly vigorous CH4 production and storage (Bastviken et al., 2004;Hughes-Allen et al., 2021). In the Yedoma landscapes of central Yakutia, lake stratification occurs during winter for lakes whose depth exceeds winter ice cover (< ~ 1.5 m) and during the summer, when the temperature gradient between surface water and bottom water is greatest (Hughes-Allen et al., 2021). CH4 production and storage in stratified lake water are collectively called 'storage flux'.

Storage flux is released by diffusion during times during times of lake overturn such as spring and fall (Bastviken et al., 2004;Hughes-Allen et al., 2021). Like diffusive flux, the CH4 can experience significant oxidation before release. CH4 can also be released through plant facilitated emission along lake shores with nascent vegetation. This flux is regulated by CH4 production and oxidation within sediments and the characteristics of the littoral vegetation (Bastviken et al., 2004).

CO2 can be released from lake waters by oxidation of CH4 during diffusive and storage flux, respiration within an oxic water column and sediments (MacIntyre et al., 2018), lateral transfer of dissolved inorganic carbon from surrounding soils by groundwater flow (Zolkos et al., 2018), and photo-oxidation of dissolved organic carbon (Ward et al., 2017). Sediment respiration is the decomposition of organic matter through various process such as oxidation, metabolism, degradation, and mineralization. In oxic conditions, organic matter is oxidized by O2 and/or sunlight, and CO2 and recycled nutrients (i.e., nitrogen (N) and phosphorous (P)) are produced (Ward et al., 2017). Given the multitude of GHG emission pathways from lakes into the atmosphere, it is difficult to quantify the contribution of waterbodies to the atmospheric carbon budget.

This thesis focuses on Central Yakutia and the deep, ice and carbon rich permafrost of this area in order to ameliorate our understanding of Yedoma permafrost and its role in local and global carbon cycles in a continued warming scenario. The objectives of this thesis are threefold: 1)

understand the developmental history of central Yakutian lakes, particularly in terms of the organic carbon accumulation and preservation, 2) quantify differences in dissolved GHG concentrations and diffusive fluxes between three different lake types and seasons in central Yakutia, 3) conduct a large-scale investigation of lake occurrence in central Yakutia with focus on lake development through time and differences in developmental history between lake types. This information will then be used to upscale in situ greenhouse gas measurements to a larger landscape unit. This thesis combines paleolimnological analysis, in situ observations, and remote sensing analysis to provide a comprehensive understanding of lake development in permafrost landscapes and the contribution of these lakes to the global carbon budget.

14,000-year carbon accumulation dynamics in a Siberian lake reveal catchment and lake productivity changes

Article

Introduction

Permafrost is a dominant landscape feature in Siberia, Alaska, and Canada, occupying more than 20 million square kilometers and representing 24% of land cover within the northern hemisphere (Brown et al., 1997). Regional climate conditions, landscape cover, and other factors control the spatial extent and thickness of permafrost, and ground ice content can vary widely across landscapes (Grosse et al., 2016;Strauss et al., 2017). The Yedoma ice complex (referred to as 'Yedoma' in the following) is a particular type of ice-rich (50-90 vol%), relatively low organic content (2-4 wt%) permafrost that can reach depths of up to 40 m (Schirrmeister et al., 2013;Hugelius et al., 2014). Permafrost landscapes of all types are sensitive to changes in temperature and other local disturbances including forest fires and forest removal for agricultural purposes as well as lake formation and development (Grosse et al., 2013;Ulrich et al., 2019), which can have widespread implications for local and regional hydrology and the global carbon cycle (Walter Anthony et al., 2016). Yedoma is particularly susceptible to localized abrupt thaw based on its high ice content (Vonk et al., 2013a). The Arctic is currently warming at a disproportionately high rate and magnitude compared to global averages, with mean annual air temperature predicted to increase by as much as 5.4°C

within the 21 st century in the absence of significant and directed global effort to reduce greenhouse gas emissions [START_REF] Ipcc ; H.-O. Pörtner | IPCC Special Report on the Ocean and Cryosphere in a Changing Climate[END_REF]. This will likely herald a period of dynamic changes within permafrost landscapes.

The surface of Yedoma landscapes in many places is covered by ponds and lakes that document thermokarst processes (Strauss et al. 2017). In addition to the dominance of thermokarst lakes in Yakutia, lake formation can occur within the dune landscape which is widespread in the lower section of the Vilyui River and the middle part of the Lena River (Pestryakova et al., 2012). These sand dunes (also called tukulans) likely originated during the early stages of the interglacial epochs and dune lakes can often form in deflation basins (Pestryakova et al., 2012). Similar to thermokarst lakes, the high heat capacity of water relative to the air causes preferential thawing of surrounding permafrost as well as settling of sand to facilitate minor deepening of the lake (Sumgin et al., 1940;Kachurin, 1961;Zhirkov, 1983). While the genesis of dune lakes clearly differs from thermokarst lakes, there is overlap in the developmental history between these two lake types. There is, however, a relative paucity of studies which examine the paleolimnological history of dune lakes located in permafrost landscapes, and in particular differences in the potential for substantial release of stored carbon upon permafrost thawing between these lakes and thermokarst lakes.

Permafrost thaw can release substantial amounts of organic and mineral matter, including carbon, to surrounding terrestrial and aquatic ecosystems (Vonk et al., 2015). Soils across northern permafrost regions could contain twice as much carbon as currently exists in the atmosphere (Schuur et al., 2015). Total global terrestrial (non-marine) stores alone are estimated to hold 1,672 Pg of carbon (PgC), with Yedoma deposits accounting for more than 500 PgC of this total (Hugelius et al., 2014). Permafrost carbon stores consist primarily of the remnants of terrestrial vegetation such as leaves and root detritus as well as microorganisms which have accumulated in the perennially frozen soil over thousands of years (Davidson and Janssens, 2006;Vonk and Gustafsson, 2013;Schuur et al., 2015). Lakes, however, store both terrestrial material (allochthonous) and also the organic matter (OM) that was produced by algae in the aquatic ecosystem (autochthonous) (Schuur et al., 2015). However, there is a relative paucity of studies which study long-term carbon storage in Siberian lakes (Mendonça et al., 2017). Lakes can also act as hotspots for greenhouse gas emissions in permafrost landscapes. These emissions originate from the mineralization of OM (both allochthonous and autochthonous) stored within lake sediments (Bouchard et al., 2015;Hughes-Allen et al., 2021;Prėskienis et al., 2021) and understanding past lake dynamics can inform predictions about future greenhouse gas emissions from permafrost landscapes.

Paleolimnological studies often rely on proxy analyses based on carbon concentrations, carbon isotopes, and nitrogen concentrations. However, the applicability of these variables for interpretation of the relationships between carbon accumulation, carbon degradation, and the permafrost catchment are not fully understood (Biskaborn et al., 2019a). Therefore, the general approach of our study was to gain insights into how these variables can contribute to a better understanding of carbon dynamics in Siberian lake systems, facilitating future paleolimnological studies. Here, we present the results from a multi-proxy analysis of an approximately 6.6 m-long sediment core, covering the last 14000 years, from Lake Malaya Chabyda within the Central Yakutia region to understand the history and processes of carbon accumulation through time within a permafrost landscape. The specific objectives of the study were: 1. present the development history of Lake Malaya Chabyda, 2. quantify and understand the accumulation of organic carbon in the lake, 3. distinguish between organic matter produced within the lake itself and within the surrounding catchment, 4. discern carbon preservation trends based on climate variability and changing lake dynamics.

Study Site

Lake Malaya Chabyda (Озеро Малая Чабыда) (61.9569 °N, 129.4091 °E) is located approximately 15 km southwest of the City of Yakutsk (Central Yakutia, Eastern Siberia).

This lake is at 188 m a.s.l, has an area of 0.24 km 2 , and a max depth of three meters (Kumke et al., 2007). During initial surveying in July 2005, Lake Malaya Chabyda had a pH of 6.71, a conductivity of 131 (µS/cm), and a temperature of 18 ºC (Pestryakova et al., 2012). The lake catchment is 10 km 2 and also includes Lake Ulakhan Chabyda (Tarasov et al., 1996), which is a lake four times larger than Lake Malaya Chabyda, approximately three km to the northwest (Figure 1). Lake Ulakhan Chabyda has an area of 2.1 km 2 , an average depth of 0.5 m and a maximum depth of 2.0 m. There are no surface inflows into Lake Ulakhan Chabyda (Pestryakova et al., 2012), but this lake does discharge into Lake Malaya Chabyda during times of high water (i.e. after spring melting). Lake Malaya Chabyda sits on massifs of spearshaped dunes which have been fixed in place by vegetation growth since the onset of the Holocene. Both lakes sit on the former Lena River erosion-accumulation plain, within the central Yakutian Depression. This plain is composed of Quaternary loams overlying Cambrian limestones (Pestryakova et al., 2012). . Fedorov pers. comm.). Yearly evaporation rates exceed total precipitation in this region (Fedorov et al., 2014b). Between 1996 and 2016, the mean annual air temperature of Central Yakutia has increased by 0.5-0.6°C per decade (Gorokhov and Fedorov, 2018).

Permafrost in this region is continuous, thick (> 500 m deep), and the upper 30-50 m (Pleistocene-age fluvial and aeolian sediments) can be extremely rich in ground ice (50-90% by volume) (Ivanov, 1984). The active layer typically reaches depths of 0.5-2.0 m, varying depending on landscape factors that include vegetation cover type, general topography, soil type, and subsurface water content (Ulrich et al., 2017b). Central Yakutia is dominated by a middle taiga landscape regime (Fedorov et al., 2014) and larch, pine, and birch forests are prevalent (Ulrich et al., 2017a). Grasslands are abundant in unforested areas, including land cleared for farming, ranching, or in the remnant depressions of old thaw lakes known as 'alases'. After the cold temperatures and low precipitation of the Younger Dryas, Central Yakutia experienced a slow, but relatively persistent increase in temperature and precipitation (Müller et al., 2009;Biskaborn et al. 2012;Nazarova et al., 2013). These conditions resulted in widespread permafrost degradation (Biskaborn et al., 2012), including the development of 'alas' depressions that are now widespread in Central Yakutia (Soloviev, 1959;Brouchkov et al., 2004). The grasslands in the alases consist of halophytic steppe-like and bog plant communities (Ulrich et al., 2017b).

Nearly half of the landscape has been affected by thermokarst since the Early Holocene, resulting in the formation of thousands of partly drained 'alas' depressions (Soloviev, 1959;Brouchkov et al., 2004). However, recent thermokarst activity related to natural landscape evolution, increasing air temperatures and/or human-induced landscape modifications (agriculture, clear-cutting, and infrastructure) is also occurring in the region, resulting in the development of numerous small, fast-developing lakes and retrogressive thaw slumps along lake shores (Fedorov et al., 2014;Séjourné et al., 2015). Dune lakes are found in the lower section of the Vilyui River and the middle section of the Lena River, where there are abundant unvegetated massifs of undulating and cross-bedded quartz loamy sand that can be up to 25 m thick. These dune features are fairly common in the Lena watershed and can account for nearly 30 percent of the landscape in some areas (Galanin et al., 2018). Central

Yakutian dunes are considered to be aeolian in origin and likely formed at the end of the last glacial epoch between 27.0-12.0 cal kBP on Lower Cambrian carbonates, although the actual abundance of CaCO3 in the dune sand is not known. Lakes are found in deflation and aeoliandammed basins (Pestryakova et al., 2012).

Methods

Field Sampling

Eight overlapping sediment cores, representing an approximately 6.6 m-long composite sequence (Figure 2 Two parallel drilling holes, approximately 1 meter apart, were used alternatingly to obtain core sections which overlap by approximately 20 cm. To penetrate ca. 1 m of lake ice cover, 250-mm-diameter holes were drilled using a hand-held Jiffy ice auger. Water depth was measured using an Echo sounder (HONDEX PS-7 LCD) and a calibrated rope for verification. Each individual core sample consisted of a 100 cm-long core collected at 2 m water depth using a Russian peat corer and supported by an UWITEC gravity coring system.

Core samples were taken alternatingly from the first drill hole and then the second until the entire length of the core was obtained. Care was taken when coring to ensure that there was 20 cm of overlap between core sections. Cores were stored in waterproof sealed, transparent PVC plastic tubes in cool and dark conditions. After the field season, the cores were transported to Potsdam, Germany and stored at 4°C in the cold rooms at AWI. The cores did not experience any visible drying or surface oxidation during storage.

Sediment core subsampling and dating

After XRF scanning in early 2014, subsampling of the cores for laboratory analyses began in November 2018 with a simple visual description and photography of the eight cores.

One-cm-thick discrete subsamples were then taken at approximately 10 cm intervals using an inox spatula. Each subsample was split into two parts containing approximately equal amounts of material (4-10 g) and weighed. One subsample was kept in the cold room for potential future analysis. The remaining subsample was used for all subsequent analysis.

Four bulk samples and three organic vegetal macro remains were extracted from the sediment cores and sent for 14 C dating to the MICADAS radiocarbon lab at AWI, Germany (Table 1). These samples were placed in glass containers and dried at 50°C and analyzed for radiocarbon dating using accelerator mass spectrometry (AMS) after using an acid treatment (method outlined in (Vyse et al., 2020). We applied Bacon in R (Blaauw and Christen, 2011) and the IntCal20 calibration curve (Reimer et al., 2020) to model the age-depth relationship.

The surface of the core represents 2013 CE (the time of core retrieval) and the linear relationship of the samples from the bottom to the top of the core shows that there is no significant reservoir effect in the lake.

X-ray fluorescence (XRF) analysis

High-resolution X-ray fluorescence (XRF) analyses were carried out with 10 mm resolution on the entire sequence using an Avaatech XRF core scanner at AWI (Bremerhaven, Germany) with a Rh X-ray tube at 10 kV (without filter, 12 s, 1.5 mA) and 30 kV (Pd-thick filer, 15 s, 1.2 mA). The sediment surface was cleaned, leveled, and covered with a 4µm ultralene foil to avoid sediment desiccation prior to XRF scanning. Individual element counts per second (CPS) were transformed using a centered log transformation (CLR) and element ratios were transformed using an additive log ratio (ALR) to account for compositional data effects and reduce effects from variations in sample density, water content, and grain size (Weltje and Tjallingii, 2008). Statistical analysis was completed using the Python programming language (Python Software Foundation, https://www.python.org/).

XRF analysis of the sequence indicated 24 detectable elements and a subset of these were selected for analysis based on low element Chi-Square ( 2 ) values.  2 values are produced by the WinAxil Software to help determine the goodness of fit of the mathematical model.

Provided that the  2 value does not exceed 3, it is considered acceptable. These selected elements include the major rock forming elements of Silicon (Si) (Chi 2 1.4), Calcium (Ca)

(Chi 2 6.3), Titanium (Ti) (Chi 2 1.3), Rubidium (Rb) (Chi 2 0.6), Strontium (Sr) (Chi 2 0.7),
Zircon (Zr) (Chi 2 0.6) and the redox sensitive, productivity indicating elements of Manganese (Mn) (Chi 2 1.3), Iron (Fe) (Chi 2 2.5), and Bromine (Br) (Chi 2 0.8).

Grain size analysis

All subsequent analyses took place after the extracted subsamples had been freezedried until completely dry (approximately 48 hours). Grain size analysis was conducted on 16 samples that were chosen to span the entire sequence at relatively regular intervals. The samples were first treated for five weeks with H2O2 (0.88 M) in order to isolate clastic material. After treatment, seven samples were eliminated from the analysis because the remaining inorganic sediment fraction was too low for detection by the laser grain size analyzer. The remaining samples were homogenized using an elution shaker for 24 h and then analyzed using a Malvern Mastersizer 3000 laser. Standard statistical parameters (mean, median, mode, sorting, skewness, and kurtosis) were determined using GRADISTAT 9.1 (Blott and Pye, 2001).

Dry bulk density, sedimentation and organic carbon accumulation calculations

Total organic carbon concentrations (see below) were used to determine the organic vs.

mineral matter content (OM vs. MM) in each sample, assuming that bulk OM contains about 50% of organic carbon (Pribyl, 2010). OM and MM concentrations were used to derive average particle densities, based on values of 1.25 g cm -3 and 2.65 g cm -3 for OM and MM, respectively (Avnimelech et al., 2001). Dry bulk density (DBD, in g cm -3 ) values were then inferred by multiplying particle densities by porosity values, which had been calculated using wet and dry weights (thus the water content before and after sediment drying by freeze drying). Sedimentation rate (SR, in cm a -1 ) was calculated using the R function "accrate.depth", which estimated mean sedimentation rate derived from the age-depth model at 0.5 cm increments downcore. All iterations at each depth from the bacon modelling output where then used in a student t-test to calculate the 95% confidence range and the p-values for SR at each 0.5 cm increment. Sediment mass accumulation rate (MAR, in g cm -2 a -1 ) was obtained by multiplying DBD by SR. Finally, the organic carbon accumulation rate (OCAR, in g OC m -2 a -1 ) was inferred as the adjusted product of MAR and the total organic carbon concentration. OCAR and MAR uncertainties were calculated from the 95% uncertainty ranges of SR.

Biogeochemical analysis

Total carbon (TC), total organic carbon (TOC), and total nitrogen (TN) analyses were completed after the freeze-dried subsamples were ground in a Pulverisette 5 (Fritsch) planetary mill at 3000 rpm for 7 minutes. TC and TN were measured in a carbon-nitrogensulphur analyzer (Vario EL III, Elementar). Five mg of sample material were encapsulated in tin (Sn) capsules together with 10 mg of tungsten-(VI)-oxide. The tungsten-(VI)-oxide ensures complete oxidation of the sample during the measurement process. Duplicate capsules were prepared and measured for each subsample. Blanks and calibration standards were placed every 15 samples to ensure analytical accuracy (< ± 0.1 wt%). Between each sample spatula was cleaned with KIMTECK fuzz-free tissues and isopropyl.

Analysis of TOC began by removing the inorganic carbon fraction by placing each subsample in a warm hydrochloric acid solution (1.3 molar) for at least three hours and then transferring the sample to a drying oven. The TC measured for each subsample in the previous analysis was used to determine the amount of sample required for the TOC analysis.

The appropriate amount of sample was weighted in a ceramic crucible and analyzed using the Vario Max C, Elementar. The TOC/TN ratio was converted to the TOC/TNatomic ratio by multiplying the TOC/TN ratio by 1.167 (atomic weight of carbon and nitrogen) (Meyers and Arbor, 2001). Total inorganic carbon (TIC) analysis was completed using a Vario SoilTOC cube elemental analyzer after combustion at 400 ºC (TOC) and 900 ºC (TIC) (Elementar Corp., Germany).

Calculation of δ 13 C was completed twice for a subset of samples using two different methodologies. The analysis completed at the AWI Potsdam ISOLAB Facility removed carbonate by treating the samples with hydrogen chloride (12 M HCl) for three hours at 97 °C, then adding purified water and decanting and washed three times. Once the chloride content was below 500 parts per million (ppm), the samples were filtered over a glass microfiber (Whatman Grade GF/B, nominal particle retention of 1.0 µm). The residual sample was dried overnight in a drying cabinet at 50 °C. The dry samples were manually ground for homogenization and weighted into tin capsules and analyzed using a

ThermoFisher Scientific Delta-V-Advantage gas mass spectrometer equipped with a FLASH elemental analyzer EA 2000 and a CONFLO IV gas mixing system. In this system, the sample is combusted at 1020 °C in O2 atmosphere so that the OC is quantitatively transferred to CO2, after which the isotope ratio is determined relative to a laboratory standard of known isotopic composition. Capsules for control and calibration were run in between. The isotope composition is given in permil (‰) relative to Vienna Pee Dee Belemnite (VPDB).

The analysis of a small subset of samples which took place at Laboratoire des sciences du climat et de l'environnement Isotopic Laboratory for methodological comparison underwent a slightly different treatment, as follows. The sediment underwent a soft leaching process to remove carbonate using pre-combusted glass beakers, HCl 0.6N at room temperature, ultra-pure water and drying at 50 °C. The samples were then crushed in a precombusted glass mortar for homogenization prior to carbon content and δ13 C analysis. The handling and chemical procedures are common precautions employed with low-carboncontent sediments. Analysis was performed online using a continuous flow EA-IRMS coupling, that is, a Fisons Instrument NA 1500 Element Analyzer coupled to a

ThermoFinigan Delta+XP Isotope-Ratio Mass Spectrometer. Two in-house standards (oxalic acid, δ13C =-19.3% and GCL, _13C =-26.7 %) were inserted every five samples. Each inhouse standard was regularly checked against international standards. The measurements were at least triplicated for representativeness. The external reproducibility of the analysis was better than 0.1 %, typically 0.06 %. Extreme values were checked twice.

Those samples for which the carbonate was leeched at the room temperature, with lower HCl concentration (0.6N), and without a filtration step (samples analyzed at Laboratoire des sciences du climat et de l'environnement Isotopic Laboratory) had δ 13 C values 0.1‰ to 1.0‰ (average 0.5‰) higher than the samples treated at the higher temperature (97.7 ºC). However, the plotted δ 13 C curve is nearly identical for the subset of samples which were subjected to both treatments (Figure 3). There is some heterogeneity in the amount of offset between the two treatment methods. This might be related to an asymmetrical distribution of hot acid-soluble organic compounds throughout the sediment core. A correction of ca. +0.5‰ was applied to the results of the high temperature treatment.

These values were then combined with the low temperature results to provide a complete dataset for the whole core. The standard deviation of the typical measurement error (1σ) is generally better than δ 13 C = ±0.15‰.

Statistical Analysis

A hierarchical agglomerative cluster analysis was performed to divide the sediment core sequence into distinct units. XRF elemental data (CLR values for Ca, Fe, Ti, Mn, Sr, Zr, Rb, Br, and Si), δ 13 C, TN (wt%), C (wt%), TOC (wt%), and TIC (wt%) were used to perform the cluster analysis.

A principal component analysis (PCA) was completed on a subset of XRF elements (Ca, Sr, Fe, Zr, and Ti), and biogeochemical parameters (TOC (wt%), C (wt%), TN (wt%), TOC/TNatomic) to explore the dimensionality of the dataset and the relationships between the included variables. All values were z-transformed. The data when considered for the whole core, do not have a Gaussian distribution, but the data do approach a normal distribution when considered by unit. Therefore, the PCA was performed individually for each unit, rather than on the whole core. Only a subset of XRF elements were included in the PCA in order to reduce dimensionality, while preserving as much variability as possible. All variables were then plotted onto an ordination plot for interpretation purposes. The clustering analysis and PCA analysis were performed using the Python programming language (Python Software Foundation, https://www.python.org/). Packages used include: numpy (Harris et al. 2020) and pandas (The pandas development team, 2020).

Results

Chronology and sedimentation rates

The results of the 14 C dating analysis of seven samples indicate that the oldest lake sediments at the base of the core were deposited around 14.1 cal kBP according to the age depth model (Figure 4). The middle horizon (dark brown clayey silt with laminated sections)

decreases in age consistently from 12.3 cal kBP at 584 cm depth to ~9.0 cal kBP at 376 cm depth. The upper horizon (homogenous, brown, clayey silt) decreased from ~9.0 cal kBP at 376 cm depth to present-day at the top of the sediment core. Mean 95% confidence over the entire core was 874 yr. 100% of the dates overlap with the age-depth model (95% ranges). It is possible that the bulk-sediment radiocarbon dates are older by multiple centuries or even millennia than recorded due to reworking of old organic material (Strunk et al., 2020).

However, based on the linear alignment of dated samples until the time of sampling, we can assume neglectable effects from the presence of old carbon in the samples.

Figure 4. (a) General stratigraphy and (b) age-depth model of the sediment core, based on radiocarbon dating of 7 samples, shown as green/yellow stars (see details in Table 1). (c) Inferred sedimentation and mass accumulation rate (SR and MAR) curves along the core, based on the age-depth model.

The age depth model indicates a higher average sedimentation rate (SR) of 0.13 cm a -1 below 290 cm, and a lower average SR of 0.07 cm a -1 in the upper horizon. Considering the dry bulk density (DBD) of each sample, instantaneous sediment mass accumulation rates (MARs) range from 0.02 g cm -2 a -1 above 300 cm depth to 0.06 g cm -2 a -1 near the bottom of the core. Considering each unit separately, average MAR values are 0.05 g cm -2 a -1 for Unit 1, 0.03 g cm -2 a -1 for Unit 2, and 0.02 g cm -2 a -1 for Unit 3. Average OCAR is similar in Unit 1 (57 g m -2 a -1 ) and Unit 3 (55 g m -2 a -1 ) and highest in Unit 2 (77 g m -2 a -1 ) (Figure 5).

Figure 5. Organic carbon accumulation rate (OCAR, in g OC m -2 a -1 ) for each biogeochemical subsample along the core (see Figure 8), inferred from the age-depth model, as well as sedimentation and mass accumulation rates (see Figure 4).

General Stratigraphy

The composite sequence is separated into three broad stratigraphic units based on the sedimentological and biogeochemical analyses (Figure 4) and the cluster analysis (Figure 6).

Due to an existing talik (area of unfrozen ground surrounded by permafrost) below Lake Malaya Chabyda (Bakulina et al., 2000), the entire sequence was unfrozen. From bottom to top, stratigraphic zones are described as follows:

• Unit 1 (663-584 cm) (14.1 cal kBP-12.3 cal kBP). The bottom 80 cm of the sequence consists of dark brown massive (i.e., not laminated) clayey silt. The lower section (663-619 cm) is notably dry and has a 'crumbly' texture. A small gypsum aggregate particle was found within this unit (623 cm) and identified using a binocular microscope. The upper section (636-584 cm) appears less dry and displays mm-scale, sections of lighter grayish brown.

• Unit 2 (584-376 cm) (12.3 cal kBP-9.0 cal kBP). This unit is composed of dark brown, laminated, clayey silt with interstratified light brown to white laminations.

These laminations are well-defined horizontal layers of non-calcareous sediment

Grain-size distribution

The grain size of the sediment core ranges from uni-, bi-, to trimodally distributed, and spans a size range from medium silt to very-fine sandy coarse silt. The mean grain size for the entire core is 29.05 m (± 6 m), with a minimum size of 14.4 m at 156 cm depth and a maximum size of 40.9 m at 624 cm depth (Figure 7). The sorting ranges from poorly to very-poorly sorted. Sorting becomes generally more poorly sorted below 400 cm depth (Figure 7). The grain size distribution of clay, silt, and sand remains fairly consistent throughout the length of the core, with silt dominating the distribution (~ 80%) throughout the entire core (Figure 7).

Biogeochemistry

From the base of the sequence to the top of Unit 1 (584 cm depth), TOC content decreases consistently from 17 to 6.5 wt% (Figure 8). There is a sharp increase in TOC content from the lowest value of 6.5 % at 584 cm depth (top of Unit 1) to ~35 wt% at just above 500 cm depth (base of Unit 2). The TOC content is relatively stable in the upper two units (Unit 2 and Unit 3) of the sediment sequence, remaining between 30-40 wt%. OCAR, controlled by both TOC and MAR, experiences two low peaks (~20 g OC m -2 a -1 ) at 600 cm depth and 500 cm depth. OCAR decreases throughout Unit 1 from approximately 100 g OC m -2 a -1 to 24 g OC m -2 a -1 . Apart from the low peak at 500 cm depth, OCAR increases relatively steadily throughout Unit 2. There is a high peak of 127 g OC m -2 a -1 at 438 cm depth. OCAR then decreases into Unit 3 where it stabilizes around 50 g OC m -2 a -1 from approximately 300 cm depth. Similar to TOC, TN is lowest in Unit 1 at the base of the core (below 1 wt%) and increases steadily throughout Unit 2 to approximately 3 wt% (Figure). TN is relatively stable throughout Unit 3 at approximately 3.5 wt%. There is less variability in TN above 150 cm depth and a steadily increasing trend to a maximum value of nearly 4 wt% at the top of the core (Figure 8). The TOC/TNatomic is highest at the bottom of the sequence (TOC/TNatomic = 20) and decreases consistently throughout Unit 1. After a transient peak at 570 cm depth (TOC/TNatomic = 18), the TOC/TNatomic ratio is very stable in the top two horizons of the sediment sequence, remaining at nearly 12 (Figure 8). The TOC/TNatomic ratio is highest for Unit 1 (20), while Unit 2 and Unit 3 have similar values of 14 and 12, respectively (Figure 9).

TIC content is low throughout the entire sequence (0.0-2.4 wt %), although Unit 1 (below 584 cm) has slightly higher average values compared to the other units (Unit 1 average: 1.9 wt % versus Unit 2 average: 0.6 wt%, Unit 3 average: 0.7 wt %) (Figure 8). δ 13 C values decrease consistently throughout Unit 1 (Figure 8). In the overlying units, δ 13 C values increase from the lowest value (-31.9‰) at 475 cm depth (middle of Unit 2) until approximately 150 cm depth (middle of Unit 1) to -20.3‰. Above 150 cm depth, δ 13 C values decrease slightly and then remain stable until the top of the core (Figure 8). Although identifying the distinct sources of organic material in a lacustrine environment is complicated by competing signals from the lake and catchment, the relationship between TOC/TNatomic and δ 13 C values can help further distinguish OM origin (Meyers and Teranes, 2006). Unit 1 plots within the C3 land plant zone, while Units 2 and 3 plot within the lacustrine algae zone (Figure 10). Again, much of Unit 3 plots above typical values for lacustrine algae.

Inorganic elemental composition

The main rock forming elements analyzed include rubidium (Rb) (Kalugin et al., 2013), zirconium (Zr) (Marshall et al., 2011), titanium (Ti) (Balascio et al., 2011), silicon (Si) (Marshall et al., 2011;Martín-Puertas et al., 2011), calcium (Ca) and strontium (Sr) (Bouchard et al., 2011) Redox sensitive elements analyzed include manganese (Mn) and iron (Fe) (Haberzettl et al., 2007;Bouchard et al., 2011). Mn values are consistent throughout the entire core, with some small changes in the variability of the values between units. Fe values in Unit 1 are high, with low variability. There is a low peak in the lower half of Unit 2, with values returning to Unit 1 levels by the top of Unit 2. Fe values fluctuate between high and low peaks from the bottom of Unit 3 until approximately 150 cm depth (4.0 cal kBP), where values decrease consistently. Fe/Ti (Figure 11), which also represents reducing conditions as well as a possible reduction in grain-size (Marshall et al., 2011;Davies et al., 2015), Elements and elemental ratios that are sensitive to increased organic content include bromine (Br) (Kalugin et al., 2007(Kalugin et al., , 2013;;Bouchard et al., 2011) and silicon/titanium (Si/Ti) (Melles et al., 2012). Br values decrease throughout Unit 1 (Figure 11). Values in Unit 2 increase from the bottom of this unit until a large peak at 520 cm depth, after which, the values decrease. Unit 3 has decreasing Br values until 250 cm depth, after which there is a sharp, but brief increase. Br values decrease from approximately 230 cm depth until the top of the sediment core. Si/Ti, which is particularly sensitive to increases in biogenic silica (notably diatoms) (Melles et al., 2012), follows a different pattern. Unit 2 has the highest Si/Ti values (Figure 11). Unit 3 has Si/Ti values which are slightly lower than Unit 1, with a slight increase that peaks around 100 cm depth (2.6 cal kBP).

decreases
Figure 11. XRF profiles of selected elements and elemental ratios along the core. The thin black lines represent continuous measurements (sampling interval of 10 mm), and the thicker red lines represent the 5-point running mean.

PCA Analysis

For Unit 1, PC1 accounted for 56.8% of explained variance, while PC 2 accounted for 22.7% variance (Figure 12a). Ca, Sr, and TOC/TNatomic are positioned together in the upper left quadrant, while Fe, C (wt%), TOC (wt%), and TN (wt%) are positioned together in the lower left quadrant. Zr is positioned in the upper right quadrant, while Ti is positioned in the lower right quadrant. δ 13 C is very well negatively correlated to PC2.

For Unit 2, PC1 accounted for 39.6% of explained variance, while PC 2 accounted for 24.4% of explained variance (Figure 12b). Ca and Sr are positioned together in the upper left quadrant, while TOC/TNatomic, Zr, and Ti are positioned together in the bottom left quadrant.

C (wt%), TOC (wt%), and TN (wt%) are very well positively correlated to PC1. Fe and δ 13 C are well negatively correlated to PC 2.

For Unit 3, PC 1 accounted for 32.8% of explained variance, while PC2 accounted for 22.4% of explained variance (Figure 12c). Sr is very well negatively correlated to PC 1, while

Discussion

Multiproxy-inferred paleolimnological history

Unit 1 (Late Pleistocene)

The observed depositional history of Lake Malaya Chabyda starts in Unit 1, which spans a time frame from approximately 14.1 -12.3 cal kBP (Late Pleistocene). The TOC/TNatomic ratio at the bottom of the core represents the maximum (20) measured in the entire recovered sequence, indicating a stronger contribution of carbon produced by vascular land plants than aquatic algae (Figure 9) (Meyers, 1994). Both, the TOC (average = 12 wt%) and TN (average = 0.68 wt%) values are lower in Unit 1 compared to Unit 2 and Unit 3, although relatively high compared to other sites in this region (Vyse et al., 2021). Decreasing TOC and MAR values, resulted in decreasing OCAR values throughout this Unit culminating in a low peak at 600 cm depth (Figure 4, Figure 5).

During the time of deposition, this region was still experiencing cold temperatures associated with the Late Pleistocene deglaciation period in the Northern Hemisphere and active microbial decomposition would have been restricted to a short period of time after spring thawing and before the onset of winter cold temperatures. Therefore, microbial decomposition would not have been particularly abundant, limiting significant OM degradation (Davidson and Janssens, 2006). As a result, the lower levels of TOC observed in Unit 1 cannot be explained solely by more active microbial decomposition compared to Unit 2 and Unit 3. It is more likely that Unit 1 experienced lower OC input compared to the upper two units. It is also possible that longer periods of ice cover compared to Unit 2 and Unit 3 restricted in situ autotrophic production (algae) in the lake. However, these TOC values are high compared to other reported values from Yedoma permafrost (Windirsch et al., 2020;Vyse et al., 2021; Figure 8). Moderately high δ 13 C values, high TIC values, and relatively low TOC (Figure 8) compared to Units 2 and 3 corroborate lower levels of bioproductivity within the lake and indicate input from sources of inorganic carbon during this time (Schirrmeister et al., 2011), as shown in other lake records from Yakutia (Biskaborn et al., 2012). These trends could also reflect an increase in accumulation of authigenic CaCO3, which was not explored in this study, but remains a possibility.

Unit 1 is massive, exhibiting no layering, and has fairly homogenous elemental composition for most of the elements examined using XRF. Analysis of detrital elements Rb, Zr, Si, and Ti indicate that Unit 1 had the highest level of terrestrial input compared to Unit 2 and Unit 3. Higher mineral detrital input coupled with low TOC, high TOC/TNatomic ratio, and moderate δ 13 C values indicate overall low productivity within the lake. The absence of soil stabilizing vegetation surrounding the lake due to prevalent cold conditions could have facilitated the input of clastic material from the catchment (Subetto et al., 2002;Nazarova et al., 2013). Si/Ti ratios and preliminary inspections of smear-slides (one slide for each unit)

suggest that there are diatoms present in all units, indicating that there was an aquatic system at the core location (Vyse et al., 2020). This matches with a continuous diatom record over 14 cal kBP from Lake Ulakhan Chabyda (Figure 1) 4 km northwest of Lake Malaya Chabyda (Pestryakova et al., 2012;Herzschuh et al., 2013). High Ti values can also indicate increased run-off from precipitation events and/or increased aeolian deposition (Davies et al., 2015).

This could have been the case at the study site, with active soil erosion and transport associated with poorly developed vegetation cover during this period of cold, Late

Pleistocene climate (Biskaborn et al., 2021b). Precipitation was generally low in the study site region during the time of deposition, making increased aeolian deposition more likely than increased run-off from precipitation events (Biskaborn et al., 2012). Relatively high Ca, Ca/Ti, and Sr values can be indicative of a persistent alkaline environment under semi-arid conditions associated to Lake Malaya Chabyda during the deposition of Unit 1. These conditions are also observed in neighboring lakes (Ulkahan Chabyda and Temje Lake) (Herzschuh et al., 2013;Nazarova et al., 2013). Both low precipitation and/or strong summer insolation would enhance evaporation. Additionally, during this period, there were open landscapes, strong winds, and, as a result, increased evaporation from the surface of the lakes.

The absence of soil cover and swamps around the lake contributed to the absence of humic acids in the surface runoff and the water in the lake was predominantly alkaline, with low bioproductivity and ultra-oligotrophic in contrast to the later Holocene period. It is also possible that changes in Ca, Ca/Ti, and Sr trends throughout the core are more related to changes in sources for these elements.

The Mn/Fe ratio for the Lake Malaya Chabyda core is generally higher and more variable in the lower sections of the core, specifically in Unit 1. Higher Mn/Fe compared to the rest of the core indicates predominately oxic conditions as Mn is more readily reduced (dissolved) under anoxic conditions compared to Fe (Davies et al., 2015;Vyse et al., 2020).

Shallow lake depth (i.e. light can penetrate to the bottom of the lake) as well as low productivity and stronger lake water mixing could be responsible for maintaining oxygenated bottom waters and a shallow lake could experience evaporation indicated by the proxies described above.

Grain-size distribution and sorting provide information about the dominant denudation, erosional and transport processes at the study site. Grain size ranges from coarse silt to fine silt and poor to very-poorly sorted sediments are present throughout the entire core (Figure 7), indicating either a relatively short sediment transport distance or the influence of a combination of different erosional and transport processes (Folk and Ward, 1957). It is also possible that the early lake had a very different shore structure than is present today and that active slumping and abrasive processes could have taken place, creating a large amount of unsorted mineral matter input to the lake. A longer transport distance or the influence of a single dominant transport process, such as aeolian activity, would likely have resulted in stronger sorting of grain sizes throughout the core (Biskaborn et al., 2013). Unit 1 sediments are slightly more 'poorly sorted' compared to Unit 2 and Unit 3. This is perhaps related to higher detrital input from the surrounding catchment. The polymodal nature of the grain size distributions suggests that the source sediments are heterogeneous and polygenetic in origin (Figure 13) (Schirrmeister et al., 2011;Wang et al., 2015;Ulrich et al., 2019). According to grain-size results, there is no sand-dominated material that could have clearly confirmed the sand-dune origin ('tukulan') of the lake, as sand percentage remained < 20% along the whole core, including at its base (Figure 7). In fact, the top of the core has the highest percentage of sand (18%). However, we expect that such sand-rich deposits lie deeper, in older sediments. The proxy analyses suggest that Unit 2, spanning the Late Pleistocene-Early Holocene transition (~ 12.3 -9.0 cal kBP), was deposited under variable limnological conditions.

Unit 2 has a consistently high level of TOC (30 -36 wt%; average 31 wt%), and its TOC/TNatomic ratio (average of 14, ranging from <10 to 18) is close to the 4 -10 range indicating production by phytoplankton or algae (Meyers, 1994;Ulrich et al., 2019) (Figure 9). There is a consistent reduction in TIC throughout Unit 2, while  13 C first decreases slightly, then increases in the upper half of Unit 2 (Figure 8).  13 C values can be varying inversely to water depth and/or in a direct relationship with primary productivity (Meyers, 2003), suggesting that Lake Malaya Chabyda could have experienced lake deepening and then shallowing and/or a decrease followed by an increase in primary productivity within the lake as supported by strongly increasing Si/Ti values (Figure 11). During the Holocene increasing thermokarst activity likely caused lake deepening, followed by lake shallowing as thermokarst activity stabilized at the lake depression began to be filled by sediments (Andreev et al., 2003). Analysis of these proxies nonetheless suggests a transition to a dominance of lacustrine versus terrestrial carbon source between 12.3 cal kBP and 9.0 cal kBP (Biskaborn et al., 2013).

Sedimentation experienced an initial decrease in the bottom half of Unit 2 (MAR decreasing from 0.03 g cm -2 a -1 bottom of Unit 2 to 0.006 g cm -2 a -1 at 500 cm depth), before gradually increasing to 0.04 g cm -2 a -1 cm at 450 cm depth (Figure 4). MAR then decreases slightly until about 300 cm depth, where it stabilizes until the top of the core. OC delivery was high during the Early Holocene (average OCAR 83 g OC m -2 a -1 above 500 cm of composite depth, i.e., after ~ 11 cal kBP until about 300 cm depth, i.e., 7.9 cal kBP). Such rates are clearly well above the reported values for high-latitude lake basins and notably higher than global modern values (Vyse et al., 2021 and references therein).

Unit 2 is the only unit in the Lake Malaya Chabyda sediment core which exhibits any clear layering or laminations (Figure 4). The subsection of Unit 2 between 535 cm depth and 500 cm depth exhibits the most pronounced layers, with mm scale alternations between light colored non calcareous layers (determined using dilute acid) and the brown, massive clayey silt layers. The presence of these layers indicates a lack of bioturbation and enhanced preservation of OM (Melles et al., 2012). Analysis of detrital elements, Rb and Zr indicate that Unit 2 experienced lower levels of detrital input compared to Unit 1 and higher rate of deposition of organic material produced within the lake, i.e., by algae. Unit 2 has the lowest Ti values compared to Unit 1 and Unit 3, further suggesting a decrease in detrital input. This signal could also be related to a decrease in the relative proportion of detrital inputs due to an increase in accumulation of lighter OM produced within the lake coupled with a decrease the detrital input (Balascio et al., 2011;Davies et al., 2015). This signal is also corroborated by a slight decrease in SR between Unit 1 and Unit 2 and a decrease in MAR between Unit 1 and Unit 2. A shift in bulk sediment deposition from mainly allochthonous in Unit 1 to mainly autochthonous sources in Unit 2 is likely indicative of an increase in biological activity within the lake contributing to an increase in the deposition of OM compared to sediment input from the surrounding catchment. In the Holocene, vegetation cover and soils are formed within the catchment area, resulting in a weakening of the processes of denudation and erosion (Andreev et al., 2003). Only dissolved substances enter the surface runoff, and groundwater and sedimentation begin to be dominated by intra-reservoir processes, particularly, the bioproductivity of the lake ecosystem (Subetto et al., 2017).

Values of Ca and Sr in Unit 2 are more variable but have slightly lower average values compared to Unit 1. There is a low peak in both elements at approximately 530 cm depth (approximately 11.3 cal kBP) suggesting that the lake did not experience any significant evaporative stages or that there were changes in the availability and input source of Ca-rich clastic material. These variabilities in Ca and Sr could also be associated with changing wind directions and reallocation of sediment sources controlled by migrating dune features, as well as changing patterns of detrital input from the surrounding catchment. These variabilities may also be due to a change in the hydrochemical characteristics of the water in the lake. Above this low peak, Ca and Sr values return to Unit 1 levels, but with higher variability, suggesting continued changes within the catchment that affected the source and/or availability of Ca and Sr source materials. Unit 2 has the highest values for Ca/Ti, which can indicate carbonate deposition (Davies et al., 2015). Ca/Ti can be controlled by in-lake (autochthonous) carbonate precipitation, biologically mediated calcite production, and/or input of additional old carbonate from catchment sources. Haberzettl et al. (2007) suggest that high Ca/Ti values indicate lower lake levels, resulting in a proliferation of biologically induced calcite production (i.e., P. lenticularis) related to warmer water temperatures and higher concentrations of nutrients. Warmer water temperatures and higher concentrations of nutrients could be related to shallow lake conditions for Unit 2. These hypotheses are corroborated by  13 C, which increase consistently throughout Unit 2, suggesting a decrease in water depth and/or an increase in biological activity (Meyers, 2003).

The Mn/Fe ratio for Unit 2 decreases consistently throughout the entire unit. Oxygendepleted conditions could be related to continued microbial activity beneath the winter-ice cover and/or higher levels of primary productivity associated with mixing of the bottom waters in the summer season (Hughes-Allen et al., 2021). This is a possible scenario as the Early Holocene experienced very cold winters (Meyer et al., 2015) and thus prolonged periods ice cover (Biskaborn et al., 2012). The Early Holocene is also associated with an increase in thermokarst lake development in central Yakutia related to warm summers. Warm summer temperatures, like prolonged periods of ice cover, can cause periods of anoxic bottom waters (Hughes-Allen et al., 2021). This is also supported by PCA biplots (Figure 12), suggesting a strong relationship between redox-sensitive elements (Fe) and organic matter proxies in the lower units 1 and 2 (Heinecke et al., 2017).

Unit 2 has the highest Si/Ti values compared to the other two units, which is also indicative of high levels of biological activity. The depositional characteristics of Unit 2 suggest that a transition occurred at approximately 12.3 cal kBP from the preferential deposition of terrestrial carbon and low lake primary productivity conditions of Unit 1 to preferential deposition of aquatic carbon and high lake primary productivity in Unit 2.

Unit 3 (Early Holocene to present)

The proxy analyses indicate that Unit 3, deposited since the Early Holocene (~ 9.0 cal kBP), was formed with some variability in lake depth, bottom water oxygen availability, and biological activity. Unit 3 is more homogenous than Unit 1 and Unit 2, both in terms of a lack of discernable laminations or other features and homogeneity within the biogeochemical metrics analyzed. Average TOC content is > 34 wt% (ranging from 29 to 39 wt%) in this unit, higher than elsewhere along the core and an order of magnitude above most of the reported values for lakes across Yedoma regions (Vyse et al., 2021). Presumably, this period heralded higher levels of bioproductivity and increased nutrient supply to an expanding lake at the study site. The mean TOC/TNatomic ratio of Unit 3 is 12, which indicates some dominance of production by phytoplankton or algae (Meyers, 1994;Ulrich et al., 2019) (Figure 9) with minor contribution from vascular plants in lake catchment area (Heinecke et al. 2017). Shallower lake waters would have amplified the influence of carbon from the lake catchment on the TOC/TNatomic ratio. This value suggests a dominance of lacustrine rather than terrestrial carbon input between 9.0 cal kBP and present day (Biskaborn et al., 2013).

OC delivery decreased after 7.9 cal kBP, as shown by an average OCAR of 55 g OC m -2 a -1 for Unit 3. This was the result of a lower sedimentation (average MAR decreasing from ~ 0.02 g cm -2 a -1 at the top of Unit 2 to ~ 0.01 g cm -2 a -1 in the upper section of Unit 3).

Still, OCAR values of 55 g OC m -2 a -1 are notably in the upper range of global modern values and significantly higher than elsewhere across high-latitude regions (Vyse et al., 2021; Figure 5). TIC remains consistently low throughout Unit 3, while  13 C increases in the lower half of the unit, before decreasing and then stabilizing in the upper half of the unit. This could indicate decreasing lake levels followed by increasing lake levels. Changes in  13 C could also be controlled by changes in primary productivity. In this case, primary productivity would have decreased briefly, then increased from the bottom to the top of Unit 3 (Meyers, 2003).

The profiles for detrital elements Rb and Zr are nearly identical to Unit 2 and offer no strong indications of changing catchment regimes or erosional transport. However, Ti values are much higher compared to Unit 2, returning to Unit 1 levels. It is possible that Ti deposition was more sensitive to increases in detrital input. Ti can be also be associated with increased run-off from rain events (Corella et al., 2012) and/or increased aeolian deposition (Bakke et al., 2009) and it is possible that the observed trends are linked to changes in these processes.

Unit 3 has the lowest values for elements Ca and Sr, suggesting a distinctly lacustrine environment that did not experience any significant evaporative events. Ca/Ti, which is also associated with evaporative events, is also consistently low and less variable throughout Unit 3. Mn/Fe values decrease from the bottom of Unit 3 until approximately 200 cm depth (5.4 cal kBP), suggesting increasing lake depth and more frequent periods of anoxic bottom water conditions (Davies et al., 2015). As lake depth increases, unfrozen lake water remains during winter and microbial activity can continue under the ice cover. The ice cover does not allow gas exchange with the atmosphere, creating an anoxic environment at the water-sediment interface (Hughes-Allen et al., 2021). Anoxic conditions can also occur in the summer related to high levels of primary productivity. Unit 3 has the lowest Si/Ti values (between 390 cm depth and 175 cm depth) compared to the rest of the core, and relatively low levels of Br, which can be related to decreases in the input of biogenic silica (Melles et al., 2012) and organic content (Kalugin et al., 2007(Kalugin et al., , 2013;;Bouchard et al., 2011), respectively. These low values for elements normally associated with high biological activity suggest that primary productivity is not high enough to account solely for an increase in anoxic bottom water conditions. Mn/Fe values increase from 200 cm depth until approximately 75 cm depth (1.9 cal kBP), where they remain stable until the top of the sediment core. Increasing Mn/Fe values suggest less frequent anoxic conditions, which might have resulted from moderate lake depth decrease, enhanced wind activity, and/or shorter periods of ice-cover allowing for more water column mixing (Hughes-Allen et al., 2021). Si/Ti values also increased during the same time period, potentially indicating a slight increase in the deposition of biogenic silica, which could be related to an increase in primary productivity. Perhaps a decrease in lake depth, that was not so drastic as to cause Lake Malaya Chabyda to freeze to the lake bottom in the winter, caused an increase in water temperature and increased nutrient concentrations during the biologically productive ice-free seasons. These analyses suggest that Unit 3 was deposited under very stable sediment deposition conditions, although there was some variability in lake depth, bottom water oxygen availability, and biological activity. Lacustrine algae were the main source of OC deposited in the lake during this period.

Lake Malaya Chabyda carbon accumulation rates

Total organic carbon concentration (TOC) is a crucial proxy for understanding the abundance of OM in sediments, including the proportion of OM that evaded remineralization during the sedimentation process. The concentration of OM in sediment is generally equivalent to twice the recorded TOC value (Meyers, 2003). Therefore, TOC values can suggest initial production of biomass as well as subsequent levels of degradation. Moreover, since TOC concentrations are expressed in % weight, therefore influenced by mineral/clastic matter inputs (artificially diluted or concentrated), it can be useful to infer organic carbon (OC) accumulation rates (expressed as mg OC cm -2 a -1 or g OC m -2 a -1 ) when reliable agedepth model and estimations of sediment dry bulk density for each sample are available (Meyers, 2003).

High TOC values, high TOC/TNatomic ratios, and low  13 C generally reflect OM which have not undergone significant decomposition under anaerobic conditions (Schirrmeister et al., 2011;Ulrich et al., 2019). Although Unit 1 showed the highest TOC/TNatomic ratio (20), it also had a lower weight percent of carbon (<20 percent C) compared to Unit 2 and Unit 3 (greater than 30 percent C for both units). Unit 1 can therefore be classified as a mineral sediment (<20 % C). Mineral sediments lose 6-13 % on average of their OC within one decade after exposure due to thawing or other processes (Schirrmeister et al., 2011).

Moreover, the inferred OCARs for Unit 1 indeed represent low values compared to Unit 2, with an average of ~100 g OC m -2 a -1 (ranging from 24 to 105 g OC m -2 a -1 ; Figure 5).

These low values are nevertheless higher than reported rates for temperate latitudes, such as the Great Lakes in North America (e.g., Meyers, 2003), and much higher than several arctic/subarctic sites, such as northern Québec (Ferland et al., 2012(Ferland et al., , 2014)), Finland (Pajunen, 2000), Greenland (Anderson et al., 2009;Sobek et al., 2014), as well as southeastern and northeastern Siberia (Martin et al., 1998;Vyse et al., 2021).

Unit 2 and Unit 3 have greater than 20 % carbon and are therefore classified generally as organic sediments. Unit 2, in particular, exhibited significant layering, which suggest a lack of bioturbation and enhanced preservation of OM and/or seasonal changes in sedimentation processes. Organic deposits, including deposits that are aquatic in origin (i.e.

fluvial, alluvial, and lacustrine), typically exhibit decade-long losses of 17% -34% of their OC after exposure by thawing or other processes (Schirrmeister et al., 2011). Some studies suggest that input of ancient carbon into aquatic systems may augment or even galvanize remineralization of modern dissolved OC (Vonk and Gustafsson, 2013;Mann et al., 2015;Strauss et al., 2017). This effect is likely due, in part, to low levels of carbon decomposition during deposition (i.e. colder conditions) and before thawing (Vonk et al., 2013a). A significant portion of the Lake Malaya Chabyda sediment core is classified as organic sediment, which is predicted to lose comparatively high percentages of their OC upon potential exposure. Jongejans et al. (2021) found that although the OC content of the Yukechi Yedoma ice complex sediments was relatively low, there was substantial greenhouse gas release upon thawing. These findings point to OM quality and decomposition history and more important drivers of greenhouse gas release than OM content alone (Jongejans et al., 2021). Although a lake currently exists, proxy evidence discussed above suggest that Lake Malaya Chabyda did experience high levels of evaporation, which might have brought the lake close to desiccation in the past. Changing temperature and precipitation regimes, lower precipitation and higher temperatures for example, might make drying out more likely in the future for this relatively small and shallow lake. In this case, it is important to consider OM quality and possible future greenhouse gas release. Furthermore, inferred OCARs for Unit 2

show a strong increase from the base (42 g OC m -2 a -1 ) to the top (76 g OC m -2 a -1 ) of this unit, in accordance with developing lacustrine conditions and enhanced biological productivity from algae (i.e., mostly autochthonous source of OM).

As lake Malaya Chabyda is located in a semiarid climate, with high summer temperatures and humid climate, the primary production in its ecosystems can be quite high, resulting in strong accumulation of organic matter in bottom sediments (Biskaborn et al. 2021b). Given the significant increase of OCAR (Figure 5) and TOC accumulation (Figure 8) around 11 cal kBP, Lake Malaya Chabyda likely transitioned to an OC sink at approximately the Pleistocene-Holocene transition. However, data about lake morphology (especially the ratio between lake area and mean depth) and C emissions, which control the net C balance (Ferland et al., 2012), are not available for this site. Ferland et al. (2012) found that sediments in large and shallow lakes, (i.e., Lake Ulakhan Chabyda), experience higher rates of decomposition compared to smaller lakes. This is due, in part, to the greater exposure of sediments in these lakes to warmer surface waters which accommodate higher rates of sediment decomposition. Additionally, small lakes, especially those with steeper bathymetry, can have high sedimentation rates and sediment focusing, which reduces the efficiency of sediment decomposition. These characteristics of small lakes lead to generally lower rates of OC degradation and enhanced long-term burial and storage of carbon. Compared to the nearby Lake Ulakhan Chabyda, notably larger but shallower (see study area above), Lake Malaya Chabyda likely acted as a more efficient C sink (per m 2 ) with its 'compact' morphology (Ferland et al., 2014). The relatively high sedimentation rate recorded within Lake Malaya Chabyda, especially between 13 cal kBP -7.9 cal kBP, likely also contributed to the lake acting as an efficient C sink, in part by reducing O2 exposure time (Ferland et al., 2014).

 13 C values can provide information about sources of OM and past productivity trends, as well as identifying changes in the availability of nutrients in surface waters (Meyers, 2003;Ulrich et al., 2019). As phytoplankton preferentially use the lighter carbon isotope ( 12 C), an increase in productivity causes an initial, but relatively brief decrease in  13 C in the lake water. However, as the phytoplankton deplete dissolved inorganic carbon stores in the lake, the  13 C values of the residual inorganic carbon in the lake water increases, causing enrichment in the  13 C values of the newly created OM (Meyers and Teranes, 2006). Reequilibration between lake water DIC and atmospheric CO2 will occur over time, resulting in the deposition of OM not enriched in  13 C in the sediment layer.  13 C values decrease slightly from the bottom of the core until a low point at approximately 450 cm depth (10.0 cal kBP) (middle of Unit 2). This  13 C trend is complicated by the high levels of TOC and a TOC/TNatomic value of 14 in Unit 2, which would normally accompany more enriched  13 C values (Meyers and Arbor, 2001). It is possible that these low  13 C values around 450 cm depth are related to changes in the  13 C signature of available DIC in the lake and/or to changes in the degradation pathways of the  13 C. It is also possible that low algae production in the lake resulted from long periods of ice cover and/or low temperatures unfavorable to primary production.

 13 C values increase consistently from 450 cm depth (10.0 cal kBP) to 150 cm depth (4.0 cal kBP), where they remain stable to the top of the core, suggesting a continuous increase in primary productivity. However, this presumption is complicated by the fact that water from Lake Ulakhan Chabyda can flow into Lake Malaya Chabyda during times of high water. As Lake Malaya Chabyda is relatively small, it is possible that this influx of water could affect the  13 C signal in the sediment core. An increase in  13 C values can also be related to an increase in nutrient delivery from soil erosion and/or active layer development.

Perhaps the continued expansion of Lake Malaya Chabyda by thermokarst processes has increased the availability of nitrates and phosphates, which could have enhanced primary productivity, resulting in an increase in  13 C values (Meyers and Teranes, 2006).

The TOC/TNatomic vs.  13 C relationship (Figure 10) suggests that the majority of carbon deposited between 14.1-12.3 cal kBP (Unit 1) was from terrestrial (C3) land plants (Meyers and Arbor, 2001;Ulrich et al., 2019). This suggests that the bioproductivity in the lake was very low at the time of deposition (Heinecke et al., 2017) and therefore not recorded in the sediment core. Also possible is that detrital input to the lake was very high during this time.

Decreasing MAR combined with low and decreasing TOC resulted in decreasing OCAR values within Unit 1 with a low peak at the top of this unit (Figure 4c; Figure 5). The dominance of vascular plant carbon also suggests higher input from the lake catchment to the lake compared to the upper two units. Thermokarst lakes usually do not have a classical littoral zone and the edges of the lakes are often overgrown with aquatic vegetation and surrounded by peat bogs. Therefore, there is no direct mineral runoff from the catchment area, as it is filtered by the peatlands. However, there can be significant input of organic matter coming from the surrounding peatlands to the lake. The OM found in Units 2 and 3 plots more toward an origin of lacustrine algae. Many of the points in Unit 3 plot above the expected values for lacustrine algae. The relatively high values (>-25‰) generally occur above 250 cm depth (6.7 cal kBP). Enriched  13 C values can occur is alkaline lacustrine systems where there is limited availability of dissolved atmospheric CO2 ( 13 C=-8.5‰), therefore HCO3-( 13 C= 1‰) becomes the dominant source of inorganic carbon for lake algae, shifting the isotopic signature to higher (more enriched) values (Meyers, 2003). This effect might have been limited in Lake Malaya Chabyda, where modern pH is slightly acidic (pH = 6.71). This change in carbon source is reflected by comparatively heavier isotopic compositions for lake algae than normally expected (Meyers, 2003). Periods of very high bioproductivity resulting in a depletion in dissolved CO2 can be responsible for a shift from CO2 to HCO3 -, among other possible explanations (Meyers, 2003).

Connections between the lake environment, permafrost dynamics, and climatic conditions

Lake development and thermokarst activity began in earnest in Eastern Siberia between 15.0-10.0 cal kBP (Nazarova et al., 2013;Ulrich et al., 2019), with the postglacial period experiencing gradual rising temperatures and permafrost thaw (Schirrmeister et al., 2011;Biskaborn et al., 2011). However, this period was also punctuated by periods of colder temperatures and low precipitation (Andreev et al., 2003;Nazarova et al., 2013). Andreev et al. (2003) recorded a transient climate amelioration during the Allerød in Eastern Siberia based on pollen spectra collected from a sediment core from Lake Ulakhan Chabyda (9.75 m depth). They suggest that mean July temperatures in this area were 1.5-2 ºC below present values and January temperatures were -2--5 ºC colder (mean annual temperatures -3--4 ºC colder than present). However, subsequent low precipitation levels (150 mm/year lower than modern conditions) and colder temperatures (Andreev et al., 2003) likely slowed the development of many nascent thermokarst lakes that were created during this relatively brief warm period before the onset of the Younger Dryas. Paleoclimatic estimates suggest that mean summer temperature decreased by 3 ºC and January mean temperature decreased by 6-7 ºC during the Younger Dryas (Andreev et al., 2003). It is possible that these environmental conditions lead to the signal, recorded in Unit 1 (14.1 cal kBP-12.3 cal kBP), of higher inputs of terrestrial plants from the surrounding catchment rather than deposition of OM produced in the lake. We also see lower values of Br, associated with organic matter production, and Si/Ti, usually related to biogenic silica production, in Unit 1 compared to the other two units.

Dry, cold conditions, associated with permafrost stability in the catchment (i.e., no or rare thermokarst activity) and persistent ice cover, were likely favorable to increased detrital input from the catchment to the lake relative to organic matter, as recorded by the XRF signals.

Environmental conditions during the time of deposition of Unit 2 (12.3 cal kBP-9.0 cal kBP) encompassed the Younger Dryas cool period. Conditions during the height of the Younger Dryas were likely colder and dryer than today and then experienced a slow, but persistent warming associated with the end of the Younger Dryas and the transition to the Early Holocene warm period (Andreev et al., 2003;Katamura et al., 2006;Biskaborn et al., 2012;Ulrich et al., 2019). Müller et al. (2009) records the peak of the episodic cooling and drying related to the Younger Dryas between 12.4-11.3 cal kBP based on pollen records from a sediment core taken from Lake Billyakh (approximately 400 km north-west of Lake Malaya Chabyda). However, the conditions were likely still warmer and wetter than the period before the Allerød climate amelioration (15-13.5 cal kBP) (Müller et al., 2009). It is difficult to see evidence of the Younger Dryas in the proxy data analyzed for the Malaya Chabyda sediment core. The transition to the Holocene climate optimum began around 11.3 cal kBP in Eastern Siberia as recorded in the pollen record, peaking between 6.7-5 cal kBP (Müller et al., 2009;Nazarova et al., 2013;Biskaborn et al., 2016). This transition, likely associated with enhanced thermokarst activity, is reflected in the increased deposition of OM produced within the lake, as well as a reduction in detrital input from the surrounding catchment (e.g., Bouchard et al., 2017, and references therein). Proxies (e.g., Mn/Fe) also indicate potentially stronger depletion in bottom-water oxygen caused by enhanced water column stratification during this time (Bouchard et al., 2011). Velichko et al. (1997) reconstructed January and July temperatures that were still 1 ºC and 0.5 ºC lower than present values combined with 25 mm/yr lower annual precipitation.

Lake development processes, including thermokarst, generally increase across the region at the time of the Early Holocene warm period, which corresponds to the lower boundary of Unit 3 (9.0 cal kBP-CE 2013) (Müller et al., 2009;Nazarova et al., 2013). Nazarova et al. (2013) hypothesize that temperatures between 8-4.8 cal kBP were likely warmer than present day and record and increase in sedimentation rate in Lake Temje consistent with increased biological production. July temperatures are recorded to be 1.5 ºC warmer than today, with the warmest temperatures occurring between 6.7-5 cal kBP. Our results show an increase in sedimentation rates at approximately 10.0 cal kBP, based on the age depth model, and higher TOC values after 9.5 cal kBP, resulting in notably higher OC accumulation rates during the first half of the Holocene. The upper half of Unit 2 seems to be distinctly lacustrine and the development of the lake basin was likely facilitated by rising precipitation levels and associated permafrost thaw (Nazarova et al., 2013;Schirrmeister et al., 2013;Ulrich et al., 2019).

Temperatures dropped below modern values in Central Yakutia after approximately 4.5 cal kBP with lowest temperatures estimated to have occurred between 3.0 and 2.0 cal kBP (Biskaborn et al., 2012;Nazarova et al., 2013). This cooler episode ('Neoglacial') has also been reported elsewhere across the circumpolar North, for example in the Canadian Arctic (Fortier et al., 2006;Bouchard et al., 2020). After this cool period, there is a general warming trend, interspersed with short-term temperature fluctuations to present-day values (Biskaborn et al., 2021a). It is possible that the enriched  13 C observed in the upper half of Unit 3 are in response to warmer temperatures and favorable growing conditions.

Conclusions

• Here we present an in-depth and high-resolution analysis of a nearly seven-meter-long sediment core which spans the Pleistocene-Holocene transition and encompasses a continuous Holocene time series. There was considerable variation in biogeochemical proxies both between and within three stratigraphic units.

• A shift occurred between 12.5 cal kBP to 11.0 cal kBP from predominately terrestrial land plant contribution to lacustrine algae contribution, which is recorded in OM deposition.

• Unit 2 and Unit 3 have high TOC values, high TOC/TNatomic ratios, and low  13 C, generally indicating OM which have not undergone significant decomposition under anaerobic conditions.

• There is high carbon content, including TOC (wt%) for this sediment core compared to other similar sites in Central Yakutia, and elsewhere across the Arctic. OCARs are above the highest reported values for temperate and high-latitude regions, for both past (Holocene and Late Pleistocene) and modern conditions. Lake Malaya Chabyda might have thus acted as an efficient OC sink since the Pleistocene-Holocene transition.

• Increases in lake depth and nutrient availability from the catchment increased bioproductivity within the lake and organic matter preservation and storage relative to decomposition. Compact lake morphology (relatively small surface-to-depth ratio), associated with higher sedimentation rates and less exposure to warmer and oxygenrich shallow waters, likely contributed to notably high OC preservation.

Seasonal patterns in greenhouse gas emissions from thermokarst lakes in Central Yakutia (Eastern Siberia)

Article

Introduction

Permafrost occupies more than 20 million square kilometers and represents 24% of land cover within the northern hemisphere [START_REF] Brown | Circum-Arctic map of permafrost and ground-ice conditions[END_REF] predicted to increase by as much as 5.4°C within the next century in the absence of significant and directed global effort to reduce greenhouse gas emissions (IPCC 2019). This will likely herald a period of dynamic changes within permafrost landscapes.

A prevalent pathway of permafrost degradation in areas of ice-rich permafrost is the initiation of thermokarst processes, which eventually result in the formation of numerous lakes in regions where topography is flat (Grosse et al. 2013). Thermokarst processes begin when disturbances such as warming or forest removal cause deepening of the active layer (the layer which thaws each summer), which results in either thaw subsidence or thermal erosion depending on the relief of the terrain and the ground ice content. In areas that are dominated by low relief terrain, such as Central Yakutia, thermokarst subsidence often causes ponding and shallow depressions (French 2017). The coalescence and expansion (both laterally and vertically) of these ponds eventually results in the development of larger and deeper lakes, with a portion of unfrozen water remaining under the ice cover in winter (Bouchard et al. 2020). Once formed, these lakes profoundly change the local ground thermal regime, sometimes increasing surrounding sediment temperatures by as much as 10°C above the mean annual air temperature (Brouchkov et al. 2004). Thermokarst processes, especially the presence of a lake or pond, increase the rate of permafrost thaw significantly compared to what would be predicted from increases in air temperature alone (Brouchkov et al. 2004;Schuur et al. 2015). Continued lateral lake expansion can cause slope slumping (so-called 'retrogressive thaw slumps') or further ground surface subsidence and erosion (Séjourné et al. 2015;Bouchard et al. 2017). Expansion and deepening continue until the lake depth surpasses the depth of ice-rich permafrost or thaw propagation beneath the lake is impeded by lake sediment insulation properties, halting further deepening. After this phase, lake size and depth are controlled mostly by the ratio between precipitation and evaporation [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]French 2017;Grosse et al. 2013). Although some of the mechanisms involved are not well understood, drainage, evaporation and infilling eventually result in lake disappearance [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]Desyatkin et al. 2009).

It is estimated that soils in northern permafrost landscapes contain twice as much carbon as currently exists in the atmosphere (Hugelius et al. 2014). Specifically, thermokarst lakes are considered as biogeochemical hotspots that play a pivotal role in the processing of permafrost organic matter and therefore potentially on the global climate (Walter Anthony et al. 2016). Anoxic sediments at the bottom of these lakes are sites of CH4 production that can be released through ebullition and diffusion. The CH4 released though ebullition in shallow lakes has little opportunity to be oxidized within the water column, so it typically represents a direct flux to the atmosphere (Bastviken et al. 2004;Bouchard et al. 2015). However, the CH4 released through diffusive flux, or through ebullition in deeper stratified lakes, remains in the water column long enough so that a significant fraction can be oxidized to CO2 by methaneoxidizing bacteria [START_REF] Vachon | Influence of water column stratification and mixing patterns on the fate of methane produced in deep sediments of a small eutrophic lake[END_REF]. In winter, the CH4 emitted through both processes can be oxidized under the ice cover [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]). Significant release of CO2 from thermokarst lakes in northern latitudes has also been documented (e.g., [START_REF] Abnizova | Small ponds with major impact: The relevance of ponds and lakes in permafrost landscapes to carbon dioxide emissions[END_REF].

There are three main pathways by which CO2 can be released from thermokarst lakes, besides when it originates from CH4 oxidation: 1) respiration within the oxic water column Little is known about the seasonality of greenhouse gas emissions from thermokarst lakes developed within permafrost. This is particularly true for winter and spring seasons that present logistical constrains [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]). Here we report on dissolved greenhouse gas concentrations, and derived diffusive fluxes, in thermokarst lakes of different development stages in Central Yakutia, Eastern Siberia, a region experiencing rapid degradation of ice-rich permafrost. The limnological properties and greenhouse gas concentrations were quantified four times throughout a full annual cycle (2018-2019). We tested the two-fold hypothesis that 1) marked local landscape heterogeneities (permafrost conditions, lake morphology and water balance) strongly affect greenhouse gas concentrations and fluxes in thermokarst lakes of Central Yakutia, and 2) these spatial patterns are further influenced by strong seasonal changes in lake water temperature, mixing regime, and dissolved oxygen.

Grasslands are abundant in unforested areas, such as land previously cleared for farming, ranching or in the remnant depressions of old thaw lakes known as 'alases' (see below). They consist of halophytic steppe-like and bog plant communities (Ulrich et al. 2017b).

Permafrost in this region is continuous (Figure 1), thick (> 500 m deep), and the upper 30-50 m (Pleistocene-age fluvial and aeolian sediments called 'Yedoma') can be extremely rich in ground ice (50-90% by volume) (Ivanov 1984). The amount of organic carbon stored in Yedoma varies widely. For example, deep cores from w Siberia and Alaska yielded organic carbon pool estimates of approximately 10 +7/-6 kg/m 3 (Strauss et al. 2013). A 22 m deep core in Central Yakutia (Yukechi) on the Abalakh terrace yielded a much lower value of organic carbon content of ~5 kg/m 3 (Windirsch et al. 2020), while another Central Yakutian study (Spasskaya Pad/Neleger site) of a shallow core (2 m) showed a notably higher organic carbon content of 19 kg/m 3 for the top two meters of larch forest-covered Yedoma deposits (Siewert et al. 2015).

Our study site is underlain by Yedoma silty loams, which are common to the Lena-Aldan interfluve, with abundant ground ice in the form of 1.5-3 m-wide ice wedges. Active layer depth in the region generally ranges between ~ 1 m below forested areas to > 2 m in exposed grassland areas (Desyatkin et al. 2009). Zones of unfrozen ground (or taliks) exist underneath major rivers and lakes deeper than the ice-cover thickness. Nearly half of the landscape has been affected by thermokarst since the early Holocene, resulting in the formation of thousands of partly drained alas depressions (Soloviev 1959;Brouchkov et al. 2004). However, recent thermokarst activity related to natural landscape evolution, increasing air temperatures and/or human-induced landscape modifications (agriculture, clear-cutting, and infrastructure) is also occurring in the region, as shown by the presence of numerous small and young, fast-developing lakes and retrogressive thaw slumps along lake shores (Fedorov et al. 2014b;Séjourné et al. 2015). The landscape is thus highly dynamic; yet the competing driving factors (climate vs. local geomorphology or vegetation development) and the timing of such changes (gradual vs. rapid/threshold) are complex to characterize and quantify at the regional scale.

Lake types

To understand the spatial and temporal heterogeneities in greenhouse gas concentrations from thaw lakes in the area, it was first necessary to classify the lakes at the study site ( Figure 2).

Differentiation was done based on field observations, past radiocarbon dating of lake sediments, geochemical signatures of lake waters, and a multiple-stage development model (Desyatkin et al. 2009;[START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]. Relevant statistics and properties for all lake types are given in Table 2.

Based on the compiled information, the lake types are:

• Unconnected alas lakes: residual water bodies located within endoreic or hydrologically closed basins (Desyatkin 2009). Most of these lakes likely formed during the transition between the Pleistocene and Holocene, approximately 8-10 thousand years ago (Ulrich et al. 2017b). However, the Holocene Thermal Maximum (~6700-5000 cal. yr BP) was also a time of prolific formation and continued development of new and existing alas lakes (Biskaborn et al. 2012;Ulrich et al. 2017b).

These lakes can be up to a few meters deep but are typically very shallow (1 meter deep or less) and are thus generally frozen to the bottom in winter. Despite a lack of liquid water in the winter, these frozen lakes have heat fluxes that can be one order of magnitude higher than those of the surrounding permafrost (Boike et al. 2015). The ancient lake depressions surrounding the small residual lakes of this type can be up to several kilometers wide and several meters deep. These alas lakes have already undergone much of the thermokarst processes and very little ground ice typically remains beneath the residual lake. Therefore, the thaw potential and resulting input of stored carbon to these lakes is low compared to recently formed thermokarst lakes (Ulrich et al. 2019). Yakutian people have used these depressions for agricultural purposes for centuries [START_REF] Crate | Permafrost livelihoods: A transdisciplinary review and analysis of thermokarst-based systems of indigenous land use[END_REF]).

• Connected alas lakes: lakes connected hydrologically to the watershed by streams or rivers. These lakes are consistently larger (several hundreds of meters across) and deeper (up to ~ 10 meters). Most of them were probably formed during the late Holocene, approximately 5-3.5 thousand years ago, although detailed chronology about their inception is still incomplete [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]Ulrich et al. 2017b). Local people currently use some of these lakes for fishing (e.g., Syrdakh Lake).

• Recent thermokarst lakes: thaw lakes formed over the last several decades mostly from human activities (e.g., forest fire and forest removal for agriculture, pipelines, or road construction) and rising temperature (Fedorov et al. 2014b). These lakes are generally small (meters to tens of meters across) and relatively shallow (generally one to two 3 Methods

Physicochemical characteristics of lake water

In 2018-2019, data were collected during four field campaigns: September 2018 (fall),

March-April 2019 (winter), May 2019 (spring), and August 2019 (summer) (Figure 1; Table 2; Figure 3). Physicochemical measurements were conducted during each of the four field campaigns using a YSI Pro DSS multiprobe sensor, including profiles in temperature (accuracy of + 0.2°C), dissolved oxygen (accuracy of ±1% of reading or 1% saturation, whichever is greater), specific conductivity (accuracy of ±1% of reading) and pH (accuracy of ±0.2 pH units). Vertical profiles of the water column were obtained from 17 lakes during the winter through a hole drilled in the ice cover using a hand-held auger, and from 3 lakes during the summer from an inflatable boat (Table 2). All other data were collected from the lake surface as discrete point measurements, mostly from the lakeshores.

Dissolved greenhouse gas measurements

In spring, summer, and autumn, the water was mostly sampled from the lakeshore at the surface. Summer samples were also collected from the pelagic section of the lake at the surface and bottom using a horizontal van Dorn sampler (integrating 15 cm depth) for recent thermokarst Lake 1 only. In winter, samples were collected from the center of each lake through the ice cover using a vertical van Dorn sampler (integrating ~ 40 cm depth), including a few supplementary samples taken from bottom waters. Water was quickly transferred into a 2 L gas exchange LDPE bottle (if not directly sampled from this bottle), and a headspace of 30 mL of atmospheric air was created (bottle with an inlet at the bottom, connected to a syringe). The bottle was shaken vigorously for three minutes, and 20 mL of the gaseous headspace was then transferred to a 12 mL pre-vacuumed gas-tight Exetainer glass vial previously flushed with nitrogen [START_REF] Hesslein | Carbon dioxide pressure in surface waters of Canadian Lakes[END_REF]. Samples were later analyzed by gas chromatography (Thermo Trace 1310, TRI-Plus Head-Space auto-sampler, HSQ 80/100 1.6 mm × 30.5 m in series with MS 5A 1.6 mm × 18.3 m, thermal conductivity and flame ionization detectors). The dissolved greenhouse gas concentration (Csur) was calculated using Henry's Law, and departure from saturation was obtained subtracting the gas concentration in the water at equilibrium with the atmosphere (Ceq) using global annual average values of 410 ppm for CO2 and 2 ppm for CH4 (IPCC 2019).

Diffusive fluxes of dissolved greenhouse gas were estimated as in Bouchard et al. (2015) based on a gas transfer coefficient (k600) standardized to a Schmidt number (Sc) of 600 [START_REF] Wanninkhof | Relationship between gas exchange and wind speed over the ocean[END_REF], and calculated with the wind-based model of [START_REF] Vachon | The ecosystem size and shape dependence of gas transfer velocity versus wind speed relationships in lakes[END_REF]: k600 = 2.51 + 1.48 U10 + 0.39 U10 log10 LA where U10 is the wind speed at 10 m above the ground from the closest meteorological station (at Yakutsk) and LA is the lake surface area (in km 2 ). The flux (in mmol m -2 d -1 ) was then calculated by applying the equation:

𝐹𝑙𝑢𝑥 = 𝑘 (𝐶 𝑠𝑢𝑟 -𝐶 𝑒𝑞 )
where k is the transfer coefficient for a given gas calculated as:

𝑘 = 𝑘 600 ( 𝑆𝑐 600 ) -0.5
Due to generally high pH (>9) and low wind turbulence during the ice-free seasons, chemical enhancement could further increase the high levels of CO2 uptake estimated for some lakes [START_REF] Wanninkhof | Chemical enhancement of CO2 exchange in natural waters[END_REF]. The exchange rate of gases from the atmosphere through the 123 boundary layer of a liquid surface is regulated by molecular diffusion and accelerated by increasing turbulence (wind). However, the exchange rate of CO2 molecules from the atmosphere through the boundary layer of a lake surface is also regulated by water temperature, pH and ionic strength. When favorable conditions are met (high pH and temperature, low winds), diffusion across the boundary layer can be enhanced by the hydration of CO2 to HCO3 -. The CO2 uptake is therefore a function of both diffusion and chemical reaction. The chemical enhancement factor for CO2 () is defined as:

 = kenh/ k
where kenh is the enhanced gas transfer velocity, and k is the gas transfer velocity for a nonreactive gas with the same diffusion coefficient as CO2 as calculated above [START_REF] Wanninkhof | Chemical enhancement of CO2 exchange in natural waters[END_REF].

Statistical analysis

The dissolved greenhouse gas data were not found to be normally distributed and did not meet the assumption of homoscedasticity based on the Shapiro-Wilk normality test.

Typical data transformations to improve normality, such as a Box-Cox and logarithmic transformations, were attempted with unsuccessful results. The small sample sizes likely impede the ability to achieve normality, and we did not expect the data to be inherently normal. Therefore, we used Kruskal-Wallis one-way analysis of variance to test the relationships between dissolved CO2 and CH4 concentrations vs. lake types and seasons.

Significance between groups was tested for dissolved CO2 and CH4 concentrations as follows: the values for dissolved greenhouse gas concentrations for each lake type were compared across seasons (e.g., dissolved CO2 concentrations for unconnected alas lakes between fall, winter, spring, and summer) and between lake types for a specific season (e.g., dissolved CO2 concentrations for unconnected alas lakes, connected alas lakes, and recent thermokarst lakes during the fall season). The null hypothesis that the medians of the compared groups are equal was rejected if p < 0.05. All of the statistical analysis was completed using the Python programming language (Python Software Foundation, https://www.python.org/).

Results

Seasonal conditions

Sampling began in September 2018, at the end of the ice-free season when the depth of the active layer is at its maximum. During this period, nighttime air temperatures drop near 0°C, initiating lake water mixing (fall overturn). Although limnological profiles were not The hydrochemical characteristics of the sampled lakes differed notably between lake types and seasons (surface values presented in Figure 4). Regardless of the season, unconnected alas lakes generally had the highest electrical conductivity. This inter-lake difference was particularly marked during the winter, when the mean conductivity value for unconnected alas lakes was 5711 μS cm -1 . The conductivity for unconnected alas lake #3 reached 8000 μS cm -1 , four times higher than the average conductivity for connected alas lakes. Connected alas lakes showed the lowest values (winter mean = 563 μS cm -1 ), with recent thermokarst lakes being intermediate (winter mean = 4228 μS cm -1 ) (Figure 4a-d). As expected, conductivity reached much lower values (averages < 1000 μS cm -1 ) for all lakes in the springtime, during and shortly after ice-cover break-up. Conductivity increased slightly during the summer, reaching mean values similar to the previous fall (~ 1000-2000 μS cm -1 )

and showing comparable values between recent thermokarst and unconnected alas lakes.

The percent saturation of dissolved oxygen for all three lake types was lowest in the wintertime, with mean values at lake surface between 0 and 20% (Figure 4e-h). Unconnected alas lakes maintained the highest levels of dissolved oxygen during fall, spring, and summer compared to the other two lake types (with more pronounced differences in spring), reaching values above 200% in summer (Figure 4h). Recent thermokarst lakes maintained the lowest concentrations of dissolved oxygen, rarely reaching saturation levels even in the springtime.

Connected alas lakes followed a similar trend as the other two lake types, but with intermediate values.

The pH values exhibited less seasonal variability compared to the other hydrochemical properties analyzed (Figure 4i-l). Unconnected alas lakes consistently had the highest pH levels relative to the other two lakes types, maintaining an average pH between 9.4 and 10 for all four seasons. Connected alas lakes and recent thermokarst lakes exhibited similar seasonal trends, with average pH values in fall, spring and summer near 8.8, and with lower values during the winter season (average 8.7).

Figure 4. Basic surface limnological properties of sampled lakes over four seasons in 2018-19, including conductivity (a-d), dissolved oxygen (e-h), and pH (i-l). Number of lakes sampled for each lake type for each season is presented in Table 2.

Seasonal Profiles.

Although temperature and dissolved oxygen profiles are not available for all lakes for each season, analysis of the profiles for selected representative lakes (unconnected alas lake #3, connected alas lake #2, recent thermokarst lake #1) (see lake locations in Figure 3) illustrates some broad seasonal trends. Profiles of temperature and dissolved oxygen show that the deeper lakes (i.e. connected alas and some recent thermokarst lakes) were stratified during winter and summer (Figure 5; Figure 6). bottom water ranging from near zero (at 2.1 m in unconnected alas #3) to 4°C (at 8.2 m in connected alas #2). During this period, dissolved oxygen saturation levels ranged between ~10 and 40% in surface waters (0-50 cm), while it decreased to values < 5% in bottom waters (Figure 5a-c-e). Conductivity and pH profiles showed smaller differences between surface and bottom waters during winter, except for recent thermokarst lake #1 showing a decrease of nearly 2000 μS cm -1 between surface and bottom (Figure 5b-d-f). These profiles illustrate the broad range of conductivity (496 to 8066 μS cm -1 ) and pH (7.4 -9.4) that can be found between individual lakes during the winter on such a lake-rich landscape. In August, available temperature and oxygen profiles for recent thermokarst lakes (lakes #1, #8, and #18) show strong stratification (Figure 6). At the time of sampling (0945-1200), the difference between surface and bottom temperature was 7.1°C (lake #1), 5.8°C

(lake #8), and 7.3°C (lake #18). Dissolved oxygen saturation levels ranged between ~28 and 54% in surface waters and between ~4 and 11% in bottom waters (Figure 6). All three recent thermokarst lakes showed similar profiles in pH and conductivity. Profiles showed decreasing pH values with depth (ranging from approximately 9 at the surface to 7.5 in bottom waters).

Conductivity profiles showed an increase with depth from approximately 2500 μS cm -1 at the surface to 6000 μS cm -1 in bottom waters for all three sampled recent thermokarst lakes. Dissolved oxygen (%) f

Dissolved greenhouse gas concentrations

Surface concentrations of dissolved greenhouse gas, expressed as departure from saturation, varied strongly between lake types and seasons (Figure 7). Broadly, unconnected alas lakes and recent thermokarst lakes had stronger seasonal differences in CO2 concentrations, while connected alas lakes showed less seasonal difference in CO2 concentration. Winter CO2 concentrations for all lake types were generally significantly different from all other seasons. Differences between unconnected alas lakes and recent thermokarst lakes were generally significant for most seasons (Table 3; Table 4).

Overall, recent thermokarst lakes consistently showed the highest CO2 concentrations between lake types, yet values varied considerably among the 15 lakes sampled and between seasons, with departure from saturation averages ranging from slight CO2 sink during spring sampling (-10 μM) to very strong CO2 source during winter sampling (+1440 μM).

Table 3. Results of the Kruskal-Wallis significance tests which compared dissolved greenhouse gas concentrations for each lake type between seasons. The results for CO2 and CH4 are grouped by lake type. The left-most and right-most columns indicate which seasons are being compared (e.g., winter CO2 was larger than fall CO2 for unconnected alas lakes). Significance (p < 0.05) is indicated by > or < signs, while = signs indicate that differences are not significant.

Unconnected alas lake Connected alas lake Recent thermokarst lake

CO2 CH4 CO2 CH4 CO2 CH4 Winter > = > > > > Fall Winter > = > = > > Spring Winter > = > = > > Summer Spring = = = > = > Fall Spring < > = = < = Summer Summer = = = = > > Fall
Table 4. Results of the Kruskal-Wallis significance tests which compared dissolved greenhouse gas concentrations of the three lake types during one season. The leftmost and rightmost columns indicate the lake types being compared. Significance (p < 0.05) is indicated by either > or < signs, depending on which lake type has the highest dissolved greenhouse gas concentration for that particular season, while = signs indicate that differences are not significant.

During the fall sampling period, unconnected alas lakes were under-saturated in CO2, and thus acting as CO2 sinks relative to the atmosphere, whereas recent thermokarst lakes were generally super-saturated and therefore sources of CO2 (Figure 7a). The CO2 concentrations in connected alas lakes were close to equilibrium with the atmosphere and thus neither a source nor a sink at time of sampling. A remarkable increase in CO2 concentrations occurred under the ice cover for all lake types. For example, recent thermokarst lake # 34 had 20 times more CO2 (1295 μM) in the winter than in the previous fall (62 μM) and connected alas lake #11 had 10 times more CO2 in the winter compared to the previous fall (190 μM in winter vs. 27 μM in fall) (Figure 7c). Bottom water concentration of dissolved CO2 was only sampled in the connected alas lake #2 during winter, showing a concentration similar than at the surface (155 μM in bottom water versus 168 μM at the surface).

During spring sampling, nearly all lakes (28 out of 33) had CO2 concentrations below saturation, with relatively low variability (Figure 7e). During late summer however (August), saturation levels at the lake surface returned to similar values as observed during the previous fall (early fall, September); unconnected alas lakes reverted to slight CO2 sinks, connected alas lakes returned to equilibrium, and recent thermokarst lakes became the largest CO2 sources again, with one particular lake (#35) showing a notably high value of nearly 500 µM (Figure 7g) (Figure 3). During summer sampling, the bottom waters of recent thermokarst lake #1 (Figure 3) showed extremely high saturation levels (> 1500 μM), nearly two orders of magnitude higher than at the surface, representing a substantial storage of CO2. Differences between surface and bottom water were not assessed in other lakes.

General trends in the concentrations of CH4 showed that recent thermokarst lakes typically had stronger seasonal differences in CH4 concentrations compared to the other two lake types. The summer sampling season showed the highest variability in CH4 concentration between lake types than the other three seasons (Table 3; Table 4).

During fall sampling, unconnected alas lakes were the strongest sources of CH4, while the other lake types were small or negligible sources (Figure 7b). During winter, as for CO2, dissolved CH4 showed a large increase in concentration for recent thermokarst and connected alas lakes (by one order of magnitude) compared to the previous fall (Figure 7d). Bottom water concentration of dissolved CH4 was only sampled for connected alas lake #2. The bottom water concentration was one order of magnitude higher than surface concentration during the winter sampling period (133 μM in bottom water versus 15 μM at the surface).

During spring sampling, all lakes (excluding a small number of outliers) were moderately saturated in CH4 compared to the values observed during winter sampling (non outliers < 60 μM), with similar saturation levels between lake types (Figure 7f). During summer sampling, observed CH4 saturation levels at the surface were lower (~ 10 μM), with unconnected alas lakes as the largest CH4 sources compared to the other lake types, which is a similar trend to the previous fall (Figure 7h). As observed for CO2 during summer, recent thermokarst lake #35 showed notably higher surface CH4 concentrations (saturation levels > 700 μM, identified by an arrow in Figure 7h) compared to all other sampled lakes. Bottom waters of recent thermokarst lake #1 (Figure 3) showed extremely high saturation levels during summer sampling (> 40 μM), nearly 130 times higher than at the surface (~ 0.4 μM), representing a substantial storage of CH4. Differences between surface and bottom water were not assessed in other lakes. 2.

Diffusive greenhouse gas fluxes

Considering only the ice-free season (from spring to fall; Table 5), the diffusive CO2 fluxes ranged from a minimum of -13.6 mmol m -2 d -1 (or approximately -600 mg C m -2 d -1 ) recorded in spring for unconnected alas lake #3 to a maximum of ~355 mmol m -2 d -1

(or ~1,200 mg C m -2 d -1 ) recorded in summer for recent thermokarst lake #18. The ice-free season median CO2 fluxes were relatively low for connected alas lakes and negative for unconnected alas lakes (1.5 and -4.8 mmol m -2 d -1 , respectively), and slightly positive for recent thermokarst lakes (1.8 mmol m -2 d -1 ). The diffusive CH4 fluxes were highest from recent thermokarst lakes in summer, reaching ~560 mmol m -2 d -1 (or 6720 mg C m -2 d -1 ). It is noteworthy that the average ice-free season CH4 flux for all studied lakes was substantially higher than the median value (14 mmol m -2 d -1 compared to 3.1 mmol m -2 d 1 ), reflecting the skewing effect of exceptionally high fluxes in summer from a few recent thermokarst lakes.

Unconnected alas lakes were regularly large CH4 sources during the ice-free season (median flux of ~7.0 mmol m -2 d -1 , compared to 1.0 mmol m -2 d -1 for recent thermokarst lakes). It should be noted that these results are estimations based on wind data from the only available meteorological station in the region (Yakutsk Airport), located ~ 120 km from our study site. Connected Alas 13 -6.9 1.5 6.4 27.6 0.2 1.4 3.9 21.0

Unconnected Alas 33 -13.6 -4.8 -0.5 98.1 0.2 6.9 10.6 56.8

All lake types 91 -13.6 -2.3 10.6 355.4 0.1 3.1 14.2 566.4

It is possible that the chemical enhancement effect contributed to increased CO2 uptake rates by some lakes during the ice-free seasons. Unconnected alas lakes consistently had the highest α values (mean value of 9.9 for unconnected alas lakes, reaching up to 27 for lake #20 in August). These large values occurred when pH was particularly high, lake water was warm and wind speeds were low. Alpha values calculated for the other two lakes presenting a CO2 uptake during August were lower (connected alas lake #11 α = 1.6; recent thermokarst lake #37 α = 5.3). During September, there were only two connected alas lakes and two recent thermokarst lakes that were acting as CO2 sinks (mean α 1.4), while the unconnected alas lake #10 presented the highest α value (3.8). The spring sampling season showed low variability in the calculated α values between lakes types, although unconnected alas lakes still had the highest mean value (2.8, compared to mean values ~ 1.0 for connected alas lakes and recent thermokarst lakes). It is important to note that during the spring season, most lakes (26 out of 31) were acting as CO2 sinks and the magnitude of CO2 uptake was highest compared to the other seasons (Figure 7). However, the highest chemical enhancement of CO2 flux was calculated for the summer season, when CO2 uptake is less prevalent and smaller in magnitude (Figure 7).

Discussion

The study area represents the typical taiga alas landscape found in Central Yakutia, a region encompassing approximately 17 000 km 2 in the Lena-Aldan interfluve (Ulrich et al. 2017a). Nearly ten percent of this landscape is comprised of lakes at varying stages of thermokarst development (Ulrich et al. 2017a). Lake morphology is directly linked to the evolution stage within the thermokarst sequence [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF], with rates of change substantially lower for older lakes (Desyatkin et al. 2009). The morphology of unconnected alas lakes has remained relatively stable during the last decades, with little lateral expansion or deepening. In contrast, recent thermokarst lakes have been notably dynamic, experiencing subsidence and lateral expansion consistent with the rapid landscape evolution observed since the mid-20 th century (Fedorov et al. 2014a;Ulrich et al. 2017a). Connected alas lakes, while they did not notably expand or deepen in the recent past, have shown some signs of thermokarst activity such as thaw slumps (Séjourné et al. 2015). This spatial heterogeneity resulting from the above chronological sequence can contribute to notable differences in lake water physicochemical characteristics, such as timing of spring ice cover breakup, water temperature and lake water mixing regime, electrical conductivity, and dissolved oxygen concentration. Depending on the lake type, we also found substantial seasonal variations in the amount of dissolved greenhouse gas in lake water, spanning two orders of magnitude throughout the year. Since most greenhouse gas studies are limited to summer quantification (e.g., [START_REF] Abnizova | Small ponds with major impact: The relevance of ponds and lakes in permafrost landscapes to carbon dioxide emissions[END_REF]Bouchard et al. 2015), such seasonal assessments are precious [START_REF] Langer | Frozen ponds: production and storage of methane during the Arctic winter in a lowland tundra landscape in northern Siberia, Lena River delta[END_REF][START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]). The wide seasonal ranges observed here underscore the need to consider such heterogeneities when attempting to estimate global greenhouse gas emissions from thermokarst lakes and permafrost landscapes.

Developmental stage as a driving factor on lake greenhouse gas concentrations and fluxes

We measured notable differences in physicochemical properties (Figure 4) and dissolved greenhouse gas concentrations (Figure 7) between thermokarst lakes of different development stages. The clear (statistically significant) difference between unconnected alas lakes and the other two types of lakes underlines the unique role played by these shallow lakes, acting as CO2 sinks most of the year and substantial CH4 sources, especially in fall.

Several studies have concluded that unconnected alas lakes are generally not sites of current thermokarst activity (e.g., Brouchkov et al. 2004;Ulrich et al. 2017a), which is consistent with our observations. Little ground ice remains below alas lake depressions, which significantly reduces the potential for input of stored carbon from surrounding permafrost (Ulrich et al. 2019). They are stable both laterally and vertically, and any variations in lake volume can be attributed to changes in the hydrological regime. As they are hydrologically isolated, the interannual water balance in these lakes is entirely dependent upon the ratio between annual precipitation and evaporation (Brouchkov et al. 2004;Ulrich et al. 2017a).

Based on these characteristics, evaporative enrichment during summer and substantial solute exclusion under winter ice cover are likely to play a role in the high concentration of nutrients and minerals observed within these lakes, such as found in Western Siberia [START_REF] Manasypov | Seasonal dynamics of organic carbon and metals in thermokarst lakes from the discontinuous permafrost zone of western Siberia[END_REF]. The notably high organic carbon concentrations (generally DOC > 100 mg L -1 except in spring) are characteristic of these lakes (Hughes-Allen et al. 2020).

In the open-water seasons, unconnected alas lakes are hotspots of biological activity, sometimes resulting in the entire surface of the lake being covered by floating aquatic plants (e.g. Lemna sp) [START_REF] Malyschez | Flora of Siberia[END_REF]. These lakes are particularly shallow (< 2 m), generally warmer than the other two lake types, and rich in nutrients. These characteristics favor high rates of photosynthetic activity, causing unconnected alas lakes to be efficient CO2 sinks (or negligible CO2 sources) during the open-water season. Negative fluxes are most prominent during the spring and fall seasons, probably linked to the more turbulent conditions over these periods leading to an efficient CO2 transfer from the atmosphere to the lake water. It should be noted that the CO2 uptake rates for these lakes might be underestimated. Chemical enhancement calculations indicate that CO2 uptake could be much higher during the summer (on average by a factor of 10) when conditions are met for this effect to be significant (high pH, high temperature, and low winds). These lakes are presenting the highest pH values (generally above 9 and up to 10.5).

Unconnected alas lakes also act as CH4 sources throughout the year, particularly in fall with flux reaching ~ 16 mmol m -2 d -1 (Table 5). These fluxes are generally much higher than what has been reported elsewhere across the Arctic (e.g., [START_REF] Abnizova | Small ponds with major impact: The relevance of ponds and lakes in permafrost landscapes to carbon dioxide emissions[END_REF]Bouchard et al. 2015;Matveev et al. 2016). High production of CH4 throughout the openwater season is likely linked to the high levels of primary productivity fueling microbes with recent, labile organic carbon to the system. Recently produced, autochthonous organic carbon may be more readily decomposed under anoxic conditions than organic carbon from a terrestrial source [START_REF] Grasset | Large but variable methane production in anoxic freshwater sediment upon addition of allochthonous and autochthonous organic matter[END_REF]. Consequently, we hypothesize that CH4 production is positively linked to the presence of autochthonous organic carbon in these lakes.

Unconnected alas lakes are shallow water bodies less strongly stratified during winter compared to the other lake types (deeper and more humic). Although limnological profiles were not collected during summer for unconnected alas lakes, it is likely that they are not as stratified as other lake types in summer as well due to their very shallow depth and susceptibility to wind-induced mixing. These conditions likely contribute to an earlier release of the greenhouse gas stored in the water column or sediments in spring and during summer.

The mean surface temperature of unconnected alas lakes in May was 11°C compared to 5°C

(connected alas lakes) and 7°C (recent thermokarst lakes), indicating that summer stratification started earlier for these lakes (Hughes-Allen et al. 2020) despite they are likely polymictic lakes. It is possible that they were presenting positive CO2 flux and higher CH4 flux earlier in spring, but since most unconnected alas lakes freeze to bottom in winter, storage flux would be mainly through greenhouse gas stored in the ice itself, or in the interstitial water of the sediment. In summer, periods of wind-induced mixing would likely generate regular venting of the greenhouse gas produced in summer by these shallow lakes, although the very high CH4 concentrations observed in fall indicate either some release of CH4 stored at the bottom or very high production rates during this period. Other studies on thermokarst lakes have documented that the timing and magnitude of greenhouse gas storage flux can vary between lakes even when they are subject to the same local climate conditions [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF].

Between the three lake types, recent thermokarst lakes generally showed the highest variability in dissolved greenhouse gas concentrations among lakes of this type (Figure 7).

This potentially results from the substantial diversity in morphology and thermokarst developmental stage for this lake type, with some lakes being much deeper and larger than others (depths from 1.2 to 4.6 m, surface areas from 0.1 ha to 2.0 ha) (Table 2). Morphology has been suggested to affect the mixing regime, and thus the seasonal patterns in greenhouse gas production versus consumption, storage and emissions [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]. For example, the spring overturn of shallower lakes may represent a more complete outgassing of all greenhouse gas stored under the ice in winter, while deeper lakes may only experience outgassing of the greenhouse gas stored at the lake surface, with deeper storage only venting out in the fall [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]). Morphology has thus been identified as a major characteristic controlling seasonal emission patterns, and deserving attention in future studies.

High concentrations of both CO2 and CH4 were measured in recent thermokarst lakes, especially under the ice cover (in winter) but also during spring and summer in some cases (Figure 7). Notably, thermokarst lake #35 (Figure 3), currently developing within a forested ground and away from any apparent human influence, had the highest surface concentrations observed in summer (> 450 μM of CO2 and ~ 1000 μM of CH4) (Figure 7g-h). Although a limnological profile is not available for this lake, strong summer stratification is expected such as is observed in other recent thermokarst lakes (Figure 6). Finding such high concentrations of CO2 and CH4 at the surface of a strongly stratified lake would indicate that greenhouse gas production was elevated in the pelagic portion of the lake, or that lateral transfer of littoral benthic production was efficient. It is also possible that this specific sampling day was following a partial mixing event bringing greenhouse gas enriched waters at the surface. The wind speed during the two days prior to greenhouse gas sampling reached sustained speeds up to 10 m s -1 , which is much higher than typical wind speeds in this area, usually near 1 m s -1 (Yakutsk Airport Metrological Station). These conditions could have contributed to partial mixing of the lake waters and the observed high concentrations of greenhouse gas in the surface waters. Overall, recent thermokarst lakes were particularly high CO2 emitters in summer and fall (Figure 7), and lower CH4 emitters than unconnected alas lakes (particularly in fall). In fact, CO2/CH4 molar ratios stayed relatively high over the summer and fall in recent thermokarst lakes (> 40) (Table 6), potentially indicating efficient methanotrophy in these systems [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]. It is also possible that such elevated ratios were caused by a lateral transfer of CO2 produced in soils under aerobic conditions (Campeau et al. 2018) or from the chemical weathering of carbonate (Zolkos et al. 2018), sources that would be particularly significant under active thermokarst erosion. The greenhouse gas production in recent thermokarst lakes is likely fueled by the thawing of ice-rich Pleistocene permafrost. Recent thermokarst lakes experience high rates of lateral and vertical expansion (Fedorov et al. 2014a;Ulrich et al. 2017a), and once initiated, this expansion occurs relatively continuously. However, the rate of expansion is strongly dependent on vegetation cover type, the presence of disturbed landscape surface (i.e., clearing of land for agriculture or infrastructure), and local permafrost conditions [START_REF] Lara | Polygonal tundra geomorphological change in response to warming alters future CO2 and CH4 flux on the Barrow Peninsula[END_REF]Ulrich et al. 2017a). The expansion of these lakes into surrounding permafrost results in the mobilization of organic and mineral matter which was sequestered for thousands of years by freezing temperatures (Hugelius et al. 2014;Strauss et al. 2017). For example, DOC values ranging from 8 to 838 mg L -1 were measured in 2018-2019(Hughes-Allen et al. 2020).

Previous work in the region suggests that old carbon is released within these lakes in the form of century-to millennia-old dissolved inorganic carbon, which is resulting from variable mixtures of late Pleistocene (20 kyr BP) and more recent or modern carbon (Ivanov 1984;Soloviev 1959;[START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]. Radiocarbon dating and stable isotopic signatures of gases, lake sediments and surrounding soils should provide more insight about the sources and pathways of the greenhouse gas emitted from recent thermokarst lakes of this region (Bouchard et al. 2015).

Seasonal variations in greenhouse gas concentrations

Seasonal differences in dissolved greenhouse gas concentrations were particularly striking at the studied site. The most notable seasonal change was the substantial increase (up to two orders of magnitude) in both CO2 and CH4 during the winter season (Figure 7).

Increasing concentrations under the ice cover have been noted in other Arctic, subarctic and boreal regions [START_REF] Langer | Frozen ponds: production and storage of methane during the Arctic winter in a lowland tundra landscape in northern Siberia, Lena River delta[END_REF][START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]. In productive systems, the isolation of the water mass under an ice cover quickly generates anoxic conditions favoring methanogenesis and suppressing methanotrophy. Our studied lakes showed median CO2/CH4 molar ratios that were much lower in winter (2.2 for all lakes; ranging from 0.2 for unconnected alas lakes to 10 for connected alas lakes) compared to summer (33 for all lakes; ranging from 4.4 for unconnected alas lakes to 56 for recent thermokarst lakes) (Table 6).

This suggests a strong production of CH4 and the predominance of hydrogenotrophy (consuming CO2 for methanogenesis) in the anoxic conditions of winter, such as proposed in other thermokarst lakes of subarctic peatlands [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF].

The greenhouse gases accumulated in the water column in late winter will be partly or completely released in spring, as soon as the ice cover degrades, a period when sampling is particularly challenging logistically. The spring overturn period has been shown to be particularly important for greenhouse gas release from lakes to the atmosphere [START_REF] Langer | Frozen ponds: production and storage of methane during the Arctic winter in a lowland tundra landscape in northern Siberia, Lena River delta[END_REF]. [START_REF] Phelps | Methane efflux from high-latitude lakes during spring ice melt[END_REF] found that greenhouse gases released during spring overturn accounted for as much as half of the total yearly emissions from Alaskan lakes. Our sampling occurred in May, likely after the major venting period, since the ice cover had already disappeared in many lakes (especially in unconnected alas lakes, which were all ice-free in early May) and surface waters had already started to warm (26 lakes were presenting temperatures above 5°C, reaching up to 17.5C). This likely explains why lake water did not present high saturations levels during this sampling period despite the high concentrations observed one month earlier under the ice cover. In fact, CO2 concentrations were mainly below saturation in May, indicating that primary producers had already started to be active.

In summer, several recent thermokarst lakes were strongly stratified (examples given in Figure 6), isolating bottom waters where oxygen depletion occurs and greenhouse gases accumulate. During August sampling, sufficient time had passed after initiation of summer stratification to observe a significant accumulation of greenhouse gases in bottom waters; for example, more than 1500 μM of CO2 and 50 μM of CH4 were measured at the bottom of recent thermokarst lake #1. This summer storage of greenhouse gas will likely be released to the atmosphere during the fall overturn period, in addition to greenhouse gas produced during the fall in sediments and water column, which can be significant considering the time lag before water and sediment temperatures cool down after air temperature has dropped [START_REF] Serikova | High carbon emissions from thermokarst lakes of Western Siberia[END_REF]. This could explain the super-saturation in CO2 (recent thermokarst lakes) and CH4 (unconnected alas lakes) measured in fall (Figure 7). Fall sampling may have been too early for the highly stratified thermokarst lakes, and thus the cumulated CH4 at their bottom would be venting later. If this is the case, surface concentration (and emissions) of CO2 would also be higher at the surface later in fall for these lakes. By contrast, connected alas lakes are relatively well mixed during the summer, as they are large water bodies influenced by wind-induced turbulence and water flows from regional streams and rivers ( Figure 6). These characteristics likely lead to the continuous release of CO2 and CH4 produced within the lake (and intermediate gas concentrations, Figure 7), rather than an isolated release of the storage flux during the fall overturn period. It is important to note that the intervals between sampling periods were not symmetrical throughout the year, and particularly short between summer (mid-late August) and fall (early September), explaining their similarities in greenhouse gas concentrations. These differing patterns underline the importance in placing automated sensors in lakes to better identify critical periods of greenhouse gas emissions [START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF].

Diffusive greenhouse gas fluxes: comparison across high-latitude regions

Diffusive greenhouse gas fluxes from lakes studied across Arctic and subarctic regions are generally limited to the ice-free season, with only a few studies providing estimations from wintertime or over a year cycle [START_REF] Langer | Frozen ponds: production and storage of methane during the Arctic winter in a lowland tundra landscape in northern Siberia, Lena River delta[END_REF][START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF].

Our estimations for diffusive greenhouse gas fluxes are based on four discrete periods over a full year cycle, including spring and fall, which have been previously identified as critical periods for greenhouse gas evasion [START_REF] Jammet | Large methane emissions from a subarctic lake during spring thaw: Mechanisms and landscape significance[END_REF][START_REF] Matveev | Winter Accumulation of Methane and its Variable Timing of Release from Thermokarst Lakes in Subarctic Peatlands[END_REF]. The diffusive CO2 fluxes estimated for this study site of Central Yakutia (range from -14 to 350 mmol m -2 d -1 ; (Table 5) are clearly outside the range reported for lakes and ponds across the Arctic. The most comparable fluxes are measured from ponds on Bylot Island (Bouchard et al. 2015) and on the Mackenzie Delta [START_REF] Tank | Northern Delta Lakes as Summertime CO2 Absorbers Within the Arctic Landscape[END_REF]. The estimated diffusive CH4 fluxes are also one order of magnitude higher (up to 560 mmol m -2 d -1 ) than published values across the Arctic (maximum values below 20 mmol m -2 d -1 for different types of lakes, including thermokarst lakes; [START_REF] Wik | Climatesensitive northern lakes and ponds are critical components of methane release[END_REF]) and almost as high as low seep ebullition fluxes from erosive margins of thermokarst lakes in Alaska and Siberia (Walter et al. 2010).

An important point to consider here is that our estimations only include diffusive fluxes, while many studies have shown that ebullition can represent a significant mode of CH4 evasion [START_REF] Bastviken | Freshwater methane emissions offset the continental carbon sink[END_REF][START_REF] Delsontro | Methane ebullition and diffusion from northern ponds and lakes regulated by the interaction between temperature and system productivity[END_REF], particularly from Arctic lakes [START_REF] Walter | Methane bubbling from northern lakes: Present and future contributions to the global methane budget[END_REF]. Some studies report that ebullition fluxes can be equivalent to diffusive fluxes, if not largely dominant in certain conditions (Walter Anthony et al. 2010;[START_REF] Wik | Climatesensitive northern lakes and ponds are critical components of methane release[END_REF], while others found that diffusion dominates in shallow thermokarst lakes (Matveev et al. 2016). The diffusive fluxes presented here already being in the upper range of published values for thermokarst lakes across the circumpolar North, they can only increase when ebullition is included, underscoring the need to consider such landscapes in global assessments.

Conclusions

Our study site in Central Yakutia is within an area of thick, ice rich permafrost, with the potential for substantial carbon mobilization in the face of continued climate change. There is considerable landscape variability and the thermokarst lakes in this region vary in age from early Holocene to the last decades. We observed temporal and spatial heterogeneity (up to two orders of magnitude) in the concentrations and fluxes of dissolved greenhouse gases from the three development stages of thermokarst lakes found in the region. There were also large differences in the relative importance of CO2 and CH4 emissions, which is fundamental for assessing their global warming potential. These differences are related in part to the mixing regime that is controlled by lake morphology and leading to variable seasonal patterns. They are also linked to variable contributions from primary producers and thermokarstic erosion shaping the light and oxygen availability. We show that diffusive fluxes of both CO2 and CH4 from thermokarst and alas lakes of Central Yakutia are among the highest presented across Arctic and subarctic regions. With their extreme bottom greenhouse gas concentrations and large but variable emissions, recent thermokarst lakes need to be closely considered as they likely involve the mineralization of ancient organic carbon that can contribute to the amplification of the greenhouse effect. On the other hand, unconnected alas lakes act as active CO2 sinks where primary production fuels methanogeny, likely dampening total emissions from the region depending on the overall greenhouse gas budget. The high levels of temporal and spatial heterogeneity between seasons and lake types illustrate the complexities of such permafrost landscapes and their responses to changes in temperature, precipitation, and anthropogenic influences. The fate of the large amounts of organic carbon stored in permafrost landscapes of Central Yakutia needs to be assessed to identify in which cases lakes contribute to its transfer to the atmosphere, especially considering the continued influence of climate change on these sensitive environments.

Automated identification of thermokarst lakes using machine learning in the permafrost landscape of Central Yakutia (Eastern

Introduction

Permafrost landscapes cover 20 million km 2 of the northern hemisphere and are particularly abundant in Siberia, Alaska, and northern Canada (Brown et al., 1997;Obu et al., 2019). Regional and local hydrological and geological factors influence its spatial distribution, thickness, and ground ice content (Grosse et al., 2013). An important feature of permafrost is its storage of enough organic carbon (OC) to significantly impact global climate if released into the atmosphere as greenhouse gas (GHG) (Strauss et al., 2017). This feature has recently propelled permafrost into the spotlight as a key component of the global cryosphere. The concern is that, as climate warming causes permafrost to thaw, the OC which has previously been sequestered by freezing temperatures will be mineralized and released as carbon dioxide (CO2) and methane (CH4) (Schuur et al., 2015;Hughes-Allen et al., 2021).

The greenhouse effects of these two compounds will amplify current warming trends, causing more permafrost thaw, subsequent OC release, and so on (Schuur et al., 2015). It is estimated that permafrost landscapes currently store approximately 1,600 billion tons of carbon, more than twice the amount that currently exists in the atmosphere today (Hugelius et al., 2014).

The idea of substantially increasing the amount of carbon which currently exists in the atmosphere if all the carbon stored in permafrost is released is certainly a concern. Climate change and other human activities have already had measurable impacts on permafrost landscapes. Researchers have recorded deepening of the active layer (the surface layer of soil on top of permafrost which freezes and thaws annually) (Park et al., 2016), increased thawing and slumping (Nitze et al., 2018), as well as increases in the number and surface area extent of thaw lakes (Nitze et al., 2017). In addition to carbon emission, permafrost thaw destabilizes infrastructure and transportation and can render farmland unusable, a concern that is likely to get even more pressing by the middle of this century, with considerable costs (Hjort et al., 2022). Areas of continuous permafrost (where permafrost underlays 90-100% of the landscape) and high ground-ice content (50-90% by volume) are particularly sensitive to changes in temperature, precipitation, and other human disturbances like forest clearing for agriculture (Grosse et al., 2013;Ulrich et al., 2017a). Not only are permafrost landscapes particularly sensitive to climate warming, but the magnitude and rate of temperature rise across the Arctic is much higher (2-3 times) compared to global averages, a phenomenon referred to as 'Arctic amplification' (Serreze and Barry, 2011). Mean annual air temperature in the Arctic is predicted to rise by as much as 5.4 ºC within the coming century in the absence of significant and directed global efforts to reduce GHG emissions (Pörtner et al., 2019).

However, like most natural Earth systems, permafrost landscapes are spatially heterogeneous and complicated. For example, landscape type (waterbody, forest, grassland, etc.) greatly affects the GHG emission characteristics of a particular area. Desyatkin et al. (2009) found large differences in CH4 emissions when comparing forest, dry grassland, wet grassland, and pond surfaces. Pond surfaces were found to have CH4 emissions more than two orders of magnitude greater than the other landscape types.

Type of waterbody and season can also cause significant differences in GHG emissions. Hughes-Allen et al. (2021) found that recent thermokarst lakes (lakes formed within the last few decades mostly from anthropogenic climate change and other human activities) released consistently higher levels of CO2 to the atmosphere in all seasons compared to the other lake types. Small, hydrologically unconnected alas lakes (residual lakes that exist in former lake depressions) acted as CO2 sinks during fall and spring, drawing CO2 out of the atmosphere, but acted as CO2 sources during summer. All lake types released CH4 to the atmosphere during all three ice-free seasons. Again, recent thermokarst lakes had higher values than the other lake types in summer, while unconnected alas lakes had higher values in fall and spring. Such striking temporal and spatial heterogeneities in GHG dynamics have also been observed elsewhere across the Arctic region, for example in Northern Canada (e.g., Bouchard et al., 2015;Prėskienis et al., 2021). The relationship between permafrost thaw and GHG emissions is complicated and nuanced by local hydrology and geomorphology.

Thermokarst lakes are often hotspots for GHG emissions from permafrost landscapes (Desyatkin, 2009;Walter Anthony et al., 2016;Hughes-Allen et al., 2021;Prėskienis et al., 2021). Thermokarst processes are generally linked with disturbances such as warming temperatures or forest removal for agriculture or by wildfires, which cause deepening of the active layer (Grosse et al., 2013). When such deepening induces melting of ground ice, which is often the case in areas of ice-rich permafrost, then thermokarst processes may start. These lakes profoundly change the local ground thermal regime, sometimes increasing surrounding sediment temperatures by as much as 10º C above the mean annual air temperature (Grosse et al., 2013). Thermokarst processes amplify permafrost thaw compared to what would be expected from rising air temperatures alone (Brouchkov et al., 2004;Schuur et al., 2015).

Lake expansion and deepening will generally continue until the accumulation of lake sediments over time creates an insulating layer between the lake water and surrounding permafrost and/or the lake becomes deeper than the layer of ice-rich permafrost (Grosse et al., 2013;French, 2017). Once the lake is no longer expanding, its size and depth are controlled by surface and subsurface inflows/outflows, as well as the balance between precipitation and evaporation. Drainage (progressive or catastrophic), evaporation, terrestrialization, and infilling will eventually result in lake disappearance (Desyatkin, 2009;Grosse et al., 2013;Bouchard et al., 2017).

Thermokarst lakes are important contributors to the global carbon cycle, however, global lake inventories used in Earth system modeling are strongly biased toward larger lakes and generally only include lakes greater than 10 ha (Verpoorter et al., 2014). The results from Hughes-Allen et al. (2021) andothers (in 't Zandt et al., 2020;Elder et al., 2021) show how important small unconnected alas lakes (mean area from = 5 ha) and recent thermokarst lakes (mean area = 0.5 ha) are to the carbon cycle in permafrost landscapes. Also important is understanding long-term changes in the number, distribution, and size of thermokarst lakes in permafrost landscapes. The limited studies which have conducted long-term analyses of thermokarst lake distribution in permafrost landscapes have found that the number and size of lakes in areas of continuous permafrost have generally increased in recent decades (Tarasenko, 2013) Some areas of continuous permafrost, including the lower Mackenzie River, Canada and northern Alaska, have experienced declines in lake number and areas (Travers-Smith et al., 2021;Chen et al., 2022). Interestingly, areas of discontinuous or isolated permafrost have seen stable or decreasing overall lake areas due to lake drainage (Grosse et al., 2008;Nitze et al., 2017. Remote sensing techniques using satellite images have become a powerful tool for analyzing lake area change in the expansive regions of continuous permafrost found in Eastern Russia.

Until recently, remote sensing studies of permafrost lakes have been limited to comparisons of imagery spanning relatively narrow time frames (Karlsson et al., 2014;Boike et al., 2016;Nitze et al., 2017). This has been due, in large part, to the lack of high-resolution imagery available at sufficient and regular time intervals. However, a combination of SPOT (Satellite pour l'Observation de la Terre) imagery and declassified American surveillance satellite imagery can provide a lengthy, high-resolution record of permafrost landscapes in Central Yakutia.

In this study, we present a long-term analysis of lake surface area change within a 4,500 km 2 area of Central Yakutia, Russia. The identified lakes are also classified based on the lake type designation developed in Hughes-Allen et al., 2021. The main objectives of this study are: 1) quantify changes in lake surface area (overall and for the three different lake types) through 1967-2020 timeframe, 2) compare changes in surface area to historical precipitation and temperature data, 3) identify trends in the spatial distribution of lake types and lake development over time, and 4) test the hypothesis that unconnected alas lakes are more susceptible to changes in precipitation than the other two lake types.

Study Site

Central Yakutia experiences extreme subarctic continental climate with long, cold, and dry winters (Coldest month being January having a mean temperature around -40°C) and warm summers (July mean temperature around +20°C), causing strong seasonal variability.

The winter season (defined by the presence of an ice cover on lake surface) usually lasts from early October until early May. The small amount of precipitation (150-250 mm) that does accumulate each year usually falls during the summer season. Average snow depth for winter months (January to April) ranges from 24 cm in January to a maximum of 30 cm in March, and then decreasing to 10 cm at the end of April (1980-2020 recorded values from Yakutsk weather station). The snow which falls in this region generally has very low water content due to cold temperatures (Fedorov et al., 2014). Yearly evaporation rates exceed total precipitation in this region (Fedorov et al., 2014). Central Yakutia, like other high latitude regions, is warming disproportionately faster than lower latitudes. Between 1996 and 2016, the mean annual air temperature of Central Yakutia has increased by 0.5-0.6°C per decade (Gorokhov and Fedorov 2018). Spring snow cover has been disappearing 3.4 days earlier per decade over the pan-Arctic terrestrial region and climate models predict decreases in snow cover duration between 10-20% by 2050 (Czerniawska and Chlachula, 2020). Changes in average snow depth and precipitation, however, are highly spatially heterogeneous with some areas of Eurasia experiencing increasing snow depth totals and precipitation (Czerniawska and Chlachula, 2020).

The study site (62.55 °N; 130.98 °E) lies approximately 130 km north-east of the city of Yakutsk on a lowland plain between the Lena River to the west and the Aldan River to the east. (Figure 1). The region is covered mostly by late Pleistocene sediments, including silty clays and sandy silts of fluvial, lacustrine or aeolian origin (Ivanov, 1984). Since the Pleistocene, numerous fluvial terraces have been formed from the activity of the Lena and Aldan rivers, and their smaller tributaries (Soloviev, 1959). Two Pleistocene-age fluvial terraces underlay this region: the Tyungyulyu terrace, which covers the western section of the study area, 50-200 m above sea level (asl), dated 14-22 kyr BP, and the higher Abalakh terrace in the eastern sector of the study area, 200-280 m asl, dated 45-56 kyr BP (Ivanov, 1984;Soloviev 1959). This region is dominated by larch, pine, and birch forests and is characterized as a middle taiga landscape regime (Fedorov et al., 2014). Grasslands are abundant in unforested areas, including land previously cleared for farming, ranching or in the remnant depressions of old thaw lakes known as 'alases'. They consist of halophytic steppe-like and bog plant communities (Ulrich et al., 2017b). Permafrost in this region is continuous (Figure 1a), thick (> 500 m deep), and the upper 30-50 m (Pleistocene-age fluvial and aeolian sediments called 'Yedoma') can be extremely rich in ground ice (50-90% by volume) (Ivanov, 1984). The amount of OC stored in Yedoma varies widely. For example, deep cores from Northern Siberia and Alaska yielded OC pool estimates of approximately 10 +7/-6 kg/m 3 (Strauss et al., 2013). A 22 m deep core in Central Yakutia (Yukechi) on the Abalakh terrace yielded a much lower value of OC content of ~5 kg/m 3 (Windirsch et al., 2020), while another Central Yakutian study (Spasskaya Pad/Neleger site) of a shallow core (2 m) showed a considerably higher OC content of 19 kg/m 3 for the top two meters of larch forest covered Yedoma deposits (Siewert et al., 2015).

Yedoma silty loams, which are common to the Lena-Aldan interfluve, underlay much of the study site, with abundant ground ice in the form of 1.5-3 m-wide ice wedges. Active layer depth in the region generally ranges between ~ 1 m below forested areas to > 2 m in exposed grassland areas (Desyatkin, 2009). Zones of unfrozen ground (or taliks) exist underneath major rivers and lakes whose depth exceeds that of the ice-cover in the winter.

Nearly half of the landscape has been affected by thermokarst since the early Holocene, resulting in the formation of thousands of partly drained alas depressions (Soloviev, 1959;Brouchkov et al., 2004). However, recent thermokarst activity related to natural landscape evolution, increasing air temperatures and/or human-induced landscape modifications (agriculture, clear-cutting, and infrastructure) is also widespread in the region. There are numerous small, recently developed, and expanding lakes and retrogressive thaw slumps along lake shores (Fedorov et al., 2014;Séjourné et al., 2015).

The thermokarst lakes in this region can be divided into three broad categories based on field observations, past radiocarbon dating of lake sediments, geochemical signatures of lake waters, morphology, and a multiple-stage development model [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]Desyatkin et al., 2009). An illustrative example from the area is presented in Figure 2. Lakes of each type have strong differences in lake physiochemistry, dissolved greenhouse gas concentration and greenhouse gas flux. Defining the nuances of lake development based on lake type in Central

Siberia is crucial to understanding potential impacts of this permafrost landscape on the regional and global carbon cycle. The three lake types are as follows (illustrated in Figure 2):

• Unconnected alas lakes: residual water bodies located within hydrologically closed basins (Desyatkin et al., 2009). Most of these lakes likely formed during the transition between the Pleistocene and Holocene, approximately 10-8 cal kBP or during the Holocene Thermal Maximum (~6.7-5 cal kBP) (Biskaborn et al., 2012;Ulrich et al., 2017b). These lakes can be up to a few meters deep but are typically very shallow (1 meter deep or less) and are thus generally frozen to the bottom in winter. The ancient lake depressions surrounding the small residual lakes of this type can be up to several kilometers wide and several meters deep and are relatively easy to distinguish on satellite images. These alas lakes have already undergone much of the thermokarst processes and very little ground ice typically remains beneath the residual lake.

Therefore, the thaw potential and resulting input of stored carbon to these lakes is low compared to recently formed thermokarst lakes (Ulrich et al., 2019). Yakutian people have used these depressions for agricultural purposes for centuries.

• Connected alas lakes: lakes connected hydrologically to the watershed by streams or rivers. These lakes are consistently larger (several hundreds of meters across) and deeper (up to ~ 10 meters). Most of them were probably formed during the mid-Holocene, approximately 5-3.5 thousand years ago, although detailed chronology about their inception is still incomplete [START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]Ulrich et al. 2017b). Local people currently use some of these lakes for fishing.

• Recent thermokarst lakes: thaw lakes formed over the last several decades mostly from human activities (e.g., forest fire and forest removal for agriculture, pipelines, or road construction) and rising temperature (Fedorov et al. 2014). These lakes are generally small (meters to tens of meters across) and relatively shallow (generally one to two meters deep) and are still expanding downwards and laterally due to active layer deepening and thermokarst processes. On satellite images, these lakes lack the higher albedo lake depression characteristic of unconnected alas lakes and appear as low

Image Data

In this study, we leveraged the entire archive of SPOT data available for the study region between 1986 and 2016 (Table 1). Developed by the Centre National d'Études Spatiales (CNES), the SPOT family includes five decommissioned satellites which operated between 1986 and 2015 (SPOT 1-5) and two operational satellites, SPOT 6 and SPOT 7, which were launched in 2012 and 2014, respectively. The SPOT images were filtered to include only months between June and October and a cloud cover less than 70 percent. In Central Yakutia, ice cover on the lakes can persist into mid-May and start as early as late October. Only icefree scenes were used in this study. Unfortunately, many of the images which were taken between June and September had very high percent cloud cover and could not be used for analysis. Frequent high percent cloud cover combined with a satellite return rate of 26 days resulted in a surprisingly small subset of workable scenes. We were able to augment this dataset slightly by including declassified military intelligence photographs (Table 1). The declassified military satellite systems code-named CORONA, ARGON, andLANYARD operated between 1960-1972 collecting photographs of the USSR and China.

These satellites carried a variety of camera systems. The early systems (KH-1, KH-2, KH-3, KH-6) carried a single panoramic camera or a single frame camera (KH-5). KH-4, KH-4A, KH-4B carried two panoramic cameras with a separation angle of 30º with one forward looking and one aft looking camera. These images are not georeferenced.

Georeferencing was completed in QGIS 3.16 (QGIS, 2022).

Defining lake boundaries and lake types

Lake boundaries were frequently defined by a clear demarcation between water and land.

In the instances where the water/land boundary was ambiguous, every effort was made to delineate all liquid water associated to the lake. However, it was sometimes unclear whether darker pixels surrounding a lake were liquid water or heavily saturated mud (Figure 3). In these cases, best judgment was used to include all pixels which correspond to lake area. Some subjectiveness is inherent in this process. Frequently, unconnected alas lakes develop into a half moon shape or a peripheral ribbon of liquid water surrounding dry ground (Figure 4). In these cases, only the liquid water was included in the lake area. In QGIS, lake types were manually assigned to the lake polygons. In the absence of field observations, past radiocarbon dating of lake sediments, or geochemical information of lake water, lake type classification was determined solely based on lake morphology.

Connected Alas (referred further as CA) lakes were easily identified based on their generally large size and the presence of inflow and/or outflow rivers and streams (Figure 5). Some CA lakes experienced such significant surface area reduction in some years that they were reclassified as UnConnected Alas lakes (referred further as UCA) for those particular scenes.

UCA lakes were identified based on their characteristic surrounding dry depression (Figure 5). The size of these depressions varies year to year and lake to lake depending on precipitation levels, the lake's phase in the multiple-stage development model, and surrounding topography. Recent Thermokarst (RT) lakes were generally small and directly surrounded by forest or other vegetation cover (Figure 5). For a small number of lakes (<10 in each scene), it was difficult to determine whether it was an UCA lake or a RT lake. In these cases, lake morphology and surrounding environment were carefully considered, and a best guess decision was made. method that is used to identify different objects in an image (i.e., pedestrians on a sidewalk, animals in a field, etc.) (He et al., 2017). Our implementation builds on top of the existing reference PyTorch implementation (Paszke et al., 2019). The backbone of the neural network is 'resnet50', a convolutional neural network that is pre-trained on the "Microsoft common objects in context (COCO)" dataset (Lin et al., 2014).

The neural network can receive an input image between 800-1333 x 800-1333 pixels.

Each satellite image, however, is approximately 30,000 x 30,000 pixels. Therefore, every satellite image was split into ~900 smaller images (depending on original image size) and the neural network treated every image separately. This resulted in lake predictions which exhibited undesirable artifacts. For example, a lake spanning two or more small images was artificially divided into smaller polygons. To alleviate this problem, each scene was split a second time into ~900 smaller images. These images were positioned in a staggered overlap of the 'base' images such that the corners of these 'overlapping' images were in the middle of each of the base images (Figure 6). Each of the image sections (base and overlapping) were treated identically by the model. This process facilitated the fusing of lake polygons which spanned multiple small images into single polygons, which occurred later in the process.

Fine tuning and training

The construction of the AI model was a three-step process of initial fine tuning using a very limited data set, a second round of fine tuning using a substantially larger dataset, and lastly, the full model training using four complete scenes. The first step was fine tuning of the model using a limited amount of data to get preliminary results. Fine tuning refers to using the saved parameters of a model which has been previously trained on a different dataset for a new deep learning process. The fine-tuning dataset was created using a SPOT 7 image polygon was predicted where a lake exists) were calculated. The false prediction rate is the sum of the false positive rate and the false negative rate. In contrast, the relative error in total predicted lake area is the difference of the false positive and negative rates, which is therefore bounded by the false prediction rate.

Surface area change analysis

The available scenes were separated into two general study areas to facilitate lake surface area change analysis: center and south (Figure 10; Table 3). Each area was approximately 40

x 40 km 2 . These study sites were drawn to maximize the temporal and spatial coverage of available scenes and to statistically sample enough of the different types of lake present in the area. All lakes with available in situ measurements (as presented by Hughes-Allen et al., 2021 andcompiled in Hughes-Allen et al., 2020) are within the center study site (Figure 10).

Lake type assignment followed the same method as described above. Version 2 polygons were manually corrected and used to calculate lake surface area for all analyzed scenes.

Overall lake surface area and lake surface area by lake type were then compared between scenes for each study site. Additionally, the complex morphology of this lake was not well predicted by the algorithm for any scene and manual correction of the lake was time consuming and laborious.

Temperature and precipitation

An exceptionally long record of temperature and precipitation data is available from the Yakutsk airport (62.0866° N, 129.7500° E) (1888-present). These data were compiled from daily records to monthly sums (precipitation) and monthly averages (temperature) from 1900-2020. To account for all precipitation that might have influenced lake surface area, the year start date was shifted to October 1 st of the previous year. For example, for a scene taken on September 1, 2000, the yearly precipitation would have included precipitation data from October 1, 1999-September 31, 2000.

These data were then compared to a 30-year moving average (minimum window 10 years). The Mann-Kendall Trend test was used to determine whether a trend existed in the temperature and precipitation data. All analyses were completed using the Python programming language (Python Software Foundation, http://www.python.org/).

Spatial distribution of lake types

In both study sites and every scene, UCA lakes were by far the dominant lake type in both count and total surface area (Figure 15). The mean surface area for UCA lakes was 3.0 ha, which is consistent with the findings of Hughes-Allen et al. ( 2021) (Table 4). UCA lakes were widely distributed throughout each scene and there was no apparent trend to their locations. RT lakes were the second most abundant lake type in terms of count, although they accounted for proportionally less of the total surface area due to their generally small size (Table 4). Although some RT lakes do appear to have formed in the absence of any human disturbance (see further in, i.e., Figure 18), most of these lakes form within 1-2 km of roads or cleared land for pastoral practices or infrastructure development. Clusters of RT lakes can be seen near the City of Borogontsy (Борогонцы) and the village of Syrdakh (Сырдах). RT lakes which formed beyond the obvious influence of human activity are likely the result of the considerable climate warming observed in central Yakutia. CA lakes were the least numerous, but second in terms of total surface area. The average surface area of CA lakes was 60 ha (median 3.3 ha), much larger than either RT lakes or UCA lakes (Table 4). Some thawing and thermokarst lake formation (Grosse et al., 2013). Direct impacts of land cover removal and infrastructure development are usually limited to within 100 m of the disturbance, but the effects can last for decades despite revegetation (Yu et al., 2015).

By 2007 (August 2), the surface area of UCA lakes had increased to ~100 ha and the number of RT lakes had increased to 12 and surface area to 5 ha. A RT lake has formed in the middle of the of the image seemingly in the absence of human interference like clearing for agriculture or a road. In 2010, UCA lake surface decreased slightly to 87 ha, while RT lake surface area remained at 5 ha. In 2011 (September 9), UCA lake surface decreased slightly to 62 ha and RT lake surface area increased to 5.2 ha. Despite the small increase in RT lake surface area, a RT lake disappears completely between 2010 and 2011 (Figure 18). It is not clear from the images what could have caused the demise of this lake. The observed trend of increasing or stable RT lake surface area compared to decreasing UCA lake surface area is perhaps related to RT lake expansion into surrounding permafrost through thermal erosion and thaw slumping, while the surface area of UCA is controlled by the balance between evaporation and precipitation rates. From 1990 onwards, there is an increase in the frequency of years which have above average yearly temperature and every year after 2009 had above average MAAT (Figure 13).

It is possible that this trend can explain the observed reduction in UCA lake surface area and the increase in number and surface area of RT lakes. Higher temperatures result in higher evaporation rates, which will cause decreases in UCA surface area while at the same time promoting RT lake expansion into surrounding permafrost.

Lake surface area change: Center study site

The 1967 (September 20) scene had substantially lower total lake surface area compared to the other five scenes (Figure 19). Many alas basins are occupied only by a very small, The unexpectedly high CA lake surface area value for the 1980 (September 20) scene is related mostly to a single large lake (indicated by the black arrow in Figure 12). By 2010, this lake had lost approximately half its surface area (becoming an UCA lake) and is relatively non-existent in subsequent scenes. This large lake is surrounded by substantial agricultural activities, and it is possible that it is used to irrigate nearby fields. The proliferation of farming in the area might have caused the drainage and eventual demise of this large lake.

There is an overall trend of increasing RT lake surface area through time. The 2016 (September 11) scene has a substantially greater number and a higher surface area value for RT lakes compared to the other four scenes. As previously stated, this is likely related to increases in temperature as well as other human disturbances.

Conclusions

• Mask R-CNN instance segmentation method is an effective and efficient way to delineate the lake polygons of large satellite images.

• Correction of the polygons generated by the Mask R-CNN was much less time consuming than manual digitization. Manual digitizing one 60 x 60 km 2 scene without the aid of the polygons generated by the neural network would take at least one full week of motivated work. Correction of the polygons generated by the neural network for one scene could be completed in two-three hours or less.

• The limited availability of clear scenes and the single band nature of the images made automatic detection of lake polygons difficult and more fine tuning is still needed before this method can be fully implemented at a large scale.

• In particular, the texture of scenes before 2005 was substantially different from scenes after 2005. The algorithm was only trained on scenes after 2005. An additional round of training is needed to be able to maximize use of the CORONA dataset.

• However, interesting trends are still visible. Overall lake surface area seems closely tied to changes in yearly precipitation and UCA lakes are particularly sensitive to these patterns. These lakes are hydrologically isolated, and their surface area is controlled only by evaporation and precipitation. RT lakes and CA lakes were less affected, and their lake levels are controlled by expansion into surrounding permafrost and connecting streams and rivers, respectively.

• Many RT lakes formed as a result of human disturbance (forest removal, road building, etc.). Some RT lakes, however, formed in the absence of any disturbance, likely because of climate warming.

Conclusions and perspectives

The work prepared and conducted during this thesis combines paleolimnological analysis, in situ observations and measurements, and remote sensing analysis to provide a comprehensive understanding of lake development in permafrost landscapes and the contribution of these lakes to the global carbon budget. Understanding the developmental history of central Yakutian lakes, particularly in terms of the OC accumulation preservation, can provide insights into the lability of stored permafrost carbon and the consequences of its liberation by warming temperatures. There are few studies which consider the seasonality and temporal variability of dissolved GHG concentrations and fluxes, and this work provides a first step toward filling this gap. The automation of lake polygon identification in Central Yakutia from satellite images will facilitate efficient and effective temporal analysis of lake surface area change in response to short-and long-term variations in temperature and precipitation and local human disturbances. The time consuming and odious task of hand digitizing lake polygons is no longer a barrier to long term temporal analysis of lake surface area change, especially for small lakes that are usually not considered in regional/global mapping efforts.

The first part of this PhD work (Chapter 2) consisted of a detailed and high-resolution analysis of an approximately 7m long sediment core from Lake Malaya Chabyda in Central Yakutia. This sediment core spans the Pleistocene-Holocene transition and several key past climactic shifts. Based on the proxy analyses, organic matter deposition in the lowest section of the core was dominated by terrestrial vegetation. This trend changes between 12.5 cal kBP to 11.0 cal kBP to organic matter deposition dominated by lacustrine algae. Proxy analysis indicated that the upper two sections of the sediment core experienced high levels of deposition of well-preserved organic matter, which had not undergone substantial decomposition before incorporation into the sediment column. In the upper two sections of the sediment core, increases in lake depth and nutrient availability from the catchment increased bioproductivity within the lake and organic matter preservation and storage relative to decomposition. This sediment core has high carbon content, including TOC (wt%), compared to other similar sites in Central Yakutia, and elsewhere across the Arctic. The OCARs are above the highest reported values for temperate and high-latitude regions, for both past (Holocene and Late Pleistocene) and modern conditions. These proxy analyses indicate that Lake Malaya Chabyda might have acted as an efficient OC sink since the Pleistocene-Holocene transition, a conclusion rarely seen in the current 'permafrost-carbon' literature.

There is still ample information in this sediment core which has not yet been investigated.

In the context of a post-doc position, laboratory analysis to determine mercury (Hg) concentrations along the length of the sediment core will be conducted and combined with the sedimentological and paleolimnological analyses conducted during this thesis. In a warming Arctic, heavy metals, including Hg, are released from permafrost soils and may be introduced to the hydrological system and food chain as temperature and other disturbances cause permafrost thawing. While Hg is a naturally occurring element, its concentration in atmospheric deposition has increased threefold since industrialization (Streets et al., 2011).

Non-vascular plants, which are abundant in permafrost landscapes, take up particularly high levels of atmospheric Hg compared to vascular plants (Olson et al., 2019). This has resulted in the significant accumulation of Hg in permafrost soils (Obrist et al., 2009).

In certain conditions such as the anoxic waters common to waterbodies found in permafrost landscapes, Hg can be transformed by bacteria to CH3Hg + , which is a quite hazardous and toxic form of organomercury. Organomercury compounds are easily transferred across biological membranes and bioaccumulated and there is thus concern that increasing concentrations of Hg in permafrost landscapes due to climate warming might have negative impacts on local communities (Ha et al., 2017). Several studies indicate the importance of surface permafrost soils in Hg release in a warming climate (Obrist et al., 2009;Schuster et al., 2011;Burke et al., 2018;Lim et al., 2019;Schaefer et al., 2020;Schirrmeister et al., 2020), but there is currently limited knowledge about Hg concentrations in deep permafrost. This future study will improve our understanding of Hg concentrations in deep permafrost.

The second axis of this thesis (Chapter 3) investigated the contribution of thermokarst lakes to GHG emissions in permafrost landscapes. These lakes are hotspots of greenhouse gas emissions, but with substantial spatial and temporal heterogeneity across the Arctic and sub-Arctic. We measured dissolved CO2 and CH4 concentrations in thermokarst lakes of Central Yakutia and their seasonal patterns over a yearly cycle. We observed temporal and spatial heterogeneity (up to two orders of magnitude) in the concentrations and fluxes of dissolved GHG from the three development stages of thermokarst lakes found in the region. There were also larger differences in CO2 and CH4 concentrations and fluxes between lake types and seasons. These differences were likely controlled by the mixing regime of the lake, which is related to lake morphology and seasonal patterns. Results from this study found diffusive fluxes of both CO2 and CH4 from thermokarst and alas lakes in the study site are among the highest presented across Arctic and subarctic regions. With their extreme bottom GHG concentrations and large but variable emissions, recent thermokarst lakes need to be closely considered as they likely involve the mineralization of ancient organic carbon that can contribute to the amplification of the greenhouse effect. On the other hand, unconnected alas lakes act as active CO2 sinks where primary production fuels methanogeny, likely dampening total emissions from the region depending on the overall greenhouse gas budget.

affected by changes in temperature and precipitation. This study indicated that unconnected alas lakes are more sensitive to annual changes in precipitation. The total surface area of unconnected alas lakes was reduced after years of low precipitation and possibly related to the trend of increasing temperatures in Central Yakutia. This trend has implications for GHG emissions from permafrost landscapes in Central Yakutia as unconnected alas lakes were identified as strong emitters of CH4. Connected alas lakes are more resilient to these changes thanks to inflow from streams and rivers. A trend of increasing temperature in Central Yakutia seemed to increase the number and total surface area of recent thermokarst lakes.

These lakes are expanding into surrounding permafrost, which increases their size vertically and laterally. Recent thermokarst lakes often developed quite quickly, over the span of onetwo years in some instances and were often found within 1 km of human disturbances like roads or farmland.

Development and refinement of this part of the thesis is ongoing and a manuscript is currently being prepared for submission to a peer reviewed journal. Implementation of the methods used in this study required many iterations and adjustments before satisfactory levels of false prediction rate were achieved. The single band nature of the satellite images and the inclusion of grainy declassified military intelligence photographs complicated the lake polygon detection process. An additional round of fine tuning of the model with an older image (ex. 1980) is likely necessary to further reduce the false prediction rates of images take before 2000. It will be interesting to compare the lake dynamics in the north study site to those of the two other study sites after the analysis of additional images. A summer 2019 scene has been ordered and will also be included in the analysis for all three study sites. The inclusion of this scene will also make it possible to extrapolate in situ GHG concentration and flux measurements to a larger area. This upscaling analysis will help improve understanding of the contribution of thermokarst lakes to the global carbon budget.

Figure 1 .

 1 Figure 1. Global permafrost distribution based on spatial extent.(https://www.eea.europa.eu/legal/copyright)

Figure 3 .

 3 Figure 3. Active layer in a forested area near Syrdakh village (Sakha Republic). Frozen ground is visible below the tape measure.

Figure 4 .

 4 Figure 4. Diagram of generalized permafrost structures. The active layer is generally deeper under non forested areas. A talik will often form beneath lakes. These taliks can be open (one side is in contact with the active layer) as illustrated or closed (surrounded on all sides by permafrost).

Figure 5

 5 Figure 5. a) Ice wedges exposed by a thermo-erosion gully near Syrdakh, Yakutia. b) Crack in the soil c) Hummocky landscape (Syrdakh, Yakutia)

Figure 6 .

 6 Figure6. Thermokarst lake formation and evolution in ice-rich, continuous permafrost. a) Thermokarst lake inception generally starts with pooling of water between troughs and melting of ice wedges. b) Coalescence of the small ponds creates small, shallow lakes. c) Thermo-erosion and thaw slumping expand the lake vertically and laterally. d) Mature lake stage is reached when the surrounding permafrost has been completely thawed. (Adapted fromBouchard et al. (2017)

Figure 7 .

 7 Figure 7. Distribution of soil organic carbon content (Richter-Menge et al., 2019)

  through at least four major pathways: ebullition flux, diffusive flux, storage flux, and flux through aquatic vegetation (plant mediated transport) (Figure 8). All four pathways are dependent on a multitude of variables including waterbody size and morphology, season, presence and type of aquatic vegetation, and local weather patterns.

Figure 8 .

 8 Figure 8. Four primary GHG emission pathways (ebullition flux, diffusive flux, storage flux, and flux through aquatic vegetation). (From Wik (2016)).

Figure 1

 1 Figure 1 a) Map showing general location of the study site (red box) within the Russia Federation (green) b) Location of study site (blue box). Note the Lena River and the city of Yakutsk (yellow star) c) Locations of Lakes Ulakhan Chabyda and Malaya Chabyda, near the city of Yakutsk. Lake catchment is outline in black. Coring site location is indicated by red circle.

Figure 2 .

 2 Figure 2. Diagram showing the original eight sections of sediment core.

Figure 3 .

 3 Figure 3. δ13C curves for each of the two treatments. The blue curve represents the high temperature treatment (90ºC). The red curve represents the low temperature treatment (50ºC).

(

  based on room temperature HCl (10%) test), each one being approximately 1 cm thick and continuous across the width of the sediment cores. The light-colored laminations are notably visible at three sections of the unit, i.e., between 550 and 505 cm depth, and near 450 and 435 cm depth. Above 400 cm depth, the core transitions to homogenous light colored clay characteristic of Unit 3. • Unit 3 (376-0 cm) (9.0 cal kBP-CE 2013). This unit is uniformly lighter brown, silty, homogenous (non-laminated) clay. Traces of oxidation were observed between ~ 265 and 225 cm depth.

Figure 6 .

 6 Figure 6. Output of the hierarchical agglomerative cluster analysis used to delineate the three core units.

Figure 7 .

 7 Figure 7. Analysis of mean and median grain size, as well as sorting and grain size distribution (% of clay, silt and sand) of 16 samples.

Figure 8 .

 8 Figure 8. Biogeochemical proxies measured of 78 discrete samples, including total organic carbon (TOC), total nitrogen (TN), TOC/TNatomic ratio, total inorganic carbon (TIC) and δ13C trends (vertical dashed lines represent the usual isotopic signature of C3 land plants; see Figure Figure 10).

Figure 9 .

 9 Figure 9. Carbon (wt%) versus nitrogen (wt%) ratio by unit.

Figure 10 .

 10 Figure 10. δ 13 C versus TOC/TNatomic ratio of 78 discrete samples. Lacustrine algae, C3 land plant, and C4 land plant designations from Meyers and Teranes, 2006.

  and associated ratios. CLR transformed values of Rb and Zr have almost identical profiles (Figure 11). Unit 1 has the highest values of Rb and Zr with low variability. Unit 2 and Unit 3 have low values and very high variability. Ti has similar valuesin Unit 1 and Unit 3, but there is a notable decrease in Unit 2, with higher variability. Si values increase slightly through Unit 1, then decrease from the bottom of Unit 2 until the middle of Unit 3. Si values increase from 250 cm depth (6.7 cal kBP) to 100 cm depth (2.7 cal kBP), where the values remain consistent until the top of the core (Figure11). Ca and Sr are high in Unit 1 with relatively low variability (Figure11). The bottom half of Unit 2 shows a decrease in Ca and Sr; values return to similar levels seen in Unit 1 in the upper half of Unit 2. These values decrease into Unit 3, where they are low with low variability. Ca/Ti, follows a similar trend to Ca. Unit 2 has the highest values, while Unit 3 has consistently low values and low variability (Figure11).

  throughout Unit 1. Fe/Ti increases in the lower half of Unit 2, reaching the highest values at approximately 550 cm depth (~11.7 cal kBP), before decreasing in the upper half of Unit 2 and stabilizing in the lower half of Unit 3. Fe/Ti values decrease consistently throughout the upper half of Unit 3. Oxygenation of the water column is represented by Mn/Fe (Melles et al., 2012). The highest values are present in the lower half of Unit 2 around 500 cm depth (~11.0 cal kBP). Values decrease above 500 cm depth until the middle of Unit 3 at around 200 cm depth (~5.4 cal kBP), where values increase again until the top of the sediment core (Figure 11).

  Ca, TOC/TNatomic, and Fe are positioned together in the lower left quadrant. C (wt%), TOC (wt%), and TN (wt%) are positioned together in the lower right quadrant. TOC (wt%) and C (wt%) are very well correlated. Ti is positioned in the upper right quadrant. Zr and δ 13 C are very well positively correlated to PC 2.

Figure 12 .

 12 Figure 12. Ordination plot of the PCA for XRF variables (Ca, Sr, Fe, Ti, Zr), carbon variables (C (wt%), TOC (wt%), TN (wt%) TOC/TN), and δ 13 C for Unit 1, Unit 2, and Unit 3.

Figure 13 .

 13 Figure 13. Sediment core profile with a subset of six grain size distribution graphs. Unit 2 (Late Pleistocene-Holocene transition) Unit 2 (Late Pleistocene-Holocene transition)

  published in Limnology and Oceanography special edition: Biogeochemistry and ecology across arctic aquatic ecosystems in the face of change. doi: 10.1002/lno.

  GREENHOUSE GAS CONCENTRATIONS ........................................................... 5.3 DIFFUSIVE GREENHOUSE GAS FLUXES: COMPARISON ACROSS HIGH-LATITUDE REGIONS ........................ 6 CONCLUSIONS ....................................................................................................................................... 7 REFERENCES ......................................................................................................................................... 8 ACKNOWLEDGMENTS ........................................................................................................................

  photo-oxidation of dissolved organic carbon(Ward et al. 2017). Lateral transport of CH4

Figure 1 .

 1 Figure 1. Study area location and context. a) and b) The study area is located within the continuous permafrost zone of Eastern Siberia, about 120 km from the city of Yakutsk (c.) Sampling sites (lakes) were visited during four different seasons (Syrdakh village appears in the background).

  meters deep), and are still expanding downwards and laterally due to active layer deepening and thermokarst processes. Their waters can be notably rich in dissolved organic carbon (DOC), with concentrations typically ranging from ~ 100 to 300 mg L - 1 (Hughes-Allen et al. 2020).

Figure 2 .

 2 Figure 2. Distribution of lake types in a subset of the study area, including unconnected alas lakes (a; blue), recent thermokarst lakes (b; red), and connected alas lakes (c; purple).

Figure 3 .

 3 Figure 3. Map showing the complete distribution of sampled lakes for all four seasons by lake type. Unconnected alas lakes are indicated by a blue point, recent thermokarst lakes are indicated by a red point, and connected alas lakes are indicated by a purple point.

  done during the fall sampling season, water surface temperatures (ranging from 4.7 to 12.5°C depending on the lake;Hughes-Allen et al. 2020) indicate that some lakes may have been sampled before the water column had completely mixed to bottom. During the winter, lakes were covered by 0.7 to 1.2 meters of ice and 30 to 50 cm of snow (measured in March-April 2019). Shallow lakes, mostly unconnected alas lakes, were frozen completely to the bottom during the winter, whereas conditions were met for inverse stratification in deeper lakes (defined as the water bodies keeping liquid water during winter, i.e. of a maximal depth > 1.2 m). Sampling occurred late in winter (< 3 April), but well before the thaw season, which usually begins in May in this region. Many lakes (larger and/or deeper) were still covered in ice at the beginning of May when the site was visited for the spring sampling, but the ice cover melted progressively until surface water reached 4°C and spring overturn started. The surface temperature of most lakes was already above 4°C at the time of sampling in spring (ranging from 2.8°C to 17.5°C; Hughes-Allen et al. 2020), with the smaller unconnected alas lakes showing consistently higher surface temperatures (mean 10.9°C) than the connected alas lakes (5.4°C) or recent thermokarst lakes (6.8°C). Therefore, spring sampling likely occurred after the initiation of the thermal stratification. The warming of surface waters during the summer initiated the thermal and chemical stratification of the water column that may have lasted throughout the summer for deeper lakes.

Figure 5 .

 5 Figure 5. Winter limnological profiles (temperature, dissolved oxygen, conductivity and pH) in recent thermokarst lake #1 (a-b), connected alas lake #2 (c-d), and unconnected alas lake #3 (e-f).

Figure 6 .

 6 Figure 6. Summer limnological profiles (temperature, dissolved oxygen, conductivity and pH) in recent thermokarst lakes #1 (a-b), #8 (c-d) and #18 (e-f).

Figure 7 .

 7 Figure 7. Dissolved greenhouse gas concentrations in lakes from surface samples, expressed as departure from saturation, during fall 2018 (a-b), winter 2019 (c-d), spring 2019 (e-f) and summer 2019 (g-h). Scales for CO2 concentrations appear on the left, and for CH4 concentrations on the right. Number of lakes sampled for each lake type for each season is presented in Table 2.

  ........................................................................................................ 3.4.1. SOUTH STUDY SITE ........................................................................................................................... 3.4.2. CENTER STUDY SITE ......................................................................................................................... 3.5. TEMPERATURE AND PRECIPITATION ....................................................................................................... CHANGES IN TEMPERATURE AND PRECIPITATION SINCE 1900 ................................................................ 4.2. SPATIAL DISTRIBUTION OF LAKE TYPES .................................................................................................. 4.3. LAKE SURFACE AREA CHANGE: SOUTH STUDY SITE ................................................................................. 4.4. LAKE SURFACE AREA CHANGE: CENTER STUDY SITE................................................................................ 5 CONCLUSIONS ....................................................................................................................................... 6 REFERENCES .........................................................................................................................................

Figure 1 .

 1 Figure 1. Study area location and context. a) and b) The study area is located within the continuous permafrost zone of Eastern Siberia, about 120 km from the city of Yakutsk (c.) Sampling sites (lakes) were visited during four different seasons (Syrdakh village appears in the background). Adapted from Hughes-Allen et al 2020

  albedo areas surrounded directly by forest. They sometimes appear adjacent to roads or newly cleared agricultural fields. Compared to the other lake types, they have notably higher concentrations (several hundreds of mg L -1 ) of dissolved OC(Hughes-Allen et al., 2020).

Figure 2 .

 2 Figure 2. Distribution of lake types in a subset of the center study area, including unconnected alas lakes (clear blue), recent thermokarst lakes (red), and connected alas lakes (purple). Adapted from Hughes-Allen et al. 2020 3 Methods

Figure 3 .

 3 Figure 3. Image showing different lake boundaries. Black arrows are pointing to fuzzy boundaries. White arrows are pointing to sharp boundaries.

Figure 4 .

 4 Figure 4. Image showing half-moon shaped unconnected alas lakes (left) and circular unconnected alas lake (right).

Figure 5 .

 5 Figure 5. Image showing the three different lake types. A connected alas lake is circled in purple. Several unconnected alas lakes are circled in blue. Several recent thermokarst lakes are circled in red.3.3. General Workflow

Figure 6 .

 6 Figure 6. Graphical representation of the individual image sections. The black outline represents a single satellite image. The no fill squares surrounded by thin and dashed black lines represent the base sections. The gray sections represent the overlap images.

Figure 7 Figure 8 .

 78 Figure 7. Comparison of all 50 checkpoints to the corrected shapefile of the 2013 (July 14)test scene. False positive rate= a lake polygon was predicted where no lake exists. False negative rate =no lake polygon was predicted where a lake exists. False prediction rate = sum of false positive rate and false negative rate. ('wrt' = with respect to)3.3.4. EnsemblingTo improve upon the error rate of the initial model, an ensembling technique was

Figure 10 .

 10 Figure 10. Layout of the two study sites (center (blue) and south (black)). The Lena River runs north to south. The Aldan River runs west to east. Tiles: Google Map

Figure 12 .

 12 Figure 12. Center study area outlined in black overlying the 1980 (September 20) scene. The City of Borogontsy is indicated by the yellow star. Syrdakh Village is indicated by the blue star. The black arrow (south of the zone) indicates the CA lake which accounts for a substantial part of the total CA lake surface area for this scene.

4

  Results and Discussion 4.1. Changes in temperature and precipitation since 1900 Temperature records from the Yakutsk airport (62.0866° N, 129.7500° E) show an increasing frequency of years with above average annual temperatures during the early 21 st century and an overall trend of increasing temperature (Mann Kendall test: trend= increasing; p= 8.83x10 -6 ) (Figure 13). The mean annual air temperature (MAAT) in 1900 was -14.6º C and the MAAT in 2019 was -5.3º C. The 30-year average fluctuated slightly between -10º C and -8º C since 1900. There were a proportionally high number of years between 1910 and 1940, and between 2000 and 2019 with above average MAAT. The temperature records from Yakutsk airport and other studies (i.e., Czerniawska and Chlachula, 2020) indicate an annual temperature increase of 0.07º C since 1900. There are several consecutive years of above or below average precipitation. For example, years between 1941 and 1950 experienced below average precipitation and years between 2005 and 2008 experienced above average precipitation. There is a trend of increasing yearly precipitation in Central Yakutia (Mann Kendall test: trend= increasing; p= 8.15x10 -7 ).

Figure 13 .

 13 Figure 13. Deviation of average annual temperature from average annual temperature between 1951-1980.

Figure 14 .

 14 Figure 14. Deviation of the annual sum of precipitation from the average annual sum of precipitation between 1951-1980.

  CA lakes did transition to UCA lakes, sometimes as the result of persistent low precipitation reducing the flow volume of their connecting rivers and streams and isolating the lakes. The distribution of CA lakes is directly related to the presence of rivers and streams and local topography. The center study site had ~50 CA lakes, while the south study site had ~10.

Figure 15 .

 15 Figure 15. Distribution of lake types in the 2016 (September 11) image (center study site). UCA lakes are outlined in blue, CA lakes in purple, and RT lakes in red. The large CA lake in the upper right corner filled by the purple hashes is not included in the surface area change analysis. The City of Borogontsy is indicated by the yellow star. Syrdakh Village is indicated by the blue star.

Figure 16 .

 16 Figure 16. Histogram of lake surface area for the south study site for a) All lake type merged, b) UCA lake types, c) RT and d) CA lake types.

Figure 17 .

 17 Figure 17. Comparison between the July 12, 1989 image, September 25, 2005 image, and August 02, 2007 image. Two large unconnected alas lakes are outlined in blue (lakes a, b) and one connected alas lake is outlined in purple (lake c). An inflow to lake c is visible in the bottom right corner of each image.

Figure 18 .

 18 Figure 18. Comparison between the July 12, 1989 image, September 25, 2005 image, August 02, 2007 image and September 11, 2011 images. UCA lakes are outlined in blue. RT lakes are outlined in red. In the 2007 image, a newly formed RT lake is indicated by the black arrow. In the 2011 image a RT lake which disappears is indicated by the white arrow.

  no lake at all (Figure20). 1967 RT and CA lake surface area values are closer to those of the other four scenes. UCA lake surface area values are, however, much lower than the other four scenes (Figure19). The 1967 scene follows five years of below average precipitation (Figure14), possibly contributing to low UCA lake levels. By the 2010 (September 23) scene, many of these alas depressions are occupied by substantial lakes (Figure20). Between 1980 and 2000, yearly precipitation fluctuates regularly above and below the 30-year average. The 2010 scene follows several years of above average precipitation, which is likely reflected in the high lake surface area values (particularly UCA lakes).

Figure 20

 20 Figure 20. a) 1967 (September 20) scene (Un-segmented). The yellow circles are indicating the areas where there are numerous completely dry alas basins. b) 2010 (September 23) scene. Lake basins within the yellow ovals are now occupied by substantial lakes.

.................................................................................................................................... 2 STUDY SITE .............................................................................................................................................. 3 METHODS ..................................................................................................................................................

  FIELD SAMPLING ...................................................................................................................................... 3.2 SEDIMENT CORE SUBSAMPLING AND DATING ........................................................................................... 3.3 X-RAY FLUORESCENCE (XRF) ANALYSIS ................................................................................................. 3.4 GRAIN SIZE ANALYSIS ............................................................................................................................... 3.5 DRY BULK DENSITY, SEDIMENTATION AND ORGANIC CARBON ACCUMULATION CALCULATIONS ............. 3.6 BIOGEOCHEMICAL ANALYSIS.................................................................................................................... 3.7 STATISTICAL ANALYSIS ............................................................................................................................

		published in Frontiers in Earth Science special edition: Yedoma permafrost
	landscapes as past archives, present and future change areas. doi: 10.3389/feart.2021.710257
	1	INTRODUCTION ..3.1

4 RESULTS ....................................................................................................................................................

  CHRONOLOGY AND SEDIMENTATION RATES ............................................................................................. 4.2 GENERAL STRATIGRAPHY ........................................................................................................................ 4.3 GRAIN-SIZE DISTRIBUTION........................................................................................................................ 4.4 BIOGEOCHEMISTRY .................................................................................................................................. 4.5 INORGANIC ELEMENTAL COMPOSITION .....................................................................................................

	4.1
	4.6	PCA ANALYSIS.........................................................................................................................................
	5	DISCUSSION ..

............................................................................................................................................

	5.1	MULTIPROXY-INFERRED PALEOLIMNOLOGICAL HISTORY.........................................................................
	5.2	LAKE MALAYA CHABYDA CARBON ACCUMULATION RATES ....................................................................
	5.3	CONNECTIONS BETWEEN THE LAKE ENVIRONMENT, PERMAFROST DYNAMICS, AND CLIMATIC CONDITIONS
		88

6 CONCLUSIONS ......................................................................................................................................... 7 REFERENCES ........................................................................................................................................... 8 ACKNOWLEDGEMENTS .....................................................................................................................

Table 1

 1

	. Radiocarbon ( 14 C) dated samples from Lake Malaya Chabyda sediment core.
	ABA: acid-base-acid.
	Bulk: sediment bulk sample.
	Wood: terrestrial wooden remain (unknown genus) from catchment vegetation.
			Composite			Methods and
	Lab ID	sample label	depth (cm)	14C age (yr)	14C Error (yr)	material
		PG2201-2_a200-
	AWI-2601.1.1	300_80-81	66.5	1802	70	ABA, wood
		PG2201-2_b280-
	AWI-2786.1.1	380_75-76	141.5	3564	38	ABA, wood
		PG2201-2_d440-
	AWI-2602.1.2	540_59.5-60	285.5	7691	32	A, bulk
		PG2201-2_f600-
	AWI-2606.1.2	700_10-10.5	396.25	8460	37	A, bulk
		PG2201-2_f600-
	AWI-2605.1.1	700_77-78	463.5	8548	106	ABA, wood
		PG2201-2_h780-
	AWI-2609.1.2	880_9.5-10	573.25	10218	39	A, bulk
		PG2201-2_h780-
	AWI-2608.1.2	880_97-97.5	661.25	12184	38	A, bulk

  LAKE TYPES ............................................................................................................................................ PHYSICOCHEMICAL CHARACTERISTICS OF LAKE WATER ........................................................................

	.....
	2.1
	3.1
	3.2 DISSOLVED 4.1 SEASONAL

3 METHODS ................................................................................................................................................

GREENHOUSE GAS MEASUREMENTS ..................................................................................... 3.3 STATISTICAL ANALYSIS .......................................................................................................................... 4 RESULTS .................................................................................................................................................. CONDITIONS .......................................................................................................................... 4.2 PHYSICOCHEMICAL CHARACTERIZATION OF LAKE WATER ..................................................................... 4.2.1 Broad trends and seasonal averages. .......................................................................................... 4.2.2 Seasonal Profiles. ......................................................................................................................... 4.2.3 Dissolved greenhouse gas concentrations ................................................................................... 4.2.4 Diffusive greenhouse gas fluxes ...................................................................................................

Table 2 .

 2 Depths and surface areas of sampled lakes in Central Yakutia (Eastern Siberia). Limnological point samples indicated by a star. Limnological profile samples indicated by anarrow. Greenhouse gas samples indicated by +. ND = no data available.

			Lake 40 62.5535 °N 130.9175 °E	Tyungyulyu	ND	0.17	*	ND	* +	* +
			Lake 41 62.5748 °N 130.8580 °E	Tyungyulyu	ND	0.37	*	ND	* +	* +
			Lake 42 62.5746 °N 130.8547 °E	Tyungyulyu	ND	0.18	*	ND	* +	* +
			Lake 43 62.5647 °N 130.9033 °E	Tyungyulyu	3.6	0.31	*	*  +	* +	* +
			Mean		2.6	0.49
	Connected	Alas	Lake 2 62.5598 °N 130.9655 °E Lake 11 62.5909 °N 131.1871 °E Lake 36 62.5631 °N 130.8998 °E Mean	Tyungyulyu Abalakh Tyungyulyu	8.8 4.9 3.6 5.7	64.64 *  + *  + 59.35 * + *  + 179.46 * + *  + 101.15	* + * + * +	*  + * + * +
			Lake 3 62.5507 °N 130.9236 °E	Tyungyulyu	2.0	1.68	* + *  +	* +	* +
			Lake 5 62.5615 °N 131.0110 °E	Tyungyulyu	ND	4.54	ND	ND	* +	* +
			Lake 10 62.5507 °N 130.9226 °E	Tyungyulyu	ND	0.88	* +	ND	* +	* +
			Lake 12 62.5844 °N 131.1910 °E	Abalakh	1.4	40.72 * + *  +	* +	* +
			Lake 19 62.5585 °N 130.9510 °E	Tyungyulyu	ND	2.22	ND	ND	* +	* +
	Alas	Lake 20 62.5702 °N 130.917 °E Lake 44 62.5660 °N 130.9012 °E	Tyungyulyu Tyungyulyu	ND 3.1	5.35 1.93	ND *	ND *  +	* + * +	* + * +
	Unconnected	Lake 45 62.5753 °N 130.9142 °E Lake 46 62.5625 °N 130.9269 °E Lake 47 62.5467 °N 130.9660 °E Lake 48 62.5663 °N 130.9401 °E	Tyungyulyu Tyungyulyu Tyungyulyu Tyungyulyu	ND ND ND ND	0.54 1.81 3.32 3.07	* ND ND ND	ND ND ND ND	* + * + * + * +	* + ND * + ND
			Lake 49 62.5478 °N 130.9587 °E	Tyungyulyu	ND	2.44	ND	ND	*	* +
			Lake 50 62.5618 °N 130.9859 °E	Tyungyulyu	ND	1.35	ND	ND	*	*
			Lake 100 62.5864 °N 131.0380 °E	Tyungyulyu	ND	4.54	ND	ND	ND	* +
			Lake 101 62.5510 °N 130.9347 °E	Tyungyulyu	ND	4.55	ND	ND	* +	* +
			Mean		2.2	4.97

  Inverse stratification was observed in winter, with

	(μS cm -1 )	Fall Unconnected alas lake Recent thermokarst lake Connected alas lake	a	Winter	b	Spring	c	Summer	d
	Conductivity
	Dissolved oxygen (%)		e		f		g		h	150
			i		j		k		l
	p H

Table 5 .

 5 Diffusive fluxes (in mmol m -2 d -1 ) of CO2 and CH4 over four sampling seasons in 2018-2019 in lakes of different developmental stages (age) in Central Yakutia(Eastern Siberia). N is the number of lakes sampled in each category.

		CO2 Flux (mmol m -2 d -1 ) CH4 Flux (mmol m -2 d -1 )
	Type	N Min Med Mean Max Min Med Mean Max
	FALL 2018
	Recent Thermokarst	6 -5.7 19.8 17.7 39.8	0.1 0.3	3.9	21.7
	Connected Alas	4 -4.0 1.5	6.7	27.6	0.2 0.4	0.6	1.4
	Unconnected Alas	4 -10.6 -7.3 -6.4	-0.7	2.6 12.0 12.4 22.9
	All lake types	14 -10.6 0.4	7.7	39.8	0.1 0.6	5.4	22.9
	SPRING 2019
	Recent Thermokarst	23 -9.4 -3.3 -1.5 21.7	0.2 6.6	12.7 74.9
	Connected Alas	5 -6.9 -0.1	0.3	8.6	4.0 4.7	8.5	21.0
	Unconnected Alas	17 -13.6 -4.9 -5.2	-0.2	0.2 11.6 15.0 56.8
	All lake types	45 -13.6 -3.6 -2.7 21.7	0.2 8.2	13.1 74.9
	SUMMER 2019
	Recent Thermokarst	16 -3.3 20.8 51.4 355.4	0.3 0.6	36.2 566.4
	Connected Alas	4 -3.7 15.5 13.6 27.0	0.3 1.1	1.5	3.2
	Unconnected Alas	12 -7.4 -3.2	8.2	98.1	0.2 3.8	3.9	9.5
	All lake types	32 -7.4 14.5 30.5 355.4	0.2 1.1	19.8 566.4
	ICE-FREE SEASONS
	(2018-19)
	Recent Thermokarst	45 -9.4 1.8	19.9 355.4	0.1 1.1	19.9 566.4

Table 6 .

 6 The CO2 to CH4 molar ratios for all lakes sampled in Central Yakutia along the seasons in 2018-2019. N is the number of lakes sampled in each category.

	CO2/CH4 ratio

  3.1. IMAGE DATA........................................................................................................................................... 3.2. DEFINING LAKE BOUNDARIES AND LAKE TYPES ..................................................................................... 3.3. GENERAL WORKFLOW ............................................................................................................................ 3.3.1.MACHINE LEARNING MODEL ............................................................................................................... 3.3.2.FINE TUNING AND TRAINING ................................................................................................................ 3.3.3.ACCURACY ASSESSMENT OF INITIAL MODEL ....................................................................................... 3.3.4. ENSEMBLING ...................................................................................................................................... 3.3.5. COMPARISON OF TOTAL SURFACE AREA FOR PREDICTION AND CORRECTED SHAPEFILES ................... 3.4. SURFACE AREA CHANGE ANALYSIS

  . Nitze et al., 2017 found that lake area in a Central Yakutian study site increased by nearly 50% between 1999 and 2014 based on Landsat analysis. Boike et al., 2016 recorded an average increase of 17.9% in the total area covered by lakes between 2002 and 2009 in the central part of the Lena River catchment in the Yakutian region of Siberia.

Table 1 .

 1 Description of satellite parameters.

	Platform Ground sample distance Swath width	Bands		Operation dates
	SPOT 1 Panchromatic: 10 m	60 km	Pan	0.51-0.73 µm	02/1986-11/2003
	Multispectral: 20 m		Green	0.50-0.59 µm
			Red	0.61-0.68 µm
			NIR	0.79-0.89 µm
	SPOT 2 Panchromatic: 10 m	60 km	Pan	0.51-0.73 µm	01/1990-07/2009
	Multispectral: 20 m		Green	0.50-0.59 µm
			Red	0.61-0.68 µm
			NIR	0.79-0.89 µm
	SPOT 3 Panchromatic: 10 m	60 km	Pan	0.51-0.73 µm	09/1993-11/1996
	Multispectral: 20 m		Green	0.50-0.59 µm
			Red	0.61-0.68 µm
			NIR	0.78-0.89 µm
	SPOT 4 Panchromatic: 10 m	60 km	M	0.61-0.68 µm	03/1998-06/2013
	Multispectral: 20 m		Green	0.50-0.59 µm
			Red	0.61-0.68 µm
			NIR	0.79-0.89 µm
			SWIR	1.58-1.75 µm
	SPOT 5 Panchromatic: 5 m	60 km	Pan	0.48-0.71 µm	05/2002-03/2015
	Multispectral: 10 m		Green	0.50-0.59 µm
			Red	0.61-0.68 µm
			NIR	1.58-1.75 µm
			SWIR	1.58-1.75 µm
	SPOT 6 Panchromatic: 2 m	10-60 km	Pan	0.45-0.52 µm	09/2012-present
	Multispectral: 8 m		Green	0.53-0.60 µm
			Red	0.62-0.69 µm
			NIR	0.76-0.89 µm
	SPOT 7 Panchromatic: 2 m	10-60 km	Pan	0.45-0.52 µm	06/2014-present
	Multispectral: 8 m		Green	0.53-0.60 µm
			Red	0.62-0.69 µm
			NIR	0.76-0.89 µm

Table 3 .

 3 Scenes used in each study site.

Table 4 .

 4 Maximum, median, and mean surface area by lake type for both study sites. Minimum area is not included as this figure is biased by very small (i.e., < 1 pixel x <1 pixel) erroneous lake polygons that are predicted by the neural network.

	Lake Type Max area (ha) Median area (ha) Mean area (ha)
	UCA	347.4	0.9	3.0
	RT	13.6	0.3	0.7
	CA	1763.4	3.3	60.1

Acknowledgements

The authors would like to thank the Melnikov Permafrost Institute (Yakutsk) for their collaboration and fieldwork assistance. We are also grateful to A. Rudolf (AWI) for his calculation of lake catchment area. D. Scheidemann (AWI) and J. Lindemann (AWI) provided laboratory support and analyses. Many thanks to G. Fougeron for his invaluable and patient help with statistical analysis and python programming. D. Subetto was supported by Ministry of Education of the Russian Federation (project No. FSZN-2020-16). L. Hughes-Allen and F. Bouchard are supported by the Agence nationale de la recherche (ANR) through the Make Our Planet Great Again initiative (Programme d'investissements d'avenirproject no. ANR-17-MPGA-0014). Acknowledgments

We sincerely thank the gracious people of Syrdakh village for welcoming us into their village and for sharing with us their traditional environmental knowledge. We are also grateful to staff at the Melnikov Permafrost Institute (Yakutsk), especially P. Konstantinov, who were invaluable in their fieldwork assistance. A. Noret (GEOPS) and S. Duval (INRS-ETE) provided laboratory support and analyses. M. Alexis, U. Christaki, G. David, L. Gandois, C. Grenier, L. Jardillier, M. Pessel, E. Pohl and A. Saintenoy provided help for field sampling. We finally thank the Associate Editor (S. Tank) and two anonymous reviewers for their insightful comments, which greatly helped to improve this manuscript. Funding was provided by ANR-MOPGA (ANR-17-MPGA-0014) through the 'Programme d'investissements d'avenir (PIA), Institute Pierre Simon Laplace (IPSL) and Université Paris Saclay. The authors declare that they have no conflict of interest.

Data Availability

All data used in this study are available at: https://doi.pangaea.de/10.1594/PANGAEA.933411.

Chapter 1

Alas lake near Syrdakh Village, Sakha Republic, Russia

Chapter 2

Sunset over a connected alas lake near Syrdakh Village, Sakha Republic, Russia.

Chapter 3

Field work on a recent thermokarst lake near Syrdakh Village, Sakha Republic, Russia

Study site

Central Yakutia is characterized by extreme subarctic continental climate with long, cold and dry winters (January mean temperature around -40°C) and warm summers (July mean temperature around +20°C), resulting in notably strong seasonal variability. The winter season (defined by the presence of an ice cover) usually lasts from late September until early May. The low annual precipitation (190-230 mm) is generally restricted to summer. Average snow depth for winter months (January to April) ranges from 24 cm in January to a maximum of 30 cm in March, and then decreasing to 10 cm at the end of April (1980-2020 recorded values from Yakutsk weather station). Yearly evaporation rates exceed total precipitation in this region (Fedorov et al. 2014b). Central Yakutia is no exception to the observed trend of high latitude regions warming disproportionately quicker than lower latitudes. Between 1996 and 2016, the mean annual air temperature of Central Yakutia has increased by 0.5-0.6°C per decade (Gorokhov and Fedorov 2018).

The study site (62.554 °N; 130.982 °E) is located near the rural village of Syrdakh, which lies on the lowland plain between the Lena River to the west and the Aldan River to the east, approximately 130 km north-east of Yakutsk (Figure 1). The region is predominantly covered by late Pleistocene sediments, including silty clays and sandy silts of fluvial, lacustrine or aeolian origin (Ivanov 1984). The activity of the Lena and Aldan rivers, as well as their smaller tributaries, has resulted in the formation of numerous fluvial terraces during the Pleistocene (Soloviev 1959). Lakes later developed within the late Pleistoceneaged fluvial terraces: the Tyungyulyu terrace, which covers the western section of the study area, 50-200 m above sea level (asl), dated 14-22 kyr BP, and the higher Abalakh terrace in the eastern sector of the study area, 200-280 m asl, dated 45-56 kyr BP (Ivanov, 1984;Soloviev 1959[START_REF] Soloviev | Thermokarst phenomena and landforms due to frost heaving in Central Yakutia[END_REF]. This region is characterized as a middle taiga landscape regime (Fedorov et al. 2014a) and is dominated by larch, pine and birch forests (Ulrich et al. 2017a).

Chapter 4

Recent thermokarst lake near Sydrakh Village (Sakha Republic, Russia)

(September 11, 2016) (see Figure 15). The original image was split into 160 smaller (1024 x 1024 pixel) images. Ten of these 160 smaller images were chosen and manually digitized.

These digitized lake polygons were then used as a preliminary fine tuning of the model. The A training dataset consisting of a total of 8,286 training samples was created from the annotated images. Each training sample consisted of a 1-megapixel image. Standard data augmentation practices were followed (random rotation, scaling, and brightness adjustment).

The Adam optimization algorithm was used to train the model (Kingma and Ba, 2014). The training of the neural network spanned 50 epochs (for each of the 8,286 small images) and checkpoints were saved at every epoch. Each checkpoint is a set of neural network parameters ('weights').

Accuracy assessment of initial model

The 50 checkpoints were compared to the corrected shapefile and the false positive rate (a lake polygon was predicted where no lake exists) and false negative rate (no lake

Comparison of total surface area for prediction and corrected shapefiles

Three polygon shapefiles were manually corrected to assess the differences between the predicted total lake surface area of each of the three saved weights to the corrected total lake surface area. The manual corrections utilized the lake polygons generated by the neural network, but rectified where needed. Manual digitizing, as was done in the first stage of finetuning the neural network, is extremely time consuming and tedious. Manual digitizing one 60 x 60 km 2 scene without the aid of the polygons generated by the neural network would take at least one full week (usually between onethree weeks) of motivated work.

Correction of the polygons generated by the neural network for one scene could be completed in two-three hours or less. Total lake surface areas were compared between the corrected shapefile and each of the three output shapefiles from the saved weights (Figure 9). Version 1 consistently overpredicted total lake surface area (2.7-8 % false prediction rate). Version 2 underpredicted the total surface area for two out of the three corrected scenes (-9-4% false prediction rate). Version 3 underpredicted two out of the three corrected scenes (0.5--12).

The false prediction rate for the 2011 (September 11) scene version 3 was very close to zero.

based on five scenes from 1989-2011. The scenes are not evenly distributed in time and there is a particularly large gap between the 1989 scene and the next scene in 2005 (See Table 3).

Center study site

Lake surface area was compared between five scenes spanning 1967-2016 in a 1150 km 2 study area (Figure 12). The scenes are not evenly distributed in time and there is a large gap

Chapter 5

Lakes near Sydrakh Village (Sakha Republic, Russia)

The contribution of permafrost landscapes to the global carbon budget as permafrost thaws in a warming climate is significant, but potentially underestimated. The challenging logistics of obtaining year-round GHG measurements from thermokarst lakes makes these records exceedingly rare. However, the high seasonal variability that we documented in this study underscores the importance of these records to enable accurate accounting of the contribution of permafrost landscapes to the global carbon budget.

It is very likely that GHG concentrations and fluxes vary annually. Yearly fluctuations in temperature and precipitation will affect the mixing regimes of the lakes and subsequently the rates of microbial activity and methanogenesis. Obviously, additional yearly records are necessary to corroborate our results and identify longer term trends in GHG emissions from thermokarst lakes in Central Yakutia. Unfortunately, in addition to the humanitarian crisis and the global upheaval that the current Russian invasion of Ukraine has caused, this war has eliminated the possibly of foreign researchers continuing this work in Yakutia. There is the possibility to shift this research to other Yedoma sites in Alaska and the Yukon. Permafrost is abundant in Alaska. 80 percent of the state is underlain by permafrost and 30 percent of this area is made up of continuous permafrost and thermokarst landforms are abundant (Jorgenson et al., 2008). Nearly 50 percent of Canada is underlain by permafrost and the continuous zone makes up a substantial portion of this area (Fritz et al., 2016). Hopefully new study sites can be established for continued research into temporal and spatial differences in GHG emissions from thermokarst lakes.

The final axis of this thesis (Chapter 4) uses machine learning techniques to automate lake polygon identification in remote sensing images of Central Yakutia. The goal of this project is to identify trends in the spatial distribution of each of the three lake types, quantify recent (since 1960) changes in lake surface area, and understand how different lake types are Permafrost landscapes are sensitive to changes in temperature and precipitation and the effects of recent human caused climate change are already evident. The interconnectedness of climate, permafrost, and the global carbon cycle make this and similar research integral to our understanding of these complex systems. I remain optimistic that we will be able to change our current climate trajectory through significant and directed policy changes and innovative developments in infrastructure and energy sectors.