Advancement toward coupling of the VAMPER permafrost model within the Earth system model i LOVECLIM ( version 1 . 0 ) : description and validation

The VU Amsterdam Permafrost (VAMPER) permafrost model has been enhanced with snow thickness and active layer calculations in preparation for coupling within the iLOVECLIM Earth system model of intermediate complexity (EMIC). In addition, maps of basal heat flux and lithology were developed within ECBilt, the atmosphere component of iLOVECLIM, so that VAMPER may use spatially varying parameters of geothermal heat flux and porosity values. The enhanced VAMPER model is validated by comparing the simulated modern-day extent of permafrost thickness with observations. To perform the simulations, the VAMPER model is forced by iLOVECLIM land surface temperatures. Results show that the simulation which did not include the snow cover option overestimated the present permafrost extent. However, when the snow component is included, the simulated permafrost extent is reduced too much. In analyzing simulated permafrost depths, it was found that most of the modeled thickness values and subsurface temperatures fall within a reasonable range of the corresponding observed values. Discrepancies between simulated and observed permafrost depth distribution are due to lack of captured effects from features such as topography and organic soil layers. In addition, some discrepancy is also due to disequilibrium with the current climate, meaning that some observed permafrost is a result of colder states and therefore cannot be reproduced accurately with constant iLOVECLIM preindustrial forcings.


Introduction
The VU Amsterdam Permafrost (VAMPER) model is a deep 1-D heat conduction model with phase change capability.It has been previously validated for single site experiments such as in Barrow, Alaska (Kitover et al., 2012).Subsequently, it has simulated both equilibrium and transient permafrost conditions at a number of arctic/subarctic locations (Kitover et al., 2012(Kitover et al., , 2013)).The VAMPER model was built with the intention of coupling it within iLOVECLIM, an Earth system model of intermediate complexity (EMIC).Using this coupling, the goal is to capture the transient nature of permafrost growth/decay over millennia as a feedback effect during major periods of climate change.To prepare for coupling, a number of enhancements have since been made to the VAMPER model.We present validations of these improvements by simulating modern-day permafrost thickness and distribution.The goal of this paper is to describe the enhancements and then analyze the validation experiments for modeling present-day permafrost, with a detailed explanation of why mismatches occur between simulated and observed data.
The first example of VAMPER as a stand-alone deep permafrost model was for Barrow, Alaska (Kitover et al., 2012), where VAMPER simulations reproduced the presentday permafrost depth using monthly averaged observation data of ground "surface" ( − 1 cm deep) temperatures.In this same study, VAMPER was also validated by comparing results against other developed deep permafrost models (also used for millennial-scale simulations) using similar forcings and parameter settings.In both Kitover et al. (2012) and Kitover et al. (2013) at selected locations (e.g., Wyoming, West Siberia, Central Siberia) were performed using the stand-alone version of the VAMPER model, forced by iLOVECLIM-generated land surface temperatures over the last 21 000 years (Roche et al., 2011).In addition, a sensitivity analysis was presented in Kitover et al. (2013), showing the range of simulated permafrost depths under different parameter settings.
Thus far, according to the work summarized above, VAM-PER has only been employed as a site-specific permafrost model.As a next step, this paper describes the necessary developments and validation to couple VAMPER with EC-Bilt, the atmospheric component of iLOVECLIM.Specifically, this presented work introduces two enhancements to the VAMPER model: (1) inclusion of snow as optional layers and (2) change in the time step.The first in particular is an issue in modeling permafrost since snow cover is a recognized influence on the ground thermal regime (Williams and Smith, 1989) and was not an available option in the previous VAMPER model version.To compensate for this, Kitover et al. (2013) had artificially introduced the effect of snow cover via a surface offset (the difference between surface air temperature and ground temperature) of +2 • C.Not only was this an assumption based on a number of previous reports and observations, but it had to be applied as an annual surface offset since the time step was 1 year.This subsequently demonstrates the need for the other enhancement, which is a sub-annual time step, where the seasonal changes in the ground thermal conditions can be captured, allowing for representation of both the snow cover effect and the active layer.In addition to the VAMPER model enhancements, two global maps were produced (geo-processed from the original maps to fit the horizontal grid of Earth system models of intermediate complexity Climate deBilt (ECBilt)) to be used as additional input parameters to the VAMPER model: geothermal heat flux and porosity.These are used when the VAMPER is run over a horizontal grid, to allow these parameters to vary spatially.
Integrating permafrost into Earth system models has attracted increased interest since research has acknowledged the effect of climate change on permafrost temperatures (Cheng and Wu, 2007), permafrost degradation (Anisimov and Nelson, 1996), and release of carbon stored within the permafrost (Davidson and Janssens, 2006).The Coupled Model Intercomparison Project phase 5 (CMIP5) (Koven et al., 2013) analyzed how different Earth system models represent the subsurface thermal dynamics and how well this class of models simulate permafrost and active layer depth.Despite the fact that there is a variety of modeling methods and configurations for the different global coupled models, the conclusion was that there is no clear ranking among the reviewed 15+ model versions.This shows that representing permafrost in Earth system models still has some challenges, which Koven et al. (2013) attribute primarily to modeling of both the atmosphere-ground energy exchange and the subsurface thermal regime.Until recently, most simu-lations of permafrost were calibrated for regional or local study such as Li and Koike (2003) on the Tibetan Plateau, Zhang et al. (2006) in Canada, andNicolsky et al. (2009) in Alaska.A growing number of studies are now modeling permafrost across the Northern Hemisphere or globally.Simulations are done using either statistical approaches like the frost index method (Anisimov and Nelson, 1996;Stendel and Christensen, 2002) or climate models such as Dankers et al. (2011) who used the Joint UK Land Environment Simulator (JULES) land surface model and Ekici et al. (2014) (GCM) to model the connection between permafrost degradation and subsequent carbon emission as a driver for the occurrence of the Palaeocene-Eocene Thermal Maximum (PETM).Modeling permafrost changes is also of interest from the hydrological perspective.Avis et al. (2011) used a version of the University of Victory (UVic) Earth System Climate Model to examine the potential decreasing areal extent of wetlands due to future permafrost thaw.
However, it should be noted that there is a difference between coupled models which actively integrate the role of permafrost (including the thermal, hydrological, and/or carbon feedbacks) (Lawrence et al., 2011), and models which look at permafrost in a post-processing perspective (e.g., Buteau et al., 2004;Ling and Zhang, 2004) meaning they are forced by the predicted temperature changes.It is the full coupling with integrated feedbacks which is a current interest of ours, where the goal is to fully couple ECBilt and VAMPERS within iLOVECLIM.The results of the work presented here serve as an important validation stage toward this goal.In the sections following, the two enhancements to the VAMPER model are explained.This includes validation of the time step change by comparing simulated annual active layer depths with empirical-based estimates.In addition, two newly developed maps of spatially varying parameters used in the VAMPER experiments are explained.For the validation, the VAMPER model is forced by ECBilt land surface temperatures, where the results are compared with a modernday map of permafrost extent in the Northern Hemisphere and observed permafrost thickness and subsurface temperatures values in boreholes.

General description
VAMPER is a 1-D permafrost model developed to estimate permafrost thickness and was designed for eventual full coupling within iLOVECLIM.Because it must fit a relatively coarse Earth system model, it is not suitable for the soil and subsurface parameters to undergo parameterization schemes.These characteristics such as soil type, organic matter, and water content are observed and vary at a much finer spatial scale than possibly represented in iLOVECLIM.VAMPER is meant rather as a generalized model to simulate conceptual permafrost thickness based on the factors which most strongly dictate the subsurface thermal regime; most notable for our purposes and discussed by Farouki (1981), these factors are mineral composition, water content, and temperature.
Other than what is specified below, construction of the VAMPER model has not changed and the methods as described in Kitover et al. (2013) still apply.In particular, these include assuming only conductive heat transfer in the subsurface and employing well-established methods for finding the temperature-dependent thermal properties of heat capacity and thermal conductivity (Farouki, 1981;Zhang et al., 2008).The subsurface is assumed to be saturated (i.e., porosity equals the water content) and there is currently no groundwater flow either horizontally or vertically between the soil layers.
The phase change process of freeze/thaw in the subsurface is handled using a modified apparent heat capacity method from Mottaghy and Rath (2006).Their method assumes that phase change occurs continuously over a temperature range, which in our case is approximately between 0 and −2 • C. The apparent heat capacity method includes an additional latent heat term in the heat diffusivity equation as a way to account for the added energy released (consumed) during freeze (thaw) of the subsurface water content.The latent heat demand during phase change, referred to as the "zero curtain effect", slows thermal diffusivity rates near the surface as the active layer freezes and thaws but also during permafrost degradation/aggradation.

VAMPER model enhancements
There are few permafrost modeling studies which have reproduced changes in permafrost thickness over geologic time periods.In these cases, a larger time step in their numerical simulations (usually 1 month or 1 year) (e.g., Osterkamp and Gosink, 1991;Lebret et al., 1994;Lunardini, 1995;Delisle, 1998) is assumed since they only need to force the models with the low-frequency changes in air temperature or ground temperature that occur over millennia.At this timescale, it is not necessary to use a sub-annual time step.In our ear-lier work with the VAMPER model (Kitover et al., 2013), we similarly used a yearly time step.However, in light of the future coupling between ECBilt and the VAMPER model, it has become clear that the VAMPER model should run on a 4 h time step.Doing this allows the VAMPER model to match the timescale of the atmosphere, the subsystem to which the VAMPER model will be coupled.Changing to a 4 h time step also reduces error in the numerical approximation, since the change in thermal properties, which are temperature-dependent, is smoother between each time step.Fortunately, being that the VAMPER model is somewhat simplified, and hence flexible, the change to a 4 h time step only required revalidating the model performance.In addition to the change in time step, we also (newly) made it possible that the VAMPER model had an overlying snowpack.Including this option is meant to simulate the effect of thermal insulation of the ground in winter.
Note that the VAMPER model with the snow enhancement is referred to as the VAMPERS model.When referring to both versions or either version separately, the "VAMPER(S)" term is used.

Time step
To illustrate the difference between applying the same annual average temperature forcing but with two different time steps (4 h vs. yearly), a sensitivity test was performed (Fig. 1a).To generate the sub-daily surface temperature forcing (4 h), a year-long temperature time series was calculated using a standard sine function with constant amplitude 20 • C and average annual temperature of −6 • C (hereafter referred to as sensitivity run 1 or "sr1"), resulting in an annual range of temperatures between −26 and 14 • C. Therefore, the case with a yearly time step, called "sr2", used −6 • C as the constant forcing.Besides the change in time step and corresponding surface temperature forcing, the thermal conductivity and heat capacity values were also allowed to differ since these variables are temperature dependent (Fig. 1b).However, heat flux and porosity parameter settings were the same in both model runs.Each experiment was run until approximate equilibrium was reached under the same constant (respective) forcing.We consider equilibrium to be when the geothermal heat flux is approximately equal to the ground heat flux.Comparing the final depth-temperature profiles between sr1 and sr2 shows a shift in the equilibrium depth-temperature profile where using an annual time step underestimates permafrost thickness by approximately 50 m (Fig. 1a).This difference is attributed to the occurrence of the thermal offset (difference between ground temperature and top of the permafrost) within the active layer in sr1 (Fig. 1b), whereas sr2 cannot exhibit such seasonal phenomena.Because VAMPER is a simple model (absence of vegetation, organics, an unsaturated subsurface, or temporally varying water content) we can attribute the thermal offset to seasonal differences in thermal conductivity, since the thermal conductivity of ice is 4 times that of unfrozen water and therefore the freezing front is propagated more effectively than the warming front.This difference causes the mean annual subsurface temperature within the active layer to be gradually colder with depth.The offset is visible in the mean annual depth-temperature profile within the top meter of Fig. 1b.

Active layer
In permafrost modeling, an active layer can only be present when the air-ground temperature forcing varies seasonally.Thus, the time step must be sub-annual.Since a 4 h time step is now implemented, the VAMPER model produces an active layer.It is necessary within the framework of model development to then check the simulation of this active layer for validation purposes.
Most dynamical permafrost models that simulate nearsurface behavior configure the parameter settings to specifically match locally observed data.Some parameterizations include organic and mineral layer thicknesses, which give soil properties such as porosity and bulk density, and unfrozen water content characteristics.Examples of these sitespecific studies include, for example, Romanovsky and Osterkamp (1995), Buteau et al. (2004), Ling, Zhang (2004), and Zhang et al. (2008), and Nicolsky et al. (2009).Since VAMPER is not parameterized to capture site-specific behavior, it is challenging to assess the ability of the model to simulate active layer dynamics.Here, we leverage the Stefan equation, used originally in engineering applications (Fox, 1992), to estimate the thickness of the active layer when the amount of energy input and thermal characteristics are known.From French (2007), the Stefan equation is defined as where AL (m) is the thickness of the active layer, σ is the cumulative thawing index (average ground surface temperature ( • C) during the thaw season times the duration of thaw season (s)), and k mw is the thermal conductivity of unfrozen soil (W(m K) −1 ).Q i (J m −3 ) is defined further as where L is the latent heat of fusion, ρ m is the dry density of the soil (kg m −3 ), W is the total moisture content, and W u is the unfrozen water content.Table 1 gives the constant variable values applied in the Stefan equation, which are the same values used in a comparable run for the VAM-PER model.Under different forcings as a function of both average annual ground surface temperature and annual amplitude, the VAMPER model's active layer thickness versus results using the Stefan equation are shown in Table 2.We suggest that when comparing the empirically based results with the series of simulations, the VAMPER model does a suitable job of reproducing annual active layer thickness.

Snowpack parameterization
An additional option to the VAMPER model is the ability to extend the heat conduction model into the snowpack when present.Prior to this, the surface offset, illustrated in Smith and Riseborough (2002), could not be produced in the VAM-PER model.
−2 20 1.9 1.9 The VAMPERS model uses snow water equivalent (SWE) values (m) with corresponding densities to compute snow thickness layers.SWE is the depth of water that would result from the complete melting of snow.The precipitation simulated in ECBilt is computed from the precipitable water of the first atmospheric layer (Goosse et al., 2010).When the air temperature is below 0 • C, the precipitation is assumed to be snow.However, this "snow" is only assumed to be frozen water, meaning it lacks any quantifiable properties besides the actual precipitation amount, and as such is directly considered the SWE value.As a result, there is an additional set of necessary functions when coupled with VAMPERS to transfer ECBilt SWE values into a snowpack thickness (Z) at time t: where ρ w is water density and ρ s snow density (Lynch-Stieglitz, 1994).The total snow density is determined as a combination of old snow (expressed as SWE t−1 from the previous time step) and freshly fallen snow at the current time step (expressed as SWE fr ): where ρ fr is the density of fresh snow (150 kg m −3 ).
There is snowpack metamorphism that occurs from a number of different processes.Notably, Dingman (2002) distinguished these as gravitational settling, destructive, constructive, and melt.However, as these different changes occur at highly varying rates and under localized conditions (aspect, slope, vegetation cover), it is difficult to incorporate such processes in an EMIC such as iLOVECLIM.On the other hand, a snowpack always undergoes densification over time and this effect should somehow be applied to the modeled snowpack.Therefore, we apply to the total snow density an empirical densification function due to mechanical compaction.The maximum allowable density is 500 kg m −3 , which cannot hold any more liquid water (Dingman, 2002).The compaction equation used (e.g., Pitman et al., 1991;Lynch-Stieglitz, 1994) is as follows: where g is gravity (9.82 m s −2 ), N (kg) is the mass of half the snowpack, T ( • C) is the temperature of the snowpack (the average temperature of the snow layer temperatures from the previous time step), and t is the time step (s).Three snow layers are then discretized from the total snow thickness, depending on whether it is above or below 0.2 m, as outlined in Lynch-Stieglitz (1994).Thermal properties are then calculated for each snow layer based on empirical formulas K s = 2.9ρ 2 s (7) (Goodrich, 1982), C s = 1.9 × 10 6 ρ s /ρ f (8) (Verseghy, 1991), where K s is the snow thermal conductivity, C s is the snow heat capacity, and ρ f is the density of ice (920 kg m −3 ).All three snow layer are subject to the same processes and simply depend on temperature, time, and thickness for their respective deformation and/or melting.
The following is a stepped description of the snow algorithm to generate a VAMPERS snowpack from ECBilt precipitation: 1. Calculate new snow density, Eqs. ( 4) and (5), using any freshly fallen snow and old snow.
4. Discretize the individual layer thicknesses based on total snow thickness.
6. Use snow thicknesses and corresponding thermal properties as additional layers in the VAMPERS model.

General description
iLOVECLIM is a "code-fork" of LOVECLIM 1.2 (Goosse et al., 2010), both which belong to a class of climate model called EMIC (Claussen et al., 2002).This type of model, as summarized by Weber (2010), "describes the dynamics of the atmosphere and/or ocean in less detail than conventional General Circulation Models".This simplification reduces computation time, thus making EMICs suitable for simulations on millennial timescales that incorporate the components with slow feedback effects, such as ice sheets, vegetation, and permafrost.Different versions of LOVECLIM have successfully simulated past climates including the Last Glacial Maximum (LGM) (Roche et al., 2007), the Holocene (Renssen et al., 2005(Renssen et al., , 2009)), and the last millennium (Goosse et al., 2005).Although there exist some different developments between iLOVECLIM and the LOVECLIM versions, both consist of the following coupled Earth system components: the atmosphere (ECBilt), the Coupled Large-scale Ice-Ocean Model (CLIO), and vegetation (Vegetation Continuous Description -VECODE) (Fig. 2).ECBilt, the atmospheric model (Opsteegh et al., 1998), consists of a dynamical core with three vertical levels at 800, 500, and 200 hPa.
It runs on a spectral grid with a triangular truncation (T21), which translates to a horizontal grid with a resolution of approximately 5.6 • lat × 5.6 • long.The CLIO module (Goosse and Fichefet, 1999) is a 3-D ocean GCM with a free surface.
It has 3 • × 3 • horizontal resolution and 20 vertical layers.VECODE, the vegetation module (Brovkin et al., 1997), is similar to VAMPER(S) in that it was particularly designed for coupling to a coarse-resolution Earth system model.It is a reduced-form dynamic global vegetation model that characterizes the land surface as either trees, grass, or no vegetation (i.e., "bare soil") and is computed at the same resolution as the ECBilt grid.The plant types may be represented fractionally within each grid cell.Each model component of iLOVECLIM was originally developed separately and the reader is referred to Goosse et al. (2010) for a detailed description of components and coupling mechanisms.Furthermore, iLOVECLIM more recently includes other optional components including the dynamical ice-sheet model GRISLI (GRenoble Ice Shelf and Land Ice model) (Roche et al., 2014) and a stable water isotopes scheme (Roche, 2013).

Proposed ECBilt-VAMPER(S) coupling description
The VAMPER(S) model will ultimately be coupled to the atmospheric component, ECBilt within iLOVECLIM.The proposed ECBilt-VAMPER(S) coupling will be done at each time step (4 h) where the land surface temperature from EC-Bilt is passed to VAMPER(S) and the ground heat flux from VAMPER(S) is returned to ECBilt.The land surface temperature is calculated within ECBilt as a function of the heat balance equation where the major heat fluxes across the air-surface interface are incorporated: sensible heat flux, latent heat flux, shortwave radiation, long-wave radiation, and ground heat flux.The land surface temperature and ground heat flux will only be communicated between components when the respective grid cell is classified as land with no overlying ice sheet (i.e., Greenland/Antarctica at present day).With this coupling, the effect of changing permafrost conditions may be reflected in the climate via changes in the surface energy balance.As permafrost degrades, the subsurface acts as a thermal sink, absorbing additional energy to accommodate latent heat demands during phase change.However, at the same time, the active layer deepens, also redistributing the (seasonal) energy distribution at the surface.Since the VAMPER(S) model ground surface temperature is taken to be the ECBilt land surface temperature, no surface offset occurs except when there is a snowpack.In this case, the snow surface temperature (i.e., the top snow layer) is assumed to be the same as the land surface temperature.This means the VAMPERS model ground temperature forcing is buffered via the three snowpack layers as discussed in Sect.2.1.2.Using the ground surface temperature forcing, the VAMPER(S) model then computes the subsurface temperature profile.This calculation, via the implicitly solved heat equation with phase change capability, is fully described in Kitover et al. (2013).As VAMPER is a 1-D model, there is no lateral energy (heat/water) transfer between adjacent grid cells in the subsurface.Permafrost thickness is determined at an annual time step using a computed average annual temperature profile, where any depth below or equal to 0 • C is considered permafrost.Although in reality there is a freezing point depression which may occur as a result of the local pressure or dissolved salts, our permafrost definition is consistent with the thermal definition of permafrost from the International Permafrost Association: "ground (soil or rock and included ice or organic material) that remains at or below 0 • C for at least two consecutive years" (IPA, 2014).
The land surface of ECBilt consists of a single "layer" which represents a volumetric storage capacity to generate surface runoff when full.This system is referred to as a bucket model in previous text (Roche et al., 2014;Roche, 2013;Goosse et al., 2010).Currently, this hydrology portion of ECBilt is not coupled to VAMPERS.However, because the active layer is a regulator of hydrology in arctic and subarctic regions (Genxu et al., 2009;Hinzman and Kane, 1992), a next step will be to expand coupling between VAMPERS and ECBilt by connecting the active layer with this bucket model.
The first phase of the coupling between VAMPERS and ECBilt will only include the land surface temperature and the ground heat flux as discussed.It should be mentioned as a caveat that additional coupling mechanisms are possible between iLOVECLIM components and VAMPER, which include hydrology and the carbon cycle, but will not be implemented for the first coupling phase.

Geothermal heat flux
The VAMPER(S) model requires a geothermal heat flux as the lower surface boundary.In Kitover et al. (2013), a sensitivity analysis was performed to look at the equilibrium permafrost thickness as a result of varying the geothermal heat flux and found that thickness can increase by about 70 m with every decrease in flux of 10 mW m −2 .To obtain the geothermal heat flux for every cell in the ECBilt grid, we used the recent publication of Davies (2013) who determined the median of heat flux estimates over a 2 • × 2 • latitude-longitude grid based on a combination of actual measurements, modeling, and correlation assumptions.Due to the mismatch of grid resolutions between Davies (2013) and ECBilt, we determined for each ECBilt grid cell, a simple area-weighted average of the Davies (2013) estimates: each of the Davies grid cells was assigned a weighing factor based on the percentage of overlap with the ECBilt cells.Below is the original map from Davies (2013) and the averaged map applied in the VAMPER(S) experiments.A sensitivity analysis with (1) the geothermal heat flux map and (2) applying the continental global average (approx.60 mW m −2 ) showed no noticeable difference in permafrost distribution.This result is different, however, than the noticeable sensitivity of geothermal heat flux on permafrost depth (Kitover et al., 2013).

Porosity
Another variable needed to run the VAMPER(S) model are depth-dependent porosity values, which in these experiments are 3000 m below the surface.In previous VAMPER studies  (Kitover et al., 2013(Kitover et al., , 2012)), it was assumed that the land subsurface was sedimentary rock, with a porosity of 0.3, 0.4, or 0.5.However, as shown in Kitover et al. (2013), the porosity, or water content, has a noticeable effect on equilibrium permafrost thickness.That sensitivity test showed about a 50 m difference in permafrost thickness when the porosity values (assuming a saturated subsurface) ranged between 0.3 and 0.5.Therefore, to both narrow our assumptions regarding the subsurface but still maintain the simplification necessary for the coarse horizontal grid, an additional lithological classification scheme was created as an additional VAMPER(S) model parameter.We reclassified the original seven categories from the Global Lithological Map Database (GLiM) from Hartmann and Moosdorf (2012) into "Bedrock (Bed)", (e.g., granitic and metamorphic rock), and "Sedimentary (Sed)" (e.g., sandstone, limestone) (Table 3, Fig. 5).In the case of Bed, the subsurface is assumed to be quite consolidated/compressed, resulting in a low water content (Almén et al., 1986;Gleeson et al., 2014).Bed was thus assigned a low porosity of 0.1, which based on sources that showed depth profiles of bedrock sites (Schild et al., 2001;Nováková et al., 2012), that stays constant with depth.On the other hand, similar to the case studies from Kitover et al. (2013), a depth porosity function from Athy (1930) was applied for the Sed class, where the surface porosity ( ) is 0.40.This porosity represents the assumed average for sandy textured soil.Similar to application of the geothermal heat flux map, a preliminary sensitivity analysis between applying the lithology map and applying a constant value (0.4) throughout the globe showed only marginal differences in permafrost distribution.This result is different, however, than the higher sensitivity of porosity on permafrost depth (Kitover et al., 2013).3 Validation of preindustrial permafrost thickness distribution

Experimental setup
The model experiments are performed over the whole globe, with the VAMPER model forced by ECBilt land surface temperatures.These values are the lower boundary layer of the atmosphere and are calculated using a surface heat budget (Goosse et al., 2010).Referring to Fig. 3a, this means that ECBilt passes temperature values to the VAMPER(S) model (right side of Fig. 3a) but no data are returned to ECBilt (left side of Fig. 3a), leaving the climate unaffected from permafrost or changes in permafrost.The model experiments also include the spatially varying parameter values of geothermal heat flux and porosity provided by the new maps (described in Sect.2.2.3 and 2.2.4).Two different model runs were performed: one without the snow enhancement or any imposed surface offset (VAMPER) and one with the snow enhancement (VAMPERS).These two are first compared in Sect.3.2.1 of the "Results and discussion" below.
Because permafrost has a very slow thermal response (Lunardini, 1995) as compared to other components in iLOVECLIM, VAMPER(S) is not forced synchronously by ECBilt.Rather, VAMPER(S) is forced continuously for 100 years and then runs offline for 900 years using the ECBilt average land surface temperature of the previous 100 years as the forcing.This asynchronous cycle is repeated for thou- Equilibrium was determined when the lower boundary heat flux approximately matches the annual average ground surface heat flux and the permafrost thickness stabilized.Although the model approaches a steady state through the subsurface depth, we acknowledge that in reality, some observed permafrost regions are not at equilibrium since they are responding to recent warming.

Results and discussion
In order to verify the performance of VAMPER(S) forced by iLOVECLIM, a series of equilibrium experiments were performed for the preindustrial (PI) climate (∼ 1750 AD).For comparative purposes, we assume the PI state of permafrost is similar enough to the current state of permafrost that we used modern-day data to validate against the PI simulations.The simulated areal extent was compared to presentday extent using the "Circum-Arctic Map of Permafrost and Ground-Ice Conditions" (Brown et al., 2014).Unlike the model validation done by Lawrence and Slater (2005), and then subsequently critiqued by Burn and Nelson (2006), our simulations attempt to capture the extent of both continuous and discontinuous permafrost.In addition, available borehole data, for sites within the arctic/subarctic, were used to evaluate the simulated thicknesses.Therefore, there are two types of validation approaches: (1) permafrost distribution and (2) permafrost depth.

Permafrost distribution validation
The first validation demonstrates the extent to which the VAMPERS model reproduces the modern-day permafrost distribution.The results can be matched against Koven et al. (2013), who simulated permafrost areas consistent with CMIP5 model output.The areal extent of permafrost distributions found in Koven et al. (2013) bracket the extent found in the present study.The maximum is reported as 28.6 × 10 6 km 2 and minimum 2.7 × 10 6 km 2 .Our simulation using VAMPERS yields approximately 20.3 × 10 6 km 2 .This is a reasonably comparable estimate considering almost 80 % (14/18) of the model area extents from Koven et al. (2013) fall within 40 % (12-28 × 10 6 km 2 ) of our model estimates.According to the discussion by Koven et al. (2013), most of the variation seen among the compared Earth system models is primarily attributed to the subsurface modeling techniques, such as water content, using a latent heat term, and differing soil thermal conductivities.Secondary causes are attributed to factors of air-ground coupling such as incorporation of  organics and a snowpack (bulk or multilayer).These conclusions are not different from our own study in that (1) snowpack plays a marked role in permafrost modeling and inclusion/exclusion will impact the results, and (2) the air-ground coupling is also a source of potential mismatch (discussed further in Sect.3.2.2).
Using the comparison shown in Fig. 7, which overlays the simulated results on the map from Brown et al. (2014), it is clear that the experiment without the snow option overestimates permafrost extent while employing the VAMPER(S) version underestimates it.This inaccuracy between both an overestimated result and an underestimated result is at least partially due to attempting to match results from a lowresolution grid to spatial coverage of much higher resolution.Because the marginal areas of permafrost extent are the most sensitive to climate, they are highly responsive to minor temperature deviations.These deviations, whether a few degrees above or below freezing, determine from a modeling point of view, whether permafrost exists or not.In the case of VAMPER, many of these marginal grid-cell average annual ground surface temperatures fall below freezing while in the case of VAMPERS, these same grid cells now fall above freezing.However, because of the coarse grid, these estimates in either case, cannot accurately represent areas which are only partially underlain by permafrost.
Inaccuracy in model results is also expected since we cannot parameterize the snowpack characteristics that alter the effect of snow on the ground thermal regime.Although we capture the role of snow cover, which imposes a reduced thermal diffusivity effect between the air and ground, there are number of snowpack characteristics that we do not include.As opposed to our generalized snowpack parameterization scheme, described in Sect.2.1.1,high-resolution snow models are fitted to observational data by analyzing, for example, the physics of accumulation, areal distribution, and snowsoil interactions.Therefore, it is arguable from this lack of detail and the results shown in Fig. 7, whether the better option is to include a snowpack in VAMPERS or not.However, we contend that the VAMPERS model is doing a reasonable job since it is producing the surface offset that would naturally occur from the snowpack (Goodrich, 1982;Smith and Riseborough, 2002).The simulated global distribution of this surface offset is shown in Fig. 8.It is determined by calculating the difference between the mean annual ground temperature (MAGT) using VAMPERS and the MAGT using VAMPER (no snow option and no imposed surface offset).Although the maximum mean annual surface offset is about 12 • C, the average among all the grid cells with snow cover is about 2.7 • C, close to our original applied surface offset of 2 • C in Kitover et al. (2013).Values between 1 and 6 • C were reported by Gold and Lachenbruch (1973).Monitoring studies of the air-ground temperature relationship also fall within this range, e.g., Beltrami and Kellman (2003), Bartlett et al. (2005), Grundstein et al. (2005), and Zhang (2005).However, larger values of 10 • C have been recorded in Alaska (Lawrence and Slater, 2010).
Further, without the snow option, changing precipitation patterns due to climate change would otherwise have no effect on the subsurface thermal conditions.In other words, the role of snow cover will be more noticeable when us-   (Brown et al., 2014).
ing the ECBilt-VAMPERS coupling in transient simulations.An example of the effect of changing snow conditions on the ground thermal regime come from Lawrence and Slater (2010), who demonstrated through experiments with the Community Land Model that (1) increased snowfall accounted for 10 to 30 % of soil warming and (2) a shortened snow season also caused soil warming due to the ground surface's increased uncovered exposure to air temperatures.From this point forward, all analysis in this study is performed on results from VAMPERS (i.e., with the snow option).
In addition to the snowpack-induced surface offset, there are a number of additional factors which have been commonly recognized in affecting the surface offset and hence should be part of the air-ground coupling.Depending on the scale of interest, the magnitude of these can vary but they include surface organic layer, vegetation, overlying water bodies, and wind.It should be recognized that within ECBilt, some of these factors are reflected in the land surface temperature (notably wind and a simplified vegetation scheme) but the others are absent.In addition, coupling the ECBilt surface hydrology to the groundwater storage would affect both the ground thermal regime and hydrological regime.In the first case, subsurface water content affects the thermal properties of the soil.In particular, the conductivity of organics have high variation seasonally.In the second instance, frozen ground is impermeable, allowing little or no subsurface water storage, in turn affecting runoff flow rates and timing.

Permafrost thickness validation
The second validation examines the simulated depth of permafrost using borehole data taken from the Global Terrestrial Network for Permafrost (GTN-P; www.gtnp.org).Figure 9 regresses observed borehole measurements mapped in Fig. 10 against the corresponding permafrost depths simulated by iLOVECLIM.It is clear that there is a larger divergence between modeled and observed depths for the deeper permafrost than for the more shallow observations, where some points are overestimated by over 300 m and some very underestimated by over 700 m.There are a number of reasons to explain the mismatch, which can manifest in both the borehole and the model data.The first reason is that borehole estimates have a given range of uncertainty since measurement techniques and subsequent interpretations are subject to error.Osterkamp and Payne (1981) described in detail potential errors associated with the freezing point depression, thermal disturbance, and lithology.
The second reason is that we assumed implicitly that the observed permafrost depths are at equilibrium with the current (or PI) climate state.This likely explains the mismatch at the central Siberian site (66 • 26 2 N, 112 • 26 5 E) (point 1, Fig. 9), where the permafrost is estimated from the borehole data to be 1000 m thick while the corresponding modeled value is only about 375 m.Like much of the Siberian permafrost, this permafrost probably developed from the preceding glacial period (Kondratjeva et al., 1993).Another example is western Siberia, (points 2 through 4, Fig. 9), which is an area documented for having relict permafrost (Zemtsov and Shamakhov, 1993;Ananjeva et al., 2003).It is also identified in the "Circumarctic Map of Permafrost and Ground-Ice Conditions" (Brown et al., 2014) and "The Last Permafrost Maximum (LPM) map of the Northern Hemisphere" (Vandenberghe et al., 2014).But it should be noted that not all the relict permafrost in western Siberia is of late Pleistocene origin and may be from earlier cold stages (Zemtsov and Shamakhov, 1993;French, 2007).
Another reason for discrepancies between modeled and observed data is that high-resolution features in the landscape and topography cannot be captured by iLOVECLIM due to the limited spatial resolution.Such factors as vegetation and the organic layer, which can vary due to local topography and micro-climatic conditions, have been shown to affect the ac- tive layer and ground thermal regime (Shur and Jorgenson, 2007;Fukui et al., 2008;Lewkowicz et al., 2011;Wang et al., 2014).Consequently, given a specific borehole site, some discrepancy in the permafrost thickness estimate will likely occur between our simplified interpretation and that which results from including more complex and local interactions.It is possible, for example, that the observed value for point 5 (720 m) is a function of higher elevation since it is from a borehole site in the Russia Highlands but this relatively local elevation effect may not be sufficiently represented in the iLOVECLIM surface temperatures, and hence is underestimated.
The other outlying points (points 6 and 7, Fig. 9) occur in Canada, but as opposed to the relict sites as mentioned above, here iLOVECLIM overestimates permafrost thickness.These discrepancies, both occurring at high latitudes of 80 and 76 • N, reveal that VAMPERS is not reproducing the subsurface temperatures well for this area.For example, a report for the specific borehole (Gemini E-10; point 6, Fig. 9) calculated the geothermal gradient to be approximately 0.04 • C m −1 (Kutasov and Eppelbaum, 2009) whereas our model result for the corresponding grid space found a gradient of approximately 0.03 • C m −1 .Although this difference is relatively small, it hints at either a necessary increase in the averaged geothermal heat flux used in the model or a change in the subsurface thermal properties (increase in thermal conductivity), which could be altered by an adjustment in the VAMPERS water content.

Climate analysis
Finally, the remaining possibility to explain inaccuracies between the modeled results and the observed results (both in reproducing spatial extent and permafrost thickness) is the iLOVECLIM climate.Results of the VAMPER(S) model, above all other parameter settings, are most dependent on the mean annual ground surface temperature, as shown in the sensitivity study from Kitover et al. (2013), so if there exists biases or discrepancies within the forcing, it will be reflected in the output.For this portion of our analysis, we took observed MAGT measurements again from the GTN-P (International Polar Year (IPY) Thermal State of Permafrost Snapshot, IPA 2010) and regressed these values against the corresponding simulated MAGT at the same approximate depth and location (Fig. 11).Figure 12 shows a map of the selected GTN-P measurements.All the temperature comparisons are within the top 30 m of the subsurface and therefore reflect recent climate as opposed to the deeper temperatures (i.e., > 150 m) that, depending on subsurface thermal diffusivity and surface temperature perturbations, can reflect historical temperatures of at least 100 years ago (Huang et al., 2000) and up to tens of thousands of years (Ter Voorde et al., 2014).
Figure 11 illustrates that VAMPERS does a reasonable job of predicting shallow subsurface temperatures (Pearson correlation = 0.64).This result supports the notion that the preindustrial climate is well represented by iLOVECLIM.Points in Kazakhstan and Mongolia, and a few others in Russia, have a warm bias in the forcing (simulated is warmer than observed), which is probably due to an inaccurate representation of elevation temperature changes in iLOVECLIM, since many of those sites are at elevations above 1000 m.Even applying the lapse rate for a standard profile (6.5 • C km −1 ; McGuffie and Henderson-Sellers, 2005) would presumably make a significant difference on the depth since earlier sensitivity tests (Kitover et al., 2013) showed an average 55 m increase in equilibrium permafrost depth for every 1 • C colder.On the other hand, many of the other points show that predicted subsurface temperatures are on average a few degrees colder than the observed, leading to the most obvious conclusion that a cold bias exists in the iLOVECLIM climate.Although the cold bias, most obvious for Canada and Alaska, is congruent to the overestimation in permafrost thickness evident from the geographic breakdown illustrated in Fig. 10, it has not previously been substantiated in former analyses of LOVECLIM or iLOVECLIM so it is more likely that such a discrepancy is due to the air-ground coupling as opposed to simply the land surface temperature forcing.Indeed, there a number of other (sub)surface processes not included in the current ECBilt-VAMPERS coupling which may reduce the apparent cold bias.These effects alter the seasonal behavior of the thermal diffusivity in the subsurface and have been well-documented in observational studies (Williams and Burn, 1996;Woo and Xia, 1996;Fukui et al., 2008).Smith and Riseborough (2002) simplified these mechanisms into the surface offset (air to ground surface) and the thermal offset (ground surface to top of the permafrost).
Overall, the average range of error between observed and predicted is about 2.6 • C. Given that the comparisons are between point-based observations and large grid cell values, meant to represent a relatively large surface area, some variability is expected to occur.

Future development
The results of this paper demonstrate the ability of VAM-PERS forced by iLOVECLIM to model current permafrost distribution and thickness.The next step is to analyze the feedback that permafrost changes have on the climate.This has been of particular interest of the last decade since it is clear that specific feedbacks exists, most notably the release of locked up carbon in the atmosphere as permafrost degrades (Anisimov, 2007).The initial method behind a full coupling would be to activate the coupling mechanisms, shown in Fig. 3, and reanalyze the equilibrium results (since a full coupling would likely lead to an altered equilibrium permafrost state).In addition, the feedback effects would be most visible during millennial-scale transient climate shifts, when major permafrost degradation and/or disappearance is likely to occur.

Conclusions
The VAMPER model has been enhanced to allow simulations of estimated present-day permafrost thickness and distributions to be made using ECBilt land surface temperatures within the iLOVECLIM equilibrated preindustrial climate as the forcing.The VAMPER time step was reduced to 4 h to match the time step of ECBilt and allow seasonal effects, notably snow cover and the active layer, to be reflected in the simulation of permafrost.The predicted annual active layer from the stand-alone VAMPER model, under different temperature forcings, compares well with results from the Stefan equation.We also describe the snow option, which introduces the thermal insulation effects and changes in the thermal properties of snow over time due to varying snow densities.In addition, we developed and applied two new maps of geothermal heat flux and porosity.Incorporating these parameters at a global scale is an important step in improving the horizontal spatial variability of permafrost thickness/distribution while also maintaining the simplicity and efficiency of ECBilt-VAMPERS.
Equilibrium experiments for the PI climate show that when the snow component is included in the VAMPER model, the permafrost extent is noticeably reduced while the average surface offset of 2.7 • C is comparable to previous reports.We then compared permafrost thickness estimates and subsurface temperatures to corresponding observed values.Considering that we are comparing point measurements to grid-cell-based values, we consider the simulations reasonable.However, reasons for the discrepancies were discussed.One is that the relatively coarse horizontal ECBilt grid will never perfectly match the sensitivity of permafrost occurrence and depth due to local factors.This is also the case in the air-land temperature coupling, where some of the local effects will simply not be present in an EMIC.Similarly, when iLOVECLIM does not accurately represent the environmental lapse rate in areas of higher elevation, the occurrence of permafrost in these areas are overlooked by the VAMPERS model.Finally, some of the observed permafrost depths are not a function of the present (PI) climate, but rather a relict presence from previous cold periods.Therefore, when comparing measured to simulated results, some underestimations occurred.It is only with millennial-scale transient iLOVECLIM (with the ECBilt-VAMPERS coupling) model runs that we can realistically simulate, for example in areas of West Siberia, how permafrost evolved over periods of major climate change.
, a number of transient simulations Published by Copernicus Publications on behalf of the European Geosciences Union.D. C. Kitover et al.: Coupling of VAMPER within iLOVECLIM

Figure 1 .
Figure 1.(a) Plot comparing VAMPER model results using different time steps (annual vs. sub-daily) but the same annual average temperature forcing of −6 • C. (b) Plot showing the sr1 average, minimum and maximum temperature-depth profiles.Also shown in (b) is the ∼ 1 m active layer, marked as diagonal lines.

Figure 3 .
Figure 3. (a) Future iLOVECLIM coupling scheme between EC-Bilt and the VAMPER(S) model showing the variables (land surface temperature, snow water equivalent (SWE), and ground heat flux) passed between the components at each time step.(b) Land surface temperature of ECBilt and ground surface temperature of VAMPER(S).

Figure 4 .
Figure 4.The original geothermal heat flux map (top) from Davies (2013) and the weighted average version (top) for use as the lower boundary value in the iLOVECLIM experiments (bottom).

Figure 5 .
Figure 5. World maps showing (a) original map from Hartmann and Moosdorf (2012) (b) map of reclassified lithology using Table 2 and (c) the version geo-processed to match the ECBilt grid resolution.

Figure 6 .
Figure 6.An illustration of asynchronous coupling between VAM-PER(S) and ECBilt.The components are run semi-coupled for 100 years while VAMPER(S) is run the entire time.This allows VAMPER(S) to equilibrate with the climate state of iLOVECLIM using less computer resources time than a synchronous version.

Figure 8 .
Figure 8. Mean annual surface offset as a result of including the snow option.

Figure 9 .
Figure 9.A 1 : 1 scatter plot comparing simulated thickness results with corresponding permafrost thickness estimates from borehole data.Points 1-7 are outliers mentioned specifically above.

Figure 10 .
Figure 10.Map of deep Global Terrestrial Network for Permafrost (GTN-P) borehole locations with the simulated permafrost thickness (with snow enhancement) and observed permafrost extent(Brown et al., 2014).

Figure 12 .
Figure 12.Map showing locations of the MAGT measurements, collected for the International Polar Year (IPY) 2010 (GTN-P), used in the comparison to corresponding iLOVECLIM simulated subsurface temperatures.
who used the Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg (JSBACH) terrestrial ecosystem model.Other examples include Lawrence and Slater (2005), who used the Community Climate System Model (CCSM) to look at future permafrost extent and associated changes in freshwater discharge to the Arctic Ocean.Schaefer et al. (2011) used a land surface model (Simple Biosphere/Carnegie-Ames-Stanford Approach -SiBCASA) to simulate reduced future permafrost coverage and subsequent magnitude of the carbon feedback.Similarly, Schneider von Deimling et al. (2012) and Koven et al. (2011) also modeled future estimates of carbon emissions due to thawing permafrost.From a paleoclimate perspective, DeConto et al. (2012) used a version of the Global Environmental and Ecological Simulation of Interactive Systems (GENESIS) global climate model

Table 1 .
Variable values applied in the Stefan equation.

Table 2 .
Calculated maximum annual active layer thickness using both the Stefan equation and the VAMPER model under different forcing scenarios.