Site-level model intercomparison of high latitude and high altitude soil thermal dynamics in tundra and barren landscapes

. Modeling soil thermal dynamics at high latitudes and altitudes requires representations of physical processes such as snow insulation, soil freezing and thawing and subsurface conditions like soil water/ice content and soil texture. We have compared six different land models: JSBACH, ORCHIDEE, JULES, COUP, HYBRID8 and LPJ-GUESS, at four different sites with distinct cold region landscape types, to identify the importance of physical processes in capturing observed temperature dynamics in soils. The sites include alpine, high Arctic, wet polygonal tundra and non-permafrost Arctic, thus showing how a range of models can represent distinct soil temperature regimes. For all sites, snow insulation is of major importance for estimating topsoil conditions. However, soil physics is essential for the subsoil temperature dynamics and thus the active layer thicknesses. This analysis shows that land models need more realistic surface processes, such as detailed snow dynamics and moss cover with changing thickness and wetness, along with better representations of subsoil thermal dynamics.


Introduction
Recent atmospheric warming trends are affecting terrestrial systems by increasing soil temperatures and causing changes in the hydrological cycle. Especially in high latitudes and altitudes, clear signs of change have been observed (Serreze et al., 2000;ACIA, 2005;IPCC AR5, 2013). These relatively colder regions are characterized by the frozen state of terrestrial water, which brings additional risks associated with shifting soils into an unfrozen state. Such changes will have broad implications for the physical (Romanovsky et al., 2010), biogeochemical (Schuur et al., 2008) and structural (Larsen et al., 2008) conditions of the local, regional and global climate system. Therefore, predicting the future state of the soil thermal regime at high latitudes and altitudes holds major importance for Earth system modeling.
There are increasing concerns as to how land models perform at capturing high latitude soil thermal dynamics, in particular in permafrost regions. Recent studies (Koven et al., 2013;Slater and Lawrence, 2013) have provided detailed assessments of commonly used earth system models (ESMs) in simulating soil temperatures of present and future state Published by Copernicus Publications on behalf of the European Geosciences Union. of the Arctic. By using the Coupled Model Intercomparison Project phase 5 -CMIP5 (Taylor et al., 2009) results, Koven et al. (2013) have shown a broad range of model outputs in simulated soil temperature. They attributed most of the inter-model discrepancies to air-land surface coupling and snow representations in the models. Similar to those findings, Slater and Lawrence (2013) confirmed the high uncertainty of CMIP5 models in predicting the permafrost state and its future trajectories. They concluded that these model versions are not appropriate for such experiments, since they lack critical processes for cold region soils. Snow insulation, land model physics and vertical model resolutions were identified as the major sources of uncertainty.
For the cold regions, one of the most important factors modifying soil temperature range is the surface snow cover. As discussed in many previous studies (Zhang, 2005;Koven et al., 2013;Scherler et al., 2013;Marmy et al., 2013;Langer et al., 2013;Boike et al., 2003;Gubler et al., 2013;Fiddes et al., 2015), snow dynamics are quite complex and the insulation effects of snow can be extremely important for the soil thermal regime. Model representations of snow cover are lacking many fine-scale processes such as snow ablation, depth hoar formation, snow metamorphism, wind effects on snow distribution and explicit heat and water transfer within snow layers. These issues bring additional uncertainties to global projections.
Current land surface schemes, and most vegetation and soil models, represent energy and mass exchange between the land surface and atmosphere in one dimension. Using a grid cell approach, such exchanges are estimated for the entire land surface or specific regions. However, comparing simulated and observed time series of states or fluxes at point scale rather than grid averaging is an important component of model evaluation, for understanding remaining limitations of models (Ekici et al., 2014;Mahecha et al., 2010). In such "site-level runs", we assume that lateral processes can be ignored and that the ground thermal dynamics are mainly controlled by vertical processes. Then, models are driven by observed climate and variables of interest can be compared to observations at different temporal scales. Even though such idealized field conditions never exist, a careful interpretation of site-level runs can identify major gaps in process representations in models.
In recent years, land models have improved their representations of the soil physical environment in cold regions. Model enhancements include the addition of soil freezing and thawing, detailed snow representations, prescribed moss cover, extended soil columns and coupling of soil heat transfer with hydrology (Ekici et al., 2014;Gouttevin et al., 2012a;Dankers et al., 2011;Lawrence et al., 2008;Wania et al., 2009a). Also active layer thickness (ALT) estimates have improved in the current model versions. Simple relationships between surface temperature and ALT were used in the early modeling studies (Lunardini, 1981;Kudryavtsev et al., 1974;Romanovsky and Osterkamp, 1997;Shiklo-manov and Nelson, 1999;Stendel et al., 2007, Anisimov et al., 1997. These approaches assume an equilibrium condition, whereas a transient numerical method is better suited within a climate change context. A good review of widely used analytical approximations and differences to numerical approaches is given by Riseborough et al. (2008). With the advanced soil physics in many models, these transient approaches are more widely used, especially in long-term simulations. Such improvements highlight the need for an updated assessment of model performances in representing high latitude/altitude soil thermal dynamics.
We have compared the performances of six different land models in simulating soil thermal dynamics at four contrasting sites. In contrast to previous work (Koven et al., 2013;Slater and Lawrence, 2013), we used advanced model versions specifically improved for cold regions and our model simulations are driven by (and evaluated with) site observations. To represent a wider range of assessment and model structures, we used both land components of ESMs (JSBACH, ORCHIDEE, JULES) and stand-alone models (COUP, HYBRID8, LPJ-GUESS), and compared them at Arctic permafrost, Alpine permafrost and Arctic nonpermafrost sites. By doing so, we aimed to quantify the importance of different processes, to determine the general shortcomings of current model versions and finally to highlight the key processes for future model developments.  (Stevens et al., 2012) and MPIOM for the ocean (Jungclaus et al., 2013). JSBACH provides the land surface boundary for the atmosphere in coupled simulations; however, it can also be used offline driven by atmospheric forcing. The current version of JSBACH (Ekici et al., 2014) employs soil heat transfer coupled to hydrology with freezing and thawing processes included. The soil model is discretized as five layers with increasing thicknesses of up to 10 m depth. There are up to five snow layers with constant density and heat transfer parameters. JSBACH also simulates a simple moss/organic matter insulation layer again with constant parameters.

ORCHIDEE
ORCHIDEE is a global land surface model, which can be used coupled to the Institut Pierre Simon Laplace (IPSL) climate model or driven offline by prescribed atmospheric forcing (Krinner et al., 2005). ORCHIDEE computes all the  (Wang et al., 2013), soil freeze-thaw processes (Gouttevin et al., 2012a), and a vertical permafrost soil thermal and carbon module (Koven et al., 2011). The soil hydrology is vertically discretized as 11 numerical nodes with 2 m depth (Gouttevin et al., 2012a), and soil thermal and carbon modules are vertically discretized as 32 layers with ∼ 47 m depth (Koven et al., 2011). A onedimensional Fourier equation was applied to calculate soil thermal dynamics, and both soil thermal conductivity and heat capacity are functions of the frozen and unfrozen soil water content and of dry and saturated soil thermal properties (Gouttevin et al., 2012b).

JULES
JULES (Joint UK Land Environment Simulator) is the landsurface scheme used in the Hadley Centre climate model Clark et al., 2011), which can also be run offline, driven by atmospheric forcing data. It is based on the Met Office Surface Exchange Scheme, MOSES (Cox et al., 1999). JULES simulates surface exchange, vegetation dynamics and soil physical processes. It can be run at a single point, or as a set of points representing a 2-D grid. In each grid cell, the surface is tiled into different surface types, and the soil is treated as a single column, discretized vertically into layers (four in the standard setup). JULES simulates fluxes of moisture and energy between the atmosphere, surface and soil, and the soil freezing and thawing. It includes a carbon cycle that can simulate carbon exchange between the atmosphere, vegetation and soil. It also includes a multi-layer snow model , with layers that have variable thickness, density and thermal properties. The snow scheme significantly improves the soil thermal regime in comparison with the old, single-layer scheme (Burke et al., 2013). The model can be run with a time step of between 30 min and 3 h, depending on user preference.

COUP
COUP is a stand-alone, one-dimensional heat and mass transfer model for the soil-snow-atmosphere system (Jansson and Karlberg, 2011) and is capable of simulating transient hydrothermal processes in the subsurface including seasonal or perennial frozen ground (see e.g., Hollesen et al., 2011;Scherler et al., 2010Scherler et al., , 2013. Two coupled partial differential equations for water and heat flow are the core of the COUP Model. They are calculated over up to 50 vertical layers of arbitrary depth. Processes that are important for permafrost simulations, such as the freezing and thawing of soil as well as the accumulation, metamorphosis and melt of snow cover are included in the model (Lundin, 1990;Gustafsson et al., 2001). Freezing processes in the soil are based on a function of freezing point depression and on an analogy of freezing-thawing and wetting-drying (Harlan, 1973;Jansson and Karlberg, 2011). Snow cover is simulated as one layer of variable height, density, and water content. The upper boundary condition is given by a surface energy balance at the soil-snow-atmosphere boundary layer, driven by climatic variables. The lower boundary condition at the bottom of the soil column is usually given by the geothermal heat flux (or zero heat flux) and a seepage flow of percolating water. Water transfer in the soil depends on texture, porosity, water, and ice content. Bypass flow through macropores, lateral runoff and rapid lateral drainage due to steep terrain can also be considered (e.g., Scherler et al., 2013). A detailed description of the model including all its equations and parameters is given in Jansson andKarlberg (2011) andJansson (2012).

HYBRID8
HYBRID8 is a stand-alone land surface model, which computes the carbon and water cycling within the biosphere and between the biosphere and atmosphere. It is driven by the daily/sub-daily climate variables above the canopy, and the atmospheric CO 2 concentration. Computations are performed on a 30 min time step for the energy fluxes and exchanges of carbon and water with the atmosphere and the soil. Litter production and soil decomposition are calculated at a daily time step. HYBRID8 uses the surface physics and the latest parameterization of turbulent surface fluxes from the GISS ModelE (Schmidt et al., 2006;Friend and Kiang, 2005), but has no representation of vegetation dynamics. The snow dynamics from ModelE are also not yet fully incorporated. Heat dynamics are described in Rosenzweig et al. (1997) and moisture dynamics in Abramopoulos et al. (1988).
In HYBRID8 the prognostic variable for the heat transfer is the heat in the different soil layers, and from that, the model evaluates the soil temperature. The processes governing this are diffusion from the surface to the sub-surface layers, and conduction and advection between the soil layers. The bottom boundary layer in HYBRID8 is impermeable, resulting in zero heat flux from the soil layers below. The version used in this project has no representation of the snow dynamics and has no insulating vegetation cover. However, the canopy provides a simple heat buffer due its separate heat capacity calculations.

LPJ-GUESS
Lund-Potsdam-Jena General Ecosystem Simulator (LPJ-GUESS) is a process-based model of vegetation dynamics and biogeochemistry optimized for regional and global applications (Smith et al., 2001). Mechanistic representations www.the-cryosphere.net/9/1343/2015/ The Cryosphere, 9, 1343-1361, 2015 of biophysical and biogeochemical processes are shared with those in the Lund-Potsdam-Jena dynamic global vegetation model LPJ-DGVM (Sitch et al., 2003;Gerten et al., 2004). However, LPJ-GUESS replaces the large area parameterization scheme in LPJ-DGVM, whereby vegetation is averaged out over a larger area, allowing several state variables to be calculated in a simpler and faster manner, with more robust and mechanistic schemes of individual-and patch-based resource competition and woody plant population dynamics. Detailed descriptions are given by Smith et al. (2001), Sitch et al. (2003), Wolf et al. (2008), Miller and Smith (2012) and Zhang et al. (2013). LPJ-GUESS has recently been updated to simulate Arctic upland and peatland ecosystems (McGuire et al., 2012;Zhang et al., 2013). It shares the numerical soil thawingfreezing processes, peatland hydrology and the model of wetland methane emission with LPJ-DGVM WHyMe, as described by Wania et al. (2009aWania et al. ( , b, 2010. To simulate soil temperatures and active layer depths, the soil column in LPJ-GUESS is divided into a single snow layer of fixed density and variable thickness, a litter layer of fixed thickness (10 cm for these simulations, except for Schilthorn where it is set to 2.5 cm), a soil column of 2 m depth (with sublayers of thickness 0.1 m, each with a prescribed fraction of mineral and organic material, but with fractions of soil water and air that are updated daily), and finally a "padding" column of depth 48 m (with thicker sublayers), to simulate soil thermal dynamics. Insulation effects of snow, phase changes in soil water, daily precipitation input and air temperature forcing are important determinants of daily soil temperature dynamics at different sublayers.

Nuuk
The Nuuk observational site is located in southwestern Greenland. The site is situated in a valley in Kobbefjord at 500 m altitude above sea level, and ambient conditions show Arctic climate properties, with a mean annual temperature of −1.5 • C in 2008 and −1.3 • C in 2009 Rasch, 2009, 2010). Vegetation types consist of Empetrum nigrum with Betula nana and Ledum groenlandicum, with a vegetation height of 3-5 cm. The study site soil lacks mineral soil horizons due to cryoturbation and lack of podsol development, as it is situated in a dry location. The soil is composed of 43 % sand, 34 % loam, 13 % clay and 10 % organic materials. 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 was very similar to the snow depth at the Climate Basic station: ± 0.1 m when the snow depth is high (near 1 m). Strong winds (> 20 m s −1 ) have a strong influence on the redistribu-tion 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 as much as 1 week, while the end of the snow cover season is similar to that at the Climate Basic station (B. U. Hansen, personal communication, 2013; ZackenbergGIS, 2012).

Schilthorn
The Schilthorn massif (Bernese Alps, Switzerland) is situated at 2970 m altitude in the northcentral part of the European Alps. Its non-vegetated lithology is dominated by deeply weathered limestone schists, forming a surface layer of mainly sandy and gravelly debris up to 5 m thick, which lies over presumably strongly jointed bedrock. Following the first indications of permafrost (ice lenses) during the construction of the summit station between 1965 and 1967, the site was chosen for long-term permafrost observation within the framework of the European PACE project and consequently integrated into the Swiss permafrost monitoring network PERMOS as one of its reference sites (PERMOS, 2013).
The measurements at the monitoring station at 2900 m altitude are located on a flat plateau on the north-facing slope and comprise a meteorological station and three boreholes (14 m vertical, 100 m vertical and 100 m inclined), with continuous ground temperature measurements since 1999 (Vonder Mühll et al., 2000;Hoelzle and Gruber, 2008;Harris et al., 2009). Borehole data indicate permafrost of at least 100 m thickness, which is characterized by ice-poor conditions close to the melting point. Maximum active-layer depths recorded since the start of measurements in 1999 are generally around 4-6 m, but during the exceptionally warm summer of the year 2003 the active-layer depth increased to 8.6 m, reflecting the potential for degradation of permafrost at this site .
The monitoring station has been complemented by soil moisture measurements since 2007 and geophysical (mainly geoelectrical) monitoring since 1999 (Hauck, 2002;Hilbich et al., 2011). The snow cover at Schilthorn can reach maximum depths of about 2-3 m and usually lasts from October through to June/July. One-dimensional soil model sensitivity studies showed that impacts of long-term atmospheric changes would be strongest in summer and autumn, due to this late snowmelt and the long decoupling of the atmosphere from the surface. So, increasing air temperatures could lead to a severe increase in active-layer thickness (Engelhardt et al., 2010;Marmy et al., 2013;Scherler et al., 2013).

Samoylov
Samoylov Island belongs to an alluvial river terrace of the Lena River delta. The island is elevated about 20 m above the normal river water level and covers an area of about 3.4 km 2 . The western part of the island constitutes The Cryosphere, 9, 1343-1361, 2015 www.the-cryosphere.net/9/1343/2015/ 900 mm 1963 mm 200 mm 400 mm Avg. length of snow cover 7 months 9.5 months 9 months 9 months Vegetation cover Tundra Barren Tundra Tundra a modern floodplain, which is lower compared with the rest of the island and is often flooded during ice break-up of the Lena River in spring. The eastern part of the island belongs to the elevated river terrace, which is mainly characterized by moss-and sedge-vegetated tundra (Kutzbach et al., 2007). In addition, several lakes and ponds occur, which make up about 25 % of the surface area of Samoylov (Muster et al., 2012). The land surface of the island is characterized by the typical micro-relief of polygonal patterned ground, caused by frost cracking and subsequent ice-wedge formation. The polygonal structures usually consist of depressed centers surrounded by elevated rims, which can be found in a partly or completely collapsed state (Kutzbach et al., 2007). The soil in the polygonal centers usually consists of water-saturated sandy peat, with the water table standing a few centimeters above or below the surface. The elevated rims are usually covered with a dry moss layer, underlain by wet sandy soils, with massive ice wedges underneath. The cryogenic soil complex of the river terrace reaches depths of 10 to 15 m and is underlain by sandy to silty river deposits. These river deposits reach depths of at least 1 km in the delta region .
There are strong spatial differences in surface energy balance due to heterogeneous surface and subsurface properties. Due to thermo-erosion, there is an ongoing expansion of thermokarst lakes and small ponds (Abnizova et al., 2012). Soil water drainage is strongly related to active layer dynamics, with lateral water flow occurring from late sum-mer to autumn (Helbig et al., 2013). Site conditions include strong snow-microtopography and snow-vegetation interactions due to wind drift .

Bayelva
The Bayelva climate and soil-monitoring site is located in the Kongsfjord region on the west coast of Svalbard Island. The North Atlantic Current warms this area to an average air temperature of about −13 • C in January and +5 • C in July, and provides about 400 mm precipitation annually, falling mostly as snow between September and May. The annual mean temperature of 1994 to 2010 in the village of Ny-Ålesund has been increasing by +1.3 K per decade (Maturilli et al., 2013). The observation site is located in the Bayelva River catchment on the Brøgger peninsula, about 3 km from Ny-Ålesund. The Bayelva catchment is bordered by two mountains, the Zeppelinfjellet and the Scheteligfjellet, between which the glacial Bayelva River originates from the two branches of the Brøggerbreen glacier moraine rubble. To the north of the study site, the terrain flattens, and after about 1 km, the Bayelva River reaches the shoreline of the Kongsfjorden (Arctic Ocean). In the catchment area, sparse vegetation alternates with exposed soil and sand and rock fields. Typical permafrost features, such as mud boils and non-sorted circles, are found in many parts of the study area. The Bayelva permafrost site itself is located at 25 m a.s.l., on top of the small Leirhaugen hill. The domiwww.the-cryosphere.net/9/1343/2015/ The Cryosphere, 9, 1343-1361, 2015  nant ground pattern at the study site consists of non-sorted soil circles. The bare soil circle centers are about 1 m in diameter and are surrounded by a vegetated rim, consisting of a mixture of low vascular plants of different species of grass and sedges (Carex spec., Deschampsia spec., Eriophorum spec., Festuca spec., Luzula spec.), catchfly, saxifrage, willow and some other local common species (Dryas octopetala, Oxyria digyna, Polegonum viviparum) and unclassified species of mosses and lichens. The vegetation cover at the measurement site was estimated to be approximately 60 %, with the remainder being bare soil with a small proportion of stones. The silty clay soil has a high mineral content, while the organic content is low, with organic fractions below 10 % . In the study period, the permafrost at Leirhaugen hill had a mean annual temperature of about −2 • C at the top of the permafrost at 1.5 m depth. Over the past decade, the Bayelva catchment has been the focus of intensive investigations into soil and permafrost conditions (Roth and Boike, 2001;Boike et al., 2007;Westermann et al., 2010Westermann et al., , 2011, the winter surface energy balance (Boike et al., 2003), and the annual balance of energy, H 2 O and CO 2 , and micrometeorological processes controlling these fluxes (Westermann et al., 2009;Lüers et al., 2014).

Intercomparison setup and simulation protocol
In order solely to compare model representations of physical processes and to eliminate any other source of uncertainty (e.g., climate forcing, spatial resolution, soil parameters), model simulations were driven by the same atmospheric forcing and soil properties at site-scale. Driving data for all site simulations were prepared and distributed uniformly. Site observations were converted into continuous time series with minor gap-filling. Where the observed variable set lacked the variable needed by the models, extended WATCH reanalysis data (Weedon et al., 2010;Beer et al., 2014) were used to complement the data sets. Soil thermal properties are based on the sand, silt, and clay fractions of the Harmonized World Soil database v1.1 (FAO et al., 2009). All model simulations were forced with these data sets. Table 3 summarizes the details of site-driving data preparation together with soil static parameters.
To bring the state variables into equilibrium with climate, models are spun up with climate forcing. Spin-up procedure is part of the model structure, in some cases a full biogeochemical and physical spin up is implemented, whereas in some models a simpler physical spin up is possible. This brings different requirements for the spin-up time length, so each model was independently spun up depending on its model formulations and discretization scheme, and the details are given in Table 4 Table A1 in Appendix A for exact soil depths used in this plot.
all model spin-up procedures was to keep the mean annual soil temperature change less than 0.01 • C in all soil layers.
Most of the analysis focuses on the upper part of the soil. The term "topsoil" is used from now on to indicate the chosen upper soil layer in each model, and the first depth of soil temperature observations. The details of layer selection are given in Table A1 of Appendix A.  Table A1 in Appendix A for exact soil depths used in this plot.

Topsoil temperature and surface insulation effects
As all our study sites are located in cold climate zones (Fig. 1), there is significant seasonality, which necessitates a separate analysis for each season. Figure 2 shows average seasonal topsoil temperature distributions (see Table A1 in the Appendix for layer depths) extracted from the six models, along with the observed values at the four different sites. In this figure ing winter and spring, Samoylov is much colder than the other sites (Fig. 2a, b). Observed summer and autumn temperatures are similar at all sites (Fig. 2c, d), with Nuuk being the warmest site in general. For the modeled values, the greatest inconsistency with observations is in matching the observed winter temperatures, especially at Samoylov and Schilthorn (Fig. 2a). The modeled temperature range increases in spring (Fig. 2b), and even though the mean modeled temperatures in summer are closer to observed means, the maximum and minimum values show a wide range during this season (Fig. 2c). Autumn shows a more uniform distribution of modeled temperatures compared with the other seasons (Fig. 2d). A proper assessment of critical processes entails examining seasonal changes in surface cover and the consequent insulation effects for the topsoil temperature. To investigate these effects, Fig. 3 shows the seasonal relations between air and topsoil temperature at each study site. Air temperature values are the same for all models, as they are driven with the same atmospheric forcing. Observations show that topsoil temperatures are warmer than the air during autumn, winter, and spring at all sites, but the summer conditions are dependent on the site (Fig. 3). In the models, winter topsoil temperatures are warmer than the air in most cases, as observed. However, the models show a wide range of values, especially at Samoylov (Fig. 3c), where the topsoil temperatures differ by up to 25 • C between models. In summer, the models do not show consistent relationships between soil and air temperatures, and the model range is highest at the Nuuk and Schilthorn sites.
To analyze the difference in modeled and observed snow isolation effect in more detail, Fig. 4 shows the changes in snow depth from observed and modeled values. Schilthorn has the highest snow depth values (> 1.5 m), while all other sites have a maximum snow height between 0.5 and 1 m (Fig. 4). Compared with observations, the models usually overestimate the snow depth at Schilthorn and Samoylov (Fig. 4b, c) and underestimate it at Nuuk and Bayelva (Fig. 4a, d).
For our study sites, the amount of modeled snow depth bias is correlated with the amount of modeled topsoil temperature bias (Fig. 5). With overestimated (underestimated) snow depth, models generally simulate warmer (colder) topsoil temperatures. As seen in Fig. 5a, almost all models underestimate the snow depth at Nuuk and Bayelva, and this creates colder topsoil temperatures. The opposite is seen for Samoylov and Schilthorn, where higher snow depth bias is accompanied by higher topsoil temperature bias (except for ORCHIDEE and LPJ-GUESS models).
As snow can be persistent over spring and summer seasons in cold regions (Fig. 4), it is worthwhile to separate snow and snow-free seasons for these comparisons. Figure 6 shows the same atmosphere-topsoil temperature comparison as in Fig. 3 but using individual (for each model and site) snow and snow-free seasons instead of conventional seasons. In this figure, all site observations show a warmer topsoil temperature than air, except for the snow-free season at Samoylov. Models, however, show different patterns at each site. For the snow season, models underestimate the observed values at Nuuk and Bayelva, whereas they overestimate it at Schilthorn and Samoylov, except for the previously mentioned ORCHIDEE and LPJ-GUESS models.  Table A1 in Appendix A for exact soil depths used in this plot.

Figure 6.
Scatter plots showing air-topsoil temperature relation from observations and models at each site for snow and snow-free seasons. Snow season is defined separately for observations and each model, by taking snow depth values over 5 cm to represent the snow-covered period. The average temperature of all snow covered (or snow-free) days of the simulation period is used in the plots. Markers distinguish snow and snow-free seasons and colors distinguish models. Gray lines represent the 1 : 1 line. See Table A1 in Appendix A for exact soil depths used in this plot.

Subsurface thermal regime
Assessing soil thermal dynamics necessitates scrutinizing subsoil temperature dynamics as well as surface conditions. Soil temperature evolutions of simulated soil layers are plotted for each model at each site in Figs. 7-10. Strong seasonal temperature changes are observed close to the surface, whereas temperature amplitudes are reduced in deeper layers and eventually a constant temperature is simulated at depths with zero annual amplitude (DZAA). Although Nuuk is a non-permafrost site, most of the models simulate subzero temperatures below 2-3 m at this site (Fig. 7). Here, only ORCHIDEE and COUP simulate a true DZAA at around 2.5-3 m, while all other models show a minor temperature change even at their deepest layers. At the high altitude Schilthorn site (Fig. 8), JSBACH and JULES simulate above 0 • C temperatures (non-permafrost conditions) in deeper layers. Compared with other models with snow representation, ORCHIDEE and LPJ-GUESS show colder subsurface temperatures at this site (Fig. 8). The simulated soil thermal regime at Samoylov reflects the colder climate at this site. All models show subzero temperatures below 1 m (Fig. 9). However, compared with other models, JULES and COUP show values much closer to 0 • C. At the high-Arctic Bayelva site, all models simulate permafrost conditions (Fig. 10). The JULES and COUP models again show warmer temperature profiles than the other models.
The soil thermal regime can also be investigated by studying the vertical temperature profiles regarding the annual means (Fig. 11) ( Fig. 12). In Figure 11, the distribution of mean values is similar to the analysis of topsoil conditions. The mean subsoil temperature is coldest at Samoylov followed by Bayelva, while Schilthorn is almost at the 0 • C boundary (no deep soil temperature data were available from Nuuk for this comparison). JSBACH, JULES, and COUP overestimate the temperatures at Schilthorn and Samoylov, but almost all models underestimate it at Bayelva. Figure 12 shows the temperature envelopes of observed and simulated values at each site. The minimum (maximum) temperature curve represents the coldest (warmest) possible conditions for the soil thermal regime at a certain depth. The models agree more on the maximum curve than the minimum curve (Fig. 12), indicating the differences in soil temperature simulation for colder periods.
The HYBRID8 model almost always shows the coldest conditions, whereas the pattern of the other models changes depending on the site. Figure 13 shows the yearly change of ALT for the three permafrost sites. Observations indicate a shallow ALT at Samoylov (Fig. 13b) and very deep ALT for Schilthorn (Fig. 13a). All models overestimate the ALT at Samoylov (Fig. 13b), but there is disagreement among models in overor underestimating the ALT at Schilthorn (Fig. 13a) and Bayelva (Fig. 13c).  Figure 2 has shown a large range among modeled temperature values, especially during winter and spring. As mentioned in the introduction, modeled mean soil temperatures are strongly related to the atmosphere-surface thermal connection, which is strongly influenced by snow cover and its properties.

Topsoil temperature and surface insulation effects
Observations show warmer topsoil temperatures than air during autumn, winter, and spring (Fig. 3). This situation indicates that soil is insulated when compared to colder air temperatures. This can be attributed to the snow cover during these seasons (Fig. 4). The insulating property of snow keeps the soil warmer than air, while not having snow can result in colder topsoil temperatures than air (as for the HYBRID8 model, cf. Fig. 3). Even though the high albedo of snow provides a cooling effect for soil, the warming due to insulation dominates during most of the year. Depending on their snow depth bias, models show different relations between air and topsoil temperature. The amount of winter warm bias from snow depth overestimation in models depends on whether the site has a "sub-or supra-critical" snow height. With supracritical conditions (e.g., at Schilthorn), the snow depth is so high that a small over-or underestimation in the model makes very little difference to the insulation. Only the timing of the snow arrival and melt-out is important. In sub-critical conditions (e.g., at Samoylov), the snow depth is so low that any overestimation leads to a strong warm bias in the simula-The Cryosphere, 9, 1343-1361, 2015 www.the-cryosphere.net/9/1343/2015/ tion e.g., for JULES/COUP. This effect is also mentioned in Zhang (2005), where it is stated that snow depths of less than 50 cm have the greatest impact on soil temperatures. However, overestimated snow depth at Samoylov and Schilthorn does not always result in warmer soil temperatures in models as expected (Fig. 3b, c). At these sites, even though JSBACH, JULES and COUP show warmer soil temperatures in parallel to their snow depth overestimations, ORCHIDEE and LPJ-GUESS show the opposite. This behavior indicates different processes working in opposite ways. Nevertheless, most of the winter, autumn and spring topsoil temperature biases can be explained by snow conditions (Fig. 5a). Figure 5b shows that snow depth bias can explain the topsoil temperature bias even when the snow-free season is considered, which is due to the long snow period at these sites (Table 2). This confirms the importance of snow representation in models for capturing topsoil temperatures at high latitudes and high altitudes. On the other hand, considering dynamic heat transfer parameters (volumetric heat capacity and heat conductivity) in snow representation seems to be of lesser importance (JS-BACH vs. other models, see Table 1). This is likely because a greater uncertainty comes from processes that are still missing in the models, such as wind drift, depth hoar formation and snow metamorphism. As an example, the landscape heterogeneity at Samoylov forms different soil thermal profiles for polygon center and rim. While the soil temperature comparisons were performed for the polygon rim, snow depth observations were taken from polygon center. Due to strong wind drift, almost all snow is removed from the rim and also limited to ca. 50cm (average polygon height) at the center . This way, models inevitably overestimate snow depth and insulation, in particular on the rim where soil temperature measurements have been taken. Hence, a resulting winter warm bias is expected (Fig. 2a, models JSBACH, JULES, COUP).
During the snow-free season, Samoylov has colder soil temperatures than air (Fig. 6c). Thicker moss cover and higher soil moisture content at Samoylov  are the reasons for cooler summer topsoil temperatures at this site. Increasing moss thickness changes the heat storage of the moss cover and it acts as a stronger insulator (Gornall et al., 2007), especially when dry (Soudzilovskaia et al., 2013). Additionally, high water content in the soil requires additional input of latent heat for thawing and there is less heat available to warm the soil.
Insulation strength during the snow-free season is related to model vegetation/litter layer representations. 10 cm fixed moss cover in JSBACH and a 10 cm litter layer in LPJ-GUESS bring similar amounts of insulation. At Samoylov, where strong vegetation cover is observed in the field, these models perform better for the snow-free season (Fig. 6c). However, at Bayelva, where vegetation effects are not that strong, 10 cm insulating layer proves to be too much and creates colder topsoil temperatures than observations (Fig. 6d). And for the bare Schilthorn site, even a thin layer of surface cover (2.5 cm litter layer) creates colder topsoil temperatures in LPJ-GUESS (Fig. 6b). At Bayelva, all models underestimate the observed topsoil temperatures all year long (Fig. 6d). With underestimated snow depth (Fig. 4d) and winter cold bias in topsoil temperature (Fig. 3d), models create a colder soil thermal profile that results in cooling of the surface from below even during the snow-free season. Furthermore, using global reanalysis products instead of site observations (Table 3) might cause biases in incoming longwave radiation, which can also affect the soil temperature calculations. In order to assess model performance in capturing observed soil temperature dynamics, it is important to drive the models with a complete set of site observations. These analyses support the need for better vegetation insulation in models during the snow-free season. The spatial heterogeneity of surface vegetation thickness remains an important source of uncertainty. More detailed moss representations were used in Porada et al. (2013) and , and such approaches can improve the snow-free season insulation in models.

Soil thermal regime
Model differences in representing subsurface temperature dynamics are related to the surface conditions (especially snow) and soil heat transfer formulations. The ideal way to assess the soil internal processes would be to use the same snow forcing or under snow temperature for all models. However most of the land models used in this study are not that modular. Hence, intertwined effects of surface and soil internal processes must be discussed together here. Figures 7-10 show the mismatch in modeled DZAA representations. Together with the soil water and ice contents, simulating DZAA is partly related to the model soil depth and some models are limited by their shallow depth representations (Fig. A1 in the Appendix, Table 1). Apart from The Cryosphere, 9, 1343-1361, 2015 www.the-cryosphere.net/9/1343/2015/ the different temperature values, models also simulate permafrost conditions very differently. As seen in Fig. 8, JS-BACH and JULES do not simulate permafrost conditions at Schilthorn. In reality, there are almost isothermal conditions of about −0.7 • C between 7 m and at least 100 m depth at this site (PERMOS, 2013), which are partly caused by the threedimensional thermal effects due to steep topography (Noetzli et al., 2008). Temperatures near the surface will not be strongly affected by three-dimensional effects, as the monitoring station is situated on a small but flat plateau , but larger depths get additional heat input from the opposite southern slope, causing slightly warmer temperatures at depth than for completely flat topography (Noetzli et al., 2008). The warm and isothermal conditions close to the freezing point at Schilthorn mean that a small temperature mismatch (on the order of 1 • C) can result in non-permafrost conditions. This kind of temperature bias would not affect the permafrost condition at colder sites (e.g., Samoylov). In addition, having low water and ice content, and a comparatively low albedo, make the Schilthorn site very sensitive to interannual variations and make it more difficult for models to capture the soil thermal dynamics . Compared to the other models with snow representation, OR-CHIDEE and LPJ-GUESS show colder subsurface temperatures at this site (Fig. 8). A thin surface litter layer (2.5 cm) in LPJ-GUESS contributes to the cooler Schilthorn soil temperatures in summer. Differences at Samoylov are more related to the snow depth biases. As previously mentioned, subcritical snow conditions at this site amplify the soil temperature overestimation coming from snow depth bias (Fig. 5). Considering their better match during snow-free season (Fig. 6c), the warmer temperatures in deeper layers of JULES and COUP can be attributed to overestimated snow depths for this site by these two models (Fig. 9). Additionally, JULES and COUP models simulate generally warmer soils conditions than the other models, because these models include heat transfer via advection in addition to heat conduction. Heat transfer by advection of water is an additional heat source for the subsurface in JULES and COUP, which can also be seen in the results for Bayelva (Fig. 10). In combination with that, COUP has a greater snow depth at Samoylov (Fig. 5), resulting in even warmer subsurface conditions than JULES. Such conditions demonstrate the importance of the combined effects of surface processes together with internal soil physics.
Due to different heat transfer rates among models, internal soil processes can impede the heat transfer and result in delayed warming or cooling of the deeper layers. JS-BACH, ORCHIDEE, JULES and COUP show a more pronounced time lag of the heat/cold penetration into the soil, while HYBRID8 and LPJ-GUESS show either a very small lag or no lag at all . This time lag is affected by the method of heat transfer (e.g., advection and conduction, see above), soil heat transfer parameters (soil heat capacity/conductivity), the amount of simulated phase change, vertical soil model resolution and internal model time step. Given that all models use some sort of heat transfer method including phase change (Table 1) and similar soil parameters (Table 3), the reason for the rapid warming/cooling at deeper layers of some models can be missing latent heat of phase change, vertical resolution or model time step. Even though the mineral (dry) heat transfer parameters are shared among models, they are modified afterwards due to the coupling of hydrology and thermal schemes. This leads to changes in the model heat conductivities depending on how much water and ice they simulate in that particular layer. Unfortunately, not all models output soil water and ice contents in a layered structure similar to soil temperature. This makes it difficult to assess the differences in modeled phase change, and the consequent changes to soil heat transfer parameters. A better quantification of heat transfer rates would require a comparison of simulated water contents and soil heat conductivities among models, which is beyond the scope of this paper.
The model biases in matching the vertical temperature curves (minimum, maximum, mean) are related to the topsoil temperature bias in each model for each site, but also the above-mentioned soil heat transfer mechanisms and bottom boundary conditions. Obviously, models without snow representation (e.g., HYBRID8) cannot match the minimum curve in Fig. 12. However, snow depth bias (Fig. 5) cannot explain the minimum curve mismatch for ORCHIDEE, COUP, and LPJ-GUESS at Schilthorn (Fig. 12b). This highlights the effects of soil heat transfer schemes once again.
In general, permafrost specific model experiments require deeper soil representation than 5-10 m. As discussed in Alexeev et al. (2007), more than 30 m soil depth is needed for capturing decadal temperature variations in permafrost soils. The improvements from having such extended soil depth are shown in Lawrence et al. (2012), when compared to their older model version with shallow soil depth (Lawrence and Slater, 2005). Additionally, soil layer discretization plays an important role for the accuracy of heat and water transfer within the soil, and hence can effect the ALT estimations. Most of the model setups in our intercomparison have less than 10 m depths, so they lack some effects of processes within deeper soil layers. However, most of the models used in global climate simulations have similar soil depth representations and the scope here is to compare models that are not only aimed to simulate site-specific permafrost conditions at high resolution but to show general guidelines for future model developments.
Adding to all these outcomes, some models match the site observations better than others at specific sites. For example, the mean annual soil thermal profiles are better captured by JSBACH at Nuuk, by JULES and COUP at Schilthorn, by ORCHIDEE at Samoylov, and by COUP at Bayelva (Fig. 11). Comparing just the topsoil conditions at the nonpermafrost Nuuk site, JSBACH better matches the observations due to its moss layer. On the other hand, by having better snow depth dynamics (Fig. 4), JULES and COUP modwww.the-cryosphere.net/9/1343/2015/ The Cryosphere, 9, 1343-1361, 2015 els are better suited for sites with deeper snow depths like Schilthorn and Bayelva. Contrarily, the wet Samoylov site is better represented by ORCHIDEE in snow season (Fig. 2a) due to lower snow depths in this model (Fig. 4) and thus colder soil temperatures. However, the snow-free season is better captured by the JSBACH model (Fig. 2c) due to its effective moss insulation and LPJ-GUESS model due to its insulating litter layer.

Active layer thickness
As seen above, surface conditions (e.g., insulation) alone are not enough to explain the soil thermal regime, as subsoil temperatures and soil water and ice contents affect the ALT as well. For Schilthorn, LPJ-GUESS generally shows shallower ALT values than other models (Fig. 13a); it also shows the largest snow depth bias (Fig. 5), excluding snow as a possible cause for this shallow ALT result. However, if snow depth bias alone could explain the ALT difference, ORCHIDEE would show different values than HYBRID8, which completely lacks any snow representation. At Schilthorn, COUP has a high snow depth bias ( Fig. 5) but still shows a very good match with the observed ALT (Fig. 13a), mainly because snow cover values at Schilthorn are very high so ALT estimations are insensitive to snow depth biases as long as modeled snow cover is still sufficiently thick to have the full insulation effect . All models overestimate the snow depth at Samoylov (Fig. 5) and most of them lack a proper moss insulation (Fig. 6c), which seems to bring deeper ALT estimates in Samoylov (Fig. 13b). However, HYBRID8 does not have snow representation, yet it shows the deepest ALT values, which means lack of snow insulation is not the reason for deeper ALT values in this model. As well as lacking any vegetation insulation, soil heat transfer is also much faster in HYBRID8 (see Sect. 3.2), which allows deeper penetration of summer warming into the soil column.
Surface conditions alone cannot describe the ALT bias in Bayelva either. LPJ-GUESS shows the lowest snow depth (Fig. 5) together with deepest ALT (Fig. 13c), while JULES shows similar snow depth bias as LPJ-GUESS but the shallowest ALT values. As seen from Fig. 10, LPJ-GUESS allows deeper heat penetration at this site. So, not only the snow conditions, but also the model's heat transfer rate is critical for correctly simulating the ALT.

Conclusions
We have evaluated different land models' soil thermal dynamics against observations using a site-level approach. The analysis of the simulated soil thermal regime clearly reveals the importance of reliable surface insulation for topsoil temperature dynamics and of reliable soil heat transfer formula-tions for subsoil temperature and permafrost conditions. Our findings include the following conclusions.
1. At high latitudes and altitudes, model snow depth bias explains most of the topsoil temperature biases.
2. The sensitivity of soil temperature to snow insulation depends on site snow conditions (sub-/supra-critical).
3. Surface vegetation cover and litter/organic layer insulation is important for topsoil temperatures in the snowfree season, therefore models need more detailed representation of moss and top organic layers.
4. Model heat transfer rates differ due to coupled heat transfer and hydrological processes. This leads to discrepancies in subsoil thermal dynamics.
5. Surface processes alone cannot explain the whole soil thermal regime; subsoil conditions and model formulations affect the soil thermal dynamics.
For permafrost and cold-region-related soil experiments, it is important for models to simulate the soil temperatures accurately, because permafrost extent, active layer thickness and permafrost soil carbon processes are strongly related to soil temperatures. There is major concern about how the soil thermal state of these areas affects the ecosystem functions, and about the mechanisms (physical/biogeochemical) relating atmosphere, oceans and soils in cold regions. With the currently changing climate, the strength of these couplings will be altered, bringing additional uncertainty into future projections.
In this paper, we have shown the current state of a selection of land models with regard to capturing surface and subsurface temperatures in different cold-region landscapes. It is evident that there is much uncertainty, both in model formulations of soil internal physics and especially in surface processes. To achieve better confidence in future simulations, model developments should include better insulation processes (for snow: compaction, metamorphism, depth hoar, wind drift; for moss: dynamic thickness and wetness). Models should also perform more detailed evaluation of their soil heat transfer rates with observed data, for example comparing simulated soil moisture and soil heat conductivities.