Simulation of permafrost and seasonal thaw depth in the JULES land surface scheme

Land surface models (LSMs) need to be able to simulate realistically the dynamics of permafrost and frozen ground. In this paper we evaluate the performance of the LSM JULES (Joint UK Land Environment Simulator), the stand-alone version of the land surface scheme used in Hadley Centre climate models, in simulating the largescale distribution of surface permafrost. In particular we look at how well the model is able to simulate the seasonal thaw depth or active layer thickness (ALT). We performed a number of experiments driven by observation-based climate datasets. Visually there is a very good agreement between areas with permafrost in JULES and known permafrost distribution in the Northern Hemisphere, and the model captures 97 % of the area where the spatial coverage of the permafrost is at least 50 %. However, the model overestimates the total extent as it also simulates permafrost where it occurs sporadically or only in isolated patches. Consistent with this we find a cold bias in the simulated soil temperatures, especially in winter. However, when compared with observations on end-of-season thaw depth from around the Arctic, the ALT in JULES is generally too deep. Additional runs at three sites in Alaska demonstrate how uncertainties in the precipitation input affect the simulation of soil temperatures by affecting the thickness of the snowpack and therefore the thermal insulation in winter. In addition, changes in soil moisture content influence the thermodynamics of soil layers close to freezing. We also present results from three experiments in which the standard model setup was modified to improve physical realism of the simulations in permafrost regions. Extending the soil column to a depth of 60 m and adjusting the soil parameters for organic content had relatively little effect on the Correspondence to: R. Dankers (rutger.dankers@metoffice.gov.uk) simulation of permafrost and ALT. A higher vertical resolution improves the simulation of ALT, although a considerable bias still remains. Future model development in JULES should focus on a dynamic coupling of soil organic carbon content and soil thermal and hydraulic properties, as well as allowing for sub-grid variability in soil types.


Introduction
The impact of climate change on permafrost in the circumpolar arctic has received much attention in recent years (Anisimov and Nelson, 1996;Stendel and Christensen, 2002;Lawrence and Slater, 2005;Lawrence et al., 2008;Zhang et al., 2008).This is partly because of the dramatic rise in temperature at northern high latitudes that is projected by many General Circulation Models (GCMs) for the coming centuries (Kattsov et al., 2005;Christensen et al., 2007).Additionally, there has been a growing concern about the amount of organic matter stored in currently frozen soils that may start to decompose once the permafrost thaws (Zimov et al., 2006;Schuur et al., 2008).Depending on whether this decomposition takes place under aerobic or anaerobic conditions, it may result in enhanced emissions of carbon dioxide (CO 2 ) or methane (CH 4 ) and other greenhouse gases, all with the ability to feed back to the regional and global climate, thereby accelerating future warming (Walter et al., 2007;Elberling et al., 2010).Although it is fair to say that some of these concerns may have been overstated (Anisomov, 2007), it is obvious that climate models need to be able to simulate realistically permafrost dynamics and its effects on the carbon balance in order to take this feedback into account.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Recent years have seen a growing number of permafrost models, often developed with a specific application in mind (Riseborough et al., 2008).Conventional permafrost modelling often assumes the soil thermal regime is in equilibrium with the climate.In reality the thermal state of the ground depends on a number of complex interactions between soil, vegetation, snow and hydrology (Sturm et al., 2001;Smith et al., 2005).When applied over large regions, permafrost modelling is further complicated by the complex patterns of vegetation and topography typical for arctic environments (Sazonova and Romanovsky, 2003;Stendel et al., 2007) and limited data availability on meteorological conditions, as well as soil and land surface properties.Anisimov et al. (2002) and Anisimov (2009) addressed part of this problem by using a stochastic modelling approach to account for the high spatial variability in large-scale permafrost models.Aside from portraying the level of uncertainty in the input parameters, the output from such a model can be used to estimate the probability of permafrost temperature and thaw depth exceeding a given threshold at a certain location.
In weather and climate models, the interaction between the atmosphere and the land surface is simulated by land surface models (LSMs).These LSMs are designed to represent the physical processes controlling the exchange of heat and moisture in order to solve the surface energy balance, typically by partitioning the available energy between evaporative, sensible and ground heat fluxes.In climate science there has been a growing recognition that different parts of the Earth system affect one another and that these feedbacks, which often involve land ecosystem-atmosphere interactions, need to be included in the models in order to achieve improved projections for the future.Consequently, LSMs have grown in complexity in an effort to include processes such as changes in vegetation cover, carbon cycling in terrestrial ecosystems and the direct effect of rising CO 2 concentrations on plant physiology.In spite of their limitations (see e.g.Slater et al., 2001;Nijssen et al., 2003), LSMs offer the opportunity to study changes in the terrestrial environment in an integrated and consistent way and to explore key interactions between the different impacts (Betts, 2007).LSMs are thus not specific permafrost models, but if claims of LSMs being physically process-based models have any validity, they ought to be capable of simulating the dynamics of permafrost and frozen ground with some realism.
The aim of this paper is to evaluate the performance of the Joint UK Land Environment Simulator (JULES) (Blyth et al., 2006), the LSM used within the Hadley Centre climate models, in simulating large-scale features of surface permafrost, and to identify areas for improvement.To achieve this we performed a number of off-line experiments with JULES driven by observation-based climate datasets.We also present results from two experiments in which the standard model setup was modified with the aim of improving the physical realism of the simulations in permafrost regions.

JULES model description
JULES is the stand-alone version of the land surface scheme used in the Hadley Centre climate models (Blyth et al., 2006).It is based on the Met Office Surface Exchange Scheme (MOSES) described by Cox et al. (1999) and Essery et al. (2001) and combines a complex energy and water balance model with a dynamic vegetation model.As a community model JULES is available for everybody to use and/or to contribute to its further development (see the JULES website http://www.jchmr.org/julesfor more details).
JULES describes the physical, biophysical and biochemical processes that control the exchange of radiation, heat, water and carbon between the land surface and the atmosphere.It can be applied as a point or a grid model.When applied in distributed mode each grid box can have several sub-grid land cover fractions or "tiles".JULES has five vegetation tiles representing different plant functional types (broad-leaf trees, needle-leaf trees, C3 (temperate) grass, C4 (tropical) grass, and shrubs) and four non-vegetated surface tiles (urban, inland water, bare soil and ice).Each tile has its own surface temperature, shortwave and longwave radiative fluxes, sensible and latent heat fluxes, ground heat flux, canopy moisture content, snow mass and snow melt.The soil underneath the land cover tiles is, however, assumed to be homogeneous across the grid box.By default, JULES does not account for lateral transfer of energy and water between grid boxes, but the model has successfully been coupled to flow routing schemes to simulate river runoff (e.g.Dadson et al., 2010).
In this paper we use simulations with JULES version 2.1.2that has a multi-layer snow scheme (described in detail by Best et al., 2011) in which the number of snow layers varies according to the depth of the snow pack.Each snow layer has a prognostic temperature, density, grain size and solid and liquid water content.The snowmelt heat flux is calculated by solving the surface energy balance.The subsurface temperatures are updated using a discretised form of the heat diffusion equation, which is coupled to the soil hydrology module through both soil water phase changes and the associated latent heat fluxes, and the soil thermal characteristics, which are dependent on frozen and unfrozen soil moisture content.The soil heat flux at the surface is calculated from the surface energy balance, while the lower boundary condition corresponds to a zero vertical gradient in soil temperature.
The soil hydrology solution is based on a finite difference approximation to the Richards' equation (Richards, 1931), using the same vertical discretisation as in the calculation of the soil thermodynamics.Like many land surface schemes JULES uses the Brooks and Corey (1964) relations (later modified by Clapp and Hornberger, 1978) to describe the soil water retention curve and deduce hydraulic conductivity and soil water suction as a function of soil moisture content.This is done for the unfrozen rather than the total soil moisture content, which is consistent with the notion that freezing of the soil lowers the hydraulic conductivity and produces a large suction by reducing the unfrozen water content (Williams and Smith, 1989).The soil hydraulic parameters are calculated using the Cosby et al. (1984) parameterisations for different soil particle size distributions.For a detailed discussion of the process descriptions in JULES, the reader is referred to Best et al. (2011).
The ability of JULES to partition incoming radiation into sensible and latent heat and how this varies on annual, seasonal and diurnal timescales has been tested at a range of FLUXNET sites by Blyth et al. (2010).The overall performance was good, but specifically in cold climates it was found that the model continues to simulate evaporation when observations indicate that transpiration is inhibited by frozen soils.JULES has also participated in various model intercomparisons, such as the forest snow model (Rutter et al., 2009) and water model (Haddeland et al., 2011) intercomparison projects.

Model setup and input data
We set up a number of model runs with JULES in order to evaluate its performance in representing the large scale characteristics of permafrost distribution.In JULES the number of soil layers can be chosen by the user.Here we applied the model in its standard configuration of four layers with a thickness of 0.10, 0.25, 0.65 and 2.0 m, giving a total depth of 3.0 m.We chose this setup to test the model in its standard configuration and also because it is consistent with how the land surface scheme is implemented in the Hadley Centre climate models.For simulating the full permafrost dynamics, such a shallow depth is obviously insufficient.The permafrost layer can extend hundreds of meters deep where it can persist for centuries.Alexeev et al. (2007) showed that the soil column needs to be at least 30 m deep in order to properly resolve the annual cycle in temperature.For decadal to century-scale changes this is arguably even deeper (Alexeev et al., 2007).However, Alexeev et al. (2007) also note that extending the soil column in a LSM may improve the simulation of soil temperatures but not necessarily of the soil hydrology, which in a model like JULES is tightly coupled to the thermodynamics.A realistic circum-arctic simulation of deep permafrost is furthermore hampered by a general lack of knowledge about the sub-surface structure.
Here we include results from an experiment with JULES where the soil column has been extended to a depth of 60 m by adding three additional layers with increasing thickness of 5, 14 and 38 m.These additional layers act as normal model layers meaning they simulate the hydrology as well as the heat exchange.The aim of this extended soil profile is however not to simulate the thermodynamics of the full permafrost layer, which would require a much more detailed setup, but to account for the "heat sink" effect of the deeper permafrost in future climate change simulations (see e.g.Dankers et al., 2010), as well as its influence on nearsurface temperatures.
Although JULES cannot be used to simulate the dynamics of the full permafrost layer, it can however be evaluated on its simulation of the freeze/thaw status of the near-surface soil.In particular we look at how well the model is able to simulate the penetration of summer warming into the frozen soil.The uppermost layer of seasonal thawing or the active layer is an important regulator of energy and mass fluxes between the surface and the atmosphere in arctic environments (Anisimov et al., 2002) and as such it is an essential feature for LSMs to capture.The ability of JULES to model this layer is also important when used within a GCM to simulate carbon release from permafrost regions.

Forcing data
A major limitation in evaluating a complex model like JULES is the availability of micrometeorological observations of sufficient quality, frequency and duration to drive the model.LSMs like JULES are designed to provide the lower boundary conditions to an atmospheric model and typically require atmospheric forcing data at high frequency, in the order of hourly to 3-hourly, without gaps in space or time.This is no problem when coupled to a GCM but may be an important constraint when testing the model offline.JULES has been benchmarked using local driving data at a range of FLUXNET sites (Blyth et al., 2010;Van den Hoof et al., 2011) but few sites with sufficient data coverage are located in the permafrost region.Before using any particular driving data set care must be taken to ensure that the appropriate corrections for precipitation, particularly winter snowfall, have been applied.To circumvent these problems, we chose to drive JULES with the following readily available gridded meteorological datasets: GSWP2: these data are provided by the Global Soil Wetness Project 2 (GSWP2, Dirmeyer et al., 2005) that aimed at comparing a broad range of LSMs under controlled conditions.The data are a hybridisation of observational and model reanalysis from the National Centers for Environmental Prediction/Department of Energy (NCEP/DOE).This means that systematic biases in the reanalysis fields were corrected by blending the 3-hourly analysis with global observation-based gridded data at a lower temporal resolution.The precipitation product in particular has been corrected to match the observed monthly precipitation from the Climate Research Unit (CRU), Global Precipitation Climatology Centre (GPCC) and Global Precipitation Climatology Project (GPCP) databases after correcting for wind-induced undercatch.Dirmeyer et al. (2005) note that the GSWP2 product tends to overcorrect precipitation, particularly in the case of snow.This suggests the GSWP2 precipitation may in places be too high.The data are provided on a 1 for the period July 1982 to December 1995 with a 3-hourly resolution.WATCH: this dataset was created for the project Water and Global Change (WATCH, Weedon et al., 2010) and is derived from the ERA-40 reanalysis product.The reanalysis fields were interpolated to half-degree resolution, corrected for elevation, and monthly adjustments were applied based on gridded observations from the CRU (temperature, diurnal temperature range, cloud-cover and number of wet days) and GPCC (precipitation) databases.An alternative precipitation product based on the CRU observations only was also created but this includes fewer stations than the GPCC data (Weedon et al., 2010).In addition corrections were applied for atmospheric aerosol loading and separate precipitation gauge corrections for rainfall and snowfall.The WATCH forcing data cover the period 1901 to 2001 but in this paper we use an earlier version starting from 1959.The data are provided with 3-hourly time step and at a half-degree resolution.The JULES simulations described in this paper that were driven by WATCH primarily use the (standard) GPCC precipitation (WATCH-GPCC) product but we did make a comparison with the CRU precipitation (WATCH-CRU) at a limited number of sites (see Sect. 6).

Model parameters
In addition to meteorological input, JULES requires information on vegetation fractions and soil parameters (see Sect. 2).
In our runs we used "standard" input datasets that are also commonly used in Hadley Centre climate model experiments.Vegetation types were derived from the International Geosphere Biosphere Programme's (IGBP) global land cover database1 that is based on satellite data from the period April 1992 through March 1993 with a resolution of 1 km.The soil data originally derive from the Wilson and Henderson-Sellers (1985) soils dataset that provides information on soil classes on a global 1-degree grid.Based on fractions of sand, silt and clay for each soil type the soil parameters used in the Clapp and Hornberger (1978) soil hydrology scheme are calculated using the equations suggested by Cosby et al. (1984).Typical parameter values obtained in this way are given in Table 1.
In its standard setup (and as used within GCMs), JULES soil parameters are constant with depth and for mineral soils only.However, as many have pointed out (e.g.Beringer et al., 2001;Nicolsky et al., 2007), organic soil behaves very differently from mineral soil.This is particularly important in the Arctic where mineral soils are often overlain by a layer of peat, mosses or lichens, providing thermal insulation to the lower mineral soil layers (Gornall et al., 2007).Although JULES can be used to simulate the soil carbon cycle, at present it does not include the influence of soil organic matter on the thermodynamic and hydrological properties of the soil.Here we include the results of an experiment in which the standard mineral soil parameters were adjusted according to soil organic content.This adjustment is similar to that used in other models (e.g.Lawrence and Slater, 2008) and is based on a linear combination of mineral and organic soil properties, according to the organic content.The properties of organic soil are based on Letts et al. (2000) and Oke (1987) and, in contrast to Lawrence and Slater (2008), were allowed to vary with depth (see Table 2; note that in Lawrence and Slater (2008) organic content itself does change with depth).
Rather than assuming a blanket layer of organic material on top of the soil column (as in e.g.Rinke et al., 2008) we obtained data on soil organic content in permafrost regions from the Northern Circumpolar Soil Carbon Database (NCSCD) (Tarnocai et al., 2009).These data were regridded to a 1degree grid and distributed over the top three model layers.
The NCSCD provides soil organic carbon content in kg m −2 over two depths, 0-30 cm and 0-100 cm.The first data were distributed over the top two layers of JULES; the remainder was assigned to the third layer.In the fourth (bottom) model layer the carbon content was assumed to be zero.A maximum value for soil carbon of 130 kg m −3 (Farouki, 1981) was used to calculate the organic fraction from the actual carbon content in each layer.This organic fraction then determines the extent that the organic soil properties are allowed to influence the final soil properties: when the organic fraction is very low, the soil properties are similar to the original mineral soil parameters.In grid cells with a high organic fraction, the parameters are close to the properties of pure organic soil given in Table 2.

Model experiments
We ran JULES for the Northern Hemisphere land mass driven by each of the meteorological forcing datasets described above (see Table 3).In the following, each simulation is designated by its forcing, i.e.GSWP2, WATCH.The GSWP2 simulations cover areas north of 25 • N while WATCH-GPCC was run for areas north of 45 • N, meaning the latter does not include the Qinghai-Tibetan Plateau.In each experiment the grid resolution was kept the same as in the driving dataset, i.e. 1 • × 1 • in the GSWP2 runs and 0.5 • × 0.5 • in WATCH-GPCC, maintaining consistency with the forcing.The internal time step in all simulations was 1 h.Before the actual simulation period the model was spun-up by running repeatedly using the meteorological input of the first year, until soil moisture and temperature in each layer reached equilibrium.
The model experiments with extended soil profile (DEEP) and with organic soil parameters (SOC) were run driven by GSWP2 forcing only.A third experiment combining the two modifications (SOC + DEEP) was also run using GSWP2.In the runs with the extended soil profile a longer spin-up period is required to bring the slow-responding lower soil layers into equilibrium with the climate.To speed up this process, the model was run with an internal time step of 3 h for 300 yr, driven by the climatological average of the GSWP2 forcing over 1983-1995, thus discarding the interannual variability in this period.
Additionally, we performed a number of point-scale simulations for comparison with observed soil temperature and moisture content at three soil climate monitoring sites in Alaska.The data was obtained from the United States Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS) 2 .These runs used site-specific soil parameters (see Table 4) but were driven by the same WATCH global meteorological dataset described above, that we assumed to be the best available, continuous and consistent meteorological forcing for these sites.Since JULES does not have a specific tundra plant functional type we assumed C3 grass to be the closest approximation of the vegetation type.At each location the model was first spun up by running repeatedly with the meteorological input of the first year (1959), and then ran continuously to 2001.However, since the soil climate research stations were only established after the mid-1990s, only the latter part of the simulation period was used for comparison with observations.Two sets of simulations were performed, the first using mineral soil parameters only, comparable to the standard setup in JULES and in the climate models.In the second set of simulations the soil parameters were adjusted for organic carbon content, derived for these locations from the NCSCD database.These runs are thus comparable to the SOC spatial model run, although the effect on the soil parameterisation can be expected to be somewhat larger as the organic fractions were not averaged out over a larger grid.

Results -standard JULES
In the following sections we evaluate the performance of JULES in simulating soil temperatures and active layer thickness (ALT) by comparing the results from the GSWP2 and WATCH-GPCC runs with observational datasets.ALT was diagnosed from the simulated soil temperatures by fitting a thermal profile by means of linear interpolation through the midpoints of each soil layer and calculating the depth at which the profile crosses the 0 • C isotherm.Below the midpoint of the bottom layer the position of the thawing front was calculated by extrapolation from above (see Riseborough, 2008).This thawing depth was calculated for each day in the simulation period and the annual maximum thaw depth was used as an indicator of the ALT.

Permafrost extent
Figure 1 shows the extent of permafrost in the JULES simulations compared to observed permafrost extent of the International Permafrost Association (IPA) (Brown et al., 2001).
In both runs there is visually a good agreement between the areas with permafrost in JULES and where permafrost is known to occur.The total permafrost area in GSWP2 (based on areas with an ALT of less than 3 m) is 22.16 million km 2 , which at first sight corresponds well with other estimates of total hemispheric permafrost extent based on the IPA map of circum-arctic permafrost (see Zhang et al., 2003).However, the IPA classifies areas with permafrost into classes with decreasing coverage: continuous permafrost (90-100 % coverage), discontinuous (50-90 %), sporadic permafrost (10-50 %), and isolated permafrost patches (less than 10 % coverage).JULES, on the other hand, does not account for subgrid variability in soil temperature, implying that it should only simulate permafrost if the areal coverage is at least 50 %.With this in mind, the model appears to overestimate the areal extent of permafrost.
The large-scale distribution of areas with permafrost is similar in both model runs, but the mean ALT is generally somewhat deeper in GSWP2 than in WATCH-GPCC.Averaged over the permafrost region the difference is around 0.1m, which may be due to differences in forcing or in parameterisation because of the different grid resolution.Table 5 shows the result of a cross-tabulation of IPA extent classes with JULES mean ALT in GSWP2.In the IPA map the total area of classes with permafrost coverage of at least 50 % is approximately 14 million km 2 .JULES captures about 97 % (13.57× 10 6 km 2 ) of this area.In other words, only about 3 % of the area underlain by continuous and discontinuous permafrost is not identified as such in the model simulation.On the other hand, JULES overestimates the total permafrost area: about 25 % (5.55 × 10 6 km 2 ) of the total area with a simulated ALT of less than 3 m is classified as isolated or sporadic permafrost in the IPA map, and an additional 14 % (3.05 × 10 6 km 2 ) has no permafrost.The latter is partly due to a mismatch in land areas: JULES simulates permafrost in about 0.75 × 10 6 km 2 that is classified as glacier, lake or ocean in the IPA map.(Brown et al., 2001) and JULES mean ALT classes in the GSWP2 run (1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995).The table summarises the total area where each IPA extent class (rows) coincides with each JULES ALT class (columns).IPA areas may differ slightly from previously reported estimates (e.g.Zhang et al., 2003)

Active layer thickness
Observations on active layer thickness have been collected by the Circumpolar Active Layer Monitoring (CALM) since the 1990s (Brown et al., 2000).The network consists of over 100 monitoring sites, most of which are located in arctic and sub-arctic lowlands.There are three primary methods of determining the ALT: by mechanical probing at rectangular grids and/or transects of various size; by employing thawtubes; and by inferring the thaw depth from ground temperature measurements.Here we compare the annual end-ofseason thaw depth from the CALM network3 with the ALT in the JULES simulations.Unfortunately there is only a limited overlap between observations and simulations: WATCH-GPCC ends in 2000, and GSWP2 in 1995, and at many sites few observations, if any, are available before then.
In Fig. 2 we compare the simulated ALT (as inferred from soil temperature) with the observed thaw depth for those locations and years where observations are available.Note that ALT was calculated by linear interpolation of the soil temperature between the layer midpoints, which may lead to considerable bias (cf.Riseborough, 2008).It should also be kept in mind that the JULES runs are relatively coarse-scale (1 • ×1 • in GSWP2 and 0.5 • ×0.5 • in WATCH-GPCC), while the observations are essentially at point-scale or representative of a much smaller area taken by a variety of methods.Also the vertical resolution of these standard simulations is relatively coarse with the bottom layer stretching over 2 m.Soil properties as well as local climatic conditions may be markedly different from the large grid-scale averages used in JULES.Nevertheless, there seems to be a general pattern of the simulated ALT being too deep compared to the observed end-of-season thaw depth (Fig. 2).Averaged over all common data pairs, the mean difference in the GSWP2 run is 0.81 ± 0.48 m, and 0.53 ± 0.50 m in WATCH-GPCC.The root mean square error (RMSE) is 0.94 and 0.73 m, respectively.

Soil temperatures
To evaluate JULES simulated soil temperatures in permafrost regions, we used the Russian Historical Soil Temperature dataset 4 .This dataset is a collection of monthly and annual average soil temperatures measured at Russian meteorological stations.Data were recovered from many sources and compiled by Zhang et al. (2001).Soil temperatures were measured at depths of 0.02 to 3.2 m using bent stem thermometers, extraction thermometers, and electrical resistance thermistors.Data coverage extends from the 1800s through 1990 but is not continuous.At many stations data collection began in the 1930s or 1950s, and not all stations continued to take measurements through to 1990.
In Fig. 3 we compare the observed annual mean soil temperatures with results from WATCH-GPCC that has a longer overlap with the observations, for those years and stations where observations are available.Simulated soil temperatures in each of the four layers are compared with observations pooled for all of the measurement depths that fall 4 Available from http://nsidc.org/data/arcss078.htmlwithin these layers.In other words, no interpolation was applied and this may explain some of the differences that can be seen.Nevertheless, it is obvious that the simulated annual mean soil temperatures are generally too low compared to the observations, and this bias is larger at colder stations.It is also larger in winter than in summer, as can be seen in Fig. 4 where simulations and observations are compared on a monthly basis.In the top soil the bias is largest in mid-winter (December-February) and smallest towards the end of summer (September-October), when simulations are generally in relatively good agreement with the observations.Deeper in the soil this minimum bias is then shifted to later in the year: at 1.6 m depth, the bias is smallest in November.The GSWP2 experiment (not shown), that has a shorter overlap with the observations, fundamentally shows a similar picture although the cold bias appears to be somewhat less than in WATCH-GPCC.
Because of the mismatch in scale in simulations (0.5 • × 0.5 • ) and observations, and the discontinuities in the observations, there is no straightforward single explanation for the cold bias in the simulated soil temperatures.Generally speaking, JULES captures the attenuation and delay of the seasonal cycle in soil temperature reasonably well.The results seem to suggest the top soil layers are cooling down too much in winter, but note that there are only limited observations available for the top soil in the winter months (see also Fig. 3).

Results -modified JULES
Figure 5 shows the difference in mean ALT between JULES in its standard setup (GSWP2, cf.Fig. 1) and the two modifications that were tested: extending the soil profile to 60 m (DEEP) and including modified soil parameters according to soil organic content (SOC).Both modifications appear to have a relatively limited, and sometimes counteracting, effect on the average depth of the ALT.Particularly towards the southern fringes of the permafrost area, including the effect of organic soil sometimes leads to a slight deepening of the active layer.In SOC the average temperature in the top soil layer (0-10 cm) is generally cooler in summer than in GSWP2, as can be expected from the insulating properties of organic material.However, the top soil temperature in SOC is warmer in winter, as winter cold can less easily penetrate into the soil column (Fig. 6).Because of the phase delay in heat transfer deeper into the soil, this leads to somewhat warmer temperatures in spring and early summer in the sub-soil layers, and in summer and early autumn in the bottom layer.In this time of the year the thaw depth is at its maximum and, particularly at lower latitudes, reaches into this deepest layer.The slightly warmer temperatures in this layer thus explain the slightly deeper ALT.In colder regions where the active layer is generally thinner, the cooling of the top soil layers in summer results in a somewhat shallower The Cryosphere, 5, 773-790, 2011 ALT compared to standard JULES, although the difference is mostly less than ∼20 cm (Fig. 5).
In the DEEP run, the bottom layer is -in places -warmer than the standard setup in spring and early summer, but cooler by the end of summer and in autumn, as more heat is being transferred to warm up the added soil layers below (Fig. 6d).In effect, the seasonal cycle is dampened because of the soil layers underneath.In places where summer thawing reaches the bottom layer, the slightly cooler temperature results in a slightly shallower ALT but the differences are generally very small (Fig. 5b).A peculiar effect that can be seen at the boundaries of the permafrost region is that the deeper soil layers allow a deeper ALT to be calculated that in places extends below the original soil profile of 3 m.Therefore, if we apply the same threshold of an ALT below 3 m to identify grid cells with permafrost the actual permafrost area in DEEP is slightly smaller than in GSWP2.
In both SOC and DEEP the differences from the standard JULES setup are, however, too small to remedy the general overestimation of the ALT compared with the observations at CALM sites that was noted earlier.The RMSE, which was 0.94 m in GSWP2, is only marginally lower: 0.93 m in the SOC experiment, 0.90 m in DEEP, and 0.91 m when both modifications are combined.

Results -site-specific runs
To better understand the performance of JULES in cold climate regions we ran a number of point-scale simulations for a number of soil climate research stations in Alaska for which observations on both soil temperature and moisture content were available.These runs were driven by the same WATCH global meteorological forcing datasets that was used in the spatial runs, but used vegetation and soil parameters that were based on the site descriptions (see Table 4).In the following we show results for three stations: Toolik (68.6 • N, 149.6 • W, altitude 759 m); Atqasuk (70.5 • N, 157.4 • W, altitude 22 m); and Barrow-1 (71.3 • N, 156.6 • W, altitude 9 m), the latter being closest to the coast.
A major limitation when running a complex model like JULES at single sites is the availability of meteorological forcing data with sufficient duration, frequency and quality.Observations at micrometeorological tower sites (such as the FLUXNET sites) often contain gaps and errors, particularly in winter.We tried to circumvent this problem by using global forcing datasets but obviously these have limitations as well.This is illustrated in Fig. 7a where the amount of snow (in mm water equivalent, SWE) at Toolik as derived from SSM/I satellite observations5 is compared with the snow mass simulated by JULES driven by the two variants of WATCH -the standard product using the GPCC precipitation data (WATCH-GPCC) and WATCH-CRU that is based on a smaller number of precipitation observations.Clearly, the latter of the two datasets considerably underestimates the amount of cold season precipitation at this location.Over the 2000-2001 winter season, the total precipitation between October and April in WATCH-CRU is only 16 mm or about 10 % of the amount in WATCH-GPCC.As a consequence, the model simulates only a very shallow snow pack that disappears too early in spring.The WATCH dataset using the GPCC precipitation, which is based on a larger number of stations than CRU, performs much better in this respect and the simulated amount of snow is much closer to the satellite-derived SWE, although the snow volume is still lower and disappears about half a month earlier.The latter is at least partly due to the WATCH forcing being somewhat warmer than observed at Toolik in the weeks preceding snowmelt, with daytime temperature frequently climbing above zero degrees Celsius from early May onwards (Fig. 7b).Note that in general there are also uncertainties associated with satellite-based estimates of SWE, especially at high values.The observed snow depth at Toolik at 2 May 2001 (Oberbauer, 2003) amounted to 0.66 m with a standard deviation of 0.09 m.Assuming a snow density in the range of 260-310 kg m −3 (e.g.Sturm and Wagner, 2010) this translates into a SWE of 172-205 ± 28 mm.The SSM/I estimate on this day is 181 mm but day-to-day variability in these data is considerable with values as low as 125 mm in the 10 days around the measurement date.JULES, using a prognostic snow density, calculates a snow depth on this day of 0.30 m and an SWE of 115 mm when using the WATCH-GPCC forcing.
These differences in snow mass have an impact on the simulated soil temperatures (Fig. 7c-e) and calculated thaw depth (Fig. 7f).Compared to WATCH-CRU, the deeper snow pack in WATCH-GPCC leads to higher temperatures in the top soil, but primarily so in autumn and winter.In spring, on the other hand, the thin snow cover and its earlier depletion in WATCH-CRU allows the top soil temperatures to rise more quickly from its colder state, while in summer both runs are in close agreement and correspond relatively well with the observed soil temperatures (Fig. 7c, d).In the subsoil, GPCC is generally warmer than CRU except for about 2 months in (late) summer (Fig. 7e).In the deepest soil layer the difference between GPCC and CRU persists throughout the year but is smallest in late autumn.Note that in WATCH-GPCC the top soil temperatures in winter are still lower than observed while the air temperature (which is the same in both WATCH variants) is mostly in good agreement with the observations at this site (Fig. 7b).
In Fig. 7f, the observed thaw depth has been inferred from observed soil temperature profiles.The ALT that is calculated in this way is deeper than the reported end-of-season thaw depth in the CALM database, which in Toolik is based on mechanical probing in a 1×1 km grid (Hinkel and Nelson, 2003).In 2001, the average ALT obtained by this method over 94 points amounted to 0.46 ± 0.14 m.The ALT www.the-cryosphere.net/5/773/2011/The Cryosphere, 5, 773-790, 2011 inferred from temperature measurements on the day of probing is ∼0.75 m.Although this is considerably deeper, it still falls within the total range of the observations (0.21-0.88 m), highlighting the high spatial variability in active layer thickness, which reflects the local influence of vegetation, substrate properties, snow cover dynamics and terrain (Hinkel and Nelson, 2003).
In spite of the much colder soil temperatures in winter, the maximum thaw depth in the WATCH-CRU simulation is considerably deeper than in WATCH-GPCC (1.04 m vs. 0.67 m, respectively).The reason for this is the much higher soil moisture content in the GPCC simulation requiring more energy and time to melt over summer.As can be seen in Fig. 7e, the third model layer (0.35-1.00 m) remains isothermal at 0 • C throughout summer in the simulation with GPCC precipitation, while in WATCH-CRU the frozen moisture fraction in this layer melts completely allowing the temperature to rise a couple of degrees above zero.This illustrates how the soil thermodynamics in JULES are coupled to the hydrodynamics, and how uncertainties in the precipitation input affect the simulation of soil temperature and consequently the calculation of the thawing depth.
In Fig. 8 we compare simulated soil temperature and moisture content with observations at Atqasuk for the year 2000.Note in these plots the meteorological forcing in the two simulations is the same (WATCH-GPCC) but the soil parameter values are different: the standard simulation assumes homogeneous mineral soil parameters, in the SOC simulation these are adjusted for organic carbon content.The JULES soil moisture shown in Fig. 8 is the unfrozen moisture fraction rather than the total soil moisture, as this was assumed to be directly comparable to the observations.The soil moisture measurements are given in water fraction by volume (wfv), meaning that at values around 0.40-0.45 the soil is fully saturated, dependent on the porosity.JULES soil moisture was converted to wfv using the porosity values of the soil type used in the simulation.For those layers where observations are available (layer 2: 10-35 cm, and layer 3: 35-100 cm) they suggest the soil is close to saturation throughout summer up to a depth of at least 50 cm (Fig. 8e, f).No deeper observations are available but the temperature measurements indicate that below this depth the soil remains frozen.In the model, the second soil layer thaws completely and the unfrozen moisture fraction in summer reflects the total moisture content.This layer is too dry in summer, especially in standard setup using mineral soil parameters.In layer 3 on the other hand, the total moisture content in both simulations is close to saturation throughout the year but most of it remains frozen.In both layers the absolute amount of soil water is larger in the SOC experiment that has a higher porosity than the mineral soil in the standard setup, affecting the simulation of soil temperatures.Because the larger amount of water requires more time to melt, the temperature in layer 2 remains isothermal at 0 • C for a longer period of time, which subsequently affects the penetration of the thawing front into the soil (Fig. 8h).However, the overall effect of the SOC parameterisation on the deeper soil temperatures, as well as the maximum thaw depth that is reached, appears to be minimal.The reported end-of-season thaw depth in the CALM database for Atqasuk in 2000 (averaged over 105 observations) amounts to 0.43 ± 0.20 m.Calculated from observed soil temperature, the maximum seasonal thaw depth is 0.45 m.In the JULES simulations, this is 0.63 m and 0.58 m with the standard and SOC setups, respectively.The JULES ALT estimates are therefore considerably deeper but still within one standard deviation from the average across the CALM observation grid.
Similar patterns can also be observed in Barrow (Fig. 9), where observations on soil moisture content are also available for the top soil.Both JULES simulations capture the temporal variation in top soil moisture reasonably well, but appear too dry compared to the measurements (Fig. 9b).An interesting feature that can be observed at both Atqasuk and Barrow is that the top three model layers are cooling down more rapidly in autumn than observed.The observations remain isothermal close to 0 • C for an extended period of time, in some years until November, suggesting the freezing of soil water is slower than simulated by the model.The opposite sometimes happens in spring, when the model remains close to melting longer than the observations, especially in the SOC runs that have higher moisture content (cf.Fig. 8c).Throughout the rest of the year the difference in soil temperature between the two model setups remains very small.In Fig. 9h the maximum thaw depth in the two simulations (standard: 0.63 m, SOC: 0.55 m) is also close to the ALT derived from soil temperature measurements (0.61 m).In the CALM database, the observed ALT at Barrrow in 1997 ranges between 0.21 and 0.75 m with an average of 0.39 ± 0.09 m.

Discussion
When evaluating JULES with respect to its ability to represent large-scale characteristics of permafrost, it is important to note that in the above runs no model calibration or "finetuning" was applied.All large-scale model experiments that have been discussed in this paper are based on physical process descriptions and soil parameters that are derived from publicly available global datasets, comparable to how JULES is used within climate model experiments.The meteorological forcing data that were used, albeit partly based on model reanalyses, arguably are the best available estimates of historical weather conditions at these large scales.Visually, the area of permafrost in JULES (cf.Fig. 1) compares very well with the IPA permafrost map; in other words, JULES simulates permafrost where it is known to occur.However, it can be argued that the model overestimates the total permafrost extent, as it simulates permafrost also in those areas where it occurs sporadically (areal coverage less than 50 %) or even only in isolated patches.In a way one could say JULES simulates the potential upper limit of large-scale permafrost occurrence rather than the actual coverage that also depends on local conditions.
Consistent with this overestimation of the total permafrost extent is a general cold bias in soil temperatures that was found when comparing JULES with observations at a large number of stations in Russia (Figs. 3 and 4) and, to a lesser extent, also in Alaska (Figs. 7-9).Especially in winter the simulated soil temperatures are considerably lower than observed, something that has also been found in other LSMs (e.g.Nicolsky et al., 2007). PaiMazumder et al. (2008) on the other hand, found that the simulated soil temperatures in the fully coupled Community Climate System Model (CCSM) were mostly too high in winter when compared with the Russian station data, which was attributed to an overestimation of winter precipitation and consequently snow depth by CCSM.One possible explanation for the cold bias found in JULES is therefore that the snow cover in the model provides too little thermal insulation against low winter temperatures.The simulations at Toolik in Alaska (Fig. 7) demonstrate that at this location the low soil temperatures in winter can at least partly be remedied with better precipitation input yielding a deeper snow pack.In this respect it is worth noting that solid winter precipitation measurements are highly uncertain (Goodison et al., 1998;Yang and Woo, 1999).The GSWP2 and WATCH-GPCC datasets used in our model runs both correct for gauge undercatch but some bias is still likely over many areas, resulting in a snow pack that is often too thin.Alternatively, the effective thermal capacity of the modelled snow layers might be too low and/or the thermal conductivity too high.Further evaluation of the snow model in JULES and its ability to represent the hydrology of high-latitude regions is therefore required (see also Haddeland et al., 2011).
Uncertainties in the precipitation input also affect soil temperatures by changing the moisture content of the soil.In JULES and similar LSMs the soil thermodynamics are tightly coupled to the hydrology, and this is especially important where phase changes are involved.Differences in soil moisture particularly influence the propagation of the thawing front in summer, and because of the relatively coarse vertical resolution of the deeper soil layers, small differences in temperature around the melting point may have a large impact on the calculated thaw depth (cf.Fig. 7).
When comparing JULES with observed thaw depths at the CALM sites, the simulated ALT was found to be generally too deep, on average by about 80 cm in the GSWP2 run and 53 cm in WATCH (Fig. 2).At first sight this finding may appear contradictory to the cold bias in simulated soil temperatures, as it suggests the soil temperatures in the model are too warm, allowing the thawing front to penetrate too much into the frozen soil.It should be kept in mind though that the simulated temperature is in much better agreement with the observations in summer and, deeper in the ground, in autumn.The depth of the active layer provides a measure of the cumulative thermal history of the ground surface during the summer thaw period and is highly sensitive to land-atmosphere exchanges.When the ALT is determined from observed temperature profiles, the difference with the simulations is more consistent with the differences in soil temperature.More so, the thaw depth in the model also depends on soil moisture content and the associated phase changes.As long as the frozen moisture fraction has not melted completely, the corresponding model layer remains isothermal around 0 • C and the thaw depth that is calculated from the soil temperature profile tends to be stable around the mid-point of that layer.This effect can be seen quite clearly in some of the station plots (cf.Figs.7f and 8h) and explains the clustering of simulated ALT around 0.67 and 2.0 m (the mid-points of the third and bottom layer, respectively) that can be observed in Fig. 2.This raises the question whether a better simulation of thaw depth and ALT can be achieved by adopting a higher vertical resolution in the model.Riseborough (2008) has shown that linear interpolation of soil temperature between layer points to estimate the position of the 0 • C isotherm may give errors in the order of 30 % of the spacing in a simulation with temperature dependent thermal properties.This suggests that in the bottom soil layer of JULES the relatively coarse vertical resolution may result in an overestimation of the ALT of around 0.4 m.To test this, we performed an additional simulation that is equivalent to the standard www.the-cryosphere.net/5/773/2011/The Cryosphere, 5, 773-790, 2011 GSWP2 run (i.e. with homogeneous mineral soil parameters) but in which the soil column was divided into 30 layers of 10 cm each.Because of the additional computation cost such a setup is unlikely to be adopted in climate model simulations but in offline experiments this approach may be feasible.The results of this experiment are summarised in Fig. 10.Compared to the standard setup of four layers, the higher vertical resolution leads to a more differentiated ALT.Deeper thaw depths are simulated mainly towards the southern fringe of the permafrost but further north and especially in large parts of Siberia the ALT becomes shallower.At the CALM sites (Fig. 10c) the average difference with the reported end-of-season thaw depth reduces from 0.81 ± 0.48 m to 0.58 ± 0.40 m (RSME from 0.94 to 0.71 m), but a considerable bias still remains.A higher vertical resolution is thus beneficial for the simulation of thaw depth in permafrost regions but cannot solve completely the general tendency to overestimate the ALT.In either case it should be kept in mind that uncertainties are very large since we are comparing point-scale observations, reflecting local scale climatological conditions and soil properties, with grid-scale simulations driven by global datasets.As noted before, the smallscale variability in ALT can be very high resulting in a wide range in the reported thaw depth at those CALM sites where observations are made in a spatial grid.High-quality and high-frequency observations (including micrometeorological measurements) from a variety of sites in permafrost areas are therefore paramount to be able to validate and improve complex land surface models like JULES and should be a key focus of ongoing and future research programmes in the Arctic.
Although JULES is not a specific permafrost model, it is important for the model to capture the freeze/thaw status of the soil and the near-surface permafrost, especially in climate model experiments aiming at including the potential feedback effects from permafrost thaw.Therefore it is important to improve the realism of the simulations in permafrost regions.Previous studies have found considerable improvements when including the effect of organic soils and a deeper soil column in LSM simulations (Nicolsky et al., 2007;Lawrence and Slater, 2008;Rinke et al., 2008).Here we found the effect on the simulation of the ALT rather limited, and in many places counteracting each other.Our approach of including organic soil is similar to that of Lawrence and Slater (2008), but different from Rinke et al. (2008) who prescribed a pure organic layer on top with a depth varying according to land surface type.Both studies found a reduction in ground temperatures especially in summer and changes in soil moisture content.In the SOC experiment we also find a lowering of the soil temperatures in the top three layers of the model in summer but mostly warmer temperatures in winter, as the lower thermal conductivity provides a better protection against the winter cold.On an annual basis, the net effect is therefore fairly small, the difference usually being less than 1 • C, with slightly colder temperatures in the top soil layers and mostly somewhat warmer conditions in the subsoil (layer 3) and bottom layer of the model.For comparison, Lawrence and Slater (2008) mention up to 2.5 • C cooler annual mean soil temperatures, Rinke et al. (2008) found a reduction in ground temperature by 0.5 to 8 • C.
The limited influence of soil organic material that we find is partly the result of regridding the original NCSCD to a rather coarse 1-degree grid, with the consequence of organic content, which locally can be very high, being averaged out over a much larger area.As a result, the organic fraction in most grid boxes is considerably less than the assumed maximum of 130 kg m −3 , especially in the deeper soil layers, and the final soil properties are still, to a large extent, influenced by the mineral soil.In our simulations the organic fraction in the top soil was nowhere higher than 0.7, and averaged over the domain covered by the NCSCD more in the order of 0.3.A better approach would therefore be to allow for sub-grid variation in soil properties, similar to the sub-grid vegetation tiles that are already used in JULES.Allowing for such subgrid variability would give the opportunity to represent soils with a higher organic content or even purely organic soils, and it would also enable JULES to estimate fractional permafrost coverage within a given grid box.Alternatively, one could opt for adopting a higher resolution in the simulations in order to better represent the spatial variability of the landscape.However, in relatively complex LSMs like JULES there is always a trade-off to be made between the level of detail that can be included in representing the land surface, and the computation time it takes to run the model.This is especially important when used as a land surface scheme in climate model simulations.
Over all, we find that JULES has some skill in simulating the large-scale characteristics of permafrost but further work is necessary to reduce the cold bias in the simulated soil temperatures, especially during winter.For a better simulation of the ALT a higher vertical resolution and extended soil profile will be beneficial.Future model development should ensure a dynamic coupling of soil organic carbon content and soil thermal and hydraulic properties (Falloon et al., 2011), as well as allowing for sub-grid variability and uncertainty in soil properties.A separate peat module for JULES is already under development in order to study peatland carbon dynamics in the boreal zone, and work is underway to couple JULES to a more advanced model of carbon and nitrogen turnover in both mineral and organic soils (Smith et al., 2010).

Conclusions
When driven by observation-based climatology, JULES is able to represent the large-scale distribution of circumpolar permafrost reasonably well.The model simulates permafrost where it is known to occur and captures more than 95 % of the continuous and discontinuous (more than 50 % spatial coverage) permafrost.However, the total extent appears to be overestimated as JULES also simulates permafrost in areas where the spatial coverage is sporadic (less than 50 %) or even in isolated patches only.Consistent with this we find a general cold bias in the simulated soil temperatures when compared with observations, especially in winter.This may partly be the result of biases in the model forcing data.Uncertainties in the precipitation input affect the simulation of soil temperatures in two ways: by affecting the thickness of the snowpack and therefore the amount of thermal insulation in winter; and by changing the amount of water in the soil which, because of the energy required for phase changes, affects the thermodynamics of soil layers close to the freezing point of water.In turn this affects the simulation of thaw depth and ALT.Generally speaking the model appears to overestimate the annual maximum thaw depth, which can at least partly be explained by the relatively coarse vertical resolution in its standard setup, with the bottom layer spanning over 2 m.However, uncertainties in the observations are large as ALT is highly variable over small scales.

Fig. 1 .
Fig. 1.(a) Observed permafrost coverage, as compiled by the International Permafrost Association (IPA) (Brown et al., 2001); (b) and (c) Areas with permafrost in the GSWP2 and WATCH runs, respectively, with mean active layer thickness (ALT) over the period indicated.

Fig. 2 .
Fig. 2. Observed and simulated (calculated from soil temperature profiles) active layer thickness (ALT) at CALM observations sites.Values are shown for the period 1990-1995 for the JULES-GSWP2 run, and 1990-2000 for the JULES-WATCH run.

Fig. 3 .
Fig. 3. Simulated and observed mean annual soil temperatures at Russian meteorological stations.Modelled soil temperatures are from the WATCH run (1958-2000) for the level indicated; observations are from depths falling within each corresponding model level.

Fig. 4 .
Fig. 4. As Fig. 3 but averaged over all months and stations to show the mean annual cycle.

Fig. 5 .
Fig. 5. Difference in mean ALT over 1983-1995 between the GSWP2 run (standard JULES) (cf.Fig. 1) and (a) SOC (JULES with soil parameters modified according to organic content); and (b) DEEP (JULES with soil profile extended to 60 m depth).Note positive values (red colours) signify the mean ALT is deeper in the modified setups than in GSWP2, negative values (blue) means it is shallower.

Fig. 7 .
Fig. 7. Simulated snow mass, air and soil temperature and thaw depth at Toolik, Alaska in 2001 using two versions of the WATCH forcing data (based on GPCC and CRU precipitation) compared with observations: (a) snow water equivalent (SWE) derived from SSM/I satellite data (available from the ArcticRIMS project); (b) air temperature at 2.0 m (WATCH forcing) and 1.5 m (observed); (c) soil temperature in the top model layer (0-10 cm), (d) second layer (10-35 cm), and (e) third layer (35-100 cm); (f) thaw depth or active layer thickness (ALT) calculated from the soil temperature profiles.Soil temperature observations are indicated by their depth of measurement.

Fig. 10 .
Fig. 10.(a) Mean ALT over 1983-1995 in a GSWP2 run with the soil column divided in 30 layers of 10 cm (cf.Fig. 1); (b) difference in ALT with the standard setup of four layers; and (c) comparison of simulated ALT from both runs with observations at CALM sites (cf.Fig. 2).Colours in (b) are as in Fig. 5 but note the scale is different.

Table 1 .
Brooks and Corey (1964)lues for fine, medium and coarse soil used in theBrooks and Corey (1964)soil hydrology scheme.F c is clay fraction; F st is silt fraction; F s is sand fraction, all dimen-

Table 2 .
Parameters for organic soil used in the SOC experiment.For explanation of the parameters and units, see Table1.

Table 3 .
JULES Northern Hemisphere grid experiments used in this paper.

Table 4 .
JULES point simulations at USDA-NRCS soil climate research sites.

Table 5 .
Cross-tabulation of permafrost extent classes in the IPA permafrost map because of regridding and the different land/sea mask in the JULES modelling grid.Values are in million km 2 .