A simplified permafrost-carbon model for long-term climate studies with the CLIMBER-2 coupled earth system model

. We present the development and validation of a simpliﬁed permafrost-carbon mechanism for use with the land surface scheme operating in the CLIMBER-2 earth system model. The simpliﬁed model estimates the permafrost fraction of each grid cell according to the balance between modelled cold (below 0 ◦ C) and warm (above 0 ◦ C) days in a year. Areas diagnosed as permafrost are assigned a reduction in soil decomposition rate, thus creating a slow accu-mulating soil carbon pool. In warming climates, permafrost extent reduces and soil decomposition rates increase, resulting in soil carbon release to the atmosphere. Four accumula-tion/decomposition rate settings are retained for experiments within the CLIMBER-2(P) model, which are tuned to agree with estimates of total land carbon stocks today and at the last glacial maximum. The distribution of this permafrost-carbon pool is in broad agreement with measurement data for soil carbon content. The level of complexity of the permafrost-carbon model is comparable to other components in the CLIMBER-2 earth system model.


Introduction
Model projections of climate response to atmospheric CO 2 increases predict that high northern latitudes experience amplified increases in mean annual temperatures compared to mid-latitudes and the tropics (Collins et al., 2013). The large carbon pool locked in permafrost soils of the high northern latitudes ) and its potential release on thaw (Schuur et al., 2008;Harden et al., 2012) make permafrost and permafrost-related carbon an important area of study. Thus far permafrost models that have been coupled within land-surface schemes have relied on thermal heat diffusion calculations from air temperatures into the ground to diagnose permafrost location and depth within soils (Koven et al., 2009;Wania et al., 2009a;Dankers et al., 2011;Ekici et al., 2014). This approach requires a good physical representation of topography, soil types, snow cover, hydrology, soil depths and geology to give a reliable output (Riseborough et al., 2008). The physically based approach lends itself to smaller grid cells and short-timescale snapshot simulations for accuracy of model output. The aim of this work is to develop a simplified permafrost-carbon mechanism that is suitable for use within the CLIMBER-2 earth system model (Petoukhov et al., 2000;Ganopolski et al., 2001), and also suitable for long-timescale experiments. The CLIMBER-2 model with a coupled permafrost-carbon mechanism, combined with proxy marine, continental and ice core data, provides a means to model the past dynamic contribution of permafrost carbon within the carbon cycle.

Physical permafrost modelling
Several land surface models diagnose permafrost and concomitant higher soil carbon concentrations (Wania et al., 2009a, b;Koven et al., 2009;Dankers et al., 2011). These models are usually driven with climatic variables output from Published by Copernicus Publications on behalf of the European Geosciences Union.
global climate models (GCMs), and grid cell sizes are the order of 2.5 • (the order of hundreds of km) for global simulations. These models use surface air temperature and thermal diffusion calculations to estimate the soil temperature at depths, and from this the depth at which water freezes in the soil. An active layer thickness (ALT) can be determined from this, and soil carbon dynamics is calculated for the unfrozen parts of the soil. These land surface models may also include a representation of peatlands (Sphagnumdominated areas, and wetlands), which store an estimated 574 GtC in northern peatlands (Yu et al., 2010), of which a large part are located within the permafrost region (Northern Circumpolar Atlas: Jones et al., 2009). The dynamic response of carbon in permafrost soils subject to (rapid) thaw is not well constrained (Schuur, 2011) and field studies and modelling studies still seek to better constrain this. Riseborough et al. (2008) reviewed advances in permafrost modelling, identifying that modelling of taliks (pockets or layers of thawed soil at depth that do not refreeze in winter) complicates physical modelling. The importance of soil depth (lower boundary conditions) was also highlighted; Alexeev et al. (2007) demonstrated that the longer the simulation, the larger is the soil column depth required in order to produce reliable thermal diffusion-based temperature calculations: a 4 m soil depth can produce reliable temperature predictions for a 2 year simulation, and for a 200 year simulation a 30 m soil depth would be required. Van Huissteden and Dolman (2012) reviewed Arctic soil carbon stocks estimates and the permafrost-carbon feedback. They note the processes by which carbon loss occurs from thawing permafrost including active layer thickening (also caused by vegetation disturbance), thermokarst formation, dissolved organic carbon (DOC) export, fire and other disturbances. Their conclusions are that "current models are insufficiently equipped to quantify the carbon release at rapid thaw of ice-rich permafrost", which within a model would require accurate representation of local topography and hydrology as well as a priori knowledge of the ice content in the soils. Koven et al. (2013) further highlighted the importance of soil depths and of soil and snow dynamics for the accuracy of permafrost extent in CMIP (coupled model intercomparison project) models. The high computing power requirements of physical models at grid sizes where output could be an acceptable confidence level makes these kinds of models currently unsuitable for long-timescale dynamically coupled modelling studies. Current CMIP model projections of future climate reported by the IPCC (Stocker et al., 2013) do not include a possible feedback mechanism from permafrost soils. There exist some studies of the possible future response of carbon in soils of the permafrost zone that do not rely on heat diffusion calculations down the soil column (Scheafer et al., 2011;Harden et al., 2012;Schneider von Diemling et al., 2012). However, these kinds of treatments are not suitable for the study of paleoclimate as they require a priori knowledge of soil organic carbon content (SOCC) of the soils at relatively high reso-lution. This is not yet feasible when considering last glacial maximum soils (for example). Zimov et al. (2009) created a physical model for carbon dynamics in permafrost soils. This one-dimensional model was intended to simulate the carbon dynamics specifically in the permafrost region. Carbon input to the soil originates from root mortality and aboveground litter transport via organic carbon leaching and mixing by bioturbation and cryoturbation. Loss of carbon from the soils occurs via decomposition. The frozen soil active layer depth also determines the maximum root depth of vegetation. Modelled soil carbon profiles were similar to those found in present-day ground data for similar conditions. Results of experiments where the temperature zone was changed linearly from Temperate to Cold, then snapped back to Temperate (mimicking a glaciation then termination in Europe), demonstrated the characteristic of slow carbon accumulation in permafrost soils, and fast carbon release on thaw. An important result of this study was that the main driver of the high carbon content in the frozen soils was the low decomposition rates, which reduce further with depth in the soil column, as a result of permafrost underlying an active layer which cycles between freezing and thawing in the year. To estimate the amounts of carbon stored on the land and the ocean at Last Glacial Maximum (LGM), Ciais et al. (2012) used δ 18 O data and carbon cycle modelling to calculate gross primary productivity (GPP) at LGM and in the present day. They estimate that the total land carbon stocks had increased by 330 GtC since LGM, but that 700 GtC less was at present stored as inert land carbon stocks compared to LGM. Zech et al. (2011), studying two permafrost-loess paleosol sequences, concluded that on glacial timescales the effect of reduced biomass productivity may be of secondary importance to the effect of permafrost preserving soil organic matter when considering total land carbon stocks. The Ciais et al. (2012) inert land carbon stock may represent this permafrost-carbon pool.

Carbon cycle responses during a deglaciation
The current leading hypothesis for the fast rise in atmospheric CO 2 in the last glacial termination (17.5 to 12 kyr BP) (Monnin et al., 2001) is that carbon was outgassed from the ocean via a reorganisation of ocean circulation that released a deep carbon store in the Southern ocean (Sigman et al., 2010;Fischer et al., 2010;Shakun et al., 2012). The Zimov et al. (2009) model, Ciais et al. (2012 and the δ 13 CO 2 record for the last termination (Lourantou et al., 2010;Schmitt et al., 2012) suggest that permafrost may have had a role to play in the dynamics of the carbon cycle during the last termination. At the start of glacial termination 1 (from the end of the last glacial period, the transition to the interglacial climate, starting at ∼ 17.5 kyr BP) a fast drop in the δ 13 CO 2 of the atmosphere was seen from ice core data. Soil carbon has a δ 13 C signature depleted by around 18 ‰ compared to the atmosphere (Maslin and Thomas, 2003); a release of carbon from thawing permafrost soils is a possible explanation for the δ 13 CO 2 record. In this study, we aim to develop a permafrost-carbon model for long-term paleoclimate studies. We present the development of the permafrost-carbon model and validate it with present-day ground measurement data for soil carbon concentrations in high northern latitude soils.

CLIMBER-2 standard model
The CLIMBER-2 model (Petoukhov et al., 2000;Ganopolski et al., 2001) consists of a statistical-dynamical atmosphere, a three-basin averaged dynamical ocean model with 21 vertical uneven layers and a dynamic global vegetation model, VECODE (Brovkin et al., 1997). The model version we use is as Bouttes et al. (2012) and Brovkin et al. (2007). The model can simulate around 20 kyr in 10 h (on a 2.5 GHz processor) and so is particularly suited to paleoclimate and long-timescale fully coupled modelling studies. The version of CLIMBER-2 we use (Bouttes et al., 2009(Bouttes et al., , 2012 is equipped with a carbon-13 tracer, ice sheets and deep sea sediments (allowing the representation of carbonate compensation) in the ocean (Brovkin et al., 2007) as well as ocean biogeochemistry. The ice sheets are determined by scaling ice sheets' size between the LGM condition from Peltier (2004) and the Pre-Industrial (PI) ice sheet using the sea level record to determine land ice volume (Bouttes et al., 2012). The dynamic vegetation model has two plant functional types (PFTs), trees and grass, plus bare ground as a dummy type. It has two soil pools, "fast" and "slow", representing litter and humus respectively. Soils have no depth, and are only represented as carbon pools. The carbon pools of the terrestrial vegetation model are recalculated once every year. The grid cell size of the atmospheric and land surface models is approximately 51 • longitude (360/7 degrees) by 10 • latitude. Given the long-timescale applications of the CLIMBER-2 model and the very large grid size for both atmosphere and land, none of the existing approaches of modelling permafrost carbon are suitable. Thermal diffusionbased physical models would produce results with unacceptable uncertainties (error bounds) compounding over long timescales. To create the permafrost model for CLIMBER-2, the driving mechanism creating high soil carbon concentrations is a reduced soil decomposition rate in the presence of permafrost, identified by Zimov et al. (2009) as the primary driver in soil carbon accumulation for these soils.

Permafrost-carbon mechanism
CLIMBER-2 grid cells for the land surface model are very large. Two options are available to diagnose permafrost location: either by creating a sub-grid within the land grid or by diagnosing a fraction of each grid cell as permafrost, which is the approach followed here. Conceptually the sub-grid model represents keeping permafrost carbon separate from other soil carbon, and the remixing model represents mixing all soil carbon in a grid cell. Figure 1 shows a schematic representation of a CLIMBER-2 grid cell, and how the permafrost fraction of the land is defined relative to other cell parameters when permafrost is diagnosed as a fraction of each cell. For the carbon cycle the calculations of carbon fluxes between atmosphere and land grid cells are for the cell mean. Each grid cell contains cell-wide soil carbon pools (fast soil or slow soil, per plant functional type), so to account for permafrost soils either a new permafrost-soil pool needs to be created for each grid cell, or permafrost soils can be mixed back into the standard soil pools at every time-step (Fig. 2a). If the land grid is downscaled a third option is available, where each sub-grid cell maintains an individual soil carbon pool (Fig. 2b). This, however, requires an increase in computational time which slows down the run speed of the model.
The soil carbon in CLIMBER-2 is built from vegetation mortality and soil carbon decomposition is dependent on surface air temperature, the total amount of carbon in the pool and the source of carbon (i.e. trees or grass). Equation (1) shows how carbon content of each pool is calculated in CLIMBER-2. The pool is denoted by C i , where pool C 1 is plant green phytomass (leaves), C 2 is plant structural biomass (stems and roots), C 3 is a soil pool made of litter and roots residue and C 4 is a soil pool made of humus and residues of woody-type stems and roots. Hereafter, the soil pools will be referred to as Soil fast for C 3 and Soil slow for C 4 .
where C is the carbon content in the pool (kgC m −2 ); k is allocation factors (0 < ki < 1); N is net primary productivity (NPP; kgC m −2 yr −1 ); m i is decomposition rates for the carbon in each pool (yr −1 ); p is the plant functional type (trees or grass).
The residence time of carbon in soil pools is 1/m; we call this τ . For residence times corresponding to decomposition rates m 3 and m 4 , τ is: where i is the soil pool, n is a multiplier dependent on the pool type, ps5 is a constant, = 0.04, T mat is mean annual temperature at the surface-air interface, • C, T ref is a reference soil temperature, fixed in CLIMBER-2 at 5 • C. The value of n is dependent on the soil carbon type, being 900 for all slow soils, 16 for fast tree PFT soil and 40 for fast grass PFT soil. The decomposition rates for organic residue in the soils are most strongly based on soil microbial activity and the relative amount of lignin in the residues (Aleksandrova, 1970;Brovkin et al., 1997). Increasing the residence time of carbon in permafrost-affected soils reduces the decomposition rates and results in higher soil carbon concentrations. We modify the residence time, τ 3,4 , in the presence of permafrost using: where a and b are tuneable dimensionless constants and F sc is frost index, a value between 0 and 1, which is a measure of the balance between cold and warm days in a year, and is shown in Eq. (4) where DDF is degree-days below 0 • C and DDT is degree-days above 0 • C in a year for daily average surface air temperature (Nelson and Outcalt, 1987). DDF and DDT have units of • C days yr −1 . Snow cover acts to insulate the ground against the coldest winter temperatures and reduces permafrost extent (Zhang, 2005;Gouttevin et al., 2012). The subscript sc in Eqs.
(3) and (4) indicates that these values are corrected for snow cover and represent the ground-snow interface conditions, not the snow surface-air interface conditions. Including the frost index as a multiplier (in Eq. 3) for the permafrost soils' carbon residence time was needed to allow the correct tuning of the model and allow for total land carbon stocks to be in agreement with data estimates. Therefore, the decomposition rates of soil carbon in permafrost-affected cells are dependent on: mean annual temperature (as with non-permafrost soils), the fractional cover of permafrost in the cell and the frost index (a measure of the severity of cold- ness in a year). This τ perma i (Eq. 3) is only applied to the soils that are diagnosed as permafrost. The remainder of the carbon dynamics in land carbon pools was unaltered from the standard model.

1-D model
We test a one-dimensional model to compare the effect the different assumptions made for the model design. The total carbon stock in a grid cell using each method (sub-grid and remixing) was compared for equilibrium soil carbon content by running the 1-D model for 100 000 simulation years. The carbon input from vegetation mortality is the same for the remixing and the sub-grid model, as is rainfall. The variables of permafrost fraction, mean annual air-surface interface temperature (MAT) and frost index are varied one at a time to compare the model outputs. The constants a and b for Eq. (3) were set to 20 for Soil fast and 2 for Soils slow (so a and b have matching values) for the permafrost soils, and as the standard model for the non-permafrost soils. These values for a and b were chosen to compare the performance of the two methods, not for accurate soil carbon concentrations. They result in total carbon in the Soil fast and the Soil slow carbon pools being approximately equal, which studies suggest is appropriate (Harden et al., 2012;Zimov et al., 2009). Figure 3 shows the output for carbon content along a permafrost gradient, taking account of the relationship between permafrost fraction, frost index and mean annual temperature. More detail on this figure is available in Appendix A. The relationship between permafrost fraction and frost index is defined as that determined in this study for the CLIMBER-2 model in Sect. 3.2. As shown in Eq. (1), NPP exerts a control on soil carbon content via input from plant material, although note that Fig. 3 shows model output for fixed NPP. For both approaches, carbon content increases non-linearly Figure 4. Comparison of NPP, which has a control on carbon input to soils, for MODIS data set (top, mean 2000-2005) and CLIMBER-2 model for PI(eq) (modelled year 1950) plotted on the same scale (gC m −2 yr −1 ). MODIS data upscaled to CLIMBER-2 grid scale are shown against equivalent points for CLIMBER-2 NPP.
along the permafrost gradient (increasing permafrost fraction of the grid cell). The remixing model shows stronger non-linear behaviour than the sub-grid model.

CLIMBER-2 modelled NPP
The comparisons of the sub-grid to remixing approaches shown in Fig via NPP in colder climates. Figure 4 shows the CLIMBER-2 modelled NPP and the MODIS 2000-2005 mean NPP product (Zhao et al., 2011) for the present day (PI, pre-industrial for CLIMBER output). The CLIMBER-2 vegetation model shows NPP patterns similar to the MODIS data set. The bo-  Figure 6. Modelled output for 1-D models along a permafrost gradient, with correction for NPP and initial value (at 0 % permafrost). Overlaid on 1 • data for SOCC binned into 0.1 permafrost fraction mean values ±1 sigma (Hugelius et al., 2013), permafrost fraction is calculated using relationship identified in Sect. 3.2. real forest belt seen at around 60 • N in the MODIS data set is not clearly seen in the CLIMBER-2 model, mainly due to the large grid cell size. In Siberia and Alaska the NPP in CLIMBER-2 is not overestimated. The reduced NPP in the coldest regions would tend to reduce soil carbon accumulation via reduced input from plant mortality. Also shown in Fig. 4 are the upscaled data points plotted against CLIMBER-2 model output. The MODIS data set represents the earth system already subject to anthropogenic forcing, where the CLIMBER-2 model output represents the natural system only. However, the use of measurement-based data to validate CLIMBER-2 NPP was preferred due to the quite large model spread seen in output for numerical global dynamic vegetation models of higher complexity than CLIMBER-2. The fact that MODIS is for the present-day "perturbed" system (due to deforestation, for example) may also explain some of the model-data mismatch, although we consider this less significant for the permafrost zone low-NPP soils in which we are interested . In order to test the applicability of the CLIMBER-2 model for the glacial climate, a comparison of NPP for the LGM with a more complex model can be done (as measurement data are not available). Figure 5 shows LGM(eq) NPP for LPX (data courtesy of M. Martin-Calvo, Prentice et al., 2011) and for CLIMBER-2 for an LGM climate. At LGM the NPP in Siberia and the coldest permafrost regions is non-zero in both models, and CLIMBER-2 follows the same general patterns as LPX predicts. CLIMBER-2 shows slightly lower NPP in the southern parts of Russia, possibly similar to the boreal forest belt that is not well represented in the PI climate background NPP due to the large grid cell size. Again, the upscaled LPX data are shown plotted against CLIMBER-2 output, showing reasonable agreement on this scale. Overall at both periods, PI and LGM, CLIMBER-2 represents NPP reasonably well. . Measurement data for active layer thickness (CALM network, Brown et al., 2003) and frost index (Zhang, 1998)  When the soil carbon content shown in Fig. 3 is adjusted to compensate for the reduction in NPP along a permafrost gradient and for the 0 % permafrost SOCC data value (by multiplying relative value by 350), the resultant outputs are shown in Fig. 6 (more details are available in Appendix A). Now the remixing model shows a slight increase in total carbon along a permafrost gradient, where the sub-grid model shows a peak value at around 80 % permafrost coverage. Figure 6 shows a comparison between these 1-D model outputs and data for SOCC. The unadjusted data are for the top 1 m of soils, whereas model output represents the full soil column. As in Sect. 4.4, the model-data comparison is carried out by assuming that 40 % of total soil carbon is located in the top 1 m for permafrost soils (and is fully described in Appendix A). From this comparison, the change in SOCC along a permafrost gradient is relatively small, due to the combined effects of reducing soil decomposition rate and reducing NPP. Here, the remixing model represents these changes quite well. It may be possible to improve the performance of the sub-grid model by, for example, downscaling the climate variables also. However, this would represent a more significant change of the land biosphere model in CLIMBER-2, and increase the complexity and therefore reduce the computational efficiency of the model. For the remixing model: at each time-step a proportion of carbon that is accumulated in the permafrost part is then sent back to decompose as standard soil. This occurs because the high-carbon permafrost soil is mixed with the lower carbon standard soil in a grid cell at each time-step. This can be seen as similar to that which occurs in the active layer. The active layer is the top layer of the soil that thaws in warm months and freezes in cold months. In warm months the carbon in this thawed layer is available to be decomposed at "standard" soils rates, determined by local temperature. In the remixing model, the relative proportion of the permafrost soil carbon that is sent to decompose as standard soil carbon reduces along a permafrost gradient. This reduction can be seen as mimicking the characteristic of a reducing active layer thickness along a permafrost gradient, which is shown in Fig. 7 for active layer thickness data upscaled to the CLIMBER-2 grid size. Here active layer thickness mean is shown plotted against mean frost index (and permafrost fraction is directly calculated from frost index in CLIMBER-2). It must be noted that on smaller spatial scales the relationship between the mean active layer thickness and the extent of permafrost in a location may be less clear. The local conditions determine both permafrost extent and active layer thickness. Our treatment for permafrost relies entirely on the relationships between climate characteristics and soil carbon contents on the CLIMBER-2 grid scale.

CLIMBER-2 permafrost-carbon model
We implemented Eq. (3) into CLIMBER-2 using the remixing model. In order to study the effect of different carbon accumulation and release rates (the permafrost-carbon dynamics) in later modelling studies, the soil carbon residence times can be tuned to distribute the carbon more into the Soil fast pool (making a quickly responding soil carbon pool) or more into the Soil slow pool (making a more slowly responding soil carbon pool). A total of four dynamic settings are retained for later coupled climate studies (described in Sect. 3.5).

Simulated climates to tune the permafrost-carbon model
Three simulated climates were used to tune and validate the permafrost-carbon model: an LGM equilibrium climate, LGM(eq); a PI equilibrium climate, PI(eq); and a PI transient climate, PI(tr) obtained at the end of a transient deglaciation from the LGM climate. These three climates allow the total soil carbon to be tuned to the estimates of Ciais et al. (2012) for the LGM and PI climate conditions; these are described in Table 1.

Calculating permafrost extent
In order to obtain a relationship between calculated frost index and the permafrost fraction of a grid cell, measurement and ground data for frost index and permafrost location were used. For present-day mean daily surface air temperatures, the freeze and thaw indices values on a 0.5 • global grid were obtained from the National Snow and Ice Data Centre (NSIDC) database (Zhang, 1998). Using these values for freeze and thaw index, a global frost-index data set on a 0.5 • grid scale was created using Eq. (4). The presentday estimates of land area underlain by permafrost are provided by Zhang et al. (2000), using the definition of zones: "continuous" as 90-100 % underlain by permafrost, "discon- tinuous" as 50-90 %, "sporadic" as 10-50 % and "isolated" as less than 10 %. Zhang et al. (2000) used these zonations to provide area estimates of the total land area underlain by permafrost. Summing the total land area that has a frost index higher than a particular value and comparing this to the Zhang et al. (2000) estimate can identify the appropriate boundary between permafrost and non-permafrost soils. Figure 8 shows the Zhang et al. (2000) permafrost areas for the high, medium and low ranges defined by the high, medium and low % estimates of permafrost zones marked as horizontal lines. The land area indicated by green squares is the total land surface in the northern hemisphere that has a frostindex value higher (where higher indicates a colder climate) than the cut-off value shown on the x axis. Here the frostindex cut-off value of 0.57 shows good agreement with the medium (mean) estimate of the Zhang et al. (2000) total area of land underlain by permafrost.

Geographic permafrost distribution for the
present-day Figure 9 shows, coloured in blue, the land grid cells with a frost-index higher than 0.57 for 0.5 • grid, with the north located at the centre of the map. Overlaid on this map area are the limits of the permafrost zones defined by the International Permafrost Association (IPA) (Jones et al., 2009). The frost-index value cut-off at 0.57 results in a southern limit of permafrost that represents approximately the middle of the discontinuous zone with some areas showing better agreement than others. Figure 10 represents the upscaling of the 0.5 • data sets for mean frost index and permafrost coverage to the CLIMBER-2 land grid scale. It shows the percentage of land in each CLIMBER-2 size grid cell defined as permafrost (according to the 0.57 frost-index cut-off value shown in Fig. 8) plotted against the mean value of frost index for the same grid cell. LGM (equilibrium) Obtained after an 80 kyr spin-up with glacial CO 2 levels of 190 ppmv, reduced ocean volume, LGM ice sheets, LGM insolation, LGM runoff. Carbonate compensation in the ocean (Brovkin et al., 2002). Sea-level effects on coast lines are not included; land area is as PI (equilibrium). The continental shelves exposed at LGM are not accounted for in this model set-up because the fate of any carbon that may have accumulated on these shelves is not well constrained. The long time of spin-up, 80 kyr, is required to allow the soil carbon pools to equilibrate.

PI (equilibrium)
Obtained after 40 kyr spin-up with pre-industrial CO 2 levels of 280 ppmv, present-day ocean volume, present-day ice sheets, insolation, and land run-off. Circled points in Fig. 9 are where the grid cell has a large fraction of ocean (more than 75 %), and the milder ocean temperatures in winter reduce the mean frost-index value of the whole grid cell. The dashed line shows a well-defined sigmoid function that relates frost index to permafrost percentage of the land. We employ this relationship to predict permafrost area in CLIMBER-2, as the frost index can be calculated within the model from modelled daily temperatures. Permafrost fraction is thus modelled as: where A and β are defined in Table 2 and the model described in Sect. 3.5. Frost index is calculated from modelled daily surface temperatures and corrected for snow cover. The snow correction in our model is achieved using a simple linear correction of surface-air temperature, using snow thickness to estimate the snow-ground interface temperature. This correction is based on data from Taras et al. (2002). The snow correction performs reasonably well in CLIMBER-2 compared to measurement data from Morse and Burn (2010) and Zhang (2005). This is because the large grid-cell size results in non-extreme snow depths and air surface temperatures. The snow correction is described in Appendix B. Equation (6) shows this linear model for snow correction, which is only applied for daily mean surface air temperatures lower than −6 • C. This snow-ground interface temperature is used to calculate the freeze index (DDF sc ) in Eq. (4).
where T g.i is ground interface temperature ( • C), T surf is surface air temperature ( • C) and SD is snow depth (cm). Overall the effect of the snow correction within the model produced a maximum decrease in permafrost area of 8 % (compared to the uncorrected version) in the most affected grid cell for the PI(eq) simulation and is therefore significant.

Permafrost extent tuning
Using the snow-corrected frost-index value, four permafrost extent models representing the range of values for permafrost area from Zhang et al. (2000) were determined. The model settings are shown in Table 2 and refer to A and β from Eq. (5). P landfraction is limited between 0 and 1, and the functions are plotted in Fig. 11. These settings were identified by adjusting the sigmoid function to obtain total permafrost area values at the PI(eq) simulation similar to the Zhang et al. (2000) areal estimates of permafrost and to maximise the difference in area between the PI(eq) and LGM(eq) simulations permafrost extent. More complex models underestimate permafrost extent at LGM (Levavasseur et al., 2011;Saito et al., 2013) quite significantly, and so by maximising the difference between PI and LGM permafrost, we reduce the underestimate as far as possible for LGM permafrost extent.  Figure 11. CLIMBER-2P model for permafrost fraction of the land in a grid cell from frost index (snow corrected). Range of areas is within the range of estimates for present-day land area underlain by permafrost by Zhang et al. (2000). Permafrost fraction is limited between 0 and 1. Zhang estimate for total permafrost area is 12.21 to 16.98 × 10 6 km 2 . Listed from HIGH to LOW, model output is 16.35, 14.87, 14.00 and 13.21 × 10 6 km 2 .

Tuning the soil carbon model
Soil carbon content is controlled by the balance between soil carbon uptake and soil carbon decomposition. There are four soil carbon pools in CLIMBER-2: Soil fast : trees derived and grass derived; Soil slow : trees derived and grass derived (Eq. 1). Soil fast have shorter carbon residence times than Soil slow , so soil decays more quickly in Soil fast pools. The tunable constants a and b (Eq. 3) are independently applied for Soil fast and Soil slow , so carbon can be placed relatively more in the Soil fast (Soil slow ) pool as required in model tuning. Carbon is lost from permafrost soils as the permafrost fraction of a grid cell reduces. If there is relatively more (less) carbon in the Soil fast pool, this results in carbon that decays more quickly (more slowly) when the permafrost thaws.
At LGM, the area of permafrost on land was larger than today (Vandenberghe et al., 2012), but not much information on soil carbon has been conserved, especially if it has long since decayed as a result of permafrost degradation during the last termination. To constrain the total carbon content in permafrost soils we use the estimates of ; for total land carbon these are 3640 ± 400 GtC at LGM and 3970 ± 325 GtC at PI, with a total change of +330 GtC between LGM and PI. The standard CLIMBER-2 model predicts total land carbon stocks of 1480 GtC at LGM and 2480 GtC at PI, showing good agreement with the active-land-carbon estimates of Ciais et al. (2012) (of 1340 ± 500 GtC LGM and 2370 ± 125 GtC PI). Any "new" soil carbon is created via the permafrost-carbon mechanism and is assumed to be equivalent to the inert land carbon pool estimates of Ciais et al. (2012). However, the dynamic behaviour of permafrost-carbon in changing climates is not well constrained and it is for this reason that a set of four dynamic settings were sought. Here the "speed" of the dy- namic setting is determined by the ratio of total Soil fast pool to Soil slow pool carbon (fp / sp), with the "slow" dynamic being fp / sp < 0.5, "medium" being fp / sp 0.5 to 1, "fast" being fp / sp 1 to 1.5 and "extra-fast" being fp / sp > 1.5 for the PI-equilibrium simulation. The variables "a" and "b" shown in Eq.
(3) were set and each setting used to run a PI(equilibrium), LGM(equilibrium) and PI(transient) simulation to identify the settings that resulted in total land carbon pools in agreement with the Ciais et al. (2012) estimates.
The LGM is conventionally defined as being the period around 21 kyr BP, when large parts of north America were underneath the Laurentide ice sheet. According to their timeto-equilibrium (the slow carbon accumulation rate), soils in this location now free of ice may not yet have reached equilibrium. Furthermore, climate has changed significantly since the LGM so permafrost soils anywhere may not be currently in equilibrium (Rodionov et al., 2007), again due to their slow carbon accumulation rates. Due to this the PI(tr) simulation model output for total land carbon was used to tune the total land carbon stocks, as it includes a receding Laurentide ice sheet. At LGM, ice sheets were at maximum extent, so the problem of land being newly exposed does not occur in the model. For this reason, the LGM(eq) simulation is used to tune total land carbon for the LGM.
Details of the tuning for total land carbon stocks are available in Appendix C. It was found that only one permafrost area setting, the LOW-MEDIUM area, provided an acceptable range of dynamic settings, as defined by the ratio of fast to slow soil carbon. The four selected dynamics settings are shown in more detail in Fig. 12 for total land carbon stock, atmospheric CO 2 and ratio of fast to slow soil carbon pool. The a and b values for these settings are shown in Table 3.
To evaluate the effect of the different dynamic settings we ran an equilibrium PI simulation for all four selected settings for 40 kyr, followed by a permafrost switch-off for a further 10 kyr. Figure 13 shows the global total land carbon stocks for this experiment. The period between 0-40 k simulation years demonstrates the transient effects of the slow accumulation rates in permafrost soils. Depending on the dynamic setting, the total land carbon takes more than 40 kyr to fully equilibrate in PI climate conditions. On permafrost switch-off, from 40 k simulation years, the soil carbon pre-Geosci. Model Dev., 7, 3111-3134, 2014 www.geosci-model-dev.net/7/3111/2014/ viously held in permafrost soils is quickly released to the atmosphere, at a rate dependent on the dynamic setting. The xfast setting entails releasing all excess carbon within hundreds of years and the slow setting around 8000 years after total permafrost disappearance. Currently, the most appropriate carbon dynamic setting is unconstrained by measurement data. It is for this reason that the permafrost-carbon dynamics settings cover a large range. They are intended to be used in transient model simulations to better constrain permafrostcarbon dynamics in changing climate. It should be noted that the PI(eq) simulation was not used to tune the model, i.e. was not used to compare model output to Ciais et al. (2012) PI total land carbon stocks. Figure 13 demonstrates only the range of dynamic response for all four settings. This PI(eq) simulation also demonstrates the difference between transient vs. equilibrium PI simulations. The slow dynamic equilibrates (after more than 40 kyr) at far higher total carbon stocks than the xfast dynamic, but for the PI(tr) simulation these two settings show very similar total land carbon stocks (we selected them for this behaviour).

Model performance
Hereafter, "CLIMBER-2P" denotes the model in which the permafrost-carbon mechanism operates fully coupled within the dynamic vegetation model. Figure 14a shows the spatial pattern of permafrost as predicted in CLIMBER-2P with the snow correction included for the LOW-MEDIUM area setting. The modelled PI(tr) permafrost extent estimates fairly well the location of the present-day southern boundary of the discontinuous permafrost zone (Jones et al., 2009), with overestimate of permafrost extent in the western Siberian grid cell, and underestimate over the Himalayan plateau. Total permafrost area extent is shown in Table 4.

Permafrost areal coverage and spatial distribution
Comparing this to performance of other models (Levavasseur et al., 2011), the PI(eq) total permafrost area is closer to Zhang et al. (2000) estimates, but it must be kept in mind that for CLIMBER-2P the area was tuned to be in agreement with a mean estimate from Zhang et al. (2000). The PI(tr) total permafrost area is higher by around 4 × 10 6 km 2 compared to the PI(eq). This is due to the North Pacific region being colder in PI(tr) than that of the PI(eq) simulation, and may be related to the land run-off, which is kept at LGM settings for the transient simulations. For the LGM period, the best PMIP2 model in the Levavasseur study (interpolated case) underestimated total permafrost area by 22 % with respect to data estimates (of 33.8 × 10 6 km 2 ), and the "worst" model by 53 %, with an all-model-median value of 47 % underestimate. The LOW-MEDIUM CLIMBER-2P setting gives an LGM total permafrost area underestimate of around 40 %, slightly better than the median for PMIP2 models' permafrost area. Figure 14b shows the LGM CLIMBER-2P permafrost extent with the reconstructed continuous and discontinuous southern boundaries (Vandenberghe et al., 2012;French and Millar, 2013) overlaid. In the LGM simulation for CLIMBER-2P, coastlines do not change so the Siberian Shelf  (Zhang et al., 2000) 33 LGM(eq) simulation for LOW-MEDIUM permafrost area. Overlaid in orange are data estimates from Circumpolar Atlas (Jones et al., 2009) for the present day, Vandenberghe et al. (2008) for LGM Eurasia, French and Millar (2013) for LGM N. America. and other exposed coastlines in the northern polar region are not included in the CLIMBER-2P permafrost area estimate. These coastal shelves cover an estimated area of 5 to 7×10 6 km 2 . Another area that is not diagnosed as permafrost in CLIMBER-2P is the Tibetan plateau, which would be an additional estimated 6 × 10 6 km 2 . If these two regions were added (totalling around 12 × 10 6 km 2 ) to the LGM area estimate it would bring the modelled permafrost area (then totalling around 33 × 10 6 km 2 ) much closer to the data estimate as reported in the Levavasseur et al. (2011)  Tibetan plateau are problematic, resulting in a possible toowarm climate (compared to the real world) in this region.

Soil carbon dynamics
Accumulation rates show general agreement with the Zimov et al. (2009) model and the Wania et al. (2009b) (LPJ) model, although the fast and xfast dynamic settings accumulate carbon faster than these comparison models. Figure 15 shows output for all permafrost dynamics for the PI (equilibrium) spin-up. The north-west Siberia site can be compared to the the Ayach-Yakha location from Wania et al. (2009b) and to the extra-cold conditions from Zimov et al. (2009). The Ayach-Yakha modelled site in Wania et al. (2009b) has a time to equilibrium of greater than 80 kyr and soil carbon content of greater than 200 kg m −2 ; the Zimov model predicts that 200 kg m −2 soil carbon content can be reached within 10 kyr in the top layer of the soil and 150 kg m −2 for the full soil column taking longer than ∼ 50 kyr to reach equilibrium. The N. Canada (Fig. 15) location takes a longer time to reach equilibrium than soils in the NW Siberia grid cell. NPP in the N Canada grid cell is less than one third of that for the NW Siberia grid cell. Due to the lower soil carbon input there is a lower range in the output between the different carbon dynamic settings for the N Canada grid cell. Northern Canada was underneath the Laurentide ice sheet at LGM. Since the demise of the Laurentide ice sheet around 13 kyr ago (Denton et al., 2010) there has not been enough time for these soils to equilibrate, which takes longer than 40 kyr according to our model. As well as this, this region has very high water contents (and islands) that are not represented in CLIMBER-2P, which may modify soil carbon concentrations. Although we do not account for water content, we can take account of the demise of the Laurentide ice sheet and the time that these soils have had to accumulate carbon. The PI climate condition and soil carbon content that we applied to tune and validate the model are the PI(tr), the transient simulation, which includes ice sheet evolution.

Soil carbon stocks
The total land carbon stocks were tuned using data from Ciais et al. (2012). An assumption made in this study is that all "extra" soil carbon, relative to the standard model, in the Arctic region is located in permafrost soils and only by the mechanism of increased soil carbon residence time in frozen soils.  Table 7). The soil types that are found in the continuous and discontinuous permafrost zone are the Cryosols (circumpolar atlas) or Gelisols (soil taxonomy). Within this group are subgroups: Turbels, which are subject to cryoturbation and characterise the continuous permafrost zone; Orthels, which are less affected by cryoturbation and are related to discontinuous permafrost; and Histels, which relate to peat growth (histosols) and have permafrost at less than 2 m depth. Histels are not directly represented in the simplified model, as they are dominated by peat growth (Sphagnum), a distinct PFT not represented in CLIMBER-2P.
The Tarnocai et al. (2009) SOCC estimates for the present day for relevant soils are shown in Table 6. Summing "All soils" with loess soils and Deltaic deposits gives the 1672 GtC estimated total SOCC for the permafrost region. The extra land carbon stocks created in our model in permafrost soils range between 1620 to 2226 GtC (Table 8) compared to  at 1672 GtC and to  at 1600 ± 300 GtC for inert land carbon for the present day. For the LGM climate, the model shows a range of 1987 to 2117 GtC for extra soil carbon compared to the Ciais estimate of 2300 ± 300 GtC for inert land carbon. The "medium" dynamic setting shows total land carbon stocks in the present day outside the range estimated by , However, during tuning (see Appendix C) this overestimate could not be improved upon.

Soil carbon contents validation
The carbon content of Orthels and Turbels decreases with depth, but high carbon contents are still found at depths of 3 m and more . For Orthels (with alluvium), around 80 % of their carbon content was found in the top 200 cm, and for Turbels 38 % of carbon content was found in the top 100 cm. Based on these values, to compare the CLIMBER-2P output with ground spatial data, it is assumed that 40 % of the modelled total soil-column carbon is located in the top 100 cm for all permafrost-affected soils.
Soil carbon data from Hugelius et al. (2013) were used to compare against the CLIMBER-2P output. The Orthels and Turbels dominate the continuous and coldest permafrost areas, with Histels and other soils becoming more dominant towards the southern parts of the permafrost region. As no peatlands or wetlands are represented in our simplified model, only Orthel and Turbel soils were used as comparison points for SOCC. SOCC data from Hugelius et al. (2013), for grid cells with 50 % or more Orthel and Turbel soils, were upscaled to the CLIMBER-2P grid. These mean SOCC data values for the top 1 m of soil were plotted against CLIMBER-2P model output for matching grid cells; this is shown in Fig. 16. Also shown in Fig. 16 is the standard model output, which has no permafrost mechanism. Two grid cells show very much higher SOCC than data suggest, with around a 3fold overestimate, and are located in Siberia. All other grid cells are within a range of ±80 %, heavily dependent on the soil carbon dynamic setting. The standard model shows progressively worse performance as mean SOCC increases in the data. The permafrost model shows an increasing SOCC trend more similar to data. The spatial location of SOCC can be compared to data using Fig. 17. The two grid cells with very high SOCC compared to data are central and eastern Siberia. These grid cells are both 100 % permafrost and have had a total of 101 kyr (80 k for LGM(eq) plus 21 k to PI(tr)) to accumulate carbon. This is in contrast to the North American continent grid cells, which were underneath the ice sheet until the deglaciation so have had less time to accumulate carbon.
The assumption that all permafrost region soil carbon acts as Turbels and Orthels has an impact on the physical location of the SOCC with respect to data. Turbels and Orthels are located in the northern parts of the permafrost zone, with Histels and other soils becoming more dominant to the south. Compared to SOCC in ground data (Fig. 17), a northern bias in SOCC is seen in model output, as expected. Histels (peatland soils) and other soil types of the permafrost zone, with  , 2009), are also not represented in our model. One example of this is the Ob River and Gulf of Ob, located in western Siberia, which, combined with dominance of Histels in this region (Hugelius et al., 2013), cause a high SOCC in data. The model does not represent well the boreal forest belt (see Fig. 4), which is also located in the southern region of the permafrost zone. This results in carbon input to soils in this region being underestimated in our model. Figure 18 shows the model outputs for the LGM climate. No soil carbon is present underneath ice sheets and the highest carbon concentrations are seen in present-day southeastern Russia and Mongolia, with quite high soil carbon concentrations in present-day northern Europe and northwestern Russia. Comparing this output to the permafrost extent model (Fig. 14), the SOCC is likely located too far north for the same reasons as the PI(tr) SOCC but also because permafrost extent is underestimated for the LGM(eq) climate. The northern China region, according to data, was continuous permafrost at LGM, as was the southwest Russia region. These regions would have higher SOCC in model output if the modelled permafrost area were closer to data estimates. The same would be true of the Siberian shelf. This means that the extra soil carbon tuned to the Ciais et al. (2012) estimate  (Table 5) is more concentrated in a central band in Eurasia than the model would predict if permafrost extent were more like the data estimate for LGM.

Applications
The simplified permafrost mechanism is intended to be used for the study of carbon-cycle dynamics on timescales of centuries/millennia and longer. It represents an improvement on the previous terrestrial carbon cycle model in CLIMBER-2, which did not include any effects of frozen soils. It is not intended for the study of carbon cycle dynamics on scales shorter than centuries due to the simplifications made and many processes not accounted for. The permafrost-carbon mechanism is dependent on the relationship between climate, soil carbon content and active layer thickness on the CLIMBER-2 grid scale. To apply this parameterisation of permafrost-carbon to other grid scales, the relationship of active layer thickness and climate variables would need to be reassessed. The relationship between permafrost fraction of a grid cell and SOCC is non-linear. The values for "a" and "b" would need to be retuned in order to output total land carbon stocks in agreement with Ciais et al. (2012) for grid scales different to the CLIMBER-2 grid.
The permafrost-carbon mechanism is fully dynamic and responds to changes in: insolation (orbit), atmospheric CO 2 Geosci. Model Dev., 7, 3111-3134, 2014 www.geosci-model-dev.net/7/3111/2014/ (via changes in NPP and climate), and land area in response to coverage by ice sheets extending or contracting. This could not be easily achieved if a box model representation of permafrost-carbon were applied as the model response to the drivers (orbit, CO 2 and ice sheet) are dependent on spatial location.

Simplifications and limitations
The permafrost model does not make any changes to soil carbon based on hydrology or ice contents. Precipitation only affects vegetation growth, not soil formation.
No account is taken of the effect of peatland soils in permafrost regions as the PFT for Sphagnum species, which accounts for most of peat soil vegetation cover, is not included in the model. The effect of frozen ground inhibiting root growth (to depth) is not accounted for, which may have an impact on the GPP and soil formation in very cold regions.
During glacial climates, no extra land is exposed as sea level drops in the CLIMBER-2P model; all the carbon used to tune the carbon dynamics for the LGM period is located on land that is at present above sea-level.
Wetlands and river deltas increase the spatial spread of the soil carbon in the real world, and these are not represented in CLIMBER-2P. Therefore, it is also not intended that the spatial location of the highest soil carbon concentrations should be used as a very good indicator of the real world case.
Slow accumulation rates in permafrost soils result in the characteristic that in the real world during thaw (or deepening of the active layer) the youngest soils would decompose first. In CLIMBER-2P all soil is mixed, so the age of carbon down the soil column cannot be represented. This age of the soils is important for the correct modelling of 14 C then seen in the atmosphere. The model has no soil "depth" (only a carbon pool) so 14 C cannot be used as a useful tracer as part of CLIMBER-2P in its current configuration. The CLIMBER-2P model does have a 13 C tracer within the carbon cycle which is intended to be used in conjunction with the permafrost model to constrain carbon cycle dynamics.
The possible impact of high dust concentrations on soil formation during glacial climates is not accounted for in the model. Loess soils, those created by wind-blown dust or alluvial soils, are not represented. For our study it is assumed that the ratio of loess to non-loess soils is the same in the present day as it was during glacial climates. This is not the case in the real world, where high dust concentrations in the dry atmosphere increased loess deposition at LGM (Frechen, 2011). However, the LGM climate is only representative of the coldest and driest part of the last glacial period. Evidence suggests that soils were productive in cold conditions in the permafrost region of the last glacial period, with loess accumulation only more widely significant towards the harsh conditions of the LGM (Elias and Crocker, 2008;Chlachula and Little, 2009;Antoine et al., 2013;Willerslev et al., 2014).
No changes were made to the vegetation model or to controls on soil input, which are only dependent on temperature and NPP; the mammoth-steppe biome is not explicitly modelled (Zimov et al., 2012).
Underneath ice sheets soil carbon is zero; as an ice sheet extends over a location with soil carbon (and vegetation), that carbon is released directly into the atmosphere. As an ice sheet retreats and exposes ground, the vegetation (and soil) can start to grow again. So, our model does not account for any carbon that may have been buried under ice sheets (Wadham et al., 2012).

Conclusions/summary
This permafrost-carbon model is a simplified representation of the general effect of frozen ground on soil carbon decomposition. In the presence of frozen ground the soil carbon decays more slowly. The method by which permafrost is diagnosed relies only on the balance between warm (above  0 • C) and cold (below 0 • C) days, which removes the problem of compounding errors in thermal diffusion calculations (for example). As such, the permafrost-carbon model would perform just as well in distant past climates as it does in pre-industrial climate. In order to account for uncertainties in carbon accumulation and release rates in frozen (and thawing) soils, a range of dynamic settings are retained that agree with total land carbon estimates of Ciais et al. (2012). Due to the slow accumulation in permafrost soils, soil carbon has a long time to equilibrium and therefore the present-day climate must be treated as a transient state, not as an equilibrium state. We showed that the model performs reasonably well at pre-industrial and present-day conditions. The permafrostcarbon model creates a mechanism that slowly accumulates soil carbon in cooling or cold climates and quickly releases this high soil carbon in warming climates, caused either by changes in insolation patterns or by global increases in temperature and climatic changes due to greenhouse gas feedbacks and ocean circulation changes. It can thus be used to quantitatively evaluate the role of permafrost dynamics on the carbon build-up and release associated with this specific physical environment, over supra-centennial to glacialinterglacial timescales. Appendix A: 1-D models Figure A1 shows the results of sensitivity experiments comparing these two approaches for one CLIMBER-2 land grid cell. Baseline settings of permafrost fraction = 0.6, frost index = 0.6 and mean annual air temperature = −10 • C have a relative soil carbon concentration of 1. The sub-grid method outputs a linear-type relationship between permafrost fraction and soil carbon stored. The remixing model outputs lower soil carbon concentration for lower fractional per- Figure A2. Relationships between frost index and MAT on the CLIMBER-2 grid scale (data from Zhang, 1998;Jones et al., 1999). Frost index determines permafrost fraction according to the model described in Sect. 3.2 (main text). NPP data for the permafrost zone from MODIS plotted against permafrost fraction (calculated from frost index values of Zhang, 1998) on the CLIMBER-2 grid scale. mafrost coverage rising quickly when permafrost fraction approaches 1. For the air temperature as variable, the two approaches show a similar response. For higher frost index the soil carbon concentration increases, with the sub-grid method showing slightly more sensitivity than the remixing model. The variables of permafrost fraction, frost index and mean annual temperature are interrelated, and covary. The relationships between these variables are shown in Fig. A2a. For permafrost fraction to frost index, the relationship is defined as that determined in the main text for the CLIMBER-2 grid scale in Sect. 3.2.
When the effect of NPP is included, the equilibrium total carbon contents are scaled according to the relationship between NPP and permafrost fraction. Figure A2b shows MODIS data for NPP plotted against frost index (calculated from data from Zhang (1998)  permafrost fraction (calculated from the frost-index value). The values are only for NPP in the high northern latitudes. To compare the model to data, it is assumed that 40 % of total soil column carbon is located in the top 1 m for permafrost soils . To convert SOCC (top 1 m) to full column, the SOCC data are multiplied by (2.5· permafrost_fraction). This soil carbon content is plotted against calculated permafrost fraction; that is, using the model from Sect. 3.2 to get permafrost fraction from frostindex data. These SOCC data are then binned into 0.1 increases in permafrost fraction and the mean value is shown with ±1 sigma in Fig. 7 (main text).

B1 Linear model
In more complex physical models, snow correction of ground temperature is achieved by modelling the thermal diffusion characteristics of the snow cover, a function of snow depth Figure B2. Model error when the linear snow correction model is used to predict temperatures at snow-depth or snow-ground interface for data from Morse and Burn (2010) (measurement data are down snow column temperatures). Positive numbers indicate that the linear model output is too warm compared to data. and snow type (for example snow density). A thermal diffusion model is used to make an estimate of the snowground interface temperature using the surface air temperature; the thermal gradient is also dependent on the initial snow-ground interface temperature. Within the CLIMBER-2 model, snow is already modelled (Petoukhov et al., 2000) as it has a significant effect on overall climate (Vavrus, 2007). Snow depth in CLIMBER-2 is available as well as snow fraction per cell, but snow type and snow density are not individually modelled. Attempting to model the thermal diffusion in the snow does not make sense for CLIMBER-2, as with permafrost location. Rather the approach is to use measurement data to create a general relationship between air temperature and snow-ground interface temperature based only on the snow depth.
The snow correction linear model is based on data from Taras et al. (2002) giving a correction for snow-ground interface temperature from snow depth and air temperature. Figure B1a shows the data from Taras et al. (2002) and the linear regressions (labelled as A, B and C) of these data re-plotted per snow depth (Fig. B1b). Equation (B1) shows this linear model for snow correction, which is only applied for surface air temperatures lower than −6 • C. This snow-ground interface temperature is used to calculate the freeze index (DDF sc ) in Eq. (4) in the main text.
where T gi is ground interface temperature ( • C), T surf is surface air temperature ( • C) and SD is snow depth (cm).

B2 Snow correction validation
This simple snow correction was tested against data from Morse and Burn (2010). Figure B2 shows the error made by the linear model when used to predict the snow-ground interface temperature (or snow depth temperature) from Morse and Burn measurement data. In the more extreme conditions, the error of the linear model is far higher, for example in deep snow and cold temperatures. Figure B3 shows the Geosci. Model Dev., 7, 3111-3134, 2014 www.geosci-model-dev.net/7/3111/2014/ Figure B3. CLIMBER-2 model output for snow depth (m) plotted against surface air temperature ( • C) for the PI(eq) (green circles) and LGM(eq) (blue squares) climates. Model output does not show extreme conditions for snow cover due to the very large grid-cell size.
outputs from CLIMBER-2 for snow depths plotted against surface air temperatures for the PI(eq) pre-industrial climate and LGM(eq) glacial climate, for all grid cells. The large CLIMBER-2 grid size means that extreme conditions are not present in the model output. Comparing Figs. B2 and B3 shows that the linear correction can provide an estimated confidence within ±8 • C for the deepest snow cover and highest temperatures of CLIMBER-2P data output, and within ±2 • C for the majority of CLIMBER-2P data outputs. A similar performance is found when comparing to snow thickness and snow-ground interface temperatures from Zhang (2005) for a site in Zyryanka, Russia. The most extreme temperatures and snow conditions produce a larger error from the linear model, but the intermediate conditions, those seen in CLIMBER-2P data points, agree better with the data. Overall the effect of the snow correction within the model produced a maximum decrease in permafrost area of 8 % (compared to the uncorrected version) in the most affected grid cell for the PI(eq) simulation and is therefore significant. Table C1 shows all the settings for "a" and "b" per soil pool (Eq. (3), main text) that were tested to obtain total soil carbon contents for the LGM and the PI simulations. Figure C1 shows the modelled total land carbon (GtC) for all simulations sorted by permafrost area function. Green dashed lines on the LOW-MEDIUM area setting indicate the dynamic settings chosen to represent the "slow", "medium", "fast" and "extra-fast" permafrost-carbon dynamic settings. The total land carbon content is clearly very sensitive to permafrost area, and despite many simulation tunings only the LOW-MEDIUM area setting provided a good enough range of dynamics that could be used to later investigate the permafrostcarbon dynamics. Within the settings chosen, the "medium" dynamic setting overestimated the present-day total land carbon estimate from Ciais et al. (2012), but further tuning experiments did not improve this overestimate.