Impact of model developments on present and future simulations of permafrost in a global land-surface model

Abstract. There is a large amount of organic carbon stored in permafrost in the northern high latitudes, which may become vulnerable to microbial decomposition under future climate warming. In order to estimate this potential carbon–climate feedback it is necessary to correctly simulate the physical dynamics of permafrost within global Earth system models (ESMs) and to determine the rate at which it will thaw. Additional new processes within JULES, the land-surface scheme of the UK ESM (UKESM), include a representation of organic soils, moss and bedrock and a modification to the snow scheme; the sensitivity of permafrost to these new developments is investigated in this study. The impact of a higher vertical soil resolution and deeper soil column is also considered. Evaluation against a large group of sites shows the annual cycle of soil temperatures is approximately 25 % too large in the standard JULES version, but this error is corrected by the model improvements, in particular by deeper soil, organic soils, moss and the modified snow scheme. A comparison with active layer monitoring sites shows that the active layer is on average just over 1 m too deep in the standard model version, and this bias is reduced by 70 cm in the improved version. Increasing the soil vertical resolution allows the full range of active layer depths to be simulated; by contrast, with a poorly resolved soil at least 50 % of the permafrost area has a maximum thaw depth at the centre of the bottom soil layer. Thus all the model modifications are seen to improve the permafrost simulations. Historical permafrost area corresponds fairly well to observations in all simulations, covering an area between 14 and 19 million km2. Simulations under two future climate scenarios show a reduced sensitivity of permafrost degradation to temperature, with the near-surface permafrost loss per degree of warming reduced from 1.5 million km2 °C−1 in the standard version of JULES to between 1.1 and 1.2 million km2 °C−1 in the new model version. However, the near-surface permafrost area is still projected to approximately half by the end of the 21st century under the RCP8.5 scenario.

Abstract. There is a large amount of organic carbon stored in permafrost in the northern high latitudes, which may become vulnerable to microbial decomposition under future climate warming. In order to estimate this potential carbon-climate feedback it is necessary to correctly simulate the physical dynamics of permafrost within global Earth system models (ESMs) and to determine the rate at which it will thaw.
Additional new processes within JULES, the land-surface scheme of the UK ESM (UKESM), include a representation of organic soils, moss and bedrock and a modification to the snow scheme; the sensitivity of permafrost to these new developments is investigated in this study. The impact of a higher vertical soil resolution and deeper soil column is also considered.
Evaluation against a large group of sites shows the annual cycle of soil temperatures is approximately 25 % too large in the standard JULES version, but this error is corrected by the model improvements, in particular by deeper soil, organic soils, moss and the modified snow scheme. A comparison with active layer monitoring sites shows that the active layer is on average just over 1 m too deep in the standard model version, and this bias is reduced by 70 cm in the improved version. Increasing the soil vertical resolution allows the full range of active layer depths to be simulated; by contrast, with a poorly resolved soil at least 50 % of the permafrost area has a maximum thaw depth at the centre of the bottom soil layer. Thus all the model modifications are seen to improve the permafrost simulations.
Historical permafrost area corresponds fairly well to observations in all simulations, covering an area between 14 and 19 million km 2 . Simulations under two future climate scenarios show a reduced sensitivity of permafrost degradation to temperature, with the near-surface permafrost loss per degree of warming reduced from 1.5 million km 2 • C −1 in the standard version of JULES to between 1.1 and 1.2 million km 2 • C −1 in the new model version. However, the near-surface permafrost area is still projected to approximately half by the end of the 21st century under the RCP8.5 scenario.
Permafrost is of interest not only because of the physical effects of permafrost thaw, but because it contains large quantities of stored organic carbon, approximately 1300-1370 Pg (Hugelius et al., 2014), which may be released to the atmosphere in a warming climate. This could have a significant climate feedback effect, which needs to be included in global Earth system models (ESMs) in order to account for the full carbon budget in the future (Koven et al., 2011;MacDougall et al., 2012;Burke et al., 2012;Schneider von Deimling et al., 2012. In order to simulate permafrost carbon feedbacks, the land-surface components of ESMs should include both an appropriate carbon cycle and a representation of the physical dynamics. The amount of carbon released from the soil is strongly dependent on the physical state of the ground -the temperature of the permafrost and the rate at which it thaws (Schuur et al., 2009;Gouttevin et al., 2012b). It is therefore important that the physical dynamics of permafrost are addressed and thoroughly evaluated in models before carbon cycle processes are considered.
There were major problems with the permafrost representation in the majority of the CMIP5 global climate model simulations . In many cases, permafrost processes were not represented, and the frozen land area in many of the climate model simulations differed drastically from the observations. Koven et al. (2012) show that there is little difference in the 0 • air temperature isotherm between models, suggesting that the differences are mainly caused by the land-surface dynamics rather than by the driving climate. Several land-surface schemes have since been modified to better represent processes that are important for permafrost, for example by including soil freezing, soil organic matter and improving the representation of snow (Beringer et al., 2001;Gouttevin et al., 2012a;Ekici et al., 2014;Paquin and Sushama, 2014). This paper demonstrates the impact of adding new permafrost-related processes into JULES (Joint UK Land Environment Simulator; Best et al., 2011;Clark et al., 2011), the land-surface scheme used in the UK Earth system model (UKESM). Although the scarcity and uncertainty of global data on permafrost limit the detail with which it can be represented in a large-scale model like JULES, it is possible to capture the broad spatial patterns of permafrost and active layer thickness (ALT) and to realistically simulate presentday conditions. Chadburn et al. (2015) describe in detail the relevant developments within JULES. These include the effects of organic matter, moss, a deeper soil column and a modification to the snow scheme. Chadburn et al. (2015) also show how these developments impact model simulations at a high Arctic tundra site. This paper now applies them to large-scale simulations, showing that they improve the model performance on a large scale, and significantly impact the simulation of permafrost under future climate scenarios.
These developments result in a more appropriate representation of the physical state of the permafrost -a necessary precursor to considering the permafrost carbon feedback.

Standard model description
JULES is the stand-alone version of the land-surface scheme in the Hadley Centre climate models Clark et al., 2011) and was originally based on the Met Office Surface Exchange Scheme (MOSES) (Cox et al., 1999;Essery et al., 2003). It combines a complex energy and water balance model with a dynamic vegetation model. JULES is a community model and is publically available from http://www. jchmr.org/jules. The work discussed here uses a JULES version 3.4.1 augmented with improved physical processes.
JULES represents 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 at a point or over a grid and requires temporally continuous atmospheric forcing data at frequencies of 3 h or greater. Each grid box can contain several different land covers or "tiles", including a number of different plant functional types as well as non-vegetated tiles (urban, water, ice and bare soil). Each tile has its own surface energy balance, but the soil underneath is treated as a single column and receives aggregated fluxes from the surface tiles.
Recently a multilayer snow scheme has been adopted in JULES (described in 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. This scheme significantly improves simulations of winter soil temperatures in the northern high latitudes (Burke et al., 2013). In the old zero-layer snow scheme, the insulation from snow was incorporated into the top layer of the soil. This scheme is currently still used when the snow depth is below 10 cm. Snow in the zero-layer scheme has a constant thermal conductivity that is added in series to the conductivity of the top layer of soil. In the multilayer snow model, thermal conductivity is parametrised as a function of snow density. Snow albedo is parametrised as a function of snow grain size .
The subsurface temperatures are modelled via a discretisation of both heat diffusion and heat advection by moisture fluxes. The soil thermal characteristics depend on the moisture content, as does the latent heat of freezing and thawing. A zero-heat-flux condition is applied at the lower boundary. The soil hydrology is based on a finite difference approximation to the Richards' equation (Richards, 1931), using the same vertical discretisation as the soil thermodynamics (Cox et al., 1999). JULES uses the Brooks and Corey (1964) relations to describe the soil water retention curve and calculate The Cryosphere, 9, 1505-1521, 2015 www.the-cryosphere.net/9/1505/2015/ hydraulic conductivity and soil water suction. The soil hydraulic parameters are calculated according to Cosby et al. (1984). The default vertical discretisation is a 3 m column modelled as four layers, with thicknesses of 0.1, 0.25, 0.65 and 2 m. The land-surface hydrology (LSH) scheme simulates a deep water store at the base of the soil column and allows subsurface flow from this layer and any other layers below the water table. Topographic index data are used to generate the wetland fraction and saturation excess runoff (Gedney and Cox, 2003).

Recent model developments
Recent developments of permafrost-related processes in JULES are described fully in Chadburn et al. (2015). This development work builds on previous studies of these processes in land-surface models (for example Beringer et al., 2001;Dankers et al., 2011;Paquin and Sushama, 2014). The implementation of these processes within JULES is briefly highlighted below.

Extended soil depth and resolution
Firstly, the number and resolution of the soil layers were increased, a functionality already available in JULES. The soil column was extended from 3 to 10 m, with 14 layers in the top 3 m and a further 14 layers in the lower 7 m, giving 28 layers in total. This is a high number compared with other models, since it was our intention to simulate a well-resolved soil (for example the maximum for the CMIP5 models is 23 layers for a 10 m soil column in GFDL-ESM2M). To make sure that all of the freeze-thaw dynamics would be captured, 10 m was chosen as a large value.
Secondly, a subroutine was added to represent bedrock. When this process is switched on in JULES, the bottom boundary of the ordinary soil column is joined on to a further column in which only thermal diffusion occurs. The heat flux across the bottom boundary of the ordinary soil column is now no longer 0, and the bedrock temperatures are modelled via a discretised heat diffusion equation. The purpose of this is partly to make a deeper soil column more computationally tractable, as hydrology and freeze-thaw dynamics form a large part of the computational load and these processes do not take place in the bedrock layers.
The number and thickness of bedrock layers are set by the user when running the model. In this study, the bedrock column was run with 100 layers of 0.5 m each, making a 50 m column, thus bringing the total soil column up to 60 m. There is a zero-heat-flux condition at the base of the bedrock column, which in future could be changed to a geothermal heat flux.
In fact, soil is often shallower than 10 m and the bedrock should start at varying depths depending on the spatial location, but initial tests showed that starting the bedrock at 3 m instead of 10 m made a very small difference compared with the impact of including a deep heat sink or not, and the same has been shown in other models, such as in .

Organic soil parametrisation
The model uses an improved implementation of the organic soil properties that were first introduced by Dankers et al. (2011). A vertical profile of soil carbon is prescribed for each grid cell (see Eq. 1) and the soil properties are calculated accordingly for each model level.
For some of the properties the organic fraction was used to provide a linear weighting of organic and mineral characteristics (as in Dankers et al., 2011). However, the saturated hydraulic conductivity, dry thermal conductivity and saturated soil water suction were calculated using an appropriate non-linear aggregation. As a result, the organic components of the dry thermal conductivity and saturated water suction have a larger effect than if they were calculated via a linear weighted average.
The Dharssi et al. (2009) parametrisation of soil thermal conductivity was extended to take account of organic soils using a modified relationship between saturated and dry thermal conductivity.

Moss layer at surface
In order to include the insulating effects of mosses, the thermal conductivity of the top soil layer was modified to account for their presence. The thermal parameters for the moss layer are based on Soudzilovskaia et al. (2013). These are also consistent with purely organic soils. It is assumed that the water potential in the moss layer is in equilibrium with that of the top soil layer. At present, hydrological processes within the moss are not explicitly represented in JULES.

Change to snow scheme
In the original multilayer snow scheme, numerical stability requires that the layered snow is only used when the snow depth is 10 cm or greater, and the old zero-layer snow scheme is used for shallower snow. The modification introduced in Chadburn et al. (2015) allows the multilayer snow scheme to run with arbitrarily thin layers, thus removing the zero-layer snow scheme from the model altogether.
In the zero-layer snow scheme the heat capacity of snow is neglected and melt water is passed directly to the soil model to be partitioned into infiltration and runoff. In the multilayer snow scheme the snow is treated as a separate layer with its own heat capacity, and a fraction of the mass in a snow layer can be retained as liquid water instead of passing straight into the soil model. This water will freeze if the layer temperature falls below 0 • C. Thus the snow mass will be slightly different in the multilayer scheme, and in general the model's behaviour is more realistic.

Applying model developments on a global scale
The following sections describe how the spatial distribution of organic matter and moss was determined for a large-scale simulation. The depth of the soil column was fixed across the globe, although there is scope for further improvement to this, for example using a spatially variable depth from soilsurface to bedrock, as in recent work by Paquin and Sushama (2014).

Global organic matter distribution
The organic fractions were calculated from a combination of the Northern Circumpolar Soil Carbon Database (NCSCD) (Hugelius et al., 2013) where available and the Harmonized World Soil Database (HWSD) (FAO/IIASA/ISRIC/ISS-CAS/JRC, 2012) for the rest of the land surface. These databases include some rather limited information about the vertical distribution of soil carbon. Using this, an approximate vertical profile of soil carbon was prescribed for each grid cell (see Eq. 1) and the soil properties calculated accordingly for each model level.
Organic carbon quantities can be obtained from both the NCSCD and HWSD data sets for the top 30 cm (C 30 ) and top 1 m (C 100 ). The profile of soil carbon was assumed to be a constant plus an exponential term. The total for the top 1 m adds up to the observed amount.
when z < 3 m, and C(z) = 0 otherwise. This form of the profile is based on the generic profiles in (Harden et al., 2012) (Fig. 2). It assumes that an exponential distribution of carbon is appropriate and that there is no carbon below 3 m. In reality some carbon will be found below 3 m, but it is not likely to have a great impact on the soil properties, which are somewhat uncertain anyway for the deeper ground. Figure 1 shows profiles generated using this method for a warm soil grid cell and a high-latitude grid cell.

Dynamic moss
There are no data sets showing the pan-Arctic distribution of mosses. In addition in a changing climate the distribution of moss may also change. Therefore, moss was implemented in JULES so that it can be run either with a static map that is input at the start of the run or with dynamic growth determined by environmental conditions in the model.
In order to determine the presence of moss in any grid cell, the model takes account of the local temperature, moisture, light, snow cover and, in some cases, wind speed (see Table 1). The moss cover is then determined by a "health" variable, whose value is updated once a day depending on the conditions over the past 24 h. Good conditions add to health and bad conditions subtract from it. It is contrained within bounds that result in maximum health within a year given optimum growing conditions. The conditions for good and poor growth are given in Table 1. The water suction is taken as the minimum of water suctions in the top soil layer and the atmosphere, the temperature, T s , is that of the top soil layer, and the light is the radiation at the bottom of the canopy. These values are chosen for being closest to the soil surface where moss is located.
The temperature, moisture and light ranges for good growth are based on the values in Proctor (1982). The light saturation and compensation curves (L sat and L comp respectively) were estimated from Proctor (1982) and are given by where light compensation is the level at which photosynthesis balances respiration, and light saturation is the level at which photosynthesis is no longer light limited (increasing radiation levels no longer increase the rate of photosynthesis). The heat and moisture conditions that cause the moss to die off are also taken from Proctor (1982). As well as dying in very hot or dry conditions, moss can suffer badly from wind damage when it is frozen (Longton, 1982). Thus the model includes a third scenario in which moss dies off: when it is cold, windy and there is no protective snow cover.
Snow protects moss from harsh conditions in winter, but of course it cannot actually grow under deep snow, so a small value is subtracted from the health under deep snow. The same occurs when it is too dark or too cold for photosynthesis to take place. Moss growth has been observed up to 2 weeks before the end of snowmelt (Collins and Callaghan, 1980), so the threshold value for growth under snow was assumed None of the above. * Must not simultaneously satisfy snow mass < 0.1 kg m −2 , wind speed > 12 m s −1 and temperature < −5 • C.
based on snow mass values in JULES 2 weeks before the end of snowmelt. The magnitude of the values added to and subtracted from the health variable were calibrated at several sites where the moss cover was known.
When the moss health is positive, it is taken to have maximum cover in the grid cell. When the health becomes negative, the percentage cover drops off linearly to 0, and the cover is 0 for the lowest quartile of health values. If there is other vegetation present, the fractional cover of moss is capped at 0.7 for those vegetation tiles. This value was pragmatically chosen given that moss can have around 50-100 % cover in forests (Beringer et al., 2001). Moss cover is assumed to be 0 for the urban or ice fraction of a grid box. The relationship of moss health to moss cover would benefit from further calibration. Figure 2 shows a moss distribution simulated by JULES in the northern high latitudes and compares it with data from the Euskirchen et al. (2007) land-cover map. In general the most densely moss-covered areas correspond to the tundra land-cover classes, which are shown in bright green. Our scheme gives some general representation of this low vegetation cover, which is otherwise missing in JULES. Moss also grows in the boreal forest (shown in purple on the Euskirchen map), but in JULES it does not grow in the deciduous needleleaf zone, which may require some investigation. Evaluating the distribution in lower latitudes is a subject for future work.

Historical meteorological forcing data
The Water and Global Change (WATCH) project produced a meteorological forcing data set (WATCH Forcing Data Era-Interim, WFDEI) for use with land-surface and hydrological models (Weedon et al., 2010(Weedon et al., , 2011Weedon, 2013). It is based on Era-Interim reanalysis data (ECMWF, 2009), with corrections generated from Climate Research Unit (Mitchell and Jones, 2005) and Global Precipitation Climatology Centre data (http://gpcc.dwd.de). It covers the time period 1979-2012 at half-degree resolution globally and at 3 hourly tem-  2007). "Shrub tundra" includes prostrate, dwarf shrub and low shrub tundra, and "boreal forest" includes boreal evergreen needleleaf and boreal broadleaf deciduous. Lower plot: mean moss cover simulated in JULES for the year 2000, from orgmossD historical simulation (see Table 2).  vided, including a combined precipitation variable rather than separate rain and snow, for two future scenarios -RCP4.5 and RCP8.5. As described in the protocol for the permafrost carbon MIP, these anomalies are combined with historical data either by addition (temperature, pressure, humidity, wind speed) or multiplication (short-wave and longwave radiation, precipitation). The anomalies for air temperature and precipitation are shown in Fig. 3. There is not a large trend in precipitation, although in RCP8.5 there is a small increase, along with an increase in variability. Air temperature, however, shows a clear increase by the end of the century, up to 10 • C in RCP8.5. The change in air temperature is much larger in winter than in summer.
For the simulations discussed in this paper, the anomalies were re-gridded to 0.5 • resolution and applied to a repeating sequence of 8 years of WFDEI reanalysis data (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). Using anomalies rather than directly using climate model data makes the climate variability in the future simulations more consistent with the historical data. The main disadvantage is that small-scale features are not captured, since submonthly variability comes from the base data set (Hempel et al., 2013).

CALM Network
The Circumpolar Active Layer Monitoring Network (CALM) (Brown et al., 2000(Brown et al., , 2003) is a network of over 100 sites at which on-going measurements of the end-of-season thaw depth (the ALT) are taken. Measurements are available from the early 1990s, when the network was formed. The data are available from http://www.gwu.edu/~calm/data/north.html.

Historical soil temperatures
The Russian historical soil temperature data set is described in Frauenfeld et al. (2004). Soil temperatures were measured at 242 stations, over different time periods starting as early as 1890. The measurements used in this paper were taken at depths of 0.2, 0.4, 0.8, 1.6 and 3.2 m using extraction thermometers, with additional measurements at 0.6, 1.2 and 2.4 m. At some of the sites the natural vegetation cover was removed and at others there is some possibility of site disturbance, however the majority of these measurement sites retained their natural vegetation and snow cover.
International Polar Year Thermal State of Permafrost (IPY-TSP) borehole inventory data were compiled in 2007-2009 from both new and existing boreholes, achieving a wide spatial coverage of soil temperature data . Data are available in the most part from 2006 to 2009 at a daily resolution, with temperatures measured at a variety of depths.

Globsnow
The European Space Agency snow water equivalent (SWE) product, GlobSnow (Takala et al., 2011) covers the years 1978-2010 and is available on an EASE-Grid at 25 km resolution. GlobSnow is produced using a combination of satellite-based microwave radiometer and ground-based weather station data. The GlobSnow documentation (Luojus et al., 2013) records an exponentially increasing bias with larger snow mass, which is particularly significant at snow mass greater than 200 kg m −2 .

Permafrost distribution data
The Circum-Arctic map of permafrost and ground-ice conditions (Brown et al., 1998) gives a historical permafrost distribution, which can be compared with permafrost area in the model. The data set contains information on the distribution and properties of permafrost and ground ice in the Northern Hemisphere (20-90 • N), with gridded data available at 12.5, 25 km and 0.5 • resolution. It records continuous, discontinuous, sporadic and isolated permafrost regions, for which the estimated permafrost area is 90-100, 50-90, 10-50 and < 10 % respectively. In this work, an estimate of the observed permafrost area is used to compare with the simulated area. For the maximum we assume that permafrost would be simulated in the whole of the continuous and discontinuous area, plus a fraction of permafrost in the areas of sporadic permafrost and isolated patches -the fractions being 0.05 and 0.3 respectively. For the minimum we assume that no permafrost would be found in the sporadic and isolated regions and that in the continuous and discontinuous zones only a fractional coverage of permafrost is found: 0.95 and 0.7 respectively for the two zones. This gives a maximum area of 17.0 million km 2 and a minimum of 13.7 million km 2 . It is important to note that there is considerable uncertainty in these data which means that the true value could fall outside of this estimate.

Simulation set-up
For the historical period, JULES was driven by the WFDEI reanalysis data at a resolution of 0.5 • (Sect. 2.4.1). JULES was driven by precipitation (sum of rain and snow in WFDEI), which was converted internally within the model to rain and snow. This maintains consistency with the future driving data. The model was spun up for 60 years by repeating the first 10 years of driving data (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988) and then run over the period 1979-2009 for the historical runs. After spin-up, the soil temperature and moisture contents were fully stabilised at the vast majority of points, i.e. there was no residual model drift. The future runs began at the start of 2006, taking their initial state from 2006 in the historical simulations and simulating both RCP4.5 and RCP8.5 scenarios until 2100.
In order to capture all the main permafrost regions in the Northern Hemisphere, the simulations were run for the region north of 25 • . The mineral soil properties, land-cover fractions and topographic index data (needed for LSH, see Sect. 2.1) were taken from HadGEM2-ES ancillary data (Collins et al., 2008). Organic soil ancillaries were generated using the method described in Chadburn et al. (2015), with the spatial distribution as described in Sect. 2.3.1.
The simulations carried out are given in Table 2. This includes the standard JULES set-up (min4l), a higherresolution soil column (min14l), a deeper soil column (minD), the effects of moss and organic soils both separately and in combination (minmossD, orgD, orgmossD) and finally the modified snow scheme (orgmossDS). When deeper soil is added in minD, this includes both the extension of the soil column to 10 m and the addition of a 50 m bedrock column. The distribution of moss was determined dynamically in the model as described in Sect. 2.3.2.

Evaluation methods
The maximum summer thaw depth or ALT was calculated by taking the unfrozen water fraction in the deepest layer that has begun to thaw and assuming that this same fraction of the soil layer has thawed. This gives significantly more precise estimates of the ALT than temperature interpolation (see Chadburn et al., 2015).
The ALT was then used to derive the near-surface permafrost extent. Our definition of near-surface permafrost is a grid cell with ALT less than 3 m for 2 or more consecutive years. There is no representation of sub-grid heterogeneity in the soils so any 0.5 • grid cell either contains 100 % nearsurface permafrost or no near-surface permafrost at all. Koven et al. (2012) calculated a range of metrics for soil temperature dynamics, which are also used here. They include the offset between the mean air, 0 and 1 m temperatures and the attenuation of the annual cycle between each of these levels. The values were calculated as in Koven et al. (2012) by first calculating the annual mean and seasonal cycle at the available depths and then interpolating these to 0 and 1 m depth. The annual cycle was interpolated between layers by assuming it falls off exponentially with depth. The model values were taken from the grid cell containing the measurement point. All IPY-TSP sites with sufficient data were used, along with the colder Russian soil temperature sites (those with a mean soil temperature below 0 • C as an indicator of permafrost). This gives a total of 86 sites with reasonable circumpolar coverage (see Sect. 2.4.4 for description of observational data).
In order to analyse future permafrost degradation, the sensitivity of near-surface permafrost area loss to climate warming was calculated via a linear regression of near-surface permafrost area against the annual mean air temperature over the historical permafrost region (region defined by the observed map in Brown et al., 1998).

Active layer and near-surface soil temperatures in historical simulation
In Figs. 4 and 5, the simulated ALT from the JULES simulations (Sect. 2.5) is compared with observations from the CALM active layer monitoring programme (Sect. 2.4.3). The low resolution in min4l significantly impairs the capacity of the model to simulate the active layer. This was seen for the point site simulation in Chadburn et al. (2015) and is even clearer in these large-scale results: the active layer in min4l has very little variability and little apparent correlation with the observations (see Fig. 5a).
Most of the new model processes reduce the active layer, bringing it into much better agreement with the observations, as shown in Fig. 4. Simulating a deeper soil column reduces the active layer mean for the CALM sites by 0.12 m (min14l to minD). The insulating effects of organic matter and moss have a greater impact, reducing the ALT by a further 0.58 m (orgmossD). The inclusion of organic matter has the single greatest effect, reducing the mean ALT by 0.44 m. The whiskers indicate the most extreme data point that is no more than 1.5 times the IQR. Outliers are not shown. Points at which either the simulated or observed active layer was very large (greater than 6 m) were removed. The model point for each CALM site is the grid box containing that site. cesses, the full range of ALT values are captured and the points fall around the one-to-one line. There is still an outlying block of points where the active layer in JULES is greater than 3.5 m and much deeper than the measurements. Many of these sites fall along the course of the Mackenzie River in northern Canada (where JULES simulates very little permafrost -see Fig. 7). Precipitation gauges are sparse in this region so there may be large uncertainties in hydrol-ogy (Weedon et al., 2011), and the observed soil temperatures vary greatly from around −1 to −7 • C as a result of various influences including land cover and snow (Burn and Kokelj, 2009). This is a subject for future investigation. The active layer thickness is determined by both the annual cycle of soil temperatures and the thermal offset between the air and the soil. Table 3 compares these dynamics in JULES with observations from the IPY-TSP data set and cold sites from the Russian soil temperature data set (see Sects. 2.4.4 and 2.6). The root mean squared error (RMSE) is calculated using the mean value of the metric for each site, so it quantifies the extent to which the variability between the sites is correctly simulated. In this table the most relevant values are the offset and attenuation of the annual cycle between the air and 1 m depth in the soil, since the soil surface is not so well defined in the observations. In the standard JULES set-up (min4l) the total offset is approximately correct, suggesting that the mean soil temperatures are simulated well. However, the annual cycle is nearly 25 % too large at 1 m depth.
The introduction of a well-resolved soil has a cooling effect of approximately 0.7 • C, and organic soils and moss have a further cooling effect of approximately 1 • C. As in Chadburn et al. (2015), the improvements to the snow scheme then compensate for the cooling from the other model changes, resulting in a 1 m ground temperature approximately the same in the final simulation (orgmossDS) as the original one (min4l). There is no significant improvement in the RMSE between the first and last simulation (quantifying the extent to which we capture the spatial variability), but the RMSE is between 2 and 2.5 • C for all simulations, suggesting that the mean soil temperatures are captured fairly accurately.
The annual cycle, however, is reduced overall by the model improvements, with the attenuation value in orgmossDS being very close to the observed value (less than 2 % different to be exact). The RMSE in these values is also reduced (by approximately 20 %) by the model developments, showing that spatial patterns are better simulated. This shows a significant improvement in the simulation. Figure 6 shows the mean annual cycle of soil temperatures at 90 cm for the sites used in Table 3. Figure 6a shows that the main effect of increasing soil resolution is to reduce the winter temperatures. This is very likely because of the way the snow is simulated in the standard snow scheme, where the insulating effect of shallow snow is incorporated into the top soil layer -this means that when the top soil layer is thinner the insulation will be less effective. The deeper soil gives a slight attenuation of the annual cycle, which is expected of an additional heat sink. Figure 6b shows that the main effect of organic soils and moss is to reduce the summer soil temperatures. Finally, the main effect of the modified snow scheme is to increase the winter temperatures. Overall, a reduction in summer temperatures and an increase in winter temperatures lead to a re-The Cryosphere, 9, 1505-1521, 2015 www.the-cryosphere.net/9/1505/2015/ Table 3. The attenuation of the annual cycle and the thermal offset in the JULES simulations, between the air and the top of the soil, and the top of the soil and 1 m depth, calculated as in Koven et al. (2012). This includes IPY-TSP and Russian data for cold sites. The RMSE (root mean squared error) values are based on the mean value of the metric at each site and thus give an indication of how well the variability between sites is captured. The bold numbers mark the observations as distinct from simulated values. duced annual cycle which matches better with the observations in Fig. 6c. There are still significant differences between the model and observations, shown by the size of the RMSE in Table 3. The error in winter snow depth when compared with the closest grid cells in the Globsnow data set (see Sect. 2.4.5) has a significant correlation of 0.3 (for about 260 points) with the error in winter soil temperatures in orgmossDS, suggesting that snow explains at least some of the remaining error. Langer et al. (2013) show that uncertainty in the simulated active layer depth comes from the uncertainty in soil composition, particularly ground-ice contents. Some variability is also expected when comparing a large grid cell with a point site, and this is difficult to quantify. While the RMSE in off-  simulations, 1979-1989. All grid cells with active layer ≤ 6 m are shown, with mean ALT indicated by colour. Bottom: observed permafrost map (Brown et al., 1998), based on maps made between approximately 1960 and 1990. On all plots the 0 • air temperature isotherm is shown in red (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)  set is reasonably small, the RMSE for the attenuation of the annual cycle is a significant fraction of the value itself. This suggests that the annual cycle is more difficult to simulate than the mean temperature.
In this section a comparison with CALM observations has shown that one essential feature for simulating the ALT is the resolution of the soil column, without which the active layer variability is not resolved (see Fig. 4). For capturing the near-surface soil temperature dynamics, the effects of organic soils and the improved snow scheme are both particularly significant in making the annual temperature cycle more realistic (Fig. 6 and Table 3). Organic soils also greatly improve the ALT (Figs. 4 and 5), so this process is a particularly useful addition to the model. The importance of organic soils has also been shown in e.g. Rinke et al. (2008); ; Koven et al. (2009).

Permafrost distribution in historical simulation
The simulated permafrost in JULES is shown in Fig. 7, along with observations from the Circum-Arctic map of permafrost and ground-ice conditions (Brown et al., 1998) (Sect. 2.4.6). The observed map shows areas with continuous, discontinuous and sporadic permafrost and isolated patches. There is no equivalent of discontinuous permafrost in JULES because each grid box has only a single soil column, so in order to compare the two maps we assume that a deeper active layer in JULES may correspond to discontinuous or sporadic permafrost. With this assumption, all the simulations match the observations fairly well in most areas. We can see that introducing the model developments brings in much more spatial variability in ALT, which generally matches with the patterns of continuous/discontinuous permafrost. The correlation between the ALT in JULES and the percentage cover of permafrost from (Brown et al., 1998) (100 % for continuous, 90 % for discontinuous, 50 % for sporadic and 10 % for isolated patches) is high, ranging between −0.37 and −0.51.
However, there are still places where continuous permafrost is observed but JULES does not simulate permafrost. Figure 8 shows that in most of these areas, JULES simulates far too much snow, which will mean too much insulation in winter leading to soils that are too warm. This is particularly noticeable in north-east Canada and two areas in north-west Russia. In north-east Canada, however, it has been shown that the GlobSnow data set underestimates the SWE (Langlois et al., 2014), so the over-estimation in JULES may not be as large as Fig. 8 suggests. However, the permafrost in this region is unstable to thawing (Thibault and Payette, 2009), so a small bias in the model could make the difference between simulating permafrost or not. For most of the remaining land surface, JULES slightly underestimates the SWE. Hancock et al. (2014) showed that JULES generally underestimates SWE when driven by reanalysis data sets.
In the observations, the edge of the permafrost-affected zone corresponds quite closely with the 0 • isotherm, al- though there is a gap in western Russia (red lines in Fig. 7). In the JULES simulations this relationship is less consistent, suggesting more spatial heterogeneity in the relationship between air and soil temperatures. For example, excessive snow cover such as that seen in north-east Canada in Fig. 8 could contribute to this effect.
It is also possible to consider the vertical distribution of permafrost. Figure 9 shows a breakdown of active layer depths for all near-surface permafrost points (a) and all points (b) in the simulations. The first thing that is clear from this plot is that the discretisation of the soil column has a very large effect on the simulation of permafrost. A series of kinks corresponding to the model discretisation is apparent in all the curves, which for the higher-resolution simulations does not significantly impact the overall shape of the curve, but for the low-resolution soil changes it almost beyond recognition. For about 50 % of the near-surface permafrost points in min4l, the active layer depth is between 1.8 and 2.2 m, where 2 m is the centre of the bottom model layer. Figure 9 shows that the vertical distribution of permafrost is affected by all the model improvements, but the most significant impact is when organic soils and moss are included. Here, the permafrost is generally found nearer to the surface.
In Fig. 9b the amount of near-surface permafrost in each simulation is apparent from the fraction of points that thaw to less than 3 m (generally about 30 % of points). This shows that some simulations have a shallower active layer (so colder soils) but less near-surface permafrost; for example, compare min14l with min4l and orgmossDS with minD. This is related to the vertical profile of soil temperatures. For example, due to the combination of processes, the maximum soil temperature in orgmossDS compared with minD is colder near the surface but warmer in the deeper soil. This is seen in the soil temperature profiles in Fig. 10. The total near-surface permafrost area in each simulation is more clearly seen by looking at the historical period in Fig. 11. For consistency with the definition of near-surface permafrost, these values include any grid cells where the ALT was less than 3 m for the past 2 years (so for example in the first year that a grid cell is frozen, it is not included in the permafrost area but after the second year it is added to the area, so the area changes from year to year). Comparing the deep-soil simulations, minD, minmossD, orgD and orgmossD, we see that adding insulation from organic soils and moss increases the near-surface permafrost area, which is consistent with the cooling effect seen in Table 3. In orgmossDS, the near-surface permafrost area is significantly reduced compared with orgmossD, which was also apparent in Fig. 9b. Finally, in the shallow (3 m) simulations (min4l and min14l), the near-surface permafrost area is smaller, but this is not really meaningful. This is because the zero-heatflux boundary condition is not correct at 3 m, which leads to "edge effects" close to the soil boundary. This is seen in Fig. 10, where the soil temperatures in the mineral soil simulations (min4l, min14l and minD) are very similar near the top of the soil, but the annual cycle in the shallow simulations (min4l and min14l) does not continue to fall off with depth, resulting in a maximum temperature that is much too high at the base of the soil. This shows that diagnosing permafrost as the area with ALT less than 3 m requires a soil column significantly deeper than 3 m.
A range for the observed permafrost area is also shown in Fig. 11. This is estimated from Brown et al. (1998) using assumptions described in Sect. 2.4.6. According to this, the simulated near-surface permafrost area in the mineral soil simulations (min4l, min14l, minD) falls inside the observed range, and the addition of organic soils and moss results in a simulated near-surface permafrost area that is somewhat too large. However, in the final simulation with all model improvements (orgmossDS), the near-surface permafrost area once again falls within the observed range.

Future permafrost degradation
Section 3.1 showed that the model improvements make the simulation of permafrost more realistic. Although further development is needed, these are important processes to include and it is worthwhile to study their impact on longterm permafrost dynamics; hence in this section we study the loss of near-surface permafrost and the active layer deepening over the next century. Figure 11 shows the time series of total near-surface permafrost area over the next century. Comparing minD and orgmossD shows that organic soils and moss reduce the interannual variability of the near-surface permafrost area. Although this variability cannot be measured on a global scale, permafrost tends to degrade or form over a number of years, so a high level of variability from year to year is probably unrealistic. In orgmossDS, although the near-surface permafrost area is significantly reduced compared to orgmossD, the interannual variability is similar. However, the shallow (3 m) simulations (min4l and min14l) have a much higher interannual variability in near-surface permafrost area than the www.the-cryosphere.net/9/1505/2015/ The Cryosphere, 9, 1505-1521, 2015  deep soil simulations, indicating the importance of the thermal inertia from a deep heat sink in the soil, which has been shown already in e.g. Stevens et al. (2007); Alexeev et al. (2007); . Table 4 shows the rate of near-surface permafrost loss per degree of warming (calculated by a linear regression between future mean air temperature, averaged over the historical permafrost region, and permafrost area). The sensitivity is reduced by the new model developments, from approximately 1.5 × 10 6 km 2 • C −1 in the standard JULES set-up (min4l) to between 1.1 × 10 6 and 1.2 × 10 6 km 2 • C −1 in orgmossDS: a reduction of about 25 %. In  the rate of permafrost loss in the Community Land Model is reduced by over 25 % by the inclusion of organic matter and a deeper soil column, which is an even larger effect than is found in JULES. In that study an even deeper soil column down to 125 m was used.
The loss of near-surface permafrost by the end of the 21st century is very large in all the JULES simulations, particularly in RCP8.5. Even in orgmossD and orgmossDS, the simulations with the lowest sensitivity to temperature, the area with near-surface permafrost has approximately halved by the end of the 21st century in RCP8.5 (see Fig. 11). In orgmossDS it decreases from approximately 16 million km 2 for the historical period down to approximately 7 million km 2 in RCP8.5 and 12 million km 2 in RCP4.5 by 2100.
In terms of area lost, the values for the standard JULES set-up are approximately the same as for HadGEM2-ES. However, the fractional loss in HadGEM2-ES is smaller, because the near-surface permafrost area itself is significantly larger (22.3 million km 2 for the historical period compared with approximately 15 million km 2 in min4l). This is predominantly because HadGEM2-ES uses the zero-layer snow scheme, leading to significantly colder soils. This suggests that the snow scheme does not have a great effect on the actual rate of permafrost degradation in JULES, which is supported by comparing orgmossD and orgmossDS in Table 4. A study of historical permafrost in JULES (Burke et al., 2013) showed a loss of 0.55-0.81 million km 2 per decade. In our historical simulations the loss rates generally fall within this range except for orgmossDS, for which the mean near-surface permafrost loss per decade is slightly lower at The Cryosphere, 9, 1505-1521, 2015 www.the-cryosphere.net/9/1505/2015/ 0.43 million km 2 . Compared with the other CMIP5 models, the rates of near-surface permafrost loss in the improved version of JULES (orgmossDS) are now lower than most of the other models Slater and Lawrence, 2013). Figure 12 shows the near-surface permafrost distribution at the end of the future simulations for the "standard" JULES set-up (min4l) and the final improved model version (orgmossDS). In RCP4.5 the near-surface permafrost retreats only from the edges of the permafrost zone, but there is some thickening of the active layer across the whole near-surface permafrost area of the order of 0.5 m, which is a significant change. In the more southern permafrost regions, the active layer deepens more in orgmossDS -as much as 1 m -which may reflect the fact that the ALT is initially shallower, so more deepening is possible.
In RCP8.5, much of the near-surface permafrost thaws in both simulations. However, in the improved model version, there is significantly more near-surface permafrost remaining at the end of the century, particularly in northern Russia, reflecting the reduced sensitivity to warming. However, the areas where near-surface permafrost remains in orgmossDS show a strong active layer deepening, significantly more than 1 m in some areas. This suggests that although less nearsurface permafrost is lost in this simulation, with a further in-crease in temperature it could disappear. Note that this refers just to permafrost at depths of up to 3 m -deeper permafrost may remain longer.

Conclusions
Large-scale simulations have shown improved physical permafrost dynamics in JULES, thanks to a deeper and betterresolved soil column, including the physical effects of moss and organic soils and an improvement to the snow scheme. The model developments reduce the simulated summer thaw depth and the amplitude of the annual cycle of soil temperatures and bring both to more realistic values. The rate of nearsurface permafrost loss under future climate warming is also reduced, as is the interannual variability of the near-surface permafrost area.
It is important to simulate a reasonable ALT before beginning to consider permafrost carbon feedback. For this we have shown that the depth and resolution of the soil column and the effects of organic soils are the most important considerations for the model. JULES is now able to simulate largescale patterns in ALT, as seen in Fig. 5, where points with shallower ALT are now generally simulated with shallower ALT in the model. www.the-cryosphere.net/9/1505/2015/ The Cryosphere, 9, 1505-1521, 2015 A well-resolved soil is absolutely essential for simulating active layer dynamics. With a poorly resolved soil it is not possible to simulate ALT variability: the thaw depth depends strongly on the model layers, Fig. 9, and there is very little spatial variability in ALT (Fig. 7), which is unrealistic (see Fig. 5). The importance of soil resolution has not often been emphasised in the literature.
We have confirmed the importance of including a deep soil column, showing that the thermal inertia from deeper soil has a significant impact on the permafrost dynamics (see e.g. Figs. 4 and 11). We also find that with a 3 m soil column it is not possible to meaningfully diagnose permafrost at 3 m depth due to the edge effects close to the bottom boundary of the soil (Fig. 10).
Organic soils and moss both act to insulate the soil. They reduce the ALT and increase the near-surface permafrost area (see for example Fig. 4 and Table 4) due to their insulating effect in summer, which helps to make the annual cycle of soil temperatures more realistic (Fig. 6b). Of these two processes, organic soils have the larger effect. Including moss has a smaller but significant physical impact, as seen for example by the increase in near-surface permafrost area of 0.7 million km 2 when it is included (Table 4). High-latitude vegetation is an important component that is currently missing in JULES. This simple insulating layer is only the first step and more work is needed to fully incorporate it into the vegetation model.
The improvement to the snow scheme has a large winter warming effect (Fig. 6c), and as a result it warms the deep soil and reduces the near-surface permafrost area, bringing it closer to the observational estimate. It also contributes to a more realistic annual cycle (Fig. 6c). However, the snow is less important for simulating the correct ALT (Fig. 5).
Snow has a very strong effect on soil temperatures Ekici et al., 2015;Langer et al., 2013), and there is certainly more scope to improve the snow model in JULES, for example with additional compaction processes and lateral redistribution. However, there are high uncertainties associated with global precipitation data in observations, reanalysis products and climate model output, particularly for the Arctic (Bosilovich et al., 2008;Kattsov et al., 2007;Hancock et al., 2014), so these limit the potential to improve the snow representation at present.
Soil moisture is also important for soil temperatures, and the two are linked in a complex manner. Water has a higher thermal conductivity than air, so in wetter soils more heat will penetrate and leave the ground. However, if there is soil freezing and thawing, the latent heat will reduce the rate of heat penetration, counteracting this effect. Furthermore, if there is soil freezing the mean temperature in the deeper soil is colder than that at the surface, since the thermal conductivity of ice is greater than that of water, so more heat moves upwards in winter than downwards in summer. The influence of soil temperature on soil moisture mainly comes from freez-ing, as this prevents moisture running out of the soil and may also hold liquid water above a permafrost layer.
Between the runs in this paper, the main differences in soil moisture come from organic soils, which increase the soil moisture content overall. Disentangling the complex impacts requires more specific experiments than the ones in this study, such as an experiment where specific influences of soil moisture on temperature are removed. Further investigation of the hydrology in JULES is vital and this is the subject of ongoing work.
Important improvements have been made in JULES and other global land-surface models (e.g. Gouttevin et al., 2012a;Ekici et al., 2014), but these models are now somewhat limited by their coarse spatial resolution, since permafrost processes are heterogeneous on small spatial scales and have non-linear effects. There is some progress being made towards upscaling small-scale processes (Muster et al., 2012;Langer et al., 2013), and large-scale data sets are improving over time, which is important for simulating realistic carbon fluxes.
In this paper we have analysed the large-scale degradation of permafrost under two future climate scenarios. This shows a significant reduction in near-surface permafrost area, with up to 1.5 million km 2 of near-surface permafrost loss per degree of high-latitude warming, although this is reduced to approximately 1.1 million km 2 in the improved model version, showing the importance of these model developments in assessing future permafrost thaw. The impact of organic matter is particularly large, as this alone reduces the sensitivity by approximately 15 % (Table 4). In RCP4.5, nearsurface permafrost is only lost from the edges of the permafrost zone by the end of the century, but in RCP8.5 the near-surface permafrost disappears entirely from some large regions, with large areas of near-surface permafrost remaining only in northern Canada and some parts of Russia. In areas where near-surface permafrost remains, there is a significant thickening of the active layer, which is relevant for consideration of the permafrost carbon feedback.