Freeze / thaw processes in complex permafrost landscapes of northern Siberia simulated using the TEM ecosystem model : impact of thermokarst ponds and lakes

Freeze/thaw (F/T) processes can be quite different under the various land surface types found in the complex tundra of the Arctic, such as polygonal tundra (wet center and dry rims), ponds, and thermokarst lakes. Proper simulation of these different processes is essential for accurate prediction of the release of greenhouse gases under a warming climate scenario. In this study we have incorporated the water layer into a dynamic organic soil version of the Terrestrial Ecosystem Model (DOS-TEM), having first verified and validated the model. Results showed that (1) the DOS-TEM was very efficient and its results compared well with analytical solutions for idealized cases, and (2) despite a number of limitations and uncertainties in the modeling, the simulations compared reasonably well with in situ measurements from polygon rims, polygon centers (with and without water), and lakes on Samoylov Island, Siberia, indicating the suitability of the DOS-TEM for simulating the various F/T processes. Sensitivity tests were performed on the effects of water depth and our results indicated that both water and snow cover are very important in the simulated thermal processes, for both polygon centers and lakes. We therefore concluded that the polygon rims and polygon centers (with various maximum water depths) should be considered separately, and that the dynamics of water depth in both polygons and lakes should be taken into account when simulating thermal processes for methane emission studies.


Introduction
The release of greenhouse gases from the large quantities of soil carbon preserved in Arctic regions constitutes an important feedback to climatic warming and the thawing of permafrost north of 45 • N (McGuire et al., 2009; Schneider von  Deimling et al., 2012).Reliable simulation of the dynamics of permafrost is therefore critical when predicting future climatic changes.The energy balance at the ground surface has an important influence on variations in permafrost.Heterogeneous ground surfaces with, for example, variable snowpack or organic layer thicknesses exert a major influence on the surface energy balance (Etzelmüller and Frauenfeld, 2009) and have in the past been integrated into both land surface models (Yi et al., 2007;Lawrence and Slater, 2008) and ecosystem models (Zhuang et al., 2001;Yi et al., 2009aYi et al., , b, 2010)).However, few of the current large-scale land surface models or ecosystem models take into account the effects that water bodies have on the dynamics of permafrost (Zhuang et al., 2006;Ringeval et al., 2012), with one exception being the model by Wania et al. (2009) which treated surface water in the same way as a litter layer.Water bodies of various sizes, ranging from those occupying polygon centers to large thermokarst lakes, are distributed across the Arctic coastal regions (French, 2007) resulting in considerable landscape heterogeneity.These water bodies have a marked effect on the surface energy balance and thermal dynamics of the surrounding permafrost soils (French, 2007).Their presence can lead to permafrost degradation, which in turn affects the terrestrial ecosystem's carbon budget.Outgassing of carbon S. Yi et al.: Freeze/thaw processes in complex permafrost landscapes dioxide from ponds and lakes was, for example, calculated to account for between 74 and 81 % of the net landscape-scale CO 2 emissions in polygonal tundra of the Lena Delta, Siberia (Abnizova et al., 2012).
There are a number of different techniques used to simulate permafrost dynamics (Riseborough et al., 2008).A wide range of numerical models exist, which are applied both in stand-alone permafrost simulations and in land surface schemes of climate models.The temperature and water content at different depths in soil or rock are calculated numerically.Numerical solutions for permafrost dynamics in large scale models are commonly obtained by solving finite difference equations.One category of numerical solution, referred to by Zhang et al. (2008) as "decoupled energy conservation parameterization", assumes that the soil water is homogeneous and freezes or thaws at exactly 0 • C. Soil temperature is calculated for each layer and if the temperature of a particular layer is then greater than 0 • C, some or all of any ice present will melt and the temperature is then recalculated, and vice versa.This is an efficient method and is commonly used in land surface models (Zhang et al., 2003;Oleson et al., 2004).However, the lower layers in land surface models are usually thick and the freezing or thawing fronts derived from soil temperature interpolation are not realistic (Yi et al., 2006).
A second category of numerical methods, referred to by Zhang et al. (2008) as "apparent heat capacity parameterization", assumes that soil water freezes or thaws over a range of temperatures below 0 • C and simulates both the unfrozen soil water content and the temperature, simultaneously.Since small changes in soil temperature within the freeze/thaw range will result in a large change in apparent heat capacity, an iterative procedure is required to ensure that only small temperature changes occur during each time interval (Nicolsky et al., 2007).This method is commonly applied in permafrost models (Goodrich, 1978;Nicolsky et al., 2007;Dall'Amico et al., 2011;Hipp et al., 2012;Langer et al., 2013) and has also recently been applied in a land surface model (Ringeval et al., 2012).Although the method is more physically realistic it requires greater computing resources, which may lead to limitations in the spatial resolution, the length of time that can be modeled, and the number of simulated land surface classes, etc.
Both categories of numerical models have their disadvantages when they are applied for regional permafrost simulation.Apart from numerical models, analytical solutions also exist that can be used for solving phase change problems.For example, exact Neumann solutions to freezing and thawing problems exist for idealized cases, such as for infinite or semi-infinite homogeneous material, steady upper boundary conditions, etc. (Lunardini, 1981).Stefan's equation, which was originally used to predict the thickness of sea ice, is widely used due to its simple form (Lunardini, 1981); an algorithm for applying Stefan's equation to a layered system (e.g., soil) was developed by Jumikis (1977) and applied in a hydrological model by Fox (1992).However, predictions from the Stefan algorithm usually overestimate the depths of freeze/thaw fronts as it neglects any heat transport beneath the front.In order to mitigate this problem of overestimation, Woo et al. (2004) developed a two-directional Stefan algorithm (TDSA).Yi et al. (2009a) integrated a TDSA within a terrestrial ecosystem model (TEM) in order to first simulate the depths of freezing or thawing fronts, and then update the soil temperatures for layers above the uppermost front, beneath the lowermost front, and between these two fronts.This is an efficient method and is able to track the positions of fronts within thick soil layers.
Although models using the above methods to simulate permafrost dynamics over large regions have been validated using in situ measurements, few of them have been verified against analytical solutions for both freeze/thaw fronts and soil temperatures at different depths, which is just as important for model validation (Romanovsky and Osterkamp, 1997).Furthermore, conceptual and numerical permafrost landscape models also require suitable methods for upscaling from local to global scales, which need to take into account field-based knowledge of the surface characteristics, of key processes, and of key parameters (Boike et al., 2012).
In this study we aimed to develop and test a model that could simulate permafrost dynamics under different types of land surface, i.e., different thicknesses of snow cover, of the organic layer, and of water cover.The model used in this study was a dynamic organic soil version of the TEM (DOS-TEM), which was originally developed for, and tested on, boreal forest sites.We first verified the DOS-TEM simulations of freeze and thaw processes with analytical solutions for idealized cases.We then modified the DOS-TEM to take into account the effects that water bodies of various sizes have on the thermal dynamics of permafrost and compared the output with in situ measurements from Samoylov Island in the Lena Delta, Siberia.Finally, we compared the simulations beneath different land surface types in order to investigate the vulnerability of permafrost to water bodies.

Site description
Samoylov Island (72 • 22 N, 126 • 30 E) is located in the southcentral part of the Lena River delta in Siberia (Fig. 1); it covers an area of about 4.3 km 2 .The average annual mean air temperature on Samoylov Island from 1998 to 2011 was −12.5 • C and the average total summer rainfall 125 mm (Boike et al., 2013).Samoylov Island contains two major geomorphological units: a floodplain, and an elevated Holocene terrace that is characterized by low-centered polygonal tundra.The elevated terrace comprises ∼ 70 % of the total area of the island and contains numerous ponds and thermokarst lakes.On average, the land surface of the terrace consists of 58 % dry tundra, 17 % wet tundra, and 25 % water surfaces (Muster et al., 2012).A 26.75 m borehole was drilled in 2006 in an area that consists of about 60 % polygon centers and 40 % polygon rims, with a negligible areal proportion of ponds (Fig. 1).Based on the dating of sediments on Samoylov (Schwamborn et al., 2002) we estimated the maximum ages of the lakes to be 2.8 kyr.Further information concerning the climate, permafrost, vegetation, and soil characteristics can be found in Boike et al. (2013).

Meteorological data processing
The collection of meteorological measurements on Samoylov Island started in 1998.The daily mean air temperature, wind speed, vapor pressure, net radiation, downward solar radiation, and total daily precipitation were calculated from hourly measurements.If more than 25 % of the measurements were missing in any one day, no value was recorded for that day.If more than 25 % of the daily values in a particular month were missing, no value was recorded for that month.We replaced the missing monthly values as follows: 1. Air temperature and precipitation (snow + rain) measurements for the same month, available from the nearby Stolb meteorological station (which has data sets from 1956, but with large gaps during the 1970s), were used to replace the missing values.
2. Long-term-mean values were used to replace some values for air temperature and precipitation that remained missing after step (1) above, as well as missing values for wind speed, radiation, and vapor pressure.We calculated the long-term monthly mean for air temperature and precipitation between 1981 and 2011 using measurements from the Stolb meteorological station, and for wind speed, downward shortwave radiation, and vapor pressure between 1998 and 2011 using measurements from the Samoylov site.
To illustrate the differences between different data sets, we compared the monthly air temperature and precipitation data sets from Samoylov Island with those from Stolb and the global reanalysis data set from the Climate Research Unit (CRU TS3.1), available from http://badc.nerc.ac.uk/ view/badc.nerc.ac.uk_ATOM_dataent_1256223773328276 (last access: 1 October 2012) (Fig. 2).

Model descriptions
The Terrestrial Ecosystem Model (TEM) family of models is designed to simulate the carbon and nitrogen pools within  vegetation and soil, and the carbon and nitrogen fluxes between vegetation, soil, and the atmosphere (McGuire et al., 1992).The most recent TEM version (i.e., the DOS-TEM) can simulate the dynamics of organic soil layers, which can be subject to fire disturbances and to ecological successions (Yi et al., 2010).The DOS-TEM consists of four modules, these being the environmental, ecological, fire disturbance, and dynamic organic soil modules.The environmental module operates on a daily time interval using mean daily air temperature, surface solar radiation, precipitation, and vapor pressure, which are downscaled from monthly input data (Yi et al., 2009a).It takes into account radiation and water fluxes between the atmosphere, canopy, snowpack, and soil.The soil moisture levels and temperatures of all soil layers are updated daily.A two-directional Stefan algorithm is used to predict the depths of freezing or thawing fronts within the soil (Woo et al., 2004); it first simulates the depth of the front in the soil column from the top downward, using soil surface temperature as the driving temperature; it then simulates the front from the bottom upward using the soil temperature at a specified depth beneath a front as the driving temperature (bottom-up forcing).If a layer contains a freezing or thawing front, this layer is then divided into two layers (Fig. A1b).The temperatures of soil layers above the uppermost freezing or thawing front and beneath the lowermost freezing or thawing front are updated separately by solving finite difference equations.The thermal properties of soil layers are affected by their water content (Yi et al., 2009a).

Model modifications
We made three modifications to the DOS-TEM in order to simulate the effects that water bodies (Fig. A1a) have on freezing or thawing processes.(1) As shown in Fig. A3, the volumetric water content of a polygon center (i.e., the water between the two innermost vertical dashed lines in the figure) is not equal to 1 due to the slope of the rims surrounding the water bodies.We took this into account by calculating the volumetric water content of different layers within water bodies of various sizes (Fig. A3); details are presented in Appendix B.
(2) When updating the thermal state of water layers they were treated in the same way as soil layers, but with different thermal properties.We followed the model of Hostetler and Bartlein (1990) to calculate the eddy diffusion coefficients for the water layers, which were then used, together with the molecular diffusion coefficient of water, to calculate the heat transfer within the water bodies and the heat exchange with the underlying sediments: details are presented in Appendix B.
(3) The original DOS-TEM only simulated bottom-up forcing for the deepest freezing or thawing front.However, taliks probably exist beneath some water bodies, and more than two freezing or thawing fronts may exist at the same time.We therefore implemented top-down and bottom-up forcing separately for each front (Appendix A).
The soil thermal conductivity in the DOS-TEM was initially calculated according to Farouki (1986).However, preliminary testing showed that the resulting soil thermal conductivities were higher than those derived from field measurements (Langer et al., 2011a and b).We therefore used the more realistic parameterization according to Johansen (1975) and Côté and Konrad (2005).Further details on the parameterization used are provided in Appendix C.

Comparisons with analytical solutions
Three different materials were tested in this study, i.e., water, mineral, and organic soil: the properties of these materials are listed in Table 1.The initial temperature of each material at different depths (up to 5000 m in the DOS-TEM) was set to −10 • C, and the temperature at the upper boundary of each material was set to 5 • C over the whole simulation period (100 years).We assumed zero heat flux conditions at the lower boundary, i.e., at 5000 m depth.The temperatures and the depth of the thawing front obtained from the DOS-TEM were compared with those from analytical solutions, and with those obtained using the one-directional Stefan equation.For the DOS-TEM, the temperature at a specific depth was calculated by linear interpolation between the temperatures of overlying and underlying layers.To test the sensitivity of the model to the depth used for the bottom-up forcing, we tried bottom-up forcing at different depths below the thawing front (i.e., at 50 cm, 1, 2, 5, and 20 m).In order to test the effects of total soil/water thickness, we also evaluated the DOS-TEM using different depths for the lower 3) 20 cm organic rich soil (wet organic, p = 0.9, vwc = 0.7) center 0.25 0 1 m organic soil (saturated organic, p = 0.9, vwc = 0.9) pond 1.1 1.05As for center lake 6 6 As for center boundary (50, 500, and 5000 m).The maximal thickness of the soil/water layer was set to 1, 10, and 100 m for runs with the lower boundary at 50, 500, and 5000 m depths, so that the total number of layers was constant for each run.

Test sites
We tested the DOS-TEM for soil or water temperatures at four different sites, i.e., on a polygon rim (rim), in a polygon center without standing water (center), in a polygon center with standing water (pond), and in a larger thermokarst lake (lake).These sites are considered to represent the most important types of land surface in the polygonal tundra landscape of Samoylov Island.The configurations of water and organic soil characteristics for the different land surface types used in the model are presented in Table 2.We used about 65 m of mineral soils (saturated sand with a porosity of 0.6) in 12 layers.The DOS-TEM assumes bedrock beneath the soil layer (Fig. A1a); in each case we used 420 m of bedrock in five layers to represent the frozen sediments on Samoylov Island.The ground heat flux at the bottom of the bedrock was set to 0.053 W m −2 (Pollack et al., 1993).The simulated soil temperatures for the four different land surface types were compared to temperature measurements from a 27 m borehole on Samoylov Island (Boike et al., 2013).

Surface temperatures
The DOS-TEM is not able to simulate the surface temperatures of water, land, or snow.We therefore used linear regression to establish the relationship between measured daily surface temperatures and air temperatures during those periods of 2011 in which both temperatures were above 0 • C, as follows: for water T surf = 0.563 T air + 4.735 (coefficient of determination (R 2 ) = 0.41, number of pairs of data (n) = 84), for land, T surf = 0.643 T air + 2.231 (R 2 = 0.54, n = 84).For snow, frozen soil and frozen water we assumed T surf = T air .A similar method has previously been used in Yi et al. (2013).

Snow cover
Wind drift is an important process that redistributes snow on the polygonal tundra landscape.Field measurements of annual maximum snow thickness usually show depths of 15-40 cm in polygon centers and much less on polygon rims (10-30 cm, with an average of about 15 cm) and frozen lakes (Boike et al., 2013).Zhang et al. (2012) introduced a snow drift factor into their NEST (Northern Ecosystem Soil Temperature) model.The factors for rims, centers, and water bodies (ponds and lakes) were 0.5, 0, and −0.25, respectively, with a positive value indicating a loss of snow due to wind drift.However, a preliminary model run indicated that the simulated snow thicknesses were overestimates, for all sites.In this study we therefore identified site-specific threshold values for maximum snow accumulation based on field observations, as follows: where D snw,max is the maximum snow thickness (m), and H is the microtopographic relief (m) (see Fig. A3).We also performed an additional simulation using observed snow thickness as a forcing at the center site.

Soil and water properties
For the soil thermal properties we used two sets of parameters, one derived from field temperature measurements (Langer et al., 2011a, b) and the other calculated from an algorithm proposed by Luo et al. (2009), details of which can be found in Appendix C (Table 3).For water, we increased the calculated value of the eddy diffusion coefficient by a factor of between 10 and 100 (following Subin et al., 2012), in order to take into account the effects of convection currents caused by complex lake topography and density instability.

Initialization
The rim, center, and pond sites were all initialized using a temperature of −10 • C for all water, soil, and bedrock layers; the lake site was initiated with −10 • C for all soil and bedrock layers and with 0 • C for water layers.For the equilibrium run, the model was forced by an average annual cycle that was generated using the monthly averages from climate data available from 1981 to 2011.When the difference in annual mean unfrozen soil thickness between two consecutive years was less than 0.01 cm the model was considered to be in a state of equilibrium.The equilibrium runtime was set to between 50 and 400 years for the rim, center and pond sites, and between 50 and 3000 years for the lake site.The period from 1981 to 2003 was used for model spin-up, and we compared the simulations with measurements collected after 2003.

Effects of (maximum) water depth
Polygon centers and lakes of various sizes and water depths are distributed across much of Samoylov Island.In order to investigate the effect that the size and water depth of polygon ponds and lakes have on the thickness of the underlying unfrozen soil we ran the DOS-TEM for a shallow, medium, and deep polygon pond (with maximum water depths of 20, 60, and 120 cm), and for a shallow, medium, and deep lake (with maximum water depths of 2, 4, and 6 m).For each polygon pond or lake the model was run with water depths of between 0 and 100 % of the maximum water depth, at intervals equal to 5 % of the maximum water depth.The annual mean unfrozen soil thickness from 2003 to 2011 was calculated for each run, for comparison purposes.

Comparisons with exact Neumann solutions and Stefan equations
The bottom-up forcing in the DOS-TEM is very important for accurate simulation of the position of the thawing front using Stefan's algorithm (Fig. 3).The root mean squared errors (RMSEs, n = 36 500) between thawing fronts simulated without bottom-up forcing and those from exact Neumann solutions for three different idealized cases were greater than 1.128 m.In contrast, the RMSEs between the thawing fronts simulated with bottom-up forcing and those from exact Neumann solutions were less than 0.047 m (Table 4).For all cases of water, mineral soil, and organic soil, the thawing fronts simulated without bottom-up forcing were very close to those calculated using the Stefan equation.The simulated water or soil temperatures and thawing fronts were not sensitive to the depth of bottom-up forcing (Fig. 3).For example, there was almost no difference between the thawing fronts simulated for bottom-up forcing at depths of between 0.5 and 20 m, in all three cases (water, mineral soil, and organic soil).The differences between thawing front simulations using bottom-up forcing and those from Neumann solutions were also very small (Fig. 3).Taking bottom-up forcing at a depth of 1 m beneath the thawing front as an example, most of the RMSEs for temperatures at depths shallower than 1 m were less than 0.01 • C, and   for depths greater than 1m they were approximately 0.1 • C (Fig. 4, Table 5).
The simulated temperatures were sensitive to the total thicknesses of the various materials, especially to the total thickness of mineral soil, which has the highest thermal conductivity and the lowest water content (Table 4).

Snow thicknesses
The simulated snow thickness from the DOS-TEM was greater than 80 cm at all sites for 2005-2006, and decreased thereafter (Fig. 5).However, measurements at the center site showed that the monthly maximum snow thickness was only 40 cm.After setting a maximum snow thickness, the differences in snow thickness between the four sites were similar to field observations, but the interannual variability was very small.Since we assumed that snow only accumulates on frozen layers of water or soil, the starting date for snow accumulation at pond and lake sites was usually later than at rim and center sites.The simulated starting dates for snow accumulation in the autumn of 2010 were about 1 month later than the observed starting dates.

Temperatures of shallow layers
For the rim site, soil temperatures for model runs that included snow drift compared well with actual measurements at depths of both 2 and 51 cm (Fig. 6).The simulated soil  temperatures at 51 cm were slightly lower than measured temperatures during the summer months.The simulated soil temperatures using the calculated thermal properties (Appendix C) were close to those simulated using the derived thermal properties at 2 cm depth but varied by about 1-3 • C at 51 cm depth.The effect of snow was very obvious: where no maximum snow thickness had been set the simulated soil temperatures were up to 10 • C warmer than the measured soil temperatures.
For the center site, the performance of the DOS-TEM was similar to the rim site during the summer seasons (Fig. 7).The DOS-TEM overestimated the soil temperatures at 40 cm depth in several of the winters.Using different soil thermal properties did not result in any obvious differences in soil temperature, and setting a maximum snow thickness had less effect than for the rim site.When observed snow thicknesses were used as a forcing, the simulated soil temperatures were closer to observed soil temperatures than other simulations for the winter seasons of 2003-2005 and 2006-2008 (Fig. 7), for which full seasonal records of observed snow thickness were available (Fig. 5).For both rim and center sites, the simulated soil temperatures fell rapidly during the fall of 2010, possibly due to the later snowfall in the simulation (Figs. 5,6,7).
For the pond site, the seasonal cycle of simulated water temperatures had a smaller amplitude than the observations (Fig. 8).For example, the simulated water temperature in the lower part of the pond site was 20 • C warmer than actual measurements from the winter of 2008-2009, and in summer it was ∼ 2 • C cooler on average.As an additional experiment we reduced the maximum snow thickness from 15 to 2 cm, which brought the simulated water temperatures in winter down to the measured temperatures.Changing the water eddy diffusion coefficient by a factor of between 10 and 100 did not result in any obvious differences between model runs.
For the lake site, the simulated water temperatures in the upper part of the lake were not as sensitive to the eddy diffusion coefficient as those in the lower part of the lake (Fig. 9).The simulation using the default water eddy diffusion coefficient considerably underestimated the water temperature (by about 10 • C) in the lower part of the lake.Increasing the eddy diffusion coefficient by a factor of either 10 or 100 improved the simulation.In the simulation with the eddy diffusion coefficient increased by a factor of 10, the eddy diffusion coefficient increased from 1.6 × 10 −4 m 2 s −1 at a depth of 0.8 cm to 5.8 × 10 −3 m 2 s −1 at a depth of 37 cm, and then decreased to 3.3 × 10 −4 m 2 s −1 at a depth of 2.2 m and to 2.4 × 10 −7 m 2 s −1 at a depth of 5.5 m during August, at the end of equilibrium run.The eddy diffusion coefficient of water under ice was about 1.6 × 10 −6 m 2 s −1 .
In the following two subsections we only analyze the freezing and thawing fronts and the deeper soil temperatures of the four sites on the basis of simulations that used a maximum snow thickness, derived soil thermal properties, and an eddy diffusion coefficient increased by a factor of 10.

Freezing and thawing fronts and unfrozen soil thicknesses
The simulated shapes of freezing and thawing fronts at the rim and center sites were similar from 2003 to 2011 (Fig. 10).The thawing fronts did not survive through the winter months and into the following year.However, multiple thawing and freezing fronts were simulated at the pond site.In an additional test performed with 2 cm maximum snow thickness, the soil temperature was colder than it was with 20 cm maximum snow thickness and the shapes of the thawing fronts were different (Fig. 10c, d).From 2003 to 2011 the average maximum depth of thawing fronts in soils under water was 0.47 for simulations with 2 cm maximum snow thicknesses and 3.86 m for those with 20 cm maximum snow thickness.The simulated thawing fronts at the lake site occurred at an average depth of 43.28 m below the lake floor.The volumetric water content (VWC) at the bottom of centers or lakes was only 0.06 less (at the most) than that at the top of water (Table A1).We investigated the effects of including a VWC calculation by comparing two sets of simulations: one taking into account the effects of water depth on the VWC, and the other using a VWC of 1.The differences of multiyear mean unfrozen soil thicknesses between two simulations were very small (less than 1 cm) for the lake and pond sites with 2 cm maximum snow thickness.However, for the pond site with 20 cm maximum snow thickness, the small difference in unfrozen soil thickness between VWC = 1 and VWC < 1 simulations was accumulated.For example, the difference in unfrozen soil thicknesses between simulations with VWC = 1 and VWC < 1 was about 0.33, 1.34 and 2.06 m for maximum equilibrium runtimes set to 50, 200, and 400 years, respectively.
In order to investigate the effects of eddy diffusivity and equilibrium runtime, we performed additional sensitivity tests on the lake site with three different eddy diffusivities (the original value of ke1, and additional values of ke10, and ke100) and three equilibrium runtimes (2000, 3000, and 4000 years).The simulation with ke1 and a 2000-year equilibrium runtime had the smallest talik thickness (31.71 m), while that with ke100 and a 4000-year equilibrium runtime had the greatest unfrozen soil thickness (53.02 m).For a particular equilibrium runtime, increasing the eddy diffusivity from ke1 to ke10 had a greater effect on unfrozen soil thickness than increasing it from ke10 to ke100.For a particular eddy diffusivity value, the effect of increasing the equilibrium runtime from 2000 to 3000 years was similar to that of increasing it from 3000 to 4000 years (Table 6).Previous versions of the DOS-TEM only considered topdown forcing for the uppermost front and bottom-up forcing for the lowest front.In this study, we implemented both topdown and bottom-up forcing for each front.There were very small differences at the rim, center and pond sites, but major differences at the lake site.For example, with eddy diffusivity increased by a factor of 10 and a 3000-year equilibrium runtime the simulated unfrozen soil thickness over the period from 2003 to 2011 was about 43.28 m using our version of the DOS-TEM, compared to 40.68 m using the previous version (Yi et al., 2009a).

Temperatures of deep layers
The averages of the modeled annual mean soil temperatures at 26.75 m depth over the period from 2007 to 2011 were approximately −10.16, −9.14, −0.99, and 1.08 • C for the rim, center, pond, and lake sites, respectively (Fig. 11).As with the unfrozen soil thicknesses, the soil temperature at 26.75 m depth at the lake site was sensitive to changes in eddy diffusivity and equilibrium runtime (Table 6), with the exception that there was almost no difference between simulations with an equilibrium runtime of 4000 years but different eddy diffusivities.

Effects of (maximum) water depth
For polygon centers with small maximum water depths (e.g., 20 and 60 cm), increasing the water depth caused a slight increase in the multiyear (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) mean unfrozen soil thickness (Fig. 12a).For deep polygon centers (e.g., up to 1.2 m deep), the unfrozen soil thickness started to increase rapidly when the water depth was about 0.8 m and reached a maximum (of about 7 m) at a water depth of around 1.05 m, before it decreased again rapidly for water depths greater than 1.05 m.When the water depth was less than 0.8 m, increasing the water depth caused a semi-linear increase in deep (27.5 m) soil temperature (Fig. 12c).When the water depth was close to the maximum water depth, the deep soil temperature was slightly reduced (Fig. 12c).This was due to a decrease in the maximum snow thickness, which facilitated the dissipation of energy out of the soil during the cold seasons.As was also the case for the pond site, an increase in the water depth caused an increase in the deep soil temperature when the water depth was less than 1 m (Fig. 12d).Further increasing the water depth increased the unfrozen sediment thickness when the lake water depth was between 1.0 and 2.0 m (Fig. 12b).Both the unfrozen sediment thickness and the deep soil temperature were relatively constant when water depth was greater than 2.0 m.This was due to the decrease in eddy diffusivity at greater water depths.For example, eddy diffusivities at water depths of 2.2 and 5.5 m were about 3.2 × 10 −4 and 2.4 × 10 −7 m 2 s −1 , respectively, during warm seasons.A further increase in the water depths did not result in any warming of sediments when using the current eddy diffusivity schemes.

Performance of the DOS-TEM
The DOS-TEM is able to simulate multiple freezing and thawing fronts simultaneously.For example, there was a year-long thawing front in the sediment and a freezing and a thawing seasonal front in the water of the lake site, and there were multiple freezing and thawing fronts at the pond site (Fig. 10).The simulated thawing fronts and soil/water temperatures at different depths compared very well with analytical Neumann solutions for all three materials, and the accuracy was not sensitive to the depth of bottom-up forcing (Figs. 3, 4).The simulated soil/water temperatures compared reasonably well with in situ measurements from Samoylov Island (Figs. 6,7,8,9).
Modeling studies usually consider thin layers within the uppermost meter of the soil column to depict daily or seasonal soil temperature dynamics, and thick layers at greater depths to improve computing efficiency (e.g., Oleson et al., 2004), which can result in problems when attempting to simulate the thermal dynamics of talik.The ground layering in the DOS-TEM (including snow, water, soil, and bedrock) follows a similar strategy.Unlike other models, the DOS-TEM can track multiple fronts, even within a single layer.The updating of freezing/thawing fronts is based on the Stefan algorithm, which is very efficient.The subsequent updates of temperature in each layer are similar to those of land surface models using the Crank-Nicolson algorithm, which is also an efficient algorithm.The DOS-TEM itself is therefore an efficient model; for example, it takes only about 10 s to simulate a period of 100 years.

Effects of different land cover types
Our investigations have largely concentrated the effects that snow thickness and water depth have on soil thermal dynamics for the different land cover types.Snow thickness has a strong impact on the soil/water temperatures during cold seasons by insulating the underlying soil or water (Stieglitz et al., 2003;Gouttevin et al., 2012).At the rim site, differences between soil temperatures simulated with and without taking snow drift into account were usually greater than 10 • C (Fig. 6).At the center site, the simulated snow thickness was greater than the measured thickness (Fig. 5) and the simulated soil temperatures were warmer than the measured temperatures, while when using measured snow thicknesses the simulated winter soil temperatures were close to the measured temperatures (Fig. 7).
The timing of snowfall is an additional factor affecting the thermal dynamics of soils.The snowfalls simulated for the fall of 2010 were later than the measured snowfalls (Fig. 5), which resulted in early decreases in soil temperatures at the rim and center sites (Figs. 6, 7).In the real world, early snowfall might be expected to melt in the unfrozen water of the pond and lake sites; in the model, snow accumulated only after the first 2 cm of water was frozen.In the fall of 2010 there was therefore no time lag between the simulated and measured water temperatures (Figs. 8,9).
Water ponding has a very important influence, increasing the underlying unfrozen soil thickness (Figs. 10, 12a, b) and the deep soil temperatures (Figs. 11,12c,d).The difference in water eddy diffusion coefficients between warm and cold seasons is responsible for the warming effects of water.Our sensitivity tests indicate that, under the present climate on Samoylov Island, the unfrozen soil thickness is very sensitive to water depths of between 1 and 2 m for thermokarst lakes, and between 0.8 and 1.2 m for polygon centers.This has significant implications for the development of talik under thermokarst lakes: following the melting of segregated ice under polygonal tundra (and associated surface subsidence) the development of a thermokarst lake or polygon center can accelerate if the water depth in the lake or pond exceeds a certain threshold.Similarly, talik beneath a thermokarst lake can disappear if the water depth in the lake falls below a certain threshold (van Huissteden et al., 2011).

Limitations and uncertainties
Because of the harsh Arctic environment some measurements of atmospheric variables were not available from Samoylov Island and the missing values were replaced with those from the nearby Stolb meteorological station.Air temperatures from the Stolb station (T _stolb) compared very well with those from Samoylov Island (T _samoylov): T _stolb = 0.97 T _samoylov + 0.65; R 2 (R: correlation coefficient) = 0.99; n = 80.The growing season precipitation at the Stolb station (P _stolb) also compared reasonably well with that for Samoylov Island (P _stolb = 0.62 P _samoylov + 8.35; R 2 = 0.53; n = 37), with averages of 26.4 and 29.3 mm month −1 , respectively.Since there were no precipitation measurements for the cold seasons on Samoylov Island, it is impossible to assess any uncertainty associated with snowfall.
Running ecosystem models for regional or global applications requires large-scale reanalysis data sets, such as the global data sets of the CRU, the European Centre for Medium-Range Weather Forecasts (ECMWF), or the National Centers for Environmental Prediction -National Center for Atmospheric Research (NCEP-NCAR).In this study we compared the air temperature and precipitation from the CRU data set with those from Samoylov Island.Air temperatures from the CRU were close to those for Samoylov Island in summer, but about 15 • C colder in January (Fig. 2).The monthly average precipitation in the growing season between 1998 and 2009 was 41.2 mm month −1 from the CRU and 29.3 mm month −1 for Samoylov Island.It is clearly important to investigate the uncertainties associated with input data when using models for large-scale cold region applications (Clein et al., 2007).
Wind drift is a common process involved in redistributing snow on the heterogeneous landscape of the Arctic tundra (Sturm et al., 2001).There are, however, no measurements of snowfall and snow cover thickness available for the various terrain units of Samoylov Island, making the parameterization of snow drift impossible.Zhang et al. (2012) used snow drift factors (scale method) and in this study we have set maximal snow thicknesses (capping method) to simulate the differences in snow thicknesses between different types of land surface.The capping method is better suited to this particular site since the snow depth is influenced by microtopography and wind drift.A snow drift factor would alter the entire snow dynamics which would yield incorrect results, especially during the early stages of winter when the snow accumulation is mainly controlled by precipitation.A snow drift factor would produce unrealistically slow snow accumulation rates.An assumption that the snow accumulates as a result of precipitation but cannot exceed a certain threshold seems to yield more realistic results for the polygonal tundra.The threshold value represents the height to which snow can accumulate without being blown away, which may represent an effective wind shadow height for the microtopographic structure on Samoylov Island.Both methods are, however, empirical.Measurements will in future need to be collected in situ in order to develop valid parameterizations for snow drift.
The surface temperatures of snow, soil, and water are critical boundary conditions for solving finite difference equations; they are dependent on atmospheric conditions as well as the snow/soil/water conditions (Yi et al., 2013).In models with hourly time steps, snow/soil/water surface temperatures are calculated by iteratively solving the surface energy balance equation for the different surfaces.This involves incoming and reflected solar radiation, incoming and outgoing longwave radiation, sensible and latent heat fluxes between the surface and the atmosphere, and ground heat flux (Oleson et al., 2004).In this study we have used a regression model to calculate surface temperatures on the basis of existing measurements.These algorithms performed better for the rim and center sites than for the pond and lake sites.
The exchange of energy in water bodies is not only a result of molecular diffusion and eddy diffusion, but also of other processes such as convection caused by water density instability and complex lake-bottom shapes, which have not been taken into account in this study.We followed Subin et al. (2012) to simulate these effects implicitly by increasing the eddy diffusion coefficient.For example, in order to agree with the dynamics of water temperatures at the bottom of the lake (6 m depth) for the lake site, the eddy diffusion coefficient had to be increased by at least a factor of 10 (Fig. 9).Extensive work is required to test this approach over other lakes in different regions.
Models are commonly run with a multiyear mean climate in order to obtain a state of equilibrium, which should only relate to the climatic and land surface characteristics and should not be affected by the amount of time used.However, the length of time used for an equilibrium run affects the simulated unfrozen soil thickness of polygon centers and lakes when the water body depth exceeds a threshold value.For example, it took less than 100 years for the DOS-TEM to reach equilibrium at the rim and center sites but it never reached equilibrium at the lake site, even after 4000 years.It is possible that the sediment of the thermokarst lake is always in a nonequilibrium state.The traditional method of determining the initial equilibrium state might not be suitable for thermokarst lakes.It is therefore important to have actual measurements of talik thicknesses beneath water bodies in order to determine the number of years required for an equilibrium run and to validate model outputs.Unfortunately, however, such information is not readily available at present.A new technology known as surface nuclear magnetic resonance has recently been used over thermokarst lakes to measure the underlying talik thickness (Parsekian et al., 2013) and promises to provide useful information on talik that can be used to improve modeling in future studies.
Lateral heat gradients clearly exist beneath the different land surface types of the polygonal tundra (Fig. 11), but since the DOS-TEM is a one-dimensional model it is unable to simulate lateral heat exchange.A two-or three-dimensional model would be better able to simulate the thermal processes in complex Arctic tundra landscapes (e.g., Ling and Zhang, 2003;Plug and West, 2009;van Huissteden et al., 2011;Kessler et al., 2012), but such models are difficult to apply over large regions.
Thermal processes vary under the different land surfaces of the heterogeneous polygonal tundra.For example, talik was present under the lake site, but not under the center or pond sites.The volumetric water content at the bottom of the water bodies varies according to the size of the water body (Table A1), which resulted in distinct differences in the simulated thicknesses of unfrozen soil when under thick snow cover (see Sect. 3.2.3).In our study we have assumed fixed shapes for polygon centers and lakes, but in future studies it would be desirable to take into account the dynamics of thermokarst lake development.

Outlook
Land surfaces are heterogeneous on different spatial scales and LSMs with low resolutions (usually hundreds of kilometers) take this heterogeneity into account using different technologies.Early LSMs only considered the major land surface type within each grid cell (Manabe, 1969); parameters of different land surface types were subsequently aggregated for each grid (Arain et al., 1999).With recent advances in computing power and remote sensing technology it has become possible to explicitly consider different types of land surface, such as those with different plant function types, urban areas, water, etc. (Oleson et al., 2004).Our study has indicated that the heterogeneity of Arctic polygonal tundra results in marked differences in soil thermal dynamics (Figs. 10,11).In order to simulate methane emissions from polygonal tundra ecosystems on a regional scale it is therefore crucial to distinguish between polygon rims, polygon centers (with varying water levels), and thermokarst lakes at different stages of development.The sensitivity analysis suggests that it is necessary to, at the very least, consider polygon rims, polygon centers (with maximum water depths of less than 0.8 m and with several water depths between 0.8 and 1.2 m), and lakes.The following regional inputs can be obtained for the above-mentioned types: 1.The proportion of surface water over regions of polygonal tundra ecosystem can be retrieved from remote sensing albedo data sets (Muster et al., 2013) and the maximal proportion of surface water over different periods.
2. The distribution of the area covered by polygon centers can be established following Cresto Aleina et al. (2013).
3. The relationship between water area and water depth can be established on the basis of in situ measurement data (Wischnewski, 2013).

Conclusions
In this study we have modified an ecosystem model to simulate thermal processes under the different land surface types of a polygonal tundra landscape on Samoylov Island, in the Lena Delta, Siberia.The simulated freeze/thaw dynamics and soil/water temperatures compared very well with analytical Neumann solutions for three different materials in idealized runs.Despite a number of limitations and uncertainties relating to model parameterization and data input, the simulated soil/water temperatures compared reasonably well with in situ measurements.The modified model is thus very efficient and suitable for large-scale regional applications.Microtopographic relief has an important effect on snow and water cover, which in turn exert important influences on the different thermal processes that operate under the various land surface types.Sensitivity tests have indicated that thermal processes are very sensitive to changes in water depth when the depth is between approximately 1 and 2 m for lakes, and between 0.8 and 1.2 m for polygon centers.The different land surface types of polygonal tundra ecosystems therefore need to be taken into account in large-scale ecosystem models, as well as water dynamics, in order to be able to accurately simulate methane emissions.eddy diffusion, and buoyant convection, among others.In the DOS-TEM we took into account molecular diffusion, eddy diffusion (which is usually 2-3 orders of magnitude greater than molecular diffusion; Subin et al., 2012), and other processes, by increasing the eddy diffusion coefficient by a factor of between 10 and 100.In cold seasons with snow and ice cover, the dissipation of energy to the atmosphere would only be realized by molecular diffusion, while in warm seasons and with open water the exchange of energy within the water would be much greater.The seasonal variation in energy exchange coefficients is therefore an important factor in the development of unfrozen soil beneath water bodies.Water layers were treated in the same way as soil and snow layers but with different thermal properties (Fig. A1a) when calculating the positions of freezing or thawing fronts and the temperatures within water bodies.Following Hostetler and Bartlein (1990), the governing equation for the onedimensional model is where T is the water/soil/snow temperature (K), t is the time (s), z is the depth from water surface (m), C is the volumetric heat capacity (J (m −3 K)), λ is the thermal conductivity of water/soil/snow (J mKs −1 ), K is the conductivity due to eddy diffusion (for water only, J (mKs) −1 ), and is a heat source term (W m −2 ).The detailed parameterization of K and can be found in Hostetler and Bartlein (1990).

Figure 1 .
Figure 1.(a) Circumpolar permafrost distribution (Brown et al., 1998) and the Lena River delta.(b) Location of the Samoylov study site within the Lena River delta, eastern Siberia (NASA, 2000), and (c) measurement locations on Samoylov Island.

Figure 2 .
Figure 2. Monthly air temperature (T a ) and precipitation (P ) measurements from Samoylov Island (Samoylov), from the Stolb meteorological station (Stolb), and from the CRU's global data set.

Figure 3 .
Figure 3. Comparisons of outputs from DOS-TEM simulations, exact Neumann solutions (Exact), and Stefan's equation (Stefan) for (a) water, (b) mineral soil,and (c) organic soil over a 100-year period.The term B50CM means simulations from the DOS-TEM with bottom-up forcing at 50 cm beneath the lowest freezing or thawing front, and likewise for other similar terms.NOBOT means no bottom-up forcing.The outputs from the DOS-TEM have been plotted for the middle of every 10th year and different cases have been started from different years in order to make the figures more readable.

Figure 4 .
Figure 4. Comparisons of outputs from DOS-TEM simulations (dashed lines) and exact Neumann solutions (solid lines) for (a) water, (b) mineral soil, and (c) organic soil over a period of 100 years, at depths from 0 to 20 m.

Figure 5 .
Figure 5. Comparisons of simulated maximum (monthly) snow thicknesses at the center, rim, pond, and lake sites, and the default values (with no maximum snow thickness set), with those from field measurements at the center site (center-obs), over the period from 2003 to 2011.

Figure 6 .
Figure6.Comparisons of monthly average soil temperatures at 2 and 51 cm depths below the rim site from simulations (mod), with and without a maximum snow thickness (b and nb), using derived and default thermal properties (der and def), with those from field measurements (obs), over the period from 2003 to 2011.

Figure 7 .
Figure7.Comparisons of monthly average soil temperatures at 1 and 40 cm depths below the center site from simulations (mod), with and without a maximum snow thickness setting (b and nb), using the observed snow thickness (snw-obs), and using derived and default thermal properties (der and def), with those from field measurements (obs), over the period from 2003 to 2011.

Figure 8 .
Figure 8. Comparisons of monthly average water temperatures at 1, 69 and 86 cm depths below the water surface at the pond site from simulations with 20 (mod-20cm) and 2 cm (mod-2cm) maximum snow thickness, over the period from 2007 to 2011.

Figure 9 .
Figure9.Comparisons of monthly average water temperatures at 2 (a), 4 (b) and 6 m (c) below the water surface at the lake site from simulations (mod) using default, 10, and 100 times the eddy diffusion coefficient (ke, ke10, and ke100), with field measurements (obs) over the period from 2009 to 2011.

Figure 10 .
Figure 10.Simulated soil freezing (blue) and thawing (red) fronts, and water freezing (green) and thawing (yellow) fronts (depths in meters) over the period from 2003 to 2011 for (a) the rim site, (b) the center site, (c) the pond site with 20 cm maximum snow thickness, (d) the pond site with 2 cm maximum snow thickness, and (e) the lake site.The surface of the soil was taken to be at 0 m depth, with the downward direction positive and the upward direction negative.

Figure 11 .
Figure 11.Comparisons between simulated (dashed lines) and measured (solid lines) values for annual mean (green), maximum (red), and minimum (blue) soil temperatures ( • C) averaged over the period from 2007 to 2011 for the (a) rim, (b) center, (c) pond, and (d) lake sites.

Figure 12 .
Figure 12.Responses of simulated multiyear (2003-2011) mean unfrozen soil thickness and soil temperature at 27.5 m depth to changes in water depths in (a) and (c) polygon centers (with maximum water depths of 0.2, 0.6 and 1.2 m) and in (b) and (d) lakes (with maximum water depths of 2.0, 4.0 and 6.0 m).

Figure A1 .
Figure A1.The ground components considered in the DOS-TEM (a), and a diagram of updating freezing and thawing fronts of ground components, together with temperatures (b).DX, TC, HC, and PCE stand for thickness, thermal conductivity, heat capacity, and energy used for phase change, respectively; there are altogether m layers: 1, 2, 3, n and m are layer indexes; DZT and DZB are the distances between the top-down and bottom-up driving depths (dashed and dotted horizontal lines) and the position of the front.Freezing and thawing fronts are indicated with blue and red solid lines, respectively.Vertical dashed lines indicate layers not shown.T tb , T bb , and T drv are the top boundary, bottom boundary and ground surface driving temperatures, respectively.The T bb at the bottom of the ground structure is determined by the temperature and thermal properties of the overlying layer and the prescribed heat flux.

Figure A2 .
Figure A2.The flowchart for the calculation of upper temperatures for updating the position of the front from above.(Note: front(s) between position and the front in the flowchart means there are front(s) between the front under consideration and the position of driving.)

Figure A3 .
Figure A3.Diagram of rim and polygon center or lake, which are separated by vertical dashed lines.r top is the radius of the top and r bot that of the bottom of a polygon center or lake; WD, WD max and H are the water depth, maximum water depth, and microtopographic relief height in a polygon center or lake, respectively.

Table 1 .
The thermal conductivity, volumetric heat capacity, volumetric water content, and porosity used in idealized runs for water, mineral soils, and organic soils.

Table 2 .
Water and organic soil configurations used in the model for the different test sites (N.A. not available).

Table 3 .
Thermal properties of different types of soil on Samoylov Island.The first is derived using soil temperature measurements; the second is calculated using the default scheme in the DOS-TEM.

Table 4 .
The root mean squared error (n = 36 500) between the thawing fronts (m) from exact Neumann solutions and simulated thawing fronts from the DOS-TEM, with different combinations of total thickness (50, 500, and 5000 m) and bottom-up forcing (b1m: bottom-up forcing at 1 m below front; nobot: no bottom-up forcing) for different materials.

Table 5 .
The root mean squared error (n = 36 500) between the temperatures ( • C) from exact Neumann solutions and simulated temperatures from the DOS-TEM for different materials, with 5000 m total thickness and bottom-up forcing at 1 m below the thawing front, at depths of between 0.05 and 20 m.

Table A1 .
The area, top radius (r top ), bottom radius (r bot ), maximum water depth (WD max ) of polygons or lakes, and their volumetric water contents (vwc) at the top and bottom.Area (m 2 ) r top (m) r bot (m) WD max (m) VWC Mean average surface area of the smallest polygon centers on Samoylov Island (Wischnewski, 2013); b mean average surface area of the largest polygon centers surveyed on Samoylov Island; c area of the large thermokarst lake on Samoylov Island. a