Simulating high-latitude permafrost regions by the JSBACH terrestrial ecosystem model

The current version of JSBACH incorporates phenomena specific to high latitudes: freeze/thaw processes, coupling thermal and hydrological processes in a layered soil scheme, defining a multilayer snow representation and an insulating moss cover. Evaluations using comprehensive Arctic data sets show comparable results at the site, basin, continental and circumarctic scales. Such comparisons highlight the need to include processes relevant to high-latitude systems in order to capture the dynamics, and therefore realistically predict the evolution of this climatically critical biome.


Introduction
The effects of global climate change are felt stronger in the northern high latitudes than elsewhere in the world (ACIA, 2005).During recent decades, polar regions have experienced an increase from around +0.5 to +1 • C in surface atmospheric temperatures, while the global mean has risen by only from +0.2 to +0.3 • C (Serreze et al., 2000).Furthermore, soil temperature in the Arctic is also undergoing warming, which is observed from borehole and active-layer measurements.After the International Polar Year (2007)(2008), these measurements were summarized to show that permafrost is warming and active-layer thickness is increasing in the Nordic regions, Russia, and North America (Christiansen et al., 2010;Romanovsky et al., 2010a;Smith et al., 2010).
Based on a simple relationship between air temperature and the permafrost probability, Gruber (2012) estimated that around 22 % (±3 %) of the Northern Hemisphere land is underlain by permafrost.During the past glacial/interglacial cycles vast amounts of organic matter have been accumulated in these soils (Zimov et al., 2006).With the abundant resources in interglacial periods, life has flourished and left huge amounts of organic matter behind; while the glacial periods created unfavorable conditions for decomposition and kept the remnants locked away in the frozen soil (DeConto et al., 2012;Schirrmeister et al., 2013).Supporting that, recent findings on the amount of soil carbon in northern circumpolar permafrost soils are larger than the previous estimates (Hugelius et al., 2010;Ping et al., 2008;Tarnocai et al., 2009;Zimov et al., 2006).According to Tarnocai et al. (2009), there are 1672 Pg of carbon stocked in the northern permafrost soils.With the current trend of increasing air temperature, this carbon rich soil is susceptible to thawing and being released to the atmosphere in the form of greenhouse gases and thus contributing to even further warming of the atmosphere (Heimann and Reichstein, 2008;Schuur et al., 2008;ACIA, 2005).Therefore, it is important to understand the underlying processes and to quantify future interactions of permafrost regions within a changing climate (Beer, 2008).
The recognition of this importance has spurred recent advances of dynamic global vegetation models and Earth system models by representing processes that are specific to Published by Copernicus Publications on behalf of the European Geosciences Union.

A. Ekici et al.: Simulating high-latitude permafrost regions
high-latitude regions.With the understanding of feedback mechanisms and recent estimates of vast amounts of soil carbon, progress has been made to address uncertainties in Arctic simulations.At present, most of the global models include common processes related to permafrost regions, e.g., latent heat release/consumption from the phase change of soil water (Riseborough et al., 2008), organic matter decomposition at freezing conditions, methanogenesis and methanerelated processes.Li et al. (2010) have shown a comprehensive review of different freezing schemes in sophisticated models.However, within the global models either an extra term of latent heat is added (e.g., Mölders et al., 2003;Takata and Kimoto, 2000) or the method of "apparent heat capacity" is incorporated into temperature calculations (e.g., Beer et al., 2007;Hinzman et al., 1998;Nicolsky et al., 2007;Oelke, 2003;Poutou et al., 2004;Schaefer et al., 2009).In either way, the models showed a significant improvement in simulating soil temperature or active-layer thickness (e.g., Dankers et al., 2011;Gouttevin et al., 2012a;Lawrence et al., 2012;Zhang et al., 2008).
Besides the freeze/thaw events, the coupling of snow and soil thermal constitutes the basis for the soil thermal profile during winter (Dutra et al., 2010;Slater et al., 2001;Stieglitz et al., 2003).Due to strong insulating properties of snow, the winter soil temperature is kept warmer than the much colder atmospheric temperature.Furthermore, the timing of snowmelt influences the duration of the growing season and the active-layer thickness, which is also related to the amount of infiltrating snowmelt water into the soil.Goodrich (1982), Kelley et al. (1968) and Groffman et al. (2006) found that snow cover strongly influences the ground thermal regime.Using the ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) model, Gouttevin et al. (2012b) showed that the snow cover and the disappearance of snow are important factors for the plant and soil metabolic activity and biogeochemical feedbacks between the soil and the atmosphere.However, in most cases snow is represented rather simply in the models.Due to the high complexity of snow types and snow processes, a simple parameterization yielding a realistic heat insulation effect was used (e.g., Beer et al., 2007;Koren et al., 1999;Verseghy, 1991).While more advanced snow schemes were developed in some modeling studies (Boone and Etchevers, 2001;Loth and Graf, 1998), it is not always practical for global modeling exercises to include such a complex approach due to its computational requirements.
Impacts of changing permafrost conditions on the climate system and vegetation activity have also been investigated.It is shown by Poutou et al. (2004) that including soil freezing in their model leads to dryer summers and warmer winters in different regions.Beer et al. (2007) have found out that with the permafrost-specific processes the high-latitude vegetation carbon stocks are better represented in a dynamic global vegetation model.In other modeling studies, future implications of possible permafrost carbon release are investigated and their effects on global climate are shown under different warming scenarios (Burke et al., 2012;Hayes et al., 2011;Koven et al., 2011;Schaefer et al., 2011;Schneider von Deimling et al., 2011;Zhuang et al., 2006).A good review of permafrost carbon cycle models is documented in McGuire et al. (2009).
Although progress has been undertaken on representing permafrost processes in land surface models, there is still a considerable uncertainty regarding the magnitude of the effects of permafrost feedbacks on climate.A consensus is not yet close to being reached regarding the timing of permafrost response to climate change and consequences of permafrost feedback mechanisms on the climate system.An intercomparison study of different land surface schemes especially with respect to cold regions' climate and hydrology revealed large differences between the models, even in case of a similar implementation of frozen ground physics (Luo et al., 2003).Due to missing processes and related deficiencies of their land surface schemes, climate models often show substantial biases in hydrological variables over high northern latitudes (Luo et al., 2003;Swenson et al., 2012).Thus, the representation of the complex dynamics of permafrostrelated processes within global models is a challenging yet essential task (Hagemann et al., 2013).To contribute to this progress, we have advanced the land surface model JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg) and we show the reliability of the new model version in multiscale evaluations.

Model description and improvements
JSBACH is the land surface component of the Max Planck Institute Earth System Model (MPI-ESM) that comprises ECHAM6 for the atmosphere (Stevens et al., 2012) and MPI-OM (Max Planck Institute Ocean Model) for the ocean (Jungclaus et al., 2012).It is designed to serve as a land surface boundary for the atmosphere in the coupled simulations; but it can also be used offline given that it is a comprehensive terrestrial ecosystem model with a process-based approach for representing key ecosystem functions.JSBACH simulates photosynthesis, phenology and land physics with hydrological and biogeochemical cycles (Raddatz et al., 2007;Brovkin et al., 2009).The photosynthesis scheme follows Farquhar et al. (1980) and Collatz et al. (1992).The BETHY (Biosphere Energy-Transfer Hydrology) model (Knorr, 2000) covers most of the fast canopy processes.The current version employs a relatively simple carbon cycle model (Raddatz et al., 2007).Vegetation carbon is classified as "green", "wood" or "reserve" carbon and these are transported into soil carbon pools via litter fluxes.The soil organic matter is stored in "fast" or "slow" soil carbon pools with different decomposition rates.All carbon pools have a constant turnover time, which is only modified by temperature and moisture in the case of soil carbon pools.
The current version of the model can be run with 30 min temporal resolution and global simulations are usually performed at 0.5 • spatial resolution; however, the 1-D point model can also be run for a single location.The grid cells are usually divided into tiles of homogeneous vegetation cover.In the version discussed here, we prescribed the vegetation cover and kept it constant over time (cf.Sects.2.2 and 2.3).
The soil is discretized into five layers with increasing thicknesses (Fig. 1).Heat conduction through the vertical soil layers is assumed to be the dominant method of heat transfer; therefore convective and radiative heat transfer processes are neglected.Surface temperature is calculated by considering incoming radiation and surface albedo, then it is used as the upper boundary forcing for the soil temperature calculations.During the snow period, the uppermost snow layer is forced by surface temperature and the bottom snow-layer temperature is used to force the soil column.In the simulations mentioned here, a constant moss layer is present over the soil.Hence the upper boundary condition for the soil temperature calculations is the moss-layer temperature, while a zero heat flux is assumed for the bottom boundary condition at 10 m depth.The one-dimensional heat transfer equation (Eq. 1) is solved for each layer.For each time step, and each soil layer, the numerical solution to heat conduction (first term on the right side of Eq. 1) gives the soil-layer temperature and then as a second step, this temperature is updated with respect to the heat used for (or gained from) phase change of soil water (second term on the right side of Eq. 1) in that layer.This routine continues from the top to the bottom to calculate all the soil-layer temperatures.
with T : soil-layer temperature (K), c: volumetric heat capacity of the soil layer (J m −3 K −1 ), λ: heat conductivity of the soil layer (W K −1 m −1 ), L f : Latent heat of fusion (J kg −1 ), ρ i : density of ice (kg m −3 ), θ i : volumetric soil ice content (m 3 ',m −3 ), t: time (s), and z: soil-layer depth (m).JSBACH mainly uses the physics package of ECHAM5 (Roeckner et al., 2003).This comprises the separation of rainfall and snowmelt into surface runoff and infiltration and the calculation of lateral drainage following the Arno scheme (Dümenil and Todini 1992).A new soil hydrology scheme (Hagemann and Stacke, 2014) has been implemented into JSBACH that uses the same five-layer structure (see Fig. 1) as the thermal module and calculates soil water transport by using the one-dimensional Richards equation (Richards, 1931) shown in Eq. ( 2).Here, the local change rate of moisture ∂θ w /∂t is related to vertical diffusion (first term on the right side of Eq. 2) and percolation by gravitational drainage of water (second term).Both processes are considered separately.Percolation is calculated following the Van Genuchten (1980) method and the diffusion is calculated using the Richtmyer and Morton (1967) diffusion scheme.For the latter, the soil water diffusivity D of each layer is parameterized following Clapp and Hornberger (1978).
The soil water content may be greater than 0 for each layer above the bedrock.There is no water available for the land surface scheme below the bedrock.Consequently, horizontal drainage (ECHAM4 formulation following Dümenil and Todini, 1992) may occur only from those layers above the bedrock.The formulation has been slightly modified as now drainage may only occur if the soil moisture is above the wilting point.Note that the previously used bucket model soil moisture now corresponds to the root zone soil moisture.The associated rooting depth determines the depth from where transpiration may occur.Bare soil evaporation is occurring only from the uppermost layer.
In the hydrology module, first the input/output terms (precipitation, snowmelt, evapotranspiration) are accumulated and infiltrated into (removed from) the soil.Then, the phase change routine updates the water and ice contents of each layer before the vertical water movement is executed.Each layer's field capacity is updated with the corresponding layer's ice content that is created or melted in the same time step.This allows for a more realistic water transport within the frozen layers.Finally the vertical water movement is performed as described above and the soil water content at each layer is updated.with θ w : volumetric soil water content (m 3 m −3 ), D: soil water diffusivity (m 2 s −1 ), K: soil hydraulic conductivity (m s −1 ), and S: source and sink terms (s −1 ).
As shown in Eq. ( 3), a supercooled water formulation is also incorporated to allow liquid water to coexist with ice under freezing temperatures.This approach follows the Niu and Yang (2006) formulation.
with θ w max : maximum supercooled water content (m), θ sat : soil porosity (m 3 m −3 ), T frz : freezing temperature of water (K), g: gravitational acceleration (m s −2 ), ψ sat : saturated soil matric potential (m), and b: Clapp and Hornberger exponent (-).Soil heat transfer is coupled with the hydrological scheme through latent heat from phase change and two parameters: the volumetric heat capacity (c) and the soil heat conductivity (λ) in Eq. ( 1).We have parameterized the heat capacity using the de Vries ( 1963) formulation (Eq.4) and the heat conductivity following Johansen's (1975)method (Eq.5).Equations (6-9) describe the terms in Eq. ( 5).With these formulations, the amounts of water and ice influence the soil thermal properties.In concert with the latent heat of fusion effect on temperature (second term on the right side of Eq. 1), a coupling of the hydrology and soil thermal dynamics is achieved.For Eq. ( 8), bulk density needs to be inserted with the given unit below.
with ρ s , ρ w , and ρ i : density of soil solids, water and ice, respectively (kg m −3 ); c s , c w , and c i : specific heat capacities of soil solids, water and ice, respectively (J kg −1 K −1 ).
with K e : Kersten number (-), λ sat : heat conductivity of the saturated soil (W K −1 m −1 ), λ dry : heat conductivity of the dry soil (W K −1 m −1 ), Sat: saturation ((θ w + θ i )/z/θ sat ), λ s , λ w , and λ i : heat conductivities of soil solids, water and ice, respectively (W K −1 m −1 ), and ρ bulk : soil bulk density (kg m −3 ).Snow is treated as external layers above the soil column.With increasing snow depth in winter, new layers are added up to maximum of five snow layers.The top four layers are always 5 cm in thickness, while the bottom layer is unlimited in size.A 5 cm snow layer is always kept in contact with the atmosphere in order to maintain the numerical stability due to rapid changes in air temperature.The uncertainty of representing 5 cm snow layers is assumed to be negligible when compared to having a nonlayered snow scheme.The snow properties are kept constant for simplicity.A snow density of 250 kg m −3 is used for the snow depth calculations and the snow heat conductivity is fixed at 0.31 WK −1 m −1 with a snow heat capacity of 522 500 Jm −3 K −1 .This simple approach is chosen to ensure the heat insulation for the soil rather than providing a complex snow model.For this reason, the snow layers are hydrologically inactive, meaning there is no water held within each snow layer, thus neither the transfer of meltwater within the snowpack nor refreezing effects are considered.Water infiltration from snowmelt into the soil is treated separately in the hydrology module.
In addition to the snow layers, the importance of moss cover in the Arctic is mentioned in several studies (Beringer et al., 2001;Rinke et al., 2008).The moss cover above soil affects the soil heat transfer through thermal and hydrological insulation depending on the thickness and wetness of the moss.Also in reality, the moss distribution shows great spatial differences.This geographic dependence of moss cover brings additional heterogeneity to the soil thermal dynamics in the Arctic.To have the first step to represent such complexity, a constant uniform moss cover without the hydrological effects is assumed for the entire domain.This moss layer has similar functions as the snow layers, i.e., not having dynamic hydrology but rather providing constant heat insulation for the underlying soil layers.For the simulations presented in this paper, a 10 cm thick moss layer is chosen for all the seasons.The heat parameters for the moss layer follow Beringer et al. (2001), with heat conductivity of 0.25 WK −1 m −1 and volumetric heat capacity of 2 500 000 Jm −3 K −1 .

Global forcing data
For the period 1901-1978, daily forcing data with 0.5 • spatial resolution from the EU project WATCH (Water and Global Change) has been used (Weedon et al., 2010(Weedon et al., , 2011)).This data is based on ERA-40 (ECMWF 40 year Re-Analysis) reanalysis results that were bias-corrected by using several observation-based data sets, such as climate grids from the Climate Research Unit, University of East Anglia (CRU).For the 1979-2010 period, ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim reanalysis data (Dee et al., 2011) has been used.This data set was downloaded at 0.5 • spatial resolution and biascorrected against the WATCH-forcing data following Piani et al. (2010).A more detailed description of the climate forcing data set can be found in Beer et al. (2014).With this approach, a consistent time series of climate data for the period 1901-2010 is ensured.The sand, silt and clay fractions from the Harmonized World Soil Database v.1.1 (FAO et al., 2009) were the basis for deriving the soil thermal properties.Up to four tiles per 0.5 • grid cell area are distinguished for vegetation-related model parameters (Raddatz et al., 2007).The coverage of these tiles has been estimated by combining the GLC2000 land cover map (GLC2000 database, 2003), the MODIS (Moderate Resolution Imaging Spectroradiometer) Vegetation Continuous Fields product (Hansen et al., 2003) and the WWF (World Wildlife Fund) biome map (Olson et al., 2001).
JSBACH was forced by global atmospheric carbon dioxide concentrations following the CMIP5 (Coupled Model Intercomparison Project Phase 5) protocol (Meinshausen et al., 2011).

Simulation setup
Site-level simulations were performed running the model at a single point, forced by meteorological site observations (see below section).Soil parameters were extracted from the above-mentioned global land surface data and given in Table 1.Using the observed meteorological data, an average seasonal cycle was prepared and repeated for 30 years to force a spin-up simulation for bringing the soil thermal and hydrological profiles to equilibrium.Then, the transient simulation for the site was conducted using multiple years of observed climate and the results were used for comparison with the soil temperature observations.The time period used for the site simulations is from August 2008 to December 2009 for Nuuk, and from July 2003 to October 2005 for Samoylov.
For the circumpolar simulations, the model was run using the previously described global daily forcing data for the grids above 50 • north.First, the model's physical state was brought into equilibrium with a 30-year run repeating an average seasonal cycle of climate variables from the period 1901-1930.Then, a climate-transient run with constant atmospheric CO 2 concentration at the 1901 value was executed for the same period.These 30-year model results were further used to force a 1000-year carbon balance model run in order to prepare equilibrated carbon pools.Finally, these carbon pools are used as the initial condition to start a fully transient run from 1901 to 2010.

Nuuk-site observations
The Nuuk observational site is on the southwestern coast of Greenland, 250 km south of the polar circle at around 64 • north and 51 • west.It is situated in the Kobbefjord at an altitude of 500 m a.s.l.(above sea level) close to the city of Nuuk.Ambient climate is arctic/polar with mean annual temperature of −1.5 • C in 2008 and −1.3 • C in 2009 (Jensen andRasch, 2009, 2010).Vegetation type consists of Empetrum nigrum with Betula nana and Ledum groenlandicum with a vegetation height of 3-5 cm.The study site's soil lacks mineral soil horizons due to cryoturbation and lack of podsol development due to its dry location.Soil type is categorized as mostly sandy soil with 10 % organic matter in the top 10 cm, no ice lenses in the profile and no permafrost.No soil ice or permafrost formations have been observed within the drainage basin.Snow cover is measured at the Climate Basic station 1.65 km from the soil station but at the same altitude.At the time of the annual Nuuk basic snow survey in mid-April, the snow depth at the soil station is much alike the snow depth at the Climate Basic station: ±0.1 m when the snow depth is high (near 1 m) and much alike if it is much lower.Strong winds (> 20 m s −1 ) have a strong influence on the redistribution of newly fallen snow especially in the beginning of the snow season, so the formation of a permanent snow cover at the soil station can be delayed by as much as one week; while the end of the snow cover season is more or less alike the date at the Climate Basic station.In some winters there is some depth hoar formation in the snowpack (B.U. Hansen, personal communication, 2013).
The meteorological (half-hourly incoming radiation, air temperature, precipitation, wind speed) and soil observations (hourly soil temperature) were downloaded from the Greenland Ecosystem Monitoring database web server (Za-ckenbergGIS).For the meteorological variables, the time period used was July 2008 to December 2010, while the

Samoylov-site observations
The Samoylov field site is located in northern Siberia (72.4 • N, 126.5 • E) at the Lena River delta.The site represents a typical lowland tundra landscape and is characterized by continuous and ice-rich permafrost, which reaches depths of about 200 m (Grigoriev et al., 1996).The local climate is Arctic-continental with a mean annual air temperature of about −13 • C. The annual temperature range spans from about −45 • C in winter to 20 • C in summer (Boike et al., 2013).The total annual precipitation is about 200 mm, of which about 25 % contributes to snowfall (Boike et al., 2008;Langer et al., 2011).The snow cover is strongly characterized by wind drift and is usually very shallow with maximum depths of about 0.5 m (Boike et al., 2013).The land surface at the field site is dominated by polygonal tundra mainly vegetated by mosses and sedges (Kutzbach et al., 2004).The tundra soil consists of water-/ice-saturated sandy peat with the water table usually close to the surface (Langer et al., 2011).
The volumetric mineral content is reported to range between 20 and 40 % while the volumetric organic content is on the order of 5-10 % (Kutzbach et al., 2004;Langer et al., 2011).
The peat soil complex reaches depths of 10-15 m and is underlain by sandy-to-silty river deposits reaching depths of at least 1 km.
Hourly values of air temperature, precipitation (not in winter), wind speed and incoming longwave radiation is provided by the site measurements.Winter precipitation and incoming shortwave radiation are complemented by WATCH reanalysis data.Altogether a continuous model forcing data set is created.Minor gap filling was needed to fill in the missing data.The time period for the prepared data set is from July 2003 to October 2005.

Circumarctic data sets
The International Permafrost Association's (IPA) permafrost map (Brown et al., 2002) was used for comparing the simulated permafrost extent with the observations.Although the IPA map has distinct permafrost classes, only the outer border of the discontinuous and sporadic zones was considered when comparing with the model's permafrost extent, which is calculated using the simulated soil temperatures from the circumarctic model simulation.Following the permafrost definition of IPA (soils under freezing temperatures for at least 2 consecutive years), the permafrost state of each grid box is determined.The permafrost condition for each grid box was calculated with regards to the soil temperature only.For all of the five soil layers the temperatures are checked if any of the layers are staying below 0 • for at least 2 years.For comparing with the IPA map, the 1980-1990 average values of the model's permafrost state were used.
The Circumpolar Active Layer Monitoring network's (CALM) data set (Brown et al., 2000) was used for evaluating the simulated active-layer thickness.The CALM network maintains active layer thickness measurements at more than 200 sites since the 1990s.We have chosen the CALM sites within the continuous permafrost zone in our simulation domain and compared them with the corresponding 0.5 • × 0.5 • grid box of the results from the simulation conducted using global climate and soil texture data as forcing.Using a linear piecewise interpolation the simulated soil temperatures in five soil layers are interpolated into 200 evenly spaced nodes and the depth of 0 • is calculated afterwards to represent the thawing depth at each time step.Then the maximum thawing depth during the summer season is taken to be the active-layer depth for comparison.If there were more than one CALM site within one model grid box, the most appropriate one is chosen for the comparison.Averaging several CALM sites within one grid box is avoided since the average value could represent a nonrealistic condition due to surface heterogeneity.We tried to select the site that is most comparable with the model assumptions (e.g., upland soils) and the soil conditions represented by the global soil map.Since not all the sites had recorded measurements during the 1990-2010 period, we have averaged the existing years of data and compared it with the averages of corresponding years from the model output.
Numerous borehole observations from circumarctic stations were gathered during the International Polar Year (IPY 2007(IPY /2008)).They include deep and shallow borehole temperature observations representing the state of the permafrost (Romanovsky et al., 2010b).These borehole measurements are available through Global Terrestrial Network for Permafrost (GTN-P).Observations from these borehole measurements were compared with the simulated temperatures.As in the CALM comparison, the corresponding grid box values of the JSBACH simulation results were used for comparison.Since there were more boreholes in most of the grid boxes and surface heterogeneity has less effect on deep soil temperatures (7-10 m depth), we have performed a grid averaging to compare with the model outputs.The time period chosen for the comparison follows the IPY period: averaging years 2007 and 2008 outputs.

Continental-scale maps
The Russian permafrost temperature map (Land Resources of Russia CD-ROM, 2002) was prepared by the Russian Academy of Sciences and The International Institute for Applied Systems Analysis (IIASA).This map is an upscaled product of several meteorological and soil station data that are gathered during the expeditions in the second half of the 20th century.The data is digitally available (Land Resources of Russia CD-ROM, 2002) and downloadable from the web server of the IIASA (Land Resources of Russia).In the map, permafrost temperature is distinguished as 9 temperature classes and the temperature ranges (range of 1 or 2 • C) show a scale from 0 to −17 • C. To prepare a map comparable with JSBACH simulation results, the mean of the observed temperature classes were used to plot the observational map in this paper.Since there was no detailed information about the depths of these observations, values are assumed to be representative of those at depth of no seasonal temperature change.Following the observational time period, mean JSBACH subsoil temperature (last soil layer, ca.7 m) of the 1960-1990 period was used for comparison.
The 0.5 • active-layer thickness map (Beer et al., 2013) from Yakutia is an upscaled digitized version of the map of landscapes and permafrost conditions in Yakutia (Fedorov et al., 1989(Fedorov et al., , 1991)).Covering most of eastern Siberia, this map is very useful to understand the permafrost conditions at a 1 : 2 500 000 spatial scale during the period 1960-1987.Maps of mean and standard deviation of active-layer thickness were prepared at 0.5 • spatial resolution based on 0.001 • raster images.Active-layer thickness values range from 0.4 m at the northern continuous permafrost zone to 2.5 m at the southern borders of permafrost where isolated patches dominate the landscape.
For comparison with the active-layer thickness map, soil temperatures simulated at a 0.5 • spatial scale during the period 1960-1990 were used to derive the model's active-layer thickness, and then the mean of all these years is used to prepare the comparison map.Observed soil temperature at 30 cm is plotted with the black line and the red line shows the JSBACH-simulated soil temperature in the second layer (ca.18 cm).

Arctic river runoff data
There are several big rivers flowing into the Arctic Ocean from Russia, Canada and Alaska; and they are all affected by the conditions of permafrost underlying their respective basins.By comparing the temporal dynamics of runoff values at the river mouths, the model performance in representing the interactions between permafrost processes and the hydrological scheme can be assessed all around the basin areas.For testing model hydrological processes, runoff data from the Lena and Yenisey rivers were compared to simulation results.The runoff observations at the river mouth stations were gathered from the R-ArcticNET database (Lammers et al., 2001).The simulated runoff values in all the grid boxes within river basins were accumulated.For the evaluation of the seasonal cycle, simulation results were shifted by 2 months accounting for the time lag between the further grid cells and the river mouth station, for the reason that JSBACH does not include a river routing scheme.

Site-level validation
By forcing JSBACH with the meteorological data from the Nuuk synoptic station, a site-level simulation was performed.JSBACH successfully captured the topsoil temperature dynamics during the simulation period (Figs. 2 and 3).Following the observations, summer 2008 topsoil temperatures gradually cool down to 0 • C. Simulated temperatures fluctuate around 0 • C from October to February, in agreement with the observed data.After June, when the simulated temperatures are above zero, it takes until mid-July to capture the observations again during the summer of 2009.The second site's simulation is performed at the Samoylov site, where the conditions are representative of wet tundra permafrost.Except for the overestimations during early summer, JSBACH results successfully captured the summer temperature dynamics at different soil depths (Figs. 4 and  5).During autumn, simulated temperatures are falling down with a similar slope to the observations, while the warming up period (May-June) displays an underestimation compared to the observed values.However, winter temperatures are not simulated as cold as the observed values.The minimum value of the JSBACH winter temperatures are 10-15 • C warmer than the observations (Figs.4-6).As in the Nuuk comparison, the zero curtain is also seen at Samoylov.The timing of the freezing is also well represented by the model.Both observed and simulated temperatures are stagnating at around 0 • C during the freezing period of September-October (Fig. 5). Figure 6 shows the temperatures at a yearlong frozen depth, where the model comparison to the observed values show similar dynamics.A more detailed analysis of these comparisons requires mentioning the effects of freeze/thaw cycles.The latent heat released (consumed) when the soil water freezes (melts) is responsible for delaying the cold (heat) penetration into the soil.The site simulation results show that the topsoil temperatures are wavering around 0 • C during the phase change event.This so-called zero-curtain effect is also visible in the observational data (Figs.2-5).This match indicated that the phase change is accurately represented by the model.
It is seen from both site-level comparisons that winter soil temperatures do not drop as low as might be expected due to atmospheric conditions alone.Even when the air temperature is minimal in high winter (ca.−20 • C for Nuuk and −40 • C for Samoylov, not shown), soil keeps a rather warm temperature profile (ca.−3 • C for Nuuk and −25 • C for Samoylov) as long as snow exists on top.
However, in reality, snow has rather complicated characteristics.Within the snowpack, metamorphism processes create various types of snow with different thermal properties (Loth and Graf, 1993).When there is new snowfall, fresh snow presses down to squeeze the air out of deeper snow layers, thus increasing the snow density.With higher density, the snow insulation effect decreases due to increased snow heat conductivity.However, depending on site-specific conditions, springtime snow insulation can be altered due to the effects of depth hoar formation, wind drift or snowmelt water.Snow properties can also be modified by rainwater percolation into the snowpack.Also, snowmelt water infiltration into the soil can change the temperature profile of the soil.Additionally snow albedo changes with these processes.Boike et al. (2013) explained the strong wind conditions at Samoylov, where the maximum snow depth does not exceed 0.5 m.However this is not the case in JSBACH simulations, so there is an overestimation in simulated snow depths (Fig. 7).Such effects are still not represented in the current version of JSBACH and they can explain the mismatch in simulated versus observed springtime soil temperature in the Samoylov site-level simulations (Figs. 4 and 5).Without dynamically changing snow properties and lack of these snowspecific processes, our model cannot correctly represent the lower spring insulation and keeps a colder soil temperature profile.Similar effects were also shown by Westermann et al. (2013).Zhang et al. (2005) and Langer et al. (2013) pointed out the importance of correct parameterization of the snow thermal properties in permafrost simulations.Further progress in resolving these issues will be shown in the next model version.

Circumarctic validation
To evaluate the model's reliability at circumarctic scale, we compared the IPA permafrost map (Brown et al., 2002) with the simulated permafrost extent.Depending on the permafrost coverage, the IPA map classifies the permafrost zones as continuous, discontinuous, sporadic permafrost and isolated patches.However, within a global model, we do not represent such classification inside a grid cell, but rather classify permafrost or non-permafrost conditions.Having this in mind, it is seen in Fig. 8 that in general the simulated permafrost extent is in good agreement with the IPA map.It covers all the continuous and discontinuous zones and extends further to include some parts of the sporadic permafrost zone and isolated patches.By definition, sporadic permafrost has 10-50 % of permafrost coverage and isolated patches have less than 10 %.Simulating permafrost in some of these regions is assumed to be realistic when the binary criterion permafrost/no permafrost is used in the model.
Another criterion for assessing the validity of our simulation results is to evaluate active-layer thickness.By definition, active-layer thickness is the maximum thawing depth in permafrost areas during any given year.It can be considered a good measure of climate state since it is affected by summer temperature, precipitation, timing of snowmelt and history of soil temperature combined.For this reason, we have compared the current state  of the simulated active-layer thickness with the CALM network data.JSBACH matched the active-layer thickness of some of the sites better than the others but in general there is an overestimation in simulated active-layer thickness (Fig. 9).
Reasons for this mismatch are mostly explained by scale issues and site-specific conditions together with the model's vertical resolution.First, the model output from a 0.5 • × 0.5 • grid box cannot be taken as equally comparable to the site's observations given that the grid-box average is not fully representative for the heterogeneous surface conditions in this area.Even though some of the CALM observations were averaged over 1 km ×1 km areas, the landscape variability still brings up a big uncertainty when compared to a model gridbox average.It is important to note that for this comparison, the model was driven by global climate and soil properties data sets (Sect.2.2) and not by specific characteristics at CALM stations.Hence, part of the scatter in Fig. 9 can be explained by the wrong representation of soil properties or local climate conditions.Therefore, the overestimation of site-level active-layer thickness should be interpreted in concert with the comparisons of spatial details of active-layer thickness (ALT) and permafrost temperature (see next section).All things considered, site-level model estimates are fairly comparable to observations (r 2 : 0.54, Fig. 9).Similar results are observed in some other modeling studies.Dankers et al. (2011) have shown a deeper simulated ALT using the JULES (Joint UK Land Environment Simulator) model.Lawrence et al. (2012) have shown that the coupled and uncoupled CLM (Community Land Model) model runs are resulting in deeper ALT in general; although the offline run from the CLM4 model version showed a more distributed result with positive and negative differences.Additionally, it is explained by Gouttevin et al. (2012a) that the freezing scheme brought a better match with the CALM observations but still with a positive bias.
Complementary to CALM comparisons, borehole temperature records from GTN-P were used to evaluate simulated subsoil temperatures (last model layer, ca.7 m).In general, the model can explain about 48 % of observed subsoil temperature variation with a tendency to a cold bias at some sites (Fig. 10).This cold bias can partly be related to the model assumption of zero heat flux at the bottom of the soil.Previously shown by Lawrence and Slater (2005), the CLM3 model (with 3.43 m soil depth and vanishing heat flux at the bottom) simulated strong permafrost degradation by 2100.Delisle (2007) responds to that by showing the importance of including bottom energy flux of the permafrost layer in the model.Delisle (2007) also suggests the necessity of representing soil heat transfer by moving groundwater while Burn and Nelson (2006) explain the CLM3 overestimation of permafrost loss by using wrong surface temperatures and lack of near-surface ground ice in their model.In a newer model version (CLM4), Lawrence et al. (2012) explain the cold bias in deep soil temperatures with the dry active layers in their model, which again brings up the importance of hydrology-heat transfer interactions.Also a deeper soil column representing up to 50 m is suggested to improve the permafrost temperature results by around 10 m for future model versions; although the effects of having this deep soil column is not clear yet.Alexeev et al. (2007) suggests using at least a 30 m soil depth to capture the seasonal temperature variability.So it seems there are a few possible reasons for the cold bias in the JSBACH deep soil temperatures (Fig. 10).Nonetheless, the borehole temperature comparison illustrates the current status of our model in representing permafrost temperatures and shows the need for future model developments for improvement.
Model results showing deeper active layers seem to disagree with colder soil temperatures at first.However, the active-layer thickness is more related to topsoil temperature, whereas borehole comparisons were used to evaluate deeper layers.The topsoil is strongly coupled to atmospheric conditions and hydrological changes.However, deep soil temperature is less influenced by variable surface conditions, but show a decadal trend that is strongly affected by long-term atmospheric changes, snow depth and vegetation cover dynamics and the boundary conditions at the bottom of the soil column.As described in Dankers et al. (2011), active-layer comparisons are mostly affected by phase change events in the upper layers, but the colder soil temperature in the deeper layers is not strongly related to these phase change effects.Similar cold biases in deep soil temperature are also documented in other modeling studies (Gouttevin et al., 2012a;Lawrence et al., 2012).

Continental-scale validation
Spatial details of modeled permafrost temperature were compared to the Russian permafrost temperature map (Land Resources of Russia CD-ROM, 2002).The simulated latitudinal temperature gradient acts in accordance with the observation-based map albeit with regional underestimation of the model output (Fig. 11).Figure 12 shows the spatial pattern of this cold bias.In general, permafrost temperature differs from −2 to −5 • C, except in northern Yakutia where the  difference can be as great as −16 • C. A cold bias in subsoil temperature was also seen in the borehole temperature comparison, supporting the fact that it is not a regional issue but rather a global deficiency of the model or the global climate forcing data set.As discussed above, one potential reason for the colder soil temperature is the bottom boundary zero heat flux assumption.This assumption is widely used in the global modeling community (Dankers et al., 2011;Lawrence et al., 2008), but evidently the soil column depth also plays an important role (Alexeev et al., 2007).
It is also important to mention the higher spatial heterogeneity of JSBACH soil temperature when compared to the observation-based map.Since the observations were gathered very sparsely (due to harsh climate conditions and remote locations in Siberia) and widely interpolated to create such a large regional product, many features from landscape heterogeneity were lost in making the Russian permafrost temperature map.On the contrary, the model simulates each grid box individually by using meteorological forcing and surface conditions specific to each grid box.This explains the longitudinal changes in the model output (Fig. 11).Also, representing a different snow depth as well as not matching the distribution of moss cover affect the amount of heat insulation for the soil and alter the whole soil temperature profile.
Another regional evaluation performed was the comparison of observed and simulated active-layer thickness maps.Figure 13 shows the comparison of the active-layer thickness map of Yakutia (Beer et al., 2013) and the spatial distribution of active-layer thickness estimated by JSBACH.As in the permafrost temperature comparison (Fig. 11), a similar latitudinal gradient is observed in both maps.Although the observation-based map shows smaller values in the northern coastal regions, the transition of values from 50 cm at the coast to 250 cm further inland is comparable to the JS-BACH map.The mismatches at the coast can be due to the thick ice overburden in those areas, which are not represented by JSBACH.The differences between the observed and simulated results (Fig. 14) show a more diverse spatial pattern than the map of temperature differences (Fig. 12).This is due to the complex nature of confounding factors of active-layer thickness i.e., soil temperature, snow-moss cover and soil moisture.In general there is an overestimation in simulated active-layer thicknesses.As seen from the CALM comparison (Fig. 9), JSBACH simulates deeper active-layer depths.However, regional differences are apparent in this comparison.Unlike the CALM comparison where all the sites were overestimated, the blue regions in the difference map (Fig. 14) show the underestimated active-layer depths from model results.These mismatches can be attributed to not representing vegetation and snow cover properly.Finally, the southern borders of the comparison map (Fig. 14) illustrate a strong positive bias (ca.+250 cm).As explained in Beer et al. (2013), isolated permafrost patches are dominant in these areas.However, the comparison is not very useful for these areas since the 0.5 • values in the observation-based map represent an average of values in permafrost islands while the model is simulating a mean soil temperature profile for the 0.5 • grid cell from which active-layer thickness is estimated.Therefore, model results of ALT are expected to be higher in these areas.

River runoff validation
To evaluate the hydrological processes, Arctic river runoff dynamics were compared to the model results.The Lena River was chosen since it has one of the basin areas least influenced by anthropogenic activities and represents a more natural pattern that is easily comparable to the model results.The current model version can simulate the annual changes (Fig. 15) and the monthly dynamics (Fig. 16) of the Lena River runoff close to the observations.Permafrost conditions allow the soil to block water infiltration during the snowmelt period leading to a dramatic runoff peak in spring.JSBACH successfully captured these effects.Similar results have been observed in other studies (Beer et al., 2007;Gouttevin et al., 2012a).
In addition, the Yenisey River was chosen as a secondary basin since it has one of the biggest basins among the Arctic rivers.In general, this comparison is similar to the Lena Basin comparison.JSBACH underestimated the annual runoff values (Fig. 17) but matched the monthly dynamics (Fig. 18).The only issue here is the low values of the annual runoff.Simulating the Yenisey Basin has a higher uncertainty, since more landscape types are involved.Nevertheless, JSBACH captured the temporal dynamics of the Yenisey River runoff values, thus supporting the validity of the permafrost-hydrology interactions within the model.Interestingly, the model fails to reproduce the runoff increase  since 1983.This could be partly due to a global dimming effect on stomatal conductance, which influences transpiration (Gedney et al., 2006).However, other effects, such as snowmelt dynamics have an impact as well.

Conclusions
In this paper, we have presented an advanced version of the process-oriented ecosystem model JSBACH that simulates cold regions through enhanced representation of snow and soil physics.By including the phase-change process, coupled thermal and hydrological processes and heat insulation from snow and moss cover, the current model version is a capable tool for simulating the physical state of high-latitude terrestrial environments.A multiscale evaluation was conducted and the results demonstrate the strength and weaknesses of the model.Site-level comparisons at both permafrost and non-permafrost sites indicate the importance of freezing and thawing together with snow insulation for  representing soil temperature dynamics.On the larger scale, permafrost extent is successfully reproduced.Comparisons with circumarctic observational data sets revealed that the model simulates deeper active-layer thicknesses accompanied with colder subsoil temperatures.These issues are explained by the differences in snow cover and moss-layer distributions that are not captured by the model, shallow depth of the soil column and the vanishing heat flux assumption at the bottom.Additionally, regional comparisons drew attention to the heterogeneous vegetation cover and the influence of topographic effects.In conclusion, this modeling study highlights the importance of the effects of latent heat and insulation from snow/moss cover in simulating the permafrost state in high-latitude regions.

Fig. 1 .
Fig. 1.Vertical soil model structure of the new JSBACH version.The numbers left of the soil column show the depths of the bottom of each layer while the numbers right of the soil column show layer thicknesses.Here snow and organic layers are simply shown to represent a multilayered snow scheme and constant moss layer described in the text.

Fig. 2 .
Fig. 2. Observed and simulated upper-layer soil temperature at the Nuuk site.Observed soil temperature at 5 cm is plotted with the black line and the red line shows the JSBACH-simulated soil temperature in the first layer (ca. 3 cm).

Fig. 3 .
Fig. 3. Observed and simulated soil temperature at the Nuuk site.Observed soil temperature at 30 cm is plotted with the black line and the red line shows the JSBACH-simulated soil temperature in the second layer (ca.18 cm).

Fig. 4 .
Fig. 4. Observed and simulated upper-layer soil temperature at the Samoylov site.Observed soil temperature at 6 cm is plotted with the black line and the red line shows the JSBACH-simulated soil temperature in the first layer (ca. 3 cm).

Fig. 5 .
Fig. 5. Observed and simulated soil temperature at the Samoylov site.Observed soil temperature at 21 cm is plotted with the black line and the red line shows the JSBACH-simulated soil temperature in the second layer (ca.18 cm).

Fig. 6 .
Fig. 6.Observed and simulated soil temperature at the Samoylov site.Observed soil temperature at 71 cm is plotted with the black line and the red line shows the JSBACH-simulated soil temperature in the third layer (ca.78 cm).

Fig. 7 .
Fig. 7. Observed and simulated snow depth at the Samoylov site.

Fig. 8 .
Fig. 8. Northern Hemisphere permafrost extent according to the International Permafrost Association's permafrost map (Brown et al., 2002).Different permafrost classes are plotted in different colors and the red line shows the border of the permafrost extent calculated from the JSBACH simulation (1980-1990 average values).

Fig. 9 .
Fig. 9. Scatter plot of the observed ALT from the CALM network(Brown et al., 2000) versus the JSBACH results.See text for further info.

Fig. 10 .
Fig. 10.Scatter plot of observed soil temperature from the GTN-P borehole temperature data set (Romanovsky et al., 2010b) versus simulated subsoil temperature (deepest soil layer, ca.7 m).See text for further info.

Fig. 11 .
Fig. 11.Comparison of Russian permafrost temperature.Observed (map a; see text for more details) (Land Resources of Russia CD-ROM, 2002) and simulated (map b) Russian permafrost temperature during the period 1960-1990.The average values in different temperature classes are plotted with the same color in both maps.

Fig. 13 .
Fig. 13.Observed (map a; see text for more details) (Beer et al., 2013) and simulated active-layer thickness (map b) in the Yakutsk area.

Fig. 15 .
Fig. 15.Simulated and observed annual Lena River runoff.Red line represents the JSBACH model version with the permafrost representation.The black line shows the observed values from the R-ArcticNet database (Lammers et al., 2001).

Fig. 16 .
Fig. 16.Simulated and observed monthly mean Lena River runoff.Line colors are the same as annual runoff plot (red: model values; black: observed values).Since the model does not use a river routing scheme, the model results are shifted 2 months to match the actual peak time in spring.

Fig. 17 .
Fig. 17.Simulated and observed annual Yenisey River runoff.Red line represents the JSBACH model version with the permafrost representation.The black line shows the observed values from the R-ArcticNet database (Lammers et al., 2001).

Fig. 18 .
Fig. 18.Simulated and observed monthly mean Yenisey River runoff.Line colors are as in Fig. 17.As for the simulated Lena River runoff in Fig. 16, the model results are shifted 2 months to match the actual peak time in spring due to the lack of river routing scheme in the model.

Table 1 .
JSBACH model parameters used in the site simulations.

www.geosci-model-dev.net/7/631/2014/ Geosci. Model Dev., 7, 631-647, 2014 soil
temperature was available from August 2008 to December 2009.The downloaded ASCII (American Standard Code for Information Interchange) files have been combined in a netCDF (network Common Data Form) format file; minorgap filling was needed to create a continuous climate forcing to force the Nuuk site-level simulations.