Improved soil physics for simulating high latitude permafrost regions by the JSBACH terrestrial cosystem model

Introduction Conclusions References


Introduction
Effects of global climate change are felt stronger in the northern high latitudes than the rest of the world (ACIA, 2005).During recent decades, polar regions have experienced around +0.5 • C to +1 • C increase in surface atmospheric temperature, while the global mean has risen by only +0.2 • C 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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).22 % of the Northern Hemisphere land is underlain by permafrost (Gruber, 2012).During the past glacial/interglacial cycles vast amounts of organic matter have been accumulated in these soils (Ciais et al., 2011).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).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 is 1672 petagrams of carbon stocked in the northern permafrost soils.With the current trend of increasing air temperature, this carbon is susceptible to thawing and being released to 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 spurned recent advancements of dynamic global vegetation models and Earth system models by representing the processes that are specific to 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 (Riseborough et al., 2008).At present, most of the global models include basic processes related to permafrost regions, e.g.latent heat release/consumption from the phase change of soil water.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 (M ölders et al., 2003;Takata and Kimoto, 2000) or the more popular method of "apparent heat capacity" is incorporated into the temperature calculations (Beer et al., 2007;Introduction Conclusions References Tables Figures
Besides the freeze/thaw events, the coupling of snow and soil thermal modules constitutes the basis for the soil thermal profile during the winter periods (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 growing season and the active layer thickness, which is also related to the amount of infiltrating snowmelt water into the soil.Shown by Gouttevin et al. (2012b), the snow cover and the disappearance of snow are important factors for the plant and soil metabolic activity and the biogeochemical feedbacks between the soil and the atmosphere.However, in most cases snow is represented rather simply in the models.Due to high complexity of snow types and snow processes, a simple parameterization yielding a realistic heat insulation effect was used (Koren et al., 1999;Verseghy, 1991).While more advanced snow schemes were developed in some modelling studies (Boone and Etchevers, 2001;Loth and Graf, 1998), it is not always practical for global modelling 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) et al., 2012;Hayes et al., 2012;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).Introduction

Conclusions References
Tables Figures

Back Close
Full Although progress has been undertaken on representing permafrost processes in global models, there is still a considerable uncertainty regarding the magnitude of the effects of permafrost feedbacks on climate and northern ecosystems.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 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 the implementation of frozen ground physics was constructed in a similar way (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 permafrost-related processes within global models is a challenging yet essential task.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 MPIOM 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 processbased 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 Introduction

Conclusions References
Tables Figures

Back Close
Full   et al. (1980) and Collatz et al. (1992).The BETHY 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 5 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 timestep, the numerical solution to heat conduction (first term on the right side of Eq. 1) gives the soil layer temperatures and then as a second step, these temperatures are 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 ), ρ: density of ice (kg m −3 ), θ i : volumetric soil ice content (m 3 m −3 ), t: time (s), z: soil layer depth (m).A new soil hydrology scheme (Hagemann et al., 2013) has been implemented into the current JSBACH version that uses the same five layer structure 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 ∂θ/∂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 whereat 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).
In the hydrology module, first the input/output terms (precipitation, snow melt, 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 and changes the field capacity of each layer with respect to the simulated ice volume.If the water and ice contents are fully occupying the field capacity, that layer is blocked for a further water transfer.Finally the vertical Introduction

Conclusions References
Tables Figures

Back Close
Full water movement is performed as described above and the soil water content at each layer is updated.
with: θ: volumetric soil water content (m 3 m −3 ) D: soil water diffusivity (m 2 s −1 ) K : soil hydraulic conductivity (m s −1 ) S: source and sink terms 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: freezing temperature of water (K) g: gravitational acceleration (m s −1 ) Ψ sat : saturated soil matric potential (m) Soil heat transfer is coupled with the hydrological scheme through 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 (1977) method (Eq.5).Equations ( 6)-( 9) describe the terms in Eq. ( 5).With these formulations, the amount of water and ice influence the soil thermal properties.In concert with the latent heat of fusion effect on Introduction

Conclusions References
Tables Figures

Back Close
Full temperature (second term on the right side of Eq. 1), a coupling of the hydrology and soil thermal dynamics is achieved.
with: ρ s , ρ w , ρ i : density of soil solids, water and ice respectively (kg m −3 ) c s , c w , c i : heat capacities of soil solids, water and ice respectively (J m −3 K −1 ) θ w , θ i : volumetric soil water and ice contents (m 3 m −3 ) 5) λ dry = 0.135ρ bulk + 64.7 2700 − 0.947ρ bulk ( 8) (always at the bottom).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 snow properties are kept constant for simplicity.A density of 250 kg m −3 is used for the snow depth calculations and the snow heat conductivity is fixed at 0.31 W K −1 m −1 with a snow heat capacity of 522 500 J m −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 transfer of melt water 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).For a simple representation, we assume a constant moss layer over the soil.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 W K and volumetric heat capacity of 2 500 000 J m −3 K −1 .

Global forcing data
For the period 1901-1978, daily forcing data with 0.5 • spatial resolution from the EU project WATCH was used (Weedon et al., 2010).This data is based on ERA-40 reanalysis results that were bias corrected by using several observation-based datasets, such as climate grids from the Climate Research Unit, University of East Anglia (CRU).
For the 1979-2010 period, ECMWF ERA-Interim reanalysis data (Dee et al., 2011) was used.This dataset was downloaded at 0.5 • spatial resolution and bias-corrected against the WATCH forcing data., following Piani et al. (2010).A more detailed description of the climate forcing dataset can be found in Beer et al. (2013b).With this approach, a consistent time series of climate data for the period 1901-2010 is ensured.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 Vegetation Continuous Fields product (Hansen et al., 2003) and the WWF biome map (Olson et al., 2001).JSBACH was forced by global atmospheric carbon dioxide concentrations following the CMIP5 protocol (Meinshausen et al., 2011).

Simulation setup
Nuuk site level simulations were performed running the model at a single point forced by hourly meteorological site observations.Soil parameters were extracted from the above-mentioned global land surface data.Using the observed meteorological data, an average seasonal cycle was prepared and repeated for 30 yr 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.
For the circumpolar simulations, the model was run using the previously described global daily forcing data for the grids above 50 • N. First, the model physical state was brought into equilibrium with a 30 yr run repeating an average seasonal cycle of climate  2009,2010).The meteorological (incoming radiation, air temperature, precipitation and wind speed) and soil observations (soil temperature) were downloaded from the Greenland Ecosystem Monitoring database web server (ZackenbergGIS).An hourly climateforcing file is created from the downloaded data and used to force the Nuuk site level simulations.

Circum-Arctic datasets
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 permafrost extent, which is calculated using the simulated soil temperatures from the circum-arctic model simulation.Following the permafrost definition of IPA (soils under freezing temperatures for at least two consecutive years), the permafrost state of each gridbox is determined.For comparing with the IPA map, the 1990 values of the model permafrost state were used.The Circumpolar Active Layer Monitoring Network's (CALM) dataset (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 • gridbox of Introduction

Conclusions References
Tables Figures

Back Close
Full the simulation results that is conducted using global climate and soil texture data as forcing.The simulated soil temperatures in 5 soil layers are interpolated into 200 evenly spaced nodes and the depth of zero degree is calculated afterwards to represent the thawing depth at each timestep.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 gridbox, the most appropriate one is chosen for the comparison.Averaging several CALM sites within one gridbox 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 1990-2010, 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 circum-arctic 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 gridbox values of the JSBACH simulation results were used for comparison.
Since there were more boreholes in most of the gridboxes 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 up-scaled product of several Introduction

Conclusions References
Tables Figures

Back Close
Full 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.Following the observational time period, mean JSBACH subsoil temperature (last soil layer, ca.7 m) during 1960-1990 was used for comparison.The 0.5 • active layer thickness map (Beer et al., 2013a) from the 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 East Siberia, this map is very useful to understand the permafrost conditions at 1 : 2 500 000 spatial scale during 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 0.5 • spatial scale during 1960-1990 were used to derive model active layer thickness, and then the mean of all these years are used to prepare the comparison map.

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, Lena and Yenisey river runoff data were compared to simulation results.The runoff observations at the river mouth stations Introduction

Conclusions References
Tables Figures

Back Close
Full were gathered from the R-ArcticNET database (Lammers et al., 2001).The simulated runoff values in all the gridboxes 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.
3 Results and discussions

Site level validation at Nuuk
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 (Fig. 2).Following the observations, sum- A more detailed analysis of this comparison 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 Nuuk 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 (Fig. 2).
In addition, winter temperatures are highly affected by the overlying snow pack.Without snow cover, cold atmospheric conditions lead to extremely cold soil temperature in winter.However, having snow above the soil can act as a blanket due to snow's physical properties.For these reasons, it is essential to include an advanced representation Introduction

Conclusions References
Tables Figures

Back Close
Full of snow cover in high latitudes where snow is present all winter.It is seen from Nuuk site level comparisons that winter soil temperature 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, not shown), soil keeps a rather warm temperature profile (ca.−3 • C, Fig. 2) as long as snow exists on top.
On the other hand, snow has rather complicated characteristics in reality.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.Effects from this change propagate into soil; hence, late winter snow has less insulating effect than early winter snowpack, allowing for stronger coupling between air and soil temperature towards the end of the season.These dynamics explain the mismatch in simulated versus observed springtime soil temperature in the site level simulations.Without dynamically changing snow properties, our model cannot correctly represent the spring insulation and keeps a colder soil temperature profile owing to constant snow heat parameters used all season.This underestimation of spring soil temperature also causes a time lag in soil warming and hence subsequent underestimation of soil temperature in early summer.Langer et al. (2013) pointed out the importance of correct parameterization of snow thermal properties in permafrost simulations.Further progress in resolving these issues will be shown in the next model version.

Circum-Arctic validation
To evaluate the model's reliability at circum-arctic 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 by 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 Introduction

Conclusions References
Tables Figures

Back Close
Full permafrost or non-permafrost conditions.Having this in mind, it is seen from Fig. 3 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. 4).
Reasons for this mismatch are mostly explained by scale issues and site-specific conditions.Firstly, the model output from 0.5 • × 0.5 • gridbox cannot be taken as equally comparable to the site observations given that the gridbox 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 2 ×1 km 2 areas, the landscape variability still brings up a big uncertainty when compared to a model gridbox average.In addition, there is a biased selection of CALM sites, favouring northern slopes of mountains and taking observations from logistically advantageous points rather than selecting the most representative point of the surrounding environment.Besides, the observational technique brings additional uncertainty of the observations (possible disturbance from the instruments) and not accounting for the natural process of soil subsidence in the model further complicate the comparison to our results.While the level of top soil changes due to melting of ice within the soil or disturbances from observational Introduction

Conclusions References
Tables Figures

Back Close
Full activities, the model keeps the same level for the soil column top for the entire simulation period and hence estimating deeper active layer thickness compared to observations.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.74, Fig. 4).Similar results are observed in some other modeling studies.Dankers et al. ( 2011) have shown a deeper simulated ALT using the JULES model.Lawrence et al. (2012) have shown that the coupled and uncoupled CLM 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 47 % of observed subsoil temperature variation with a tendency to a cold bias at some sites (Fig. 5).This cold bias can partly be related to the model assumption of zero heat flux at the bottom of the soil.A deeper soil column representing up to 50 m is suggested to improve the permafrost temperature results around 10 m for future model versions.In addition, the initialization of the global model run was achieved by forcing the model with 1901-1930 average climate conditions, hence neglecting true climatic conditions during previous time periods.Nonetheless, the borehole temperature comparison illustrates the general reliability of our model in representing permafrost temperature.
Model results showing deeper active layers seem to disagree with colder soil temperature 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.On the other hand, deep soil temperature is less influenced by variable surface conditions, but Introduction

Conclusions References
Tables Figures

Back Close
Full show a decadal trend that is strongly affected by long-term atmospheric changes, snow 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 is 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 a slight underestimation of the model output (Fig. 6). Figure 7 shows the spatial pattern of this cold bias.In general, permafrost temperature differs by −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.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 modelling community (Dankers et al., 2011;Lawrence et al., 2008), but evidently the soil column depth also plays an important role (Alexeev et al., 2007).We anticipate that performing a simulation with deeper soil column can overcome these cold biases in the model results.
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 gridbox individually by using Introduction

Conclusions References
Tables Figures

Back Close
Full meteorological forcing and surface conditions specific to each gridbox.This explains the longitudinal changes in the model output.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 8 shows the comparison of the active layer thickness map of Yakutia (Beer et al., 2013a) and the spatial distribution of active layer thickness estimated by JSBACH.As in the permafrost temperature comparison (Fig. 6), 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 0.5 m at the coast to 2.5 m further inland is comparable to the JSBACH map.The differences between the observed and simulated results (Fig. 9) show a more diverse spatial pattern than the map of temperature differences (Fig. 7).This is due to the complex nature of confounding factors of active layer thickness, i.e. soil temperature, snow-moss cover and soil moisture.There is an overestimation at the coastal areas.Regions with positive bias (ca.+50 cm) are also seen in eastern Siberia and mid-western parts of the map.These areas coincide with high altitudes: Verkhoyansk Range (1500 m) in the east and Central Siberian Plateau (1000 m) in the west.Mountainous regions have more complex hydrology (slope effects) and the vegetation cover changes remarkably with height.Relating to that, a 10 cm thick moss layer over the entire model domain produces further uncertainty, in particular at high altitudes.An overestimation of moss layer insulation leads to warmer topsoil temperature and deeper active layer thickness in mountainous areas.On the other hand, the model shows a negative bias (ca.−30 cm) in the lowlands (blue regions in Fig. 9).This could be the opposite consequence of not representing vegetation and snow cover heterogeneity properly: having less insulation leading to colder topsoil temperatures in lowlands.Finally, the southern borders of the comparison map (Fig. 9) illustrate a strong positive bias (ca.+120 cm).
As explained in Beer et al. (2013a), isolated permafrost patches are dominant in these areas.However, the comparison is not very useful for these areas since the 0.5 • values

Conclusions References
Tables Figures

Back Close
Full 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.Overall, the model performance in representing active layer thickness in Yakutia seems reasonable and in most of the areas the error ranges from 20 cm to 80 cm, which is quite common for global model simulations.

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 least influenced basin areas from anthropogenic activities and represents a more natural pattern that is easily comparable to the model results.The improved model version can simulate the annual changes (Fig. 10) and the monthly dynamics (Fig. 11) 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 matched the average annual runoff values (Fig. 12) and the monthly dynamics (Fig. 13).The only issue here is that the model underestimated the annual runoff, which can be explained by the bigger size of the Yenisey river basin.Having a 2 580 000 km 2 basin area, 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 the 1980.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.

GMDD Introduction
Full

Conclusions
In this paper we have presented an advanced version of the process-oriented ecosystem model JSBACH that ameliorates the simulations of cold regions through enhanced representation of soil physics.By including the phase change process, coupled thermal and hydrological processes and heat insulation from snow and moss cover, the new model version is a capable tool for simulating the physical state of high latitude terrestrial environments.A multi-scale evaluation was conducted and the results demonstrate the strength and weaknesses of the model.Site-level comparisons 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 circum-arctic observational datasets 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.

GMDD Introduction
Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | that soil freezing 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 studies, future implications of possible permafrost carbon release are investigated and their effects on global climate are shown under different warming scenarios (Burke Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the saturated soil (W K −1 m −1 ) λ dry : heat conductivity of the dry soil (W K −1 m −1 ) S: saturation (water + ice volume/layer depth) λ s , λ w , λ i : heat conductivities of soil solids, water and ice respectively (W K −1 m −1 ) ρ 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 5 snow layers.The first 4 layers are always 5 cm in thickness (always at the top), while the last layer is unlimited in size Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | variables from the period 1901-1930.Then, a climate-transient run with constant atmospheric CO 2 concentration at 1901 value was executed for the same period.These 30 yr model results were further used to force a 1000 yr carbon balance model run in order to prepare equilibrated carbon pools.Finally, these carbon pools are used as initial condition to start a fully transient run from 1901 to 2010.The Nuuk observational site is on the southwestern coast of Greenland, 250 km south of the polar circle at around 64 • N and 51 • W. It is situated in the Kobbefjord with an altitude of 500 m a.s.l.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 and Rasch, Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 • C or 2 • C) show a scale from 0 • C to −17 • C.
Discussion Paper | Discussion Paper | Discussion Paper | mer 2008 topsoil temperatures gradually cool down to 0 • C. Simulated temperatures fluctuate around 0 • C from October to late February, in agreement with the observed data.However, JSBACH soil temperatures are below zero until June while the observational records show a warming-up period to 0 • C from April to May 2009 and keep zero values until June.When the simulated temperatures are over zero after June, it takes until mid-July to capture the observations again during the summer of 2009.
Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 : 8 Fig. 1 .
Figure 1: Vertical soil model structure of the new JSBACH version.The numbers left to the soil 3 column show the depths of the bottom of each layer while the numbers right to the soil column 4 show layer thicknesses.Here snow and organic layer are simply shown to represent multi layered 5 snow scheme and constant moss layer described in the text.6 7 8

Figure 2 :
Figure 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.3cm).