Simulating the thermal regime and thaw processes of ice-rich permafrost ground with the land-surface model CryoGrid 3

. Thawing of permafrost in a warming climate is governed by a complex interplay of different processes of which only conductive heat transfer is taken into account in most model studies. However, observations in many permafrost landscapes demonstrate that lateral and vertical movement of water can have a pronounced inﬂuence on the thaw trajectories, creating distinct landforms, such as thermokarst ponds and lakes, even in areas where per-mafrost is otherwise thermally stable. Novel process param-eterizations are required to include such phenomena in future projections of permafrost thaw and subsequent climatic-triggered feedbacks. In this study, we present a new land-surface scheme designed for permafrost applications, Cryo-Grid 3, which constitutes a ﬂexible platform to explore new parameterizations for a range of permafrost processes. We document the model physics and employed parameterizations for the basis module CryoGrid 3, and compare model results with in situ observations of surface energy balance, surface temperatures, and ground thermal regime from the Samoylov permafrost observatory in NE Siberia. The comparison suggests that CryoGrid 3 can not only model the evolution of the ground thermal regime in the last decade, but also consistently reproduce the chain of energy transfer processes from the atmosphere to the ground. In addition, we demonstrate a simple 1-D parameterization for thaw processes in permafrost areas rich in ground ice, which can phe-nomenologically reproduce both formation of thermokarst ponds and subsidence of the ground following thawing of ice-rich subsurface layers. Long-term simulation from 1901 to 2100 driven by reanalysis data and climate model output demonstrate that the hydrological regime can both accelerate and delay permafrost thawing. If meltwater from thawed ice-rich layers can drain, the ground subsides, as well as the formation of a talik, are delayed. If the meltwa-ter pools at the surface, a pond is formed that enhances heat transfer in the ground and leads to the formation of a talik. The model results suggest that the trajectories of future per-mafrost thaw are strongly inﬂuenced by the cryostratigraphy, as determined by the late Quaternary history of a site.


Introduction
In the past decade, the fate of permafrost areas in a changing climate has rapidly moved into the focus of climate researchers.This has led to significant improvements in the representation of permafrost and frozen ground processes in the land-surface schemes employed in Earth system modeling schemes (e.g., Lawrence et al., 2008).These efforts have, to a large part, been motivated by the potential release of the greenhouse gases (GHGs) carbon dioxide and methane from thawing carbon-rich permafrost soils (Walter et al., 2006;Zimov et al., 2006).However, including this possibly significant climatic feedback (Schaefer et al., 2011) in future climate predictions requires a correct description of both the climate forcing of the global permafrost area and the physical processes governing the transformation of the local permafrost systems.However, there is a striking gap between Published by Copernicus Publications on behalf of the European Geosciences Union.
the temporal and spatial scales on which the forcing and the impacts occur (e.g., Lee et al., 2014).This implies that the actual mobilization of carbon is not directly accessible through process modeling on large grid cells, as it is generally done in atmospheric modeling.As an example, thermokarst and active layer detachments can trigger strong mass wasting of ice-and organic-rich soils in time periods of days to years at spatial scales of meters to hundreds of meters.Carbon fluxes from hotspots can significantly exceed the background fluxes (e.g., Walter et al., 2006;Abnizova et al., 2012;Langer et al., 2015).Integrating such highly localized and non-linear processes in large-scale grid-based models requires the development of up-and downscaling strategies, which, despite the emerging stronger collaboration between field scientists and modelers, is still in its infancy.
Many permafrost areas are dominated by ice-rich ground (Brown et al., 1997).Due to ice wedges, segregation ice, or ice lenses, the volumetric ice content can exceed the natural pore volume of the soil (i.e., excess ground ice), which in case of thawing will lead to a settling and consolidation of the ground material.The released meltwater is either lost as surface or subsurface runoff leading to an effective subsidence of the surface, or it pools up, thus forming a thermokarst pond or lake.In either case, the ground and surface properties are modified by thawing of excess ground ice, which can in turn influence the trajectories of subsequent permafrost thaw.The vast lowland permafrost areas in Siberia, Alaska, and Canada dominated by thermokarst lakes and depressions are ample evidence that thawing of excess ground ice has played a major role in the landscape development throughout the Holocene.To predict future landscape development in a warming climate that is intimately linked to the carbon cycle of these areas, new and robust model approaches are needed, which can potentially be implemented in the landsurface schemes of Earth system models (ESMs).
At the moment, there are two main directions in permafrost modeling on larger spatial scales: dedicated ground thermal models focus on the physical variables characterizing the state of the permafrost, and have been successfully applied to map permafrost distribution (e.g., Gisnås et al., 2013;Zhang, 2013;Zhang et al., 2013;Westermann et al., 2013;Fiddes et al., 2015) and derive the ground temperature evolution for past and future climate conditions (e.g., Zhang et al., 2008;Jafarov et al., 2012).They are characterized by vertical modeling domains of one hundred to several hundred meters with a vertical resolution of centimeters within the active layer.As a consequence, they provide a sufficient accuracy to resolve the annual dynamics of active layer thawing and refreezing, as well as the evolution of ground temperatures in deeper layers.For validation purposes, they can thus be directly compared to field measurements, which can usually be reproduced to within 0.2 m (active layer thickness) or 1.5 K (ground temperature), depending on the type of input data (e.g., Langer et al., 2013;Westermann et al., 2013;Jafarov et al., 2012).
On the other hand, ESM land-surface schemes are directed towards adequately providing the lower boundary for the atmospheric transfer scheme, i.e., the surface fluxes of radiation, momentum, and sensible and latent heat.There is a huge spread in the modeled present permafrost area in the CMIP5 results, with values between 1.4 (factor 10 underestimation) and 27.3 (factor 2 overestimation) million km 2 (Koven et al., 2013).While latent heat effects due to soil freezing and thawing have been included in a number of ESMs (e.g., Gouttevin et al., 2012), the conductive heat transfer in the ground remains poorly represented due to the limited number of soil grid cells and insufficient depth of the soil column.Therefore, projections on timescales of many decades over which permafrost thaw will occur can be strongly biased, which may at least partly explain the large inter-model spread in the sensitivity of permafrost extent, ground temperatures, and active layer depth towards climate change (Koven et al., 2013).
In this study, we describe a new permafrost modeling scheme, CryoGrid 3, which aims to narrow the gap between traditional transient permafrost models and landsurface schemes employed in Earth system modeling.Cryo-Grid 3 adds land-atmosphere coupling to the subsurface heat transfer model CryoGrid 2 (Westermann et al., 2013), which has been employed to simulate ground temperatures in Norway (Westermann et al., 2013) and NE Greenland (Westermann et al., 2015).In this way, it is possible to run CryoGrid 3 with the same forcing data as ESM land-surface schemes and compare results of present and future simulations of the ground thermal state.Moreover, the numerical structure is designed to be as simple as possible, which eases the requirements for adding and investigating parameterizations for permafrost processes.Following a similar philosophy as Essery (2015), CryoGrid 3 is dedicated to active users with the programming skills to modify and extend the code, rather than being a black box tool with a predefined and thus limited set of options.To illustrate the potential of the model scheme, we implement simple parameterizations for ground subsidence and initial thermokarst formation and demonstrate the effect on future simulations of the ground thermal regime for a permafrost site in NE Siberia.

CryoGrid 3 -model description
CryoGrid 3 builds upon the subsurface thermal model Cryo-Grid 2 (Westermann et al., 2013), supplemented by the calculation of the surface energy balance as upper boundary condition, as well as modifications of the snow scheme.As in CryoGrid 2, the energy transfer in the ground is governed by heat conduction and freeze-thaw processes, while the soil water balance is not explicitly accounted for, i.e., the sums of soil water and ice contents are constant in time for each soil layer.An overview of the processes implemented in the different models is given in Table 1.In the following, we pro-  vide a detailed description of the CryoGrid 3 model physics and the numerical scheme.

Driving variables
CryoGrid 3 is driven by time series of (i) air temperature T air [

The surface energy balance
The energy input to the uppermost grid cell is derived from the surface energy balance, i.e., the fluxes of shortwave radiation (S in , S out ) and long-wave radiation (L in , L out ), as well as the sensible, latent, and ground heat fluxes, Q h , Q e and Q g (see Sect. 2.4 for the latter), respectively.In the following, superscript numbers denote the number of a grid cell (increasing with depth, with grid cell 1 located at the surface).The equations follow Foken (2008), with variables and constants explained in Appendix A. The energy input to the uppermost grid is computed as While incoming shortwave and long-wave radiation are directly provided (Sect.2.1), the other fluxes are parameterized based on the driving data sets.With the albedo α, the outgoing shortwave is calculated as while the outgoing long-wave radiation L out is derived from Kirchhoff's and Stefan-Boltzmann law as The turbulent fluxes of sensible and latent heat, Q h and Q e , are parameterized in terms of gradients of air temperature and absolute humidity (see Appendix B) between the air at height h aboveground and at the surface following Monin-Obukhov similarity theory (Monin and Obukhov, 1954): The resistance terms r constitute the aerodynamic resistance, r a , characterizing the strength of the turbulent exchange, and the surface resistance against evapotranspiration, r s .The latter is an empirical parameter that can be adjusted to account for the fact that the water vapor pressure above soil surfaces is lower than the saturation vapor pressure above a water surface (Eq.B3) for non-saturated surface soils.
The aerodynamic resistances r a for sensible heat H and latent heat W , account for the atmospheric stability by including the integrated atmospheric stability functions H,W , which describe deviations from the logarithmic profile of the neutral atmospheric stratification (Appendix C).They are inversely proportional to the friction velocity which depends on the wind speed u and the integrated stability function for momentum, M (Appendix C).The variable L * is the Obukhov length, which in itself is a function of the turbulent fluxes Q h and Q e as well as the friction velocity u * : The surface emissivity ε, the roughness length z 0 and the surface resistance against evapotranspiration r s are assigned different values, depending on whether the surface is snowfree or covered by snow.The same applies to the albedo α for which a constant value is used for snow-free surfaces, while for snow-covered surfaces the albedo changes with the age of the surface snow (Sect.2.3).

Surface energy balance of snow
For snow-covered surfaces, the governing equations for the surface energy balance remain widely unchanged, with a few notable exceptions: first, penetration of shortwave radiation in the snow is accounted for by assuming exponential damping with depth (e.g., O'Neill and Gray, 1972).The shortwave radiation flux at depth d below the snow surface is given by with a bulk value for the entire spectral range for the SWradiation extinction coefficient β SW [m −1 ] assumed.At the bottom of the snow layer, the remaining shortwave radiation flux is assigned as energy input to the uppermost soil grid cell.
The albedo of the snow is calculated as a function of the difference between time t and the time t SF when the last significant snow fall occurred, following the formulations used in the ERA-Interim reanalysis (Dee et al., 2011;ECMWF, 2006).The time of the last significant snow fall t SF is set to t if the snowfall in the period [t − 1 d, t] exceeds a snow water equivalent of SWE SF .Furthermore, the snow albedo is set to the albedo of fresh snow, α max , from which it decreases with time towards the albedo of old snow, α min .For non-melting conditions, a linear decrease is assumed, while an exponential decrease is assumed for melting snow, The values for the time constants τ α,f and τ α,m [s −1 ] (Table 4) are selected according to ECMWF (2006) and correspond to an albedo change of 0.1 in 12.5 days for nonmelting conditions and an e-folding time of about 4 days for melting conditions.According to Eq. ( 1), positive temperatures can occur for snow grid cells, which will trigger melting of the snow.In this case, the temperature of the concerned grid cell is set to T = 0 • C, and the energy input ∂E i /∂t is partitioned into energy leading to a temperature change and energy leading to snowmelt, as detailed in Sect.2.5.For the uppermost grid cell, all components of the surface energy balance are calculated according to Eqs. ( 3)-( 8) with a fixed temperature

Subsurface heat transfer
The subsurface thermal scheme of CryoGrid 3 is based on conductive heat transfer as given by Fourier's law, where the generation and consumption of latent heat due to the phase change of soil water is taken into account through an effective volumetric heat capacity c eff (z, T ) [Jm −3 K −1 ] featuring a sharp peak in the freezing range of the soil water (Westermann et al., 2013).Both the effective heat capacity and the thermal conductivity k(z, T ) [Wm −1 K −1 ] of the soil are functions of the volumetric fractions of the soil constituents water, ice, air, mineral, and organic (Cosenza et al., 2003;Westermann et al., 2013).At the lower boundary of the soil domain, a geothermal heat flux F lb is assigned, which is usually assumed constant over time.
Vertical or lateral movement of water is not considered in the here described basis version of CryoGrid 3 so that the sum of the volumetric contents of ice and water remains constant (for simplicity, the densities of ice and water are assumed equal).The phase change of the soil water is determined by the freezing characteristic θ w (z, T ) (see Dall'Amico et al., 2011) depending on the soil type (three classes, sand, silt and clay), which can change with depth.The employed parameterizations for c eff (z, T ), k(z, T ), and θ w (z, T ) are identical to the ones employed in CryoGrid 2, documented in detail in Westermann et al. (2013).The subsurface thermal properties are specified as a stratigraphy of the volumetric contents of soil mineral and organic material, air and the sum of ice and water (as in Tables 2 and 3).

Energy transfer and mass balance of the snowpack
In the current version of CryoGrid 3, a constant density ρ snow is assumed for the entire snowpack, with the exception of density modifications through refreezing rain and meltwater (see below).A change in the snow water equivalent SWE due to snowfall, sublimation, or melt (Sect.2.3) leads to a change in snow depth d snow as To account for the buildup and disappearance of the snow cover, the position of the upper boundary is allowed to change dynamically by adding or removing grid cells.The initial temperature of a newly added grid cell is set to the air temperature T air .
As in soil, a thermal model based on heat conduction is employed for the snowpack, in conjunction with a 1-D coldhydrology scheme for percolation of melt and rain water in snow, following Westermann et al. (2011).However, only a limited number of natural processes occurring in the snowpack are accounted for, and CryoGrid 3 is therefore not comparable to sophisticated snow schemes, such as SNOWPACK   (Bartelt and Lehning, 2002;Lehning et al., 2002a, b) or CROCUS (Vionnet et al., 2012).In contrast to the soil domain, where latent heat effects are accounted for by the effective heat capacity, two state variables are required to describe the energy content of a snow grid cell, its temperature T and the volumetric water content θ w .Similar to the soil domain, the principal means of energy transfer within the snowpack is conductive heat transfer.We define two auxiliary variables T * and θ * w , which are equal to T and θ w at the beginning of each time step (see Sect. 2.6).We subsequently update T * according to the heat conduction equation with the snow heat capacity c snow and thermal conductivity k snow being functions of the snow density ρ snow as and (Yen, 1981) In this first step, θ * w is unchanged.Eq. ( 14) allows both for temperatures T * above the melting point of ice, T m = 0 • C, and for T * < T m for non-zero water content θ * w = 0.The temperature T and water content θ w of a cell are thus adjusted according to the energy content Hereby, E = 0 corresponds to ice (i.e., θ w = 0) at temperature T m .Since a decrease in the water content, θ w < θ  responds to refreezing of water, the ice content θ i is increased by θ * w − θ w .The water flux at the upper boundary is given by the rainfall rate, which is added to the water content of the first grid cell.If the water content in a grid cell exceeds the field capacity of the snow, θ fc w (i.e., the maximum volumetric water content the snow can hold) following melt or rainfall, the water is infiltrated in the snow cover in a routing scheme similar to the one employed in Westermann et al. (2011).Hence, a snow grid cell with temperature T can hold the field capacity, plus the amount of water required to increase its temperature to T m , i.e., θ fc w − (T − T m ) c snow /L sl .Excess water is routed to the next cell, until the cell at the bottom of the snowpack is reached, where infiltration is allowed until saturation.The water then starts to pool up from bottom to top, until it reaches the top of the snowpack, where excess water is removed from the system as runoff.The infiltration routine updates the water contents θ w for each grid cell.In a final step, the temperature T is recalculated by applying Eq. ( 17) to account for the refreezing of the infiltrated water.

Time integration
In CryoGrid 3, only the physical process of heat conduction and the associated surface energy balance are integrated in time, while all other processes, such as the infiltration and refreezing of meltwater (Sect.2.5), occur instantly at the end of each time step in the form of a physical consistency check.
To make the code intuitive to use and modify, the simplest possible time integration scheme is employed, first-order forward Euler.Since it is neither computationally efficient nor numerically stable, an adaptive integration time step is used, which prevents too large a change in the internal energy of a grid cell and thus allows a certain degree of control over the accuracy of the solution.
As a first step, the method of lines (Schiesser, 1991) is used to discretize the spatial derivatives in Eqs. ( 12) and ( 14).The change of the internal energy E i of grid cell i (i = 1, N ) over time is given by where z i is the thickness of the i-th grid cell, while z i−1,i is the spacing between and k i−1,i the thermal conductivity assigned to the midpoint of the i −1-th and the i-th grid cells.
For the lowermost grid cell, i = N, the geothermal heat flux is used so that For the uppermost grid cell, i = 1, the discretized form of Eq. ( 1) is used, i.e., Using the first-order forward Euler scheme to integrate Eqs. ( 18)-(20) in time (as given by the first term in Eqs. 12 and 14), the temperature of the grid cell i at time t + t is obtained as The integration time step t is adapted so that change of the internal energy of a grid cell between time steps does not exceed a threshold E max as To avoid the computational costs of solving the coupled equation system for the surface energy balance (Eqs.1-8), the Obukhov length L * of the time step t is employed so that the friction velocity u * and the fluxes Q h and Q e can be computed for the time step t + t following Eqs.( 1)-( 7).The Obukhov length is subsequently updated with the new values of u * , Q h , and Q e of time step t + t.While this procedure introduces numerical errors and potential instabilities, it is adequate for short time steps t during which the components of the energy balance and the associated Obukhov lengths do not change strongly.

CryoGrid 3 Xice -a simple representation of permafrost thaw processes in ice-rich permafrost
Excess ground-ice melt: in many permafrost landscapes, ground layers super-saturated with ice occur, i.e., the volumetric ice content exceeds the porosity of the matrix material when thawed.Such layers are stable when frozen, but become unstable upon thawing, with the water from the melted ground ice becoming mobile.If the meltwater can drain laterally, e.g., in sloping terrain, this will result in subsidence of the surface, while the water can pool up at the surface and create a thermokarst pond in areas with poor drainage.
Either process has the potential to drastically modify the surface energy balance and ground heat transfer and thus change the trajectories of permafrost thaw considerably.CryoGrid 3 Xice facilitates modeling of such phenomena by introducing vertical fluxes of water when a grid cell supersaturated with ice thaws for the first time (Fig. 1).Hereby, we follow the physically based subsidence scheme proposed by Lee et al. ( 2014), but extend it further with a simple representation of energy transfer within shallow water bodies.We emphasize that CryoGrid 3 Xice can only account for the melting of excess ground ice, but not to its aggregation, e.g., due to longterm sedimentation processes, growth of peat, or formation of ice wedges, ice lens, and segregation ice.
In the following, we demonstrate the scheme for the case of grid cells with equal spacing (Fig. 1).For each grid cell, a so-called natural porosity θ nat (e.g., Hopkins, 1949;Nicou et al., 1993) is assigned, which is the porosity of the consolidated soil after excess ground ice has melted and the meltwater has been removed.For saturated soils without excess ground ice, the natural porosity is equal to the ice/water content θ i,w .If a previously frozen grid cell with index i and θ nat < θ i +θ w thaws for the first time (i.e., it reaches a temperature of more than 0 • C and all ice has melted; Fig. 1 b is moved upwards to grid cell i − 1.To conserve volume, a column of mineral and organic ground material of is moved downwards from grid cell i − 1 to i, which ensures that the porosity of grid cell i is now equal to the natural porosity (Fig. 1b 2 /c 2 ).The process is recursively repeated from the initial grid cell upwards, i.e., the increase of water and decrease of mineral/organic contents in grid cell i − 1 will again trigger the condition θ nat < θ i + θ w (Fig. 1b 2 /c 2 ).Hence, water is moved further upwards, while mineral and organic contents move down.When an unsaturated layer is reached, the air is displaced and moved upwards by the water added to the grid cell, similar to Eq. ( 23).If sufficient air is present in the soil column, repeated application of the excess ice routine can lead to an uppermost grid cell, which consists entirely of air.In this case, the grid cell is removed from the soil domain (similar to the procedure during snowmelt, Sect.2.5) and the ground subsides by one grid cell.After all soil air is removed by this mechanism, ongoing melting of excess ice layers triggers the water to pool up at the surface (i.e., the uppermost grid cell consists of 100 % water, Fig. 1c 7 ).If this water is removed through lateral surface or subsurface runoff, the ground subsides (Fig. 1c S 7 ), while a thermokarst pond will form in flat landscapes where lateral fluxes of water are small (Fig. 1c 7 ).A combination of both processes is possible in landscapes with pronounced microtopography, in that thawing ice-rich ground initially subsides until it reaches the level of a surrounding wetland, after which a thermokarst pond is formed.To simulate such effects, a water table depth can be defined in CryoGrid 3 Xice above which all surface water is removed, followed by the formation of a pond at the given water table depth.
When lateral water fluxes lead to ground subsidence, it is evident that this will also change the soil moisture regime in the active layer.However, in the absence of an explicit treatment of the water balance, we cannot simulate such effects in CryoGrid 3 Xice.For simplicity, we assume that the ground remains fully saturated after the excess ice meltwater reaches the surface, i.e., only surface water is removed.We assume the surface cover to be unaffected by excess ice meltwater, so that the parameters related to the surface energy balance, such as the albedo α and the roughness length z 0 , remain unchanged.
Energy transfer in water bodies: if water pools up and forms a thermokarst pond, the surface energy balance and the heat transfer in the ground are modified systematically

Thawed
Frozen Water Ice Mineral+organic Figure 1.The scheme for thawing excess ground ice for an example with saturated ground conditions and equally spaced grid cells (with indicated fractions of ice, water and mineral): when a grid cell with excess ground ice thaws from time step t 1 to t 2 (a→b 1 ), excess water is successively moved upwards from grid cell, while the respective amounts of mineral and organics are moved downwards (b 1 →b 6 ).Following repeated thawing of excess ice grid cells, the uppermost grid cell consists entirely of water (c 7 ) and an initial pond has formed.In case of the subsidence scenario, such grid cells are removed (c S 8 ) and the ground surface is lowered.(Boike et al., 2015).If the uppermost grid cell consists entirely of water, we therefore reduce the surface albedo and set the surface resistance to evapotranspiration to zero.Shallow water bodies such as thermokarst ponds and lakes are generally well mixed during summer (Boike et al., 2015), which strongly increases the energy transfer from the surface to the bottom of the water body.We phenomenologically simulate this process by assigning a uniform temperature to all unfrozen water layers when the uppermost grid cell is unfrozen.In a first step, the temperatures of all grid cells are advanced to the next time step by the normal heat transfer scheme of CryoGrid 3. Second, the weighted average of their temperatures is assigned to all unfrozen water grid cells, thus achieving a uniform profile while conserving the internal energy within the water column.When the surface freezes, a stable stratification forms in the water column, which significantly reduces the heat transfer (Boike et al., 2015).If the uppermost grid cell is frozen, the normal heat transfer scheme of CryoGrid 3 is employed without modifications (i.e., without applying any mixing), which assigns the relatively low thermal conductivity of immobile water (0.57Wm −1 K −1 ) to all water layers.In the model, melting of the winter ice cover in spring occurs from top to bottom, i.e., the ice does not float to the surface as observed in nature.However, since the mixing is applied to all unfrozen grid cells above and below the remaining ice layer, heat is transferred to the bottom of the ice layer and leads to a rapid melting from both sides.
In essence, the normal heat conduction scheme of Cryo-Grid 3 is left unchanged for water bodies, with the only modification being the mixing of water temperatures for summer conditions.While this simple scheme can account for the strong asymmetry of the heat transfer into the water body between winter and summer, it excludes potentially important processes such as radiative heat transfer and buoyancy-driven thermal stratifications in the water column.In Sect.5.1, the scheme is validated against in situ measurements of lake bottom temperatures from Boike et al. (2015).

The Lena River delta
At its mouth in the Laptev Sea, the Lena River forms one of the largest deltas in the Arctic, covering around 25 000 km 2 .Continuous cold permafrost with a mean annual near-surface ground temperature on the order of −10 • C underlies the Lena Delta to about 400-600 m below the surface (Yershov et al., 1991).The landscape has been shaped by the river through erosion and sedimentation (Fedorova et al., 2015), but also permafrost-related thermokarst processes, which have created a large number of water bodies.Within the modern delta, three main geomorphological units can be distinguished according to their late Quaternary genesis, which has led to distinctly different subsurface (cryo-)stratigraphies.
The first river terrace, which covers large parts of the eastern and central delta, is shaped by Holocene river erosion and sedimentation processes.The sediment mainly consists of silty sands, while organic matter has accumulated in thick peat layers (Schwamborn et al., 2002), which feature volumetric organic contents of 5-10 % and mineral contents of 20-40 % (Kutzbach et al., 2004;Zubrzycki et al., 2012).Such sediments are in general super-saturated with ground ice, with average volumetric ice contents reaching 60-80 % and massive ice wedges extending up to 9 m deep in the ground (Grigoriev et al., 1996;Schwamborn et al., 2002;Langer et al., 2013).
The second river terrace in the NW part of the delta was created in a continental setting before and during the transition to the Holocene when the sea level was lower and the gradient of the Lena River steeper than today.It therefore features coarse-grained sandy deposits, which are characterized by low ground ice and organic contents (Schwamborn et al., 2002).Near the surface, an organic upper horizon is often lacking and about 50 % of the area is free of vegetation (Rachold and Grigoriev, 1999;Ulrich et al., 2009;Schneider et al., 2009).The second terrace is often referred to as Arga complex, after Arga Island in the NW part of the delta.
The third river terrace is the oldest geomorphological unit of the Lena Delta, shaped by the extremely cold climate of the last glacial period when the area was not glaciated.Often referred to as ice complex or Yedoma, the ground in the third terrace is extremely rich in ground ice and organic matter, which has not been eroded by the Lena River in the Holocene (Grigoriev, 1993;Zubrzycki et al., 2012).Massive ice wedges reach depths of 20-25 m (Schirrmeister et al., 2003;Grigoriev, 1993;Schwamborn et al., 2002) and in many parts, the ground ice can reach volumetric contents of 90 % (Schirrmeister et al., 2011).The upper soil horizon features a thick organic layer.
The three geomorphological units offer the opportunity to conduct a sensitivity study along a strong gradient of groundice contents, with little to no excess ground ice in the second terrace, while the first and third terrace feature medium and high excess ground-ice contents, respectively.
Samoylov Island provides the unique opportunity to force CryoGrid 3 with a long-term in situ record, while at the same time providing comprehensive high-quality observations for validation.An automatic weather station (AWS) has been operational since 1998 (Boike et al., 2008), measuring air temperature, net radiation, humidity, and wind speed at hourly intervals.Since 2009, all components of the radiation balance are measured, i.e., incoming and outgoing short and long-wave radiation.Eddy covariance measurements of sensible and latent heat fluxes have been ongoing since 2007 in close proximity of the AWS (Langer et al., 2011a, b), and gaps in the time series of the climate station have, when pos-sible, been filled with data from the eddy covariance stations.Furthermore, records of temporary climate stations operated in different time periods on Samoylov Island have been used to fill gaps in the AWS record (for details, see Boike et al., 2013).
The spatial distribution and the onset and termination of the snow cover were recorded with an automated camera.Systematic measurements of thaw depth have been conducted since 2002 on a 150 point grid, using a steel rod pushed vertically into the soil to the depth at which icebonded soil provides firm resistance.The grid points on the polygonal tundra include land surfaces classified as dry and wet tundra and overgrown water (Fig. 12 in Boike et al., 2013).The maximum thaw depth at the end of the summer shows a large interannual variability, with an overall increase of thaw depths since the start of the measurements in 2002.Active layer temperatures were measured since 2002 in three profiles installed across the low centered ice wedge polygon in the center, rim, and slope (the area with a very small topographic gradient between center and rim).The temperatures were recorded using thermistors (107, Campbell Scientific Ltd.), which are calibrated at 0 • C so that the absolute error was less than 0.1 • C over the temperature range ±30 • C. The temperature data used for validation of CryoGrid 3 are from the deepest sensor at 0.40 m depth below surface in the center of the polygon, which is always water-saturated and mostly comprised of peat.Additional detailed information concerning the climate, permafrost, land cover, vegetation, and soil characteristics on Samoylov Island can be found in Boike et al. (2013).

CryoGrid 3 runs for the Lena River delta
To evaluate the performance of CryoGrid 3 and demonstrate the excess ice scheme, two model runs were performed: validation runs for the period 1979-2014 were forced by measured data from the Samoylov AWS and downscaled ERA-Interim reanalysis data to benchmark the model performance.
To assess the potential impact of excess ground ice on the trajectories of permafrost thaw, long-term thaw susceptibility runs from 1901 to 2100 were conducted based on reanalysis data and climate model output.
A first validation run is conducted for a subsurface stratigraphy typical for Samoylov Island on the first river terrace so that it can be compared to in situ observations.Since the assumed snow density has a pronounced influence on ground temperatures (see also Langer et al., 2013), we conduct two runs each for confining values of 200 kg m −3 and 250 kg m −3 , the range obtained from in situ measurements (Boike et al., 2013).An additional run is conducted to benchmark the performance of the simple water body scheme (Sect.2.7) against measurements of lake bottom temperatures on Samoylov Island (Boike et al., 2015).The model setup is unchanged for this lake validation run, but the sub- surface stratigraphy of a well-developed thermokarst lake is applied.
The long-term thaw susceptibility runs investigate the influence of cryostratigraphy on thaw trajectories, using the three geomorphological units in the Lena Delta to perform a sensitivity analysis guided by real-world conditions.Hereby, two scenarios for future warming (RCP 4.5 and 8.5; see below) and two scenarios for ground drainage are considered: in the subsidence scenario, excess water pooling at the surface is completely removed (corresponding to good drainage), while no water is removed for the thermokarst scenario so that a pond forms following thawing of excess ground ice (corresponding to poor drainage).For these runs, a snow density of 225kg m −3 is employed.

Surface, snow, and ground properties
Subsurface properties: to investigate the effect of different (excess) ground-ice contents in the long-term thaw susceptibility runs, a typical stratigraphy is compiled for each of the three river terraces, oriented at available observations (Sect.3.2).In the following, these are referred to as Samoylov (first terrace), Arga (second terrace), and ice complex (third terrace) stratigraphy.
The three stratigraphies are documented in Table 2.We assume that the terraces differ in their near-surface composition whereas in greater depth (> 20 m) the ground is assumed uniform with sandy sediments lacking excess ground ice (Schirrmeister et al., 2011;Schwamborn et al., 2002).For Arga, sandy mineral soil saturated with ice/water (but without excess ground ice) is assumed close to the surface, corresponding to the lack of an upper organic horizon and the general sparsity of vegetation (Sect.3.1).For Samoylov, we assign a typical stratigraphy for a medium wet polygon center, in which soil temperature measurements employed for validation are conducted (Sect.3.2).On top, a medium dry upper organic layer is assumed, followed by a saturated layer with a high water/ice content which extends to a depth of 0.9 m, well below the modeled and measured active layer thicknesses in 2010-2015 (Sect.5).In this layer, we assume that the soil matrix structure of mineral and organic material can support the high soil porosity, so that thawing of ground ice does not lead to vertical or lateral water movement (natural porosity is equal to the true porosity, Table 2).This is changed in the layer below, extending to 9 m, a typical depth of ice wedges (Sect.3.1), in which we assume the soil matrix can only support a porosity of 40 %, so that the excess water becomes mobile upon thawing.The ice complex stratigraphy is similar to Samoylov, but features a drier upper organic layer followed by a saturated layer with high porosity which can yet be supported by the soil matrix.The super-saturated layer below features very high ice contents and extends to a typical depth of ice wedges in the third terrace (Sect.3.1).For the Samoylov and ice complex stratigraphies, the depth below, which excess ground ice is assumed, has been adjusted so that it is below the modeled maximum thaw depth for the years 1901-present (Sect.5) and thermokarst or subsidence do not occur before present-day.
For the validation runs (Table 3), we employ the Samoylov stratigraphy as described above, as well as an estimated stratigraphy for the thermokarst lake Sa_Lake_1 (Boike et al., 2015) on Samoylov Island.For the latter, we prescribe the true water depth of 6 m (Table 3), slightly deeper than the thermokarst lake that would form after the melting of all the excess ground ice in the Samoylov stratigraphy (approx.3.75 m depth).Furthermore, we assume that layers of excess ground ice are no longer present under the lake so that the lake depth remains constant over time.
As far as possible, surface properties derived from measurements on Samoylov Island are employed (Table 4), as published in Langer et al. (2011a, b) and Boike et al. (2013).In addition, the maximum snow depth in CryoGrid 3 is restricted to 0.45 m, the approximate height of the polygon rims above the centers.Snow depth observations of Boike et al. (2013) show that the maximum snow height in polygon centers rarely exceeds this values, since further snowfall is largely blown away by consistently strong winds.

Model forcing data
Validation run: the forcing data for the validation run are based on the available in situ record of incoming shortwave and long-wave radiation S in and L in , air temperature T air , absolute humidity q (computed from relative humidity, RH; Appendix B) and wind speed u from the AWS on Samoylov Island.To fill remaining gaps in the measurement time series and to generate a record for model spin-up prior to the measurement period, surface fields of the ERA-Interim reanalysis (Simmons et al., 2007;Dee et al., 2011) were statistically downscaled for Samoylov Island.For each day of the year, a linear regression was fitted between the ERA-Interim variable and the available measurements for a period of 1 month around the respective date.Using the regression results for each day of year (if a measurement was not available), a time series from 1 June, 1979, to 31 December, 2014, was compiled for the abovementioned variables.The rates of snowfall were used directly from the reanalysis, since reliable measurements are lacking.However, the snow depth is limited to the approximate height of the polygons (Sect.4.1), which ensures realistic maximum annual snow depths independent of the snowfall rates.The largest uncertainties can hence be expected during the buildup phase of the snow cover, before the threshold height is reached.
CryoGrid 3 was initialized to steady-state conditions of the forcing for the first 10 model years by re-running the model with the 1979-1989 forcing until the change of temperatures in the uppermost 100 m was less than 0.05 • C between subsequent 10-year runs.The first initial condition was determined as described in Westermann et al. (2013).
Geosci.Model Dev., 9, 523-546, 2016 www.geosci-model-dev.net/9/523/2016/Long-term thaw susceptibility run: in order to force the model from 1901 to 2100 with a consistent data set, model-derived forcing data building on an anomaly approach were employed which were subsequently downscaled for Samoylov Island similar to the validation run.Past and present-day forcing data were taken from the CRU-NCEP data set (http://dods.extra.cea.fr/data/p529viov/cruncep/), which combines high-frequency variability from the National Centers for Environmental Prediction/National Center for Atmospheric Research reanalysis (Kalnay et al., 1996) with the monthly mean climatologies from the Climatic Research Unit temperature and precipitation data sets (Harris et al., 2014).The meteorological forcing for future climate change experiments was constructed following the Representative Concentration Pathway (RCP) 4.5 and 8.5 climate change projections obtained with the CCSM4 coupled climate model (Meehl et al., 2012).Monthly climate anomalies (for the meteorological state variables: air temperature, specific humidity, surface pressure, and wind speed) and scale factors (for the meteorological flux variables: precipitation, incoming shortwave and long-wave radiation) were extracted from the "Permafrost Carbon RCN forcing data" (ESFG, 2015).The high-frequency variability in the future meteorological forcing was obtained by looping over the 1996-2005 CRU-NCEP forcing data.The base climatological period, which is used to calculate the monthly anomalies/scale factors, is also 1996-2015.These forcing data were also used in Koven et al. (2015) and in an ongoing international intercomparison exercise within the framework of the Permafrost Carbon Network (http://www.permafrostcarbon.org/;workinprogress).
As for the validation run, the data from the closest model grid cell were downscaled for Samoylov Island using the in situ record from the AWS.Since the match with the measured data was significantly worse compared to the ERA-Interim reanalysis (in particular for the CCSM4 data, but also for CRU-NCEP), the entire data set for the overlap periods 2010-2012 for CRU-NCEP and 2012-2014 for CCSM4 was used to determine a linear regression for L in , T air , and q.For the incoming shortwave radiation S in , which follows a defined diurnal pattern, the procedure was conducted separately for each of the time step 00:00, 06:00, 12:00, and 18:00.Values with S in = 0 were excluded from the regression.Downscaling of wind speeds was challenging, since there was no statistically significant correlation to measured data.Therefore, the distribution of wind speeds from CRU-NCEP/CCSM4 was scaled to fit the distribution obtained from the in situ record.For this purpose, the wind speeds were sorted in ascending order in bins containing 20 % of the available data each.The CRU-NCEP/CCSM4 data were then scaled linearly so that the average wind speed of each bin matched the average wind speed of the in situ data in the respective bin.As for the validation run, the snowfall data from CRU-NCEP/CCSM4 were used without modification, but as in the validation run, the maximum height of the snow cover is limited to the approximate height of the polgon rims.The model initialization procedure was identical to the validation run.

Comparison to in situ data for Samoylov Island
Polygon center stratigraphy: the model results of the validation run for the polygon center stratigraphy (Table 3) are compared to published in situ measurements of the surface energy balance and the ground thermal regime from Samoylov Island (Langer et al., 2011a, b;Boike et al., 2013).For the purpose of ground thermal modeling, CryoGrid 3 can reproduce measured surface temperatures very well, except for a slight cold-bias for temperatures colder than −30 • C (Fig. 2).The spread of the data points around the 1 : 1 line could to some extent be explained by small temporal offsets between the real measured data and the ERA-Interim reanalysis upon which part of forcing data are based (Sect.4.2).Furthermore, the Bowen ratio (Q h /Q e ) during the summer months agrees well with published values from Langer et al. tioning at the surface, although the modeled values in winter are slightly smaller in magnitude.The satisfactory agreement of surface energy balance and surface temperatures indicates that CryoGrid 3 can consistently reproduce the landatmosphere interactions on Samoylov Island, despite the simple treatment of summer evapotranspiration with a constant value for the surface resistance against evapotranspiration r s .Boike et al. (2013) provide snow start and end dates, as determined visually from images taken by an automatic camera system.A comparison to the corresponding dates modeled in CryoGrid 3 is not entirely straight-forward particularly for the snowmelt, since the dates obtained from the camera images represent spatial averages, while the model calculates the snowmelt for only a single realization of snow depth.On the average of the years 2001-2010, CryoGrid 3 captures the start and the end of the snow season quite well: the measured and modeled onset date of the snow cover deviate by +3.8 ± 19.2 days, while the snow end dates differ by +1.7 ± 12.2 days.However, the large standard deviations are evidence of significant deviations in some years, which can amount to up to 2 weeks for the snowmelt and more than 1 month for the onset of the snow cover.The latter is completely determined by the forcing data of snowfall (which are obtained from the reanalysis), so that the bias is explained by shortcomings in distinguishing snow-and rainfall in the reanalysis in fall.Significant biases in the snowmelt date can be explained by two reasons: first, the modeled maximum snow depth may be biased due to an erroneous amount of snowfall (which is again determined from the reanalysis).Second, the melting of the snow cover usually occurs in episodes, which are interrupted by cold periods, during which the melting ceases or even new snow falls.Even a small bias in the melt rates can therefore cause a considerable deviation in the termination of the snowmelt, for instance when the snow in the model erroneously persists until the beginning of a cold spell and only melts in the next warm period.
CryoGrid 3 results of the validation run are further compared to measurements of ground temperatures and thaw depths, with a particular focus on the annual dynamics and interannual differences.Fig. 4 displays the comparison of modeled ground temperatures at a depth of 0.4 m in a wet polygon center.In most years, the simulations can successfully reproduce the annual temperature range of about 20 to 25 • C, while the length of the zero curtain, when ground temperatures are confined to 0 • C, is in some years underestimated by up to 1 week.While the reason for this is not entirely clear, we note that the uncertainty of the snow-forcing The measurements correspond to the average of 150 locations on Samoylov Island (Boike et al., 2013).The average standard deviation of the measurements (i.e., the spatial variability of thaw depths) is 0.06 m.The blue area depicts the spread between model runs with snow densities of 200 and 250 kg m .
data is largest during the onset and buildup of the annual snow cover (Sect.4.2), which could at least partly explain biases of modeled ground temperatures during fall freezing.
An example is the abnormally early end of the zero curtain in fall 2010 in the modeled active layer temperatures (Fig. 4), which coincides with a modeled onset of the snow cover 27 days later than observed by Boike et al. (2013).
In the majority of the considered years, CryoGrid 3 results successfully capture both the absolute values of thaw depths and the timing of the thaw progression (Fig. 5).With the exception of the year 2002, the agreement between measurements and model results is excellent, and even small interannual variations on the order of a 5 cm are reproduced (e.g., 2012 vs. 2014).There is no significant difference in modeled thaw depths between model runs for different snow densities, which is in agreement with a study of model sensitivity for Samoylov Island (Langer et al., 2013).We emphasize that the average of 150 point measurements of thaw depth displayed in Fig. 5 masks considerable spatial variability (average standard deviation of 0.06 m) due to variable ground and surface properties.
We conclude that the good agreement between the modeled and measured ground thermal regime can be linked to the surface energy balance, surface temperatures, and the snow dynamics, all of which are adequately reproduced in the model.
Thermokarst lake stratigraphy: for the thermal regime and the melting of excess ground-ice layers below water bodies, the temperature at the lake bottom is a crucial variable.For validation of the simple water body scheme (Sect.2.7), we compare the 3-year in situ record of bottom temperatures from the 6 m deep Sa_Lake_1 on Samoylov Island presented in Boike et al. (2015) to lake bottom temperatures modeled with CryoGrid Xice (see Sect. 4.1 for details).ular the strong asymmetry between winter and summer temperatures.The model can furthermore reproduce the characteristic spike-like dip in temperatures visible in the measurements each fall: at this time, the surface is not frozen yet and persistent well-mixed conditions lead to a cooling of the entire water column towards 0 • C. When the surface freezes, the heat transfer in the water column switches to comparatively inefficient conductive heat transfer, so that the heat loss towards the surface is strongly reduced.The lake bottom is then warmed by the heat stored in the sediments below, before it slowly cools again in the course of winter.During this time, CryoGrid 3 Xice slightly underestimates the measured temperatures, most likely because the energy transfer from the lake sediment is underestimated.During ice break-up in spring, CryoGrid 3 Xice shows a bottom cooling to 0 • C, which is not observed in the measurements, but the timing of the subsequent strong warming at the beginning of summer is once again captured by the model.
While the model slightly underestimates the true annual average temperatures, the comparison to measured lake bottom temperatures suggests that the simple water body scheme of CryoGrid 3 Xice can in general reproduce the thermal regime at the bottom of Sa_Lake_1.We conclude that minor modifications of the model physics (namely the mixing of temperatures in the water column during summer) and changes to model parameters (e.g., the albedo) facilitate reproducing the thermal regime both in cold permafrost ground (Fig. 4, mean annual temperature about −9 • C) and at the bottom of a thermokarst lake (Fig. 6, mean annual temperature about +3.

Long-term thaw susceptibility runs with CryoGrid 3 Xice
The CRU-NCEP/CCSM4 data set facilitates simulations of the ground thermal state over a period of 200 years from 1901 to 2100 for different future warming scenarios.Since these forcing data are derived from a coupled climate model, they feature considerable biases for some forcing variables, although partly corrected as described in Sect.4.2.Therefore, the ground thermal regime modeled for Samoylov Island deviates to a certain extent from the observations.Using the stratigraphy for Samoylov Island (Table 2), the modeled maximum thaw depths for the period 2002-2014 is on average 0.66 m, approximately 0.1 m more than for the more realistic forcing of the validation run (Fig. 5).
The stratigraphy of the ground has a pronounced influence in particular on the thaw depth: for the ice complex stratigraphy, a significantly lower average maximum thaw depth of 0.43 m was obtained, which is a result of the higher ice/water contents underneath the insulating unsaturated top layer (Table 2).For the mineral soils of the Arga stratigraphy, a much deeper active layer with maximum thaw depth of on average 1.29 m was modeled.While systematic measurements of thaw depths are lacking for these two geomorphological units, the modeled values are in general agreement with the few existing observations: on the ice complex on Kurunghak Island approx.10 km SW of Samoylov Island, thaw depths have been determined to be 0.2-0.3m in August, 2013, which is significantly shallower than at the same time in the first terrace on Samoylov Island (Fig. 5).For the second terrace, a point observation on Arga Island (73 • 29 39.2 N, 124 • 22 33.1 E) about 150 km NW of Samoylov Island from 11 August, 2010, yielded a thaw depth of just below 1 m, almost double compared to Samoylov Island on the first terrace.We emphasize that all CryoGrid 3 model runs are based on the forcing data set downscaled for Samoylov Island, so that the results cannot be strictly compared to observations on the second and third terrace.However, the differences in both modeled and measured thaw depths between the three geomorphological units are systematic and significantly larger than the interannual variability on Samoylov Island, which suggests that CryoGrid 3 can reproduce existing patterns of thaw depths for the three geomorphological units in the Lena Delta.
Figures 7 and 8 display the long-term evolution of the annual freeze-thaw state for the three different stratigraphies and two future GHG emission scenarios.An important feature common to all is that stable permafrost is modeled for a period of more than a century, with the active layer undergoing annual freeze-thaw cycles on top of frozen ground.
For the Arga stratigraphy lacking excess ground ice, neither ground subsidence nor thermokarst can occur and the simulations yield identical results for the two cases.For RCP 4.5, permafrost is thermally stable until 2100, with a gradual increase of maximum thaw depth and a ground warming in 10 m depth of around 5 • C (Fig. 9).For RCP 8.5, the warming of the ground accelerates after 2050, which results in the formation of a talik towards the end of the model period, with temperatures at 10 m depth approaching 0 • C.This remarkable ground warming can be linked to a strong increase of incoming long-wave radiation and air temperatures in the CCSM4 data, in particular in fall and early winter.With air temperatures increasing by about 1 • C decade −1 for RCP 8.5, permafrost (which features modeled 10 m temperatures between −8 and −5 • C in 2010) is no longer thermally stable in 2100.The ground warming is amplified by a particularly pronounced temperature increase in fall and early winter, before the annual snow cover has reached its full thickness and thus insulating effect.In the month of November, the warm-ing is almost twice compared to the annual average, with average temperatures increasing by approximately 15 • C from 2010 and 2100, which leads to delayed freezing of the active layer and eventually the formation of a talik.
If layers with excess ground ice are present, the modifications of the model physics implemented in CryoGrid 3 Xice become relevant, as soon as the maximum thaw extends into these layers.For the Samoylov stratigraphy, this is the case around 2060, while thawing of ground layers supersaturated with ice occurs considerably earlier for ice complex.Note that for ice complex the effect occurs a few years earlier in the RCP 4.5 scenario as summer temperatures are slightly warmer than in RCP 8.5 in the 2020s in the applied model forcing data.For the subsidence scenario and RCP 4.5 (Fig. 7), the ground subsides on the order of 1 m in less than 10 years, but subsequently stabilizes since the substantial fractions of mineral and organic material in the ground layers super-saturated with ground ice (Table 2) consolidate and thus form a thicker buffer over the ice-rich layers, so that the annual ground thawing does not reach into the zone with excess ground ice any more.For the ice complex stratigraphy, on the other hand, ground subsidence continues at a reduced rate after the strong initial thaw, since the formation  of the buffer layer is slower than the increase in annual thaw depths due to the low contents of mineral and organic material in the excess ice layers.For the future warming scenario RCP 8.5 (Fig. 8), the ground warming is strong enough that a talik can form within the buffer layer at the end of the simulation period for the Samoylov stratigraphy.For ice complex, however, the simulations do not display talik formation until 2100, despite of a significant increase in thaw depths.The energy input to the ground, which for the Arga and Samoylov stratigraphies leads to talik formation, is to a large extent absorbed by thawing of ice-rich ground in relative proximity to the surface.This process, on the other hand, leads to a drastic subsidence of more than 5 m, accompanied by a warming of mean annual ground temperatures in the permanently frozen ground layers to close to 0 • C. For the subsidence scenario, ground temperatures at 10 m depth below the subsided surface remain below 0 • C in all four cases.We conclude that hydrological conditions leading to drainage of meltwater from excess ground ice delay or even prevent the onset of talik formation compared to a ground stratigraphy without excess ground ice.
For the thermokarst scenario, CryoGrid 3 Xice delivers entirely different trajectories of permafrost thaw.As a consequence of thawing of ground layers with excess ice, a thermokarst pond forms, which rapidly reaches a depth of 1 m or more.As observed in nature (Boike et al., 2012), the pond initially freezes to its bottom each winter for a certain period in the RCP 4.5 scenario (Fig. 7).Upon further warming, the pond does not freeze to its bottom anymore, and a talik forms at the bottom of the pond and the underlying sediment layer.The maximum modeled thickness of the winter ice layer on the pond is on the order of 1-2 m, which is in qualitative agreement with maximum ice thicknesses of around 2 m observed for current climate conditions on thermokarst lakes in Samoylov Island (Boike et al., 2012(Boike et al., , 2015)).For the Samoylov stratigraphy, a significantly shallower pond forms compared to the ice complex stratigraphy, but the thickness of the unfrozen sediments underneath the pond is larger due to higher mineral and organic contents in the excess ground-ice layers.For RCP 8.5 and Samoylov, the entire excess ice layer extending to a depth of 9 m is thawed by around 2090, so that the lake depth does not increase further.For ice complex the excess ice layers extend deeper and the thermokarst lake continues to deepen until the end of the simulation period, with a final depth of around 10 m.Following the initial development of a pond, ground temperatures at 10 m depth rapidly increase (Fig. 9) and eventually become positive, except for the Samoylov stratigraphy and RCP4.5.
The model results suggest that hydrological conditions preventing the drainage of meltwater from excess ground ice lead to the formation of thermokarst ponds and finally talik development even for the moderate RCP 4.5 scenario.Compared to the Arga stratigraphy without layers with excess ground ice, permafrost thaw is considerably accelerated and taliks occur earlier and expand more rapidly in deeper ground layers.
6 Discussion and outlook

CryoGrid 3 -an experimental platform for permafrost process parameterizations
The structure of CryoGrid 3 is designed to be similar to ESM land-surface schemes, so that algorithms and process parameterizations can easily be transferred.While an adaptive time step is employed in the time integration procedure (Sect.2.6), constant time steps as typical in large-scale modeling approaches are feasible, but the numerical stability cannot be guaranteed with the simple forward Euler scheme in this case.Nevertheless, this choice guarantees maximum simplicity and readability of the CryoGrid 3 code, so that the model can be easily modified to test the effect of adding different process parameterizations to the basis module.As such, the goal of CryoGrid 3 is different from existing land-surface schemes, such as GEOtop (Endrizzi et al., 2014) (Jansson and Karlberg, 2004) or SURFEX (Masson et al., 2013), which are aimed at users running the models with a predefined set of model options, largely without modifying the program code itself.In contrast, CryoGrid 3 is optimized to encourage users with the necessary programming skills to modify and extend the program code, thus providing an open platform to explore the effect of different process parameterizations on model results.
In this study, CryoGrid 3 has been applied to cold, continuous permafrost in Siberia.While the satisfactory agreement of model results and in situ observations is encouraging, it is likely that the model performance would be less satisfactory for some environmental conditions encountered in permafrost areas.We see three major shortcomings of the current version of CryoGrid 3, which should be addressed prior to application in a wider range of permafrost environments: 1. Snow physics: CryoGrid 3 employs a constant density for freshly fallen snow and does not account for snow microstructural processes leading to changes of density and thermal properties.While this simple scheme seems adequate for the conditions on Samoylov Island with relatively low snow depths and spread of snow densities (Boike et al., 2013), a more sophisticated treatment of snow processes is warranted if CryoGrid 3 shall be applied, e.g., to mountain environments characterized by a much higher snow cover.Specifically, advanced snowpack models developed for avalanche forecasting, such as CROCUS (Vionnet et al., 2012) or SNOW-PACK (Bartelt and Lehning, 2002), could significantly improve the capacity of CryoGrid 3 to simulate the ground thermal regime, and more flexible, modular approaches like the Factorial Snowpack Model (Essery, 2015) could also be incorporated.
2. Vegetation cover: while CryoGrid 3 can successfully reproduce energy transfer and ground thermal regime in the low-vegetated study area in NE Siberia, it is unlikely that it can account for the near-surface energy transfer in areas with higher vegetation, in particular forests.Many existing schemes simply add a vegetation/canopy layer with defined exchange coefficients for the fluxes of mass and energy on top of the existing soil grid, which is feasible in various levels of complexity (e.g., Zhang et al., 2003;ECMWF, 2006).
3. Water balance: in the study area, seasonal changes of the water level are generally low, in agreement with estimates of the summer water balance showing that evapotranspiration on average more or less balances precipitation (Boike et al., 2013).In many permafrost systems; however, seasonal changes in soil moisture play a pronounced role for energy transfer in the ground (e.g., Hollesen et al., 2011).While other land-surface schemes dedicated to modeling permafrost temperatures, such as GEOtop (Endrizzi et al., 2014), explicitly account for the dynamics of soil water transport by solving Richard's equation, a more simple bucket approach (e.g., Budyko, 1961Budyko, , 1974;;Zaitarian, 1992) may be sufficient in permafrost regions.

Comparison to field and model studies on ground subsidence
The CryoGrid 3 Xice scheme provides a simple way of incorporating ground subsidence and thermokarst in ESM landsurface schemes.With respect to ground subsidence, it is to a large degree similar to the ground subsidence scheme CLM EXICE (Lee et al., 2014), which was employed to model ground subsidence in Arctic permafrost areas until 2100.While Lee et al. ( 2014) found a maximum ground subsidence on the order of 0.5 m in NE Siberia, significantly larger values are obtained in this study, which can be explained by the fact that much lower total excess ice contents were assigned in the large-scale study (maximum 0.5 m water equivalent opposed to 2 m for Samoylov, and almost 10 m for ice complex).A critical issue is the evolution of the soil water content following excess ice thaw, which influences the development of the active layer thickness and thus the subsequent subsidence rates.Furthermore, the thickness and properties of organic surface layers (which have a crucial influence on thaw depths) may change as a result of consolidation or decomposition when excess ground ice below thaws and the soil moisture regime changes.The simple subsidence scheme of Cryo-Grid 3 Xice, which only removes excess water at the surface (Sect.2.7), does not account for these processes, but explicit modeling of the water balance (as discussed in Sect.6.1) could be included to improve the representation of the soil moisture regime and subsidence rates.
A validation of the modeled future subsidence rates is not possible in a strict sense, as ground subsidence does not occur at or before present-day in the long-term thaw susceptibility runs for Samoylov Island, in agreement with in situ observations (Boike et al., 2013).However, the study of Günther et al. (2015) conducted for Muostakh Island off the coast of the Lena River Delta (150 km SE of Samoylov Island, dominated by Pleistocene sediments similar to the ice complex stratigraphy), where strong ground subsidence has occurred in the last 50 years, facilitates comparison of measured and modeled properties of ground subsidence in a qualitative way.The site has experienced significant warming most likely related to changes in the sea ice conditions in the past 50 years, and Günther et al. (2015) document ground subsidence rates of up to 5 m/50 years which is comparable in magnitude to the simulations for the ice complex stratigraphy for the next 50 years.Furthermore, the study clearly documents an inverse relationship between active layer thickness and subsidence rates (Fig. 10, Günther et al., 2015), which they relate to the differences in the ground-ice content below the active layer.They document a ground subsidence of on average 2 m/50 years at a thaw depth of 0.4 m, www.geosci-model-dev.net/9/523/2016/Geosci.Model Dev., 9, 523-546, 2016 about half this value for a thaw depth of 0.65 m, and close to no subsidence for active layer depths exceeding 0.8 m.In our simulations for the three stratigraphies, we can qualitatively reproduce this relationship.We model a current active layer thickness of about 0.4 m with future ground subsidence rates 6-8 m/50 years (ice complex stratigraphy) and 0.65 m with 4 m/50 years subsidence for the Samoylov stratigraphy.For the Arga stratigraphy lacking excess ground ice, the modeled active layer thickness exceeds 1 m, while no subsidence occurs.The exact subsidence rates necessarily depend on the applied forcing (which is certainly different for the past 50 years on Muostakh Island than for the next 50 years on Samoylov Island), but the qualitative agreement with the long-term observations suggests that CryoGrid 3 Xice can capture crucial dependencies on subsurface parameters that have been observed in nature.

Towards modeling permafrost landscape dynamics in regions with ice-rich ground
CryoGrid 3 Xice is based on a 1-D-approximation for processes, which in nature are strongly influenced by lateral fluxes of energy and water.A thermokarst pond can only extend with depth in the model, while in nature it will also expand laterally due to erosion at the shores.Conversely, cooling from surrounding thermally stable permafrost could delay or stop further thawing of the ground underneath an established water body.Such processes have been described in dedicated thermokarst lake models (Plug and West, 2009;Kessler et al., 2012), but the coupled 3-D framework of heat conduction and water and mass movement requiring high spatial resolutions of around 10 m precludes application in ESM land-surface schemes.The algorithm describing excess ice melt in CryoGrid 3 Xice is, on the other hand, simple and highly practical for implementation in large-scale models: if a layer with excess ground ice thaws for the first time, the properties of the above-lying grid cells are rearranged while at the same time conserving mass.As demonstrated by Lee et al. (2014) for CLM, such a step could easily implemented in an ESM land-surface scheme.When a thermokarst pond or lake has formed, heat conduction (as already implemented in the land-surface scheme) continues to be the only means of energy transfer, while mixing of the water column during summer is simulated by assigning an average temperature for the water body after each time step.In this way, wellmixed conditions typical during summer (Boike et al., 2015) are simulated, while the relatively low thermal conductivity of liquid water limits the heat transfer for the stratified water column under an ice cover (Boike et al., 2015).However, it is not clear how this will affect the seasonal dynamics and the long-term evolution of the water body, which should be clarified in a comparison with dedicated lake schemes, such as FLake (Mironov et al., 2010).Furthermore, it must be emphasized that the satisfactory performance of the CryoGrid 3 Xice water body scheme was only confirmed for a well-developed thermokarst lake, which represents the endpoint of the thermokarst simulations for the ice complex and Samoylov stratigraphies (Figs. 7,8).It is not ultimately clear how the scheme would perform for the wide range of smaller and medium-sized ponds on Samoylov Island (Muster et al., 2012;Langer et al., 2015), which necessarily represent initial stages for the development of larger thermokarst lakes.Therefore, we consider the modeled dynamics of the transition from terrestrial permafrost to a deeper lake to be highly uncertain and recommend more research on the factors governing the thermal regime of (in particular) smaller water bodies.
Another critical missing step to describing landscape evolution is to understand the 3-D variability of near-surface ground-ice content and active layer thickness, as well as the resulting lateral fluxes of heat and water.As demonstrated in Sect.5.2, thermokarst/subsidence processes will preferably occur in areas where layers with high excess content exist close to the base of the active layer.In such areas, a small increase of active layer thickness could trigger melting of excess ice and thus initiate transformative landscape changes.

Threshold behavior and metastable states in permafrost systems?
The simulations presented here suggest that ground thawing could display a pronounced threshold behavior in ice-rich permafrost landscapes.Once the active layer extends into layers with excess ground ice, irreversible ice melting with pronounced impact on the ground surface is initiated: when the meltwater is removed, the ground subsides, bringing initially deeper layers with excess ground ice in proximity to the thaw front.The CryoGrid 3 Xice simulations suggest that this process can remove several meters of ground ice within less than a decade (ice complex stratigraphy; Figs. 7, 8).Such depressed areas would become preferential flow paths for surface waters and eventually develop into incised drainage networks (e.g., Fortier et al., 2007).On the other hand, accumulation of sediments from the melted out ice layers will eventually create a protective layer without excess ground ice on top of the remaining layers, so that the permafrost system can return to a thermally stable state with increased active layer thickness, as for the RCP 4.5 simulation for the Samoylov stratigraphy (Fig. 7).A similar metastable state could be reached by a thermokarst pond in which the entire water and ground column refreeze during winter, as typical for small ponds on Samoylov Island (Langer et al., 2015;Boike et al., 2015).In the CryoGrid 3 Xice simulations, the ponds generally form taliks within 10-15 years of their creation (Figs. 7,8), but that occurs in conjunction with ongoing warming.It should therefore be investigated whether CryoGrid 3 Xice could also yield such metastable states for present-day climate forcing on Samoylov Island.The present-day landscape in NE Siberia is evidence that melting of ground ice has shaped many permafrost land-Geosci.Model Dev., 9, 523-546, 2016 www.geosci-model-dev.net/9/523/2016/scapes since the last glacial period (e.g., Schirrmeister et al., 2008).First results provide encouragement that CryoGrid 3 Xice may be developed further into a modeling tool with which such observed landscape changes can be linked to the development of the climate forcing.The presented simulations with CryoGrid 3 Xice indicate that the trajectories of permafrost thaw are crucially influenced by the vertical and lateral distribution of ground ice and processes following excess ground-ice thaw.The stark contrast between the different thaw scenarios for the subsidence and thermokarst cases suggests that simulations of permafrost thaw based exclusively on conductive heat transfer, as employed in the vast majority of studies, are highly questionable in regions with ice-rich ground.

Conclusions
In this study, we present a novel land-surface scheme designed to reproduce the ground thermal regime in permafrost regions.The numerical structure of CryoGrid 3 is simple, using a first-order forward Euler time integration scheme.This relative simplicity makes the new model scheme a flexible platform to explore new parameterizations for a range of permafrost processes.At the upper boundary, CryoGrid 3 is forced by the surface energy balance, while conductive heat transfer is the principal means of energy transfer in the ground.Model results forced with measured time series of meteorological variables from the Samoylov permafrost observatory agree well with in situ records of ground temperatures and active layer thickness, while at the same time consistently reproducing measured surface temperatures and components of the surface energy balance.
In a second step, we introduce a simple parameterization for thawing of excess ground ice, which redistributes the generated meltwater within the 1-D soil column.In a simplified way, the scheme is capable of simulating both ground subsidence and formation of thermokarst ponds.In a long-term simulation from 1901 to 2100, we retrace the rapid transition of a permafrost ground column that is thermally stable throughout the 20th century to a thermokarst pond which eventually develops an underlying talik in the course of the 21st century warming.These initial model results suggest that the cryostratigraphy of the ground and the hydrological regime of a site have a pronounced impact on the permafrost thaw trajectory.In general, ground subsidence due to drainage of meltwater appears to delay the onset of talik formation, while the formation of a thermokarst pond has the opposite effect.In the simulations, the differences between the subsidence and the thermokarst scenario are on the order of the differences between the two investigated future anthropogenic emission scenarios, RCP 4.5 and 8.5.The results support earlier findings that melting of excess ground ice and resulting hydrological processes deserve critical attention in modeling efforts on the future evolution of the worldwide permafrost state.
The Supplement related to this article is available online at doi:10.5194/gmd-9-523-2016-supplement.

Figure 2 .
Figure 2. Modeled vs. measured daily average surface temperatures for in total 2066 days between 2002 and 2009.The measured surface temperatures are obtained from measurements of outgoing long-wave radiation at the Samoylov climate station(Langer et al.,  2011a, b;Boike et al., 2013).The black line represents the 1:1 line.

Figure 3 .
Figure 3. Modeled and measured average surface energy balance for (a) the summer period 7 Jun to 30 Aug 2008 and (b) the winter period 1 December 2007 to 30 Jan 2008, as documented in Langer et al. (2011a, b).Shown are the net radiation (Q net ), sensible, latent, and ground heat flux (Q h , Q e , Q g ), and the residual of the surface energy balance C (0 in the model).The modeled fluxes are taken from the validation run with low snow density (200 kg m −3 ), the values for the high snow density agree within 2 Wm −2 .The measured fluxes correspond to Fig. 9 in Langer et al. (2011b).

Figure 4 .
Figure 4. Modeled and measured ground temperatures at a depth of 0.4 m at a wet polygon center on Samoylov Island (Boike et al., 2013).The blue area depicts the spread between model runs with snow densities of 200 and 250 kg m −3 .

Figure 5 .
Figure 5. Modeled and measured thaw depths on Samoylov Island.The measurements correspond to the average of 150 locations on Samoylov Island(Boike et al., 2013).The average standard deviation of the measurements (i.e., the spatial variability of thaw depths) is 0.06 m.The blue area depicts the spread between model runs with snow densities of 200 and 250 kg m .
presents clear evidence that the model can reproduce the main features of the annual temperature regime, in partic-

Figure 6 .
Figure 6.Modeled and measured temperature at the bottom of a 6 m deep thermokarst lake on Samoylov Island for the years 2009 to 2012.The measurements are described in detail in Boike et al. (2015), the model results are derived from validation runs for a lake stratigraphy with snow densities of 200 kg m −3 and of 250 kg m −3 (see Sect. 5.1).

Figure 7 .
Figure 7. Annual freeze-thaw state for three different cryostratigraphies (top to bottom: Arga, Samoylov, ice complex) and the drainage scenarios subsidence (left) and thermokarst (right), as simulated by CryoGrid 3 Xice forced by CRU-NCEP/CCSM4 and the moderate future warming scenario RCP 4.5.

Figure 8 .
Figure 8. Annual freeze-thaw state for three different cryostratigraphies (top to bottom: Arga, Samoylov, ice complex) and the drainage scenarios subsidence (left) and thermokarst (right), as simulated by CryoGrid 3 Xice forced by CRU-NCEP/CCSM4 and the future warming scenario RCP 8.5.

Figure 9 .
Figure 9. Annual average temperatures at 10 m below the ground surface (i.e., the actual ground surface in case of subsidence) for the years 1901-2100; two future warming scenarios and three different (cryo-)stratigraphies of the ground.

Table 1 .
Westermann et al., 2013)t processes that are accounted for in the models CryoGrid (abbreviated CG) 2 (seeWestermann et al., 2013)and 3, as well as in the excess ice scheme Xice.Simplified parameterizations in brackets.

Table 2 .
Sediment stratigraphies of the long-term susceptibility runs, with volumetric fractions of the soil constituents and soil type for each layer given.Layers super-saturated with ice in bold.

Table 3 .
Sediment stratigraphies of the validation runs, with volumetric fractions of the soil constituents and soil type for each layer given.