An improved representation of physical permafrost dynamics in the JULES land-surface model

. It is important to correctly simulate permafrost in global climate models, since the stored carbon represents the source of a potentially important climate feedback. This carbon feedback depends on the physical state of the permafrost. We have therefore included improved physical permafrost processes in JULES (Joint UK Land Environment Simula-tor), which is the land-surface scheme used in the Hadley Centre climate models. The thermal and hydraulic properties of the soil were modiﬁed to account for the presence of organic matter, and the insulating effects of a surface layer of moss were


Introduction
The northern high latitudes (NHLs) are an important region in terms of the changing global climate. Both observations and future projections of warming are amplified in this region (Overland et al., 2004;Bekryaev et al., 2010;Stocker et al., 2013). At the land-surface scale, significant thawing of permafrost has already been observed in many areas (Camill, make future climate projections and inform emissions targets (Stocker et al., 2013).
In order to include permafrost carbon feedbacks in landsurface models, the first requirement is that the physics is simulated correctly. This includes thaw depth and rate of thaw, hydrological processes and soil temperature dynamics, which all affect soil carbon stocks and decomposition rate (Gouttevin et al., 2012b;Exbrayat et al., 2013).
While permafrost-specific models have made progress towards correctly simulating permafrost dynamics (Riseborough et al., 2008;Jafarov et al., 2012;Westermann et al., 2015), in global land-surface models the Arctic has often been neglected, leading to the large discrepancies between models and reality seen in Koven et al. (2012). One reason that the NHLs are poorly represented in global models is the difficulty of obtaining observations with which to drive and evaluate the models. Harsh conditions in the Arctic mean that much of the land area is difficult to access, and detailed simulations are only possible on small scales. However, the use of small-scale simulations where observations are available can help to improve the large-scale dynamics. Several global land-surface models have already improved their representation of permafrost physics (Beringer et al., 2001;Gouttevin et al., 2012a;Ekici et al., 2014a).
In this paper we add new permafrost-relevant processes into JULES (Joint UK Land Environment Simulator), which is the land-surface scheme in the Hadley Centre climate models and will be used in the first UK Earth system model Clark et al., 2011), improving on the past implementation of these processes (Christensen and Cox, 1995;Cox et al., 1999). We evaluate the model at a site level, where it is reasonable to compare the model directly with observational data, and a large quantity of data is available. Being able to simulate realistically at a site level shows that the physics of the model is correct, which is a prerequisite for trusting large-scale simulations. These developments are included in large-scale simulations in Chadburn et al. (2015).
JULES already includes some of the processes that are important for permafrost: the effects of soil freezing and thawing on the energy budget and, more recently, a multilayer snow scheme, which significantly improves model performance (Burke et al., 2013). However, systematic differences between JULES simulations and reality have been identified. When compared with observations of active layer thickness (ALT) (thickness of seasonally frozen layer), the simulated active layer in JULES is consistently too deep. This is seen, for example, in Dankers et al. (2011), where the simulated active layer was compared with observations from over 100 sites in the CALM (Circumpolar Active Layer Monitoring Network) active layer monitoring programme (Brown et al., 2000). This bias in ALT indicates that the soil may warm too quickly in summer, which would lead to an amplification of the annual cycle of soil temperatures. This amplification is indeed observed in JULES (Burke et al., 2013). This suggests that the model undergoes an accelerated soil warming in summer, meaning either that too much heat enters the soil or too much of that heat accumulates near the surface.
There are two controls on the amount of heat entering and leaving the soil: the land-cover above the soil and the thermal properties of the soil itself. In particular, soil organic matter and the moss layer that is often present in the low Arctic can greatly influence the ALT and summer soil temperatures (Dyrness, 1982). This is because moss and organic matter have insulating properties and can also hold more water than mineral soils. The importance of accounting for organic matter in land-surface models has been discussed in e.g. Rinke et al. (2008), , and Koven et al. (2009). Snow also insulates the soil in winter and has a very large effect on the soil temperatures and permafrost dynamics Langer et al., 2013;Ekici et al., 2014b). Thus, in this model development work, we consider implementing the physical effects of moss and organic matter, and further improving the snow scheme in JULES.
An accumulation of heat near the surface in the model can be related to the heat sink of the deeper part of the soil: if the model does not simulate a deep soil column this heat sink is missing. Several studies have shown that a shallow soil column does not give realistic temperature dynamics (Stevens et al., 2007;Alexeev et al., 2007). Finally, the resolution of the soil column affects the numerical accuracy of the simulation and also the precision to which the ALT can be resolved. The default configuration for JULES represents only the top 3 m of soil with four layers. Therefore, in this work the depth and resolution of the soil column is increased, including a thermal "bedrock" column at the base.
The impact of soil hydrology is also considered, showing that if the soil moisture were simulated correctly the simulations of soil temperature could be further improved. Soil temperatures are affected by the water content of the soil not only through its thermal properties but also via the latent heat of freezing, which slows down the rate of temperature change.
Simulations are performed of the Samoylov Island site in Siberia, adding each model development in turn. This shows the impact of the new processes and significant improvements to model performance and the representation of permafrost in JULES. Areas for future development are also clearly identified.

Model description (standard version)
JULES is a stand-alone land-surface model, which is also used in the Hadley Centre coupled climate models Clark et al., 2011), and was originally based on the MOSES (Met Office Surface Exchange Scheme) landsurface scheme (Cox et al., 1999;Essery et al., 2003). It combines a sophisticated energy and water balance model with a dynamic vegetation model. JULES is a community model JULES simulates 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 a continuous time series of atmospheric forcing data at a frequency of 3 h or greater. Each grid box can contain several different land covers or "tiles", including a number of different plant functional types (PFTs) 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.
JULES uses a multilayer snow scheme (described in Best et al., 2011) in which the number of snow layers varies according to the depth of the snowpack. Each snow layer has a prognostic temperature, density, grain size and water content. 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.
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 of Richards' equation (Richards, 1931), with 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 hydraulic conductivity and soil water suction. Soil hydraulic and thermal parameters are input to the model via an ancillary file. 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 scheme (LSH) simulates a deep water store at the base of the soil column and allows for subsurface flow from this layer and any other layers below the water table. Topographic index data is used to generate the wetland fraction and saturation excess runoff (Gedney and Cox, 2003).
JULES also includes a dynamic vegetation model, TRIF-FID, which simulates vegetation competition to determine the grid-box fraction assigned to each PFT (Cox, 2001). JULES may also be run with TRIFFID switched off and a fixed vegetation fraction, which was the case for the simulations in this paper, where the focus is on the physical processes.

Permafrost model developments
Model developments include the thermal effects of a surface moss layer, the thermal and hydrological effects of soil organic matter, a thermal "bedrock" column beneath the ordi-nary soil, and an improvement of the multilayer snow scheme to allow for arbitrarily thin layers. The resolution and depth of the soil column is also increased. These improvements are described in detail in the following subsections.

Moss
The characteristics of moss will vary between different species and ecosystems, but all mosses will insulate the soil. Therefore, the thermal conductivity of the soil was modified to represent this insulating layer. Its purpose in these simulations is to give a somewhat generic representation of the thin layer of moss-rich vegetation which is abundant in the Arctic. Although any vegetation layer in JULES has an insulating effect thanks to the canopy heat capacity , this new type is necessary because the current PFT's are not appropriate for Arctic tundra.
The thermal conductivity of moss depends on its water content. For simplicity we assume that the moss layer coincides with the top layer of the soil, and thus has the same hydraulic suction. The water content is then calculated from the suction using the Brooks and Corey (1964) where θ is the volumetric water content, ψ is the soil water suction, b is the exponent and the subscript sat refers to the values at saturation. The following hydraulic parameters were used for moss (Beringer et al., 2001): b moss = 1, ψ sat,moss = 0.12 m, θ sat,moss = 0.9. The dependence of moss thermal conductivity on water content was measured by Soudzilovskaia et al. (2013). We choose the representative values for the saturated conductivity of 0.5 Wm −1 K −1 and for dry conductivity 0.06 Wm −1 K −1 , and linearly interpolate between the two depending on the moisture content (Eq. 1). These are also consistent with the values given for organic soils in Williams and Smith (1991).
The user can choose the thickness of the moss layer; the default value is 5 cm. The thermal conductivity of the top 5 cm of soil is then modified according to the parameters above. This is applied to a fraction of the grid-box depending on a variable representing the percentage cover of moss.

Organic soils
Organic soils were previously considered in JULES by Dankers et al. (2011), who concluded that their effects were small. In this paper, however, we use an improved implementation of their impact.
As in Dankers et al. (2011) the volumetric fraction of organic soil, f org , was used to modify the soil properties to include the effects of organic matter. f org was estimated as a vertical profile using observations of organic carbon at different depths. Soil carbon observations were avail-able in kilograms per cubic metre, which were converted to a volumetric fraction using literature values for the density (800 kg m −3 ) and porosity (Dankers et al., 2011, Table 2) of organic matter.
For some of the soil properties, the organic fraction was used to provide a linear weighting of organic and mineral characteristics (Appendix A, Eqs. A1, A4, A7), as in Dankers et al. (2011) and other similar work. However, the saturated hydraulic conductivity, dry thermal conductivity and saturated soil water suction were calculated using a more appropriate non-linear aggregation (Appendix A, Eqs. A2, A3, A8). 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. See Appendix A for details. New soil properties, which include the effects of organic matter, may now be input to JULES via an updated soil ancillary file.
The current parameterisation of saturated thermal conductivity in JULES (Dharssi et al., 2009) restricts the values to those appropriate for mineral soils. Organic soils have a much lower saturated thermal conductivity, so it was necessary to modify the Dharssi parameterisation to take account of this.
The thermal conductivity of dry soil (λ dry ) is input to JULES via the ancillaries, and the saturated thermal conductivity is calculated in the model, depending on the fraction of the soil moisture that is currently frozen. The actual value of thermal conductivity is then calculated by interpolating between the dry and saturated conductivities depending on the water content. The literature values used in JULES for the dry thermal conductivity are 0.25 Wm −1 K −1 for clay soils and 0.3 Wm −1 K −1 for sandy soils, and saturated conductivity of 1.58 and 2.2 Wm −1 K −1 , respectively (for unfrozen soils) (Williams and Smith, 1991, Table 4.1).
For organic soils, the dry conductivity is approximately 0.06 Wm −1 K −1 and the saturated conductivity 0.5 Wm −1 K −1 (Williams and Smith, 1991). However, using Dharssi's method the minimum value for saturated conductivity is 1.58 Wm −1 K −1 . It was therefore necessary to implement a parameterisation of saturated thermal conductivity that extends to the appropriate values, for which a smooth curve was fitted to the data (Appendix A, Eq. A11). The two curves are shown in Fig. 1. The conductivities for mineral soils will be slightly different in the new formulation, but this difference will be small and well within the uncertainty of the literature values.
Note that the same thermal conductivity values are used for both moss and organic soil. This is consistent with the fact that, for example in peat soils, the layer of living moss can be almost indistinguishable from the surface organic layer. One good reason for treating them separately, however, is that moss can also grow in places without a pronounced organic layer. Figure 1. New method of calculating saturated thermal conductivity (λ sat0 ) from dry thermal conductivity (λ dry ), compared with the standard method. See Appendix A, Eq. (A10) for the Dharssi parameterisation and Appendix A, Eq. (A11) for the new method, which is modified to include organic soils.

Increased soil resolution and depth
The hydrologically active soil column is run in three different configurations, beginning with the standard four-layer configuration, with layer thicknesses of 0.1, 0.25, 0.65 and 2.0 m. The second configuration increases the soil resolution without increasing the depth, having 14 layers in 3 m of soil. The layer thicknesses start at 0.05 m and increase with depth according to the function dz n = 0.05n 0.75 . Finally, the highresolution column is extended to 10 m with a total of 28 layers, following the same function for dz n .
In this last case, soil column depth is increased even further by adding an extra column to the base of the hydrologically active column, to represent bedrock. This bedrock column adds another 50 m, bringing the total soil column to 60 m. See Sect. 2.2.4 for details.

Bedrock thermal dynamics
A bedrock column was added to JULES, starting at the base of the hydrologically active soil column. We assume that the heat transfer in this deep column is not influenced by hydrological processes. This allows for the representation of a deep soil column without a large computational load. Heat transfer is simulated by thermal diffusion: where T s,deep is the temperature in the deep soil column, t is time and z is vertical depth. This is discretised to first order as follows: where i indexes the timesteps and n indexes the vertical layers. This uses a constant heat capacity, C deep , and thermal conductivity, λ deep , which may be set by the user. The default values are C deep = 2.1 × 10 6 JK −1 m −3 and λ deep = 8.6 Wm −1 K −1 (the properties of the soil solids in sand from Beringer et al. (2001), and very close to the values for quartz in Williams and Smith, 1991). By default, the vertical layer thickness is dz deep = 0.5 m, with 100 layers, resulting in an extra 50 m soil column, but the user can also set these values. In most models the deep soil is not so finely resolved -in fact it is often represented as a single thick layer, but since the heat diffusion is so computationally light, there is no reason not to resolve the dynamics more accurately.
In the hydrologically active soil column an implicit solution is used for the temperature increments, but for bedrock the explicit solution is sufficient since temperature changes are slow and there are no freeze-thaw processes to consider. The heat flux across the boundary with the base of the hydrologically active soil column is where the thermal conductivity, λ base , is an interpolation between the bottom layer of the hydrological column and the top layer of the bedrock column. Here N is the number of soil hydrological layers, which interface with the bedrock column. The heat flux at the base of the bedrock column is set to zero by default, but could be set to the geothermal heat flux in future versions.

Improved snow scheme
The original release of JULES included the same simple snow model as in the MOSES land-surface scheme (Cox et al., 1999) (Burke et al., 2013), but the old snow model was retained for shallow snow of less than 10 cm depth to avoid numerical instabilities. For this study, a modification has been implemented that allows shallow snow to be represented by a distinct model layer. This is done by calculating the heat flux into the snow or soil surface according to the temperature gradient between the surface and a fixed depth below

Figure 2. Map showing location of Samoylov Island and Northern
Hemisphere permafrost distribution (Brown et al., 1998).
the surface. The snow layer temperature is stepped forward in time using the backward Euler method, which remains stable for an infinitesimal layer thickness.

Samoylov Island site information
Point-scale simulations were carried out using data from the Samoylov Island field site in the Lena River delta, Siberia. Figure 2 shows the location of Samoylov Island in the context of the whole Arctic permafrost region. There is a large quantity of data available from this site, making it a good site for detailed process evaluation (e.g. Yi et al., 2014). The landscape is formed of ice-wedge polygonal tundra with ponds and thermokarst lakes. There is an abundance of mosses and organic soil, so including the model developments described above has a notable impact on the JULES simulations. A typical soil profile is shown in Fig. 3b, highlighting the moss and organic layers. Figure 3a shows an aerial view of the monitoring site, including the meteorological station, soil temperature monitoring and active layer monitoring grid (only polygon centre points are highlighted as these data are used for evaluation). A detailed description of the site may be found in Boike et al. (2013).

Forcing data
The meteorological driving data were prepared using observations from the site combined with reanalysis data for the grid cell containing the site. For the period 1901-1979, water and global change forcing data (WFD) were used (Weedon et al., 2010(Weedon et al., , 2011. This is a meteorological forcing data set based on ERA-40 reanalysis (ECMWF, 2006), with corrections generated from Climate Research Unit (CRU) (Mitchell and Jones, 2005) and Global Precipitation Climatology Centre (GPCC) data (http://gpcc.dwd.de). Data is provided at half-degree resolution for the whole globe at 3-hourly time resolution from 1902 to 2001, with years prior to 1958 based on random years from ERA-40 but corrected with observations from the relevant time period. For the period 1979-2010, WATCH Forcing Data Era-Interim (WFDEI) was used (Weedon, 2013). This is produced using the same techniques as the WFD but is instead based on the ERA-Interim reanalysis data (ECMWF, 2009) and covers the period 1979-2012. For the time periods where observed data were available from Samoylov, correction factors were generated by calculating monthly biases relative to the WFDEI data. These corrections were then applied to the time series from 1979 to 2010 of the WFDEI data. The WFD before 1979 was then corrected to match this data and the two data sets were joined at 1979 to provide gap-free 3-hourly forcing from 1901 to 2010. Meteorological station observations were used for all variables except snowfall, which was estimated from the observed snow depth by treating increases in snow depth as snowfall events with an assumed snow density of 180 kg m −3 . Snow depth observations are available daily from 2002 to 2013, although with some missing years. These reconstructions were then used to provide correction factors to WFDEI and WFD. This leads to a more realistic snow depth in the model than using direct precipitation measurements, due to wind effects and the difficulty of accurately measuring snowfall.

Soil and land-cover characteristics
The land characteristics were chosen to represent a depressed polygon centre, and the evaluation data (soil temperatures, moisture, etc.) were also taken from polygon centre measurements (see Fig. 3a).
The mineral soil is a sandy loam and was assumed to have 50 % silt, 45 % sand, 5 % clay, which is consistent with the information in Boike et al. (2013). The soil properties were calculated using the Cosby et al. (1984) relations. Site-specific organic carbon quantities are given in Zubrzycki et al. (2013), but there is significant heterogeneity, with values for polygon centres ranging between 3 and 85 kg m −3 . The mean values of 25 kg m −3 of organic carbon above 30 cm and 35 kg m −3 from 30 cm to 1 m were used, giving a volumetric fraction f org between 0.4 and 0.6. Following the model set-up used in Langer et al. (2013), organic carbon below 1 m was taken as zero. The transition between carbon quantities above and below 30 cm was smoothed into a curve. Organic properties were then combined with the mineral properties as in Sect. 2.2.2.
To verify this parameterisation of organic soil properties in JULES, we compare the resulting thermal properties with those in Langer et al. (2011a, b). We compare saturated values in JULES with values for saturated peat. In JULES the thermal conductivity is consistent with the Langer values, lying between 0.7 and 0.9 Wm −1 K −1 when thawed and between 1.9 and 2.1 Wm −1 K −1 when frozen. The values from Langer et al. (2011a, b) are 0.72 ± 0.08 (thawed) and 1.92±0.19 Wm −1 K −1 (frozen). The heat capacity in JULES is 3.5-3.8 (thawed) and 2.2-2.3 MJ m −3 K −1 (frozen), which is again close to the Langer values of 3.8 ± 0.2 (thawed) and 2.0 ± 0.05 MJ m −3 K −1 (frozen); although the heat capacity when frozen is a little too high in JULES, this is a reasonable level of consistency given the high spatial variability in soil properties.
The vegetation at Samoylov is composed predominantly of mosses, along with grasses and small shrubs with about 10 % coverage. The land cover in JULES was taken as 10 % grass with a height of 10 cm. Moss cover was set to 90 % (or 90 % bare soil in simulations without moss).
For simulations with higher-resolution and deeper soil, the set-up is described in Sect. 2.2.3. The new bedrock routine was also used (Sect. 2.2.4), adding a further 50 m heat sink to the base of the soil. Samoylov Island sits above a deep river deposit, so the deep soil is composed of silt deposits. The estimated parameters for the deep soil are approximately C deep = 2.1 MJ m −3 K −1 and λ deep = 2 Wm −1 K −1 , so these values were used for the bedrock column in JULES (Boike et al., 2013).
The improved snow scheme was included, along with a change of the fresh snow density in all simulations from the default value of 100 to 130 kg m −3 , to better match the observed snow density and depth. The fresh snow density applies when the snow first reaches the ground, after which it undergoes standard compaction processes (see Best et al., 2011, Eq. 21), meaning that a higher fresh snow density will lead to a higher snow density year-round. In test runs the mean simulated density over all snow-covered periods was around 190 kg m −3 (compared with 165 kg m −3 with the default fresh snow density). It is possible that this estimate is Geosci. Model Dev., 8, 1493-1508, 2015 www.geosci-model-dev.net/8/1493/2015/ still too low, since the observed density for a polygon centre is around 230 kg m −3 . However, 190 kg m −3 is close to the spatial average given in Boike et al. (2013) and this is also approximately consistent with the assumed value of 180 kg m −3 used to create the driving data. This is considered further in Sect. 3.2. The LSH scheme was also switched on (see Sect. 2.1). This scheme adds a deep water store at the base of the soil and thus improves the water-holding capacity of the soil.

Simulation set-up
Simulations were performed first for the standard version of JULES using just mineral soil (min4l). The developments of increased soil discretisation (min14l), deeper soil (minD), organic soil properties (orgD), moss insulation (orgmossD) and the improved snow scheme (orgmossDS) were then systematically introduced (see Table 1), with the final simulation containing all of the model improvements (orgmossDS). The simulations labelled ρ fresh = 170 and Saturated in Table 1 are discussed in Sect. 3.2.
The simulations were spun-up for 200 years using the first 10 years of driving data (starting at 2 January 1901), by which point the soil temperatures and water contents were stable. They were then run from 1901 until the end of 2010.

Calculating active layer thickness (ALT)
Commonly used methods of calculating ALT in land-surface models make use of the soil temperatures, either by taking the depth of the deepest layer that is above 0 • C or an interpolation of soil temperatures to find the depth of 0 • C; see for example Koven et al. (2012) and Lawrence et al. (2012). However, this method is limited by the vertical discretisation. In JULES, when a given layer is freezing or thawing, the temperature of the layer remains very close to 0 • C for the duration of freeze-thaw, with the consequence that any interpolation puts the thaw depth very close to the centre of the layer. However, more information may be extracted from JULES by outputting the frozen and unfrozen water contents in the layer. In this paper, the ALT is 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 is represented by the following equation: where θ f and θ u are frozen and unfrozen water content as a fraction of saturation, and n is the deepest layer that has completely thawed (θ f,n = 0). This gives significantly more precise estimates than the usual temperature interpolation. Figure 4 shows an example of the thawing period in 2006 for one of the JULES simulations (orgmossDS , Table 1), where the thaw begins too early but the maximum depth is well simulated. The temperature interpolation method uses a linear interpolation to find the depth of 0 • C. It is clear that this method produces thaw depth in a series of steps corresponding to the JULES layers. The new method based on fraction of unfrozen soil moisture gives a much smoother curve, which corresponds better to the observations.  Figure 5 shows the simulated ALT at Samoylov over the 11year period 2000-2010. It is clear that many of the new processes in the model reduce the ALT, with the final simulation including deep soil, moss, organic properties and the new snow scheme (orgmossDS) bringing the simulated ALT very close to the observations. In orgmossDS the ALT falls within the range of observations for every year where measurements are available. A significant bias of over 1 m for the standard JULES set-up (min4l) has been removed by including these model developments.

Soil temperatures and ALT
Comparing the first two simulations, min4l and min14l, the mean ALT is reduced by 0.2 m when soil resolution is increased. The base layer in min4l is 2 m thick, and the thaw depth consistently reaches almost to the centre of this layer, although in some years earlier in the simulation (not shown) the thaw reaches only the third model layer and thus the ALT changes by approximately 0.5 m in 1 year, which is unrealistic behaviour.
A comparison of shallow and deep simulations, both with high resolution (min14l and minD in Fig. 5) demonstrates that there is a small but significant reduction in bias when the depth of the soil column is increased to 10 m and bedrock is added. In this case the mean ALT is reduced by 0.13 m.
The addition of the organic soil parameterisation has the single biggest impact on the ALT in these simulations. The mean ALT in minD is 1.03 m, and in orgD it is 0.49 m, a reduction of over half a metre. Moss on its own also has a large impact, reducing the mean ALT by 0.35 m from minD to minmossD. However, the effects are non-linear: comparing orgD and orgmossD, the mean ALT is reduced by only 0.17 m by the addition of moss.
The ALT depends on the maximum soil temperatures in summer, but it is important to simulate the correct soil temperatures for the whole year and the whole soil column. Table 2 shows some key performance metrics for the soil temperatures at different depths, most of which are significantly improved in the final simulation (orgmossDS). This shows that although the changes to the snow scheme have very little effect on the ALT, they reduce the root mean square error (RMSE) in upper soil temperatures from 4.1 to 3.4 • C, which is a significant improvement.
The simulated active layer temperatures in the mineral soil simulations (min4l, min14l and minD), shown in Table 2, are too warm (≈ 2 • C) and the annual cycle is much too large. This is consistent with the large-scale biases in JULES discussed in Sect. 1. The addition of organic soils and moss reduces the mean active layer temperature and annual cycle to more realistic values. The changes to the snow scheme then increase the active layer temperature by 1.2 • C, but the temperature bias is still less than 1 • C and the annual cycle bias is reduced from 30 % to less than 10 % in orgmossDS.
Considering the deep soil temperatures in Table 2, organic soils and moss have a cooling effect, which is offset by a warming due to the change in snow scheme (this is also true at 32 cm depth). The temperature biases are more negative (or less positive) deeper in the soil than they are at the surface, and the annual cycles also have biases that are inconsistent at different depths. This suggests that the profile of soil temperatures is not entirely realistic in the simulations. Further measurements of deep soil properties at the site would allow for a more detailed analysis of this. Figure 6a and b show the active layer soil temperatures in more detail, showing the combined effect of the model improvements. Temperatures are shown for both the whole active layer (Fig. 6a) and for a single slice at 32 cm depth (Fig. 6b). The limitation of the lower-resolution standard JULES set-up (min4l) is clear in Fig. 6a, where the soil temperature changes in a series of steps, whereas it is much smoother in both the observations and the other simulations. Figure 6b shows that the improved model matches the observed soil temperatures much better in summer, and somewhat better in the shoulder seasons (spring and autumn). Comparing minD with orgmossD shows that organic soils and moss have the main impact on summer soil temperatures. Comparing orgmossD and orgmossDS shows that the snow scheme has the greatest effect during the shoulder seasons.
At 32 cm depth, the RMSE in the warmest months (August-September) is reduced from 4.0 • C in the minD Geosci. Model Dev., 8, 1493-1508, 2015 www.geosci-model-dev.net/8/1493/2015/ −9.9 −8.6 −8.9 23 1.5 0.14 - Figure 6. (a) Soil temperatures in active layer, simulated (top three plots) and observed (lower plot). The simulations are, from top, the standard four-layer JULES set-up (min4l), deeper and better-resolved soil (minD), and adding to this organic soils, moss, and the improved snow scheme (orgmossDS). Observations are for a polygon centre (see Fig. 3). (b) Active layer soil temperatures at 32 cm depth, simulated and observed. The lines represent horizonal slices through the contour plots in Fig. 6a. Additionally, the simulation orgmossD is shown, which includes organic soils and moss but not the new snow scheme.
simulation to just 0.7 • C in orgmossDS. This suggests that the most important processes for the summer have been identified and included, namely the insulating effects of moss and organic soils. However, the temperatures in snow-covered seasons are much more difficult to simulate, with the RMSE for the other months reduced from 5.3 • C in minD to 3.9 • C in orgmossDS, which is a significant reduction but not nearly as large as for the summer. One reason for this is that snow varies dynamically on short timescales, which strongly affects the energy balance. In contrast, processes that affect the summer temperatures are relatively static -for example, the organic content of the soil will change very slowly (peat growth of around 2 mm per year is observed at the site). Snow will be considered further in Sect. 3.2.

Snow and soil moisture
The largest remaining errors in soil temperatures in the final simulation (orgmossDS) occur during the winter and shoulder seasons (see Fig. 6b). Figure 7 shows the observed and simulated snow depth over the same time period as Fig. 6b. It is clear that in winter 2003-2004, when the mid-winter Figure 7. Simulated and observed snow depth at Samoylov over the same years as soil temperatures (Fig. 6b). The simulation orgmossDS includes all model improvements (see Table 1).
soil temperatures are simulated fairly accurately, the snow depth is below that observed, whilst in winter 2004-2005 the snow depth is close to the observations but the soil temperatures are too warm. This suggests that the simulated snow density is too low. The snow density determines the thermal conductivity which, combined with the snow depth, is used to calculate the heat flow between air and soil. A further simulation was performed, increasing the fresh snow density even more from 130 to 170 kg m −3 (see Table 1). This increased the mean snow density that was simulated in JULES from around 190 to 220 kg m −3 , which matches more closely with the observational estimate specifically for polygon centres, which is about 230 kg m −3 (Boike et al., 2013). Figure 8 shows the effect of increasing snow density. The soil is now too cold in winter 2003-2004, which is consistent with there being too little snow. In winter 2004-2005, where snow depths are more realistic, the soil temperatures match better with those observed. During the coldest months (January-March) there is a strong correlation of approximately 0.85 between the error in snow depth and the error in soil temperature, for both simulations. However, the linear regression line crosses a long way above the origin in orgmossDS (4.3 • C), whereas when the fresh snow density is higher it passes closer to the origin (1.8 • C) -see Fig. 8. For these months, using ρ fresh = 170 kg m −3 reduces the RMSE in soil temperature from 3.9 to 2.4 • C. However, the wholeyear RMSE in soil temperature is increased from 3.4 to 3.7 • C, mainly because of differences in temperatures in the shoulder seasons, in particular during the freeze-up period in autumn, when the simulated zero-curtain length is too short (the zero curtain is the period for which the soil remains at or close to 0 • C during freeze or thaw). The end of the freezeup happens on average 30 days too early in orgmossDS, and when the snow density is increased it is even earlier, on average 42 days before the observed freeze-up date.
The zero-curtain duration is determined by the latent heat associated with freeze-thaw. In reality, polygon centres tend to be saturated (Boike et al., 2013). If there is not enough soil moisture, some latent heat will be missing, reducing the zero-curtain length. Figure 9 compares the volumetric soil moisture content in the observations and simulations. It is clearly improved in the organic soil simulations (orgmossD, orgmossDS) compared to the mineral soil simulations (min4l, minD), but there is still too little soil moisture, partly because the porosity is too low and partly because the soil does not always stay saturated. The offset timings of freeze and thaw are clearly seen, showing that the timing of thaw is greatly improved in orgmossD and orgmossDS, but there is little effect on the time of the freeze. Note that the unfrozen soil moisture content in winter is higher than the observations, which suggests the hydraulic parameters may need some refinement.
In order to investigate the soil moisture effect, a simulation was performed in which the soil was kept saturated all yearround, and the organic matter content in the upper soil layers was increased to 45 kg m −3 to increase the porosity to match better with the observations (Saturated in Table 1). As a result, the thermal conductivity and heat capacity generally fall within the uncertainty of the values in Langer et al. (2011a, b), although the frozen heat capacity is increased with a maximum value of 2.4 MJ m −3 K −1 , which is 20 % greater than the values in Langer et al. (2011b). The simulation results are shown in Fig. 10. Figure 10 shows a great improvement in the zero-curtain length, with the mean difference between the simulated and observed freeze-up now being only 13 days (instead of 30 in orgmossDS). The overall RMSE for the upper soil temperatures is reduced further from 3.4 to 2.7 • C, and the deep soil temperatures are also improved (see Saturated in Table 2). The summer soil temperatures are actually a little warmer and in fact this reduces the RMSE for August-September temperatures (at 32 cm) slightly from 0.7 to 0.6 • C. Figure 10a shows that the time series of unfrozen soil moisture is also greatly improved. These improvements highlight the need for more work on soil hydrology.
The simulation where the soil is held saturated is now a great improvement on the original mineral soil simulation, with RMSE at 32 cm reduced by almost half. However, the zero-curtain still falls short by nearly 2 weeks, and the mid-winter temperatures still differ significantly from observations, especially for winter [2003][2004]. The difference of winter temperatures is likely due to the driving data, which do not always result in the correct snow depth in the model nor the timing of snowfall and snowmelt (see Fig. 7). In reality, snow depth at Samoylov does not correspond very closely to snowfall, as it depends strongly on wind redistribution. This is a difficult problem to solve for a 1-D model such as JULES. The problem of zero-curtain duration, however, may also be related to the snow density, and there is scope to improve this in the model. According to the parameterisation used in the Crocus snowpack model, which depends on temperature and wind speed (Vionnet et al., 2012), the fresh snow density should be much lower than 170 kg m −3 for this Geosci. Model Dev., 8, 1493-1508, 2015 www.geosci-model-dev.net/8/1493/2015/  Table 1). The lower plot compares the error in soil temperatures and snow depths for the coldest months only (January-March) using daily values. site. This would give more snow insulation during the freezeup period, but the simulated mid-winter snow density in JULES would then be too low. This could be addressed by including compaction processes in the model that are currently not represented, such as wind compaction and temperaturegradient metamorphosis, both of which are potentially important (Sturm and Holmgren, 1998;Vionnet et al., 2012).

Conclusions and future work
Improvements have been made to the physical representation of permafrost in the JULES land-surface model. Additional processes represented include an insulating moss layer, the physical properties of organic soil, and a bedrock column. In addition, the representation of snow and discretisation of the soil have been modified. These developments are extremely relevant for the Arctic in general, since soils in the continuous permafrost zone are often organic-rich and covered by moss, which is certainly the case at Samoylov Island, where we run the model simulations. It is therefore important to include these processes in global land-surface models.
In the simulations, soil temperatures and ALT are significantly improved by the model developments. Firstly, increasing the model depth and resolution is necessary to correctly simulate the physical processes. It has been shown that a shallow soil column cannot give realistic permafrost dynamics (see e.g.  and a high enough resolution is required to correctly solve the physical equations. Once this basic function of the model has been improved, including the new, permafrost-relevant processes of organic soils and moss leads to a great improvement in summer soil temperatures. The RMSE in summer soil temperatures decreases from 4.0 to 0.7 • C, and the ALT reduces by 0.7 m to fall within within 0.1 m of the observations. This suggests that the most important processes for the summer have been included. In the shoulder seasons, the zero-curtain duration is strongly related to soil moisture. This requires further work in JULES, as the model does not obtain the saturated conditions observed in the field. The relevance of this is seen by fixing the soil moisture in the saturated simulation, which alters the timing of freeze-up from 30 days to only 13 days too early. Snow is the most important process for winter soil temperatures, which can be seen here by the high correlation (0.85) between soil temperature error and snow depth error in the winter months. Soil temperatures are particularly sensitive to shallow snow; hence, our improvement to the snow model is essential for simulating soil temperatures in the shoulder seasons. The snow on Samoylov Island is shallow and highly wind-blown, which is typical of these low-lying tundra regions. We find that the fresh snow density required to obtain the correct mid-winter snow density in JULES is too high, indicating a need for further work, potentially to include more snow compaction processes.
Another area for future development is the vegetation, since there are currently no specific high-latitude PFT's in JULES. The moss cover represented here is a first step towards simulating tundra vegetation; however, this represents only the physical effects of a constant layer of moss, leaving more work to be done, for example on growth, carbon cycling, and on other types of vegetation.
We believe that we have significantly improved the representation of permafrost processes in JULES, providing generic model improvements that could be adopted in other GCM land-surface schemes. However, this is still a work in progress for the whole community. Even if a model simulates the right processes in a 1-D column, scaling these up to represent subgrid heterogeneity in a large grid box is still an open problem (Muster et al., 2012;Langer et al., 2013). In most global land-surface models, only vertical processes are simulated, meaning the lateral flow of heat and water, and blowing snow are all omitted. Techniques to include these processes are currently under development (e.g. Tian et al., 2012;Essery and Pomeroy, 2004;Yi et al., 2014). Of course, on the large scale, models are still heavily constrained by the availability and uncertainty of observational data.