Effect of soil property uncertainties on permafrost thaw projections : A calibration-constrained analysis

Effect of soil property uncertainties on permafrost thaw projections: A calibration-constrained analysis Dylan R. Harp1, Adam L. Atchley1, Scott L. Painter2, Ethan T. Coon1, Cathy J. Wilson1, Vladimir E.Romanovsky3, and Joel C. Rowland1 1Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA 2Climate Change Science Institute, Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA 3Geophysical Institute, University of Alaska Fairbanks, USA Correspondence to: Dylan R. Harp (dharp@lanl.gov)


Introduction
Increasing Arctic air and permafrost temperatures (Serreze et al., 2000;Jones and Moberg, 2003;Hinzman et al., 2002;Romanovsky et al., 2007), the resulting increase in the thickness of soil that thaws on an annual basis (Romanovsky and Osterkamp, 1995), and the potential for greenhouse gas release due to the ensuing decomposition of previously frozen organic carbon (Koven et al., 2011;Schaefer et al., 2011) provide motivation for developing robust numerical projections of the thermal hydrological trajectory of Arctic tundra in a warming climate.Projections of permafrost thaw and the associated potential for greenhouse gas release from the accelerated decomposition of previously frozen carbon are subject to several sources of uncertainty, including (but not limited to) structural uncertainties in the climate models; uncertainty about the model forcings/inputs in the future (scenario uncertainty in the typology of Walker et al., 2003); parametric uncertainties in soil and surface properties that control the downward propagation of thaw fronts; and structural uncer-Published by Copernicus Publications on behalf of the European Geosciences Union.

D. R. Harp et al.: Effect of soil property uncertainties on permafrost thaw projections
tainties in the surface and subsurface thermal hydrological models.
Previous efforts to characterize uncertainty in permafrost thaw projections have mostly focused on climate model structural uncertainties and climate model uncertainties, presumably because of an implicit assumption that those two sources of uncertainty overwhelm the other sources.However, recent large-scale model comparisons suggest that a substantial portion of projected permafrost uncertainties is a result of structural model differences in land surface/subsurface schemes (Slater and Lawrence, 2013;Koven et al., 2013), particularly how subsurface thermal hydrologic processes are represented (Koven et al., 2013) rather than simply climate variation.Although those studies focused on structural uncertainty in surface and subsurface models and not on soil property uncertainty, the reported sensitivity to the subsurface model suggests that uncertainty in soil properties may also contribute significantly to overall uncertainty in thaw projections.
The bulk hydrothermal properties of soil that control the active layer thickness (ALT, i.e., the depth of soil that thaws on an annual basis) (Neumann, 1860;Stefan, 1891;Romanovsky and Osterkamp, 1997;Peters-Lidard et al., 1998;Kurylyk et al., 2014) vary among sites and locally within a single site, in particular being sensitive to the local organic matter content and bulk porosity (Letts et al., 2000;Price et al., 2008;O'Donnell et al., 2009;Hinzman et al., 1991;Chadburn et al., 2015a).Langer et al. (2013) identify the soil composition uncertainties, particularly the soil ice/water content, to have the largest effect on ALT.Intermediate to large-scale thermal simulations of ALT are known to be sensitive to soil properties (Hinzman et al., 1998;Rawlins et al., 2013).Because of this sensitivity, large-scale Earth system models (ESMs) were recently updated to include layers of moss and peat in order to better represent subsurface thermal conditions (Beringer et al., 2001;Lawrence and Slater, 2008;Wania et al., 2009;Subin et al., 2012;Ekici et al., 2014;Chadburn et al., 2015b).Despite the recognition of soil property uncertainty and heterogeneity as important contributors to uncertainties in permafrost conditions and extent, global and regional studies that address permafrost future conditions and extent typically apply broad soil texture classifications, such as those defined by Clapp and Hornberger (1978) and Cosby et al. (1984), to parameterize soil properties (Lawrence and Slater, 2008), usually without consideration of soil property uncertainty (Lawrence and Slater, 2005;Hinzman et al., 1998;Shiklomanov et al., 2007;Koven et al., 2013;Rinke et al., 2008).
Soil property uncertainty is different from many other sources of projection uncertainty (e.g., climate model uncertainty) in that uncertainties in soil properties may be reduced by a combination of site characterization (Hinzman et al., 1998) and model calibration (Romanovsky and Osterkamp, 1997;Nicolsky et al., 2009;Jiang et al., 2012;Atchley et al., 2015).Initial steps in that direction have been taken.For ex-ample, Romanovsky and Osterkamp (1997) calibrate thermal soil properties using a purely conductive thermal model using measured temperatures at several sites and Nicolsky et al. (2009) perform a sensitivity analysis of a calibration (data assimilation) approach to identify its ability to recover thermal soil properties using a 1-D thermal model and apply the calibration approach to several sites.Atchley et al. (2015) recently demonstrated an iterative approach for using site characterization data to simultaneously refine thermal hydrology model structure and estimate model parameters.Their approach was applied to the Barrow Environmental Observatory, but could be used at other sites to improve model structure and parameter assignments in the regional or global context.
Recognizing that permafrost projections are sensitive to subsurface model representations and that soil property uncertainties may be reduced through characterization and parameter estimation, a natural next step is to quantify how such activities will impact overall uncertainties in permafrost thaw projections in comparison to other sources of uncertainty.Here we address that question.Specifically, we consider how uncertainties in soil hydrothermal properties propagate to uncertainties in numerical projections of permafrost thaw at a well-characterized site.We go beyond a traditional unconstrained uncertainty quantification and focus on the residual uncertainties that remain after soil parameters have been carefully calibrated to borehole temperature data.The intent of the current work is to develop initial insights into how effective site characterization activities might be at reducing uncertainties associated with soil parameters.We show that with the future climate specified and with the advantage of calibration targets from a well-characterized site, significant uncertainties remain in projected ALT and other metrics important for carbon decomposition in the future climate.We evaluate both predictive uncertainty and interannual climate variability.We show that this residual uncertainty is significant, albeit less than that associated with uncertainties in the future climate.
We focus on temperature data as they are one of the easiest and most common types of soil data to collect at field sites and are often used for early site characterization.While many sites may have other types of measurements available, such as water and ice content measurements, many of these are more difficult to obtain at regular temporal intervals for extended periods of time.The incorporation of other types of data, such as water and ice content measurements, would be expected to reduce soil property uncertainty; however, this is not investigated here.
The arctic site in this investigation is the polygonal tundra within the Barrow Environmental Observatory (BEO).In particular, we focus on NGEE-Arctic site "area C", which contains degraded permafrost characterized by ∼50 cm deep troughs and shallow low centers.The polygonal tundra of the BEO is classified as a lowland, cold continuous permafrost system with a range of polygonal types and states, which includes intact low center polygons to degraded ice wedges and associated high center polygons.Much of the polygonal tundra contains an organic rich surface layer of peat overlaying a silty loam soil.Due to a low evaporative demand soils remain moist, despite relatively low annual precipitation, of which the bulk falls in the summer months (Liljedahl et al., 2011).The snowpack over the microtopography at the site is redistributed to a relatively level surface by strong winds, resulting in the deepest snowpack over troughs.Snow depth measurements collected around the site on 2 May 2013 were between 20 and 40 cm for centers, 10 and 20 cm for rims, and 40 and 60 cm for troughs, while the average snow density was 326 kg m −3 (Atchley et al., 2015).While our investigation focuses on the polygonal tundra within the BEO, other arctic landscape types are also prevalent (hillslopes, lakes, pingos).The importance of soil properties and the dominant influence of particular soil properties may change in landscapes other than polygonal tundra.

Model
We use the Arctic Terrestrial Simulator (ATS) to numerically solve the coupled groundwater flow, thermal, and surface energy balance equations.ATS is an integrated thermal hydrological code developed specifically for Arctic permafrost applications.It implements the modeling strategy outlined by Painter et al. (2013) using the multiphysics framework Arcos (Coon et al., 2016) to manage model complexity in process rich simulations such as these.Various components of ATS have already been described elsewhere; therefore, only a brief summary is provided here.
In the subsurface, the ATS solves nonlinear conservation equations for water and energy, using a three-phase (airwater-ice), single-component representation (Karra et al., 2014), which is a simplification of a more general twocomponent (water and representative gas phase) model (Painter, 2011).A recently developed constitutive model (Painter and Karra, 2014) is used to partition water between ice and liquid phases in unsaturated or saturated conditions.The partitioning model relates unfrozen water content below the nominal freezing point to the unfrozen soil moisture characteristic curve, thus avoiding empirical freezing curves.The model has been successfully compared to a variety of laboratory experiments on freezing soils (Painter and Karra, 2014;Karra et al., 2014;Painter, 2011).The material component model defines thermal conductivities and is described in detail in Appendix A of Atchley et al. (2015).Surface boundary conditions use a "fill and spill approximation", where we allow up to 4 cm of water to pond on the surface; all additional ponded water may run off the domain.The surface and subsurface thermal hydrology systems are coupled using continuity of pressure, mass flux, temperature, and energy flux, in a thermal extension of the coupling strategy presented in (Coon et al., 2015a).Additionally, we use a surface energy balance (Hinzman et al., 1998;Ling and Zhang, 2004;Atchley et al., 2015) in which surface latent and sensible heat, incoming and outgoing radiation, and conducted heat terms, along with incoming precipitation and outgoing evaporation are tracked.Finally, a dynamic, snow model is incorporated for tracking snow aging and consolidation, with resulting effects on albedo and melt (Atchley et al., 2015).As described in Sect.4.4 of Atchley et al. (2015), the snow model accounts for snow redistribution over the microtopography of the site and depth hoar formation.Additional details about the snow model are described in detail in Appendix B of Atchley et al. (2015).Not represented within this system are carbon cycle and vegetation processes, including long-term changes of peat composition, variability in peat thickness, and evolving microtopography due to degradation of ice wedges.
The subsurface domain is represented by a 2 cm layer of moss, followed by a 10 cm layer of peat, and an approximately 50 m mineral soil layer.The mesh is discretized in an increasing fashion from 1 cm at the surface to 2 m at the bottom (∼50 m).We performed a mesh discretization analysis, presented in Fig. S1 in the Supplement, to determine that the discretization was adequate.The required climate forcings for the ATS models are precipitation of snow and rain, air temperature, wind speed, relative humidity, and incoming short and longwave radiation.The lower boundary is set to a constant temperature of −9.7 • C.

Previous calibration from Atchley et al. (2015)
The uncertainty quantification is performed around a previous calibration by Atchley et al. (2015).They used 1-D column models representing the dominant microtopographical features (center, rim, and trough of polygonal ground) to calibrate hydro-thermal soil parameters using soil temperatures at the BEO measured by the Next-Generation Ecosystem Experiments Arctic (NGEE-Arctic) team during calendar year 2013.Initial conditions for the models were generated by completely freezing the fully saturated model from below and then allowing the initial conditions to emerge over www.the-cryosphere.net/10/341/2016/The Cryosphere, 10, 341-358, 2016 a 10-year spin-up simulation using daily air temperatures averaged from 10 years of data as the top boundary condition.This process allowed a shallow vadose zone to develop consistent with field observations.The calibration considered temperatures measured at 9 depths from 10 to 150 cm.The calibration was performed in a coupled fashion where each "model run" of the calibration consisted of simulating center, rim, and trough column models with the same soil parameter values for peat and mineral soil.This coupled calibration identifies soil parameters that provide a generalized fit, compromising in a least squares sense to match the data from all three models.An implicit assumption of the coupled calibration is that the soil properties may be adequately represented as independent of the microtopography.Atchley et al. (2015) first calibrated subsurface properties using 2 cm deep temperatures measured in 2013 as Dirichlet boundary conditions and temperatures measured at the considered depths as calibration targets.Then Atchley et al. (2015) performed an additional surface/subsurface calibration to verify that the surface energy balance model is capable of producing surface temperatures consistent with measurements.The coupled surface/subsurface model allows the use of future climate models as model forcings to drive hydro-thermal permafrost projections.
The calibration data period is limited to calendar year 2013 since at the time of calibration, this was the only full year of high-resolution borehole temperatures available at the site (Atchley et al., 2015).Subsequently, year 2014 data has become available.To verify that the calibration has extracted the hydrothermal properties of the system independent of the climatic conditions during the calibration, we evaluated the ability of the calibrated parameters to produce forward simulations that are consistent with 2014 data.This evaluation is presented in the results section.

Soil property uncertainty quantification
We generated an ensemble of 1,153 calibration-constrained parameter combinations by the null-space Monte Carlo (NSMC) method (Doherty, 2004).The NSMC approach samples from insensitive regions of the parameter space (i.e., the null space) determined by an eigenanalysis of parameter sensitivities calculated at the calibration point.Based on analysis of ensemble forward simulations of the calibration year (2013) and a convergence analysis of the 95 % confidence band of simulated temperatures, we consider all parameter combinations in the ensemble calibrated and equally consistent with measured temperatures.

Permafrost projections through 2100
In order to make projections of hydro-thermal permafrost conditions, we use the surface/subsurface model described in Sect.2.1.We use the Community Earth System Model (CESM) (Gent et al., 2011) driven by the Representative Concentration Pathway 8.5 (RCP8.5)greenhouse gas concentration trajectory (Moss et al., 2008) from year 2006 to 2100 as atmospheric forcings for the surface energy balance of the model.The CESM output was adjusted to fit current climate at the BEO.In this way, we hold the climate scenario constant to isolate the effect of soil property uncertainty.RCP8.5 corresponds to a business as usual warming scenario with 8.5 Wm −2 forcing by 2100.
The projection simulations took on the order of several hours (∼2-4 h) on a Linux cluster with 3.2 GHz processors.We used the Model Analysis ToolKit (MATK) Python module (http://matk.lanl.gov) to facilitate the concurrent execution of the ensemble of ATS models on Los Alamos National Laboratory high performance computing clusters.

Permafrost metrics
Predictive uncertainty of projections is determined by comparison of permafrost metrics at year 2006 and for the last decade of the projections (2091 through 2100).The metrics include (1) ALT, (2) annual thaw depth-duration (D), (3) annual mean liquid saturation (S l ), and (4) a modified Stefan number (S T ) and are described below.

Active layer thickness (ALT)
In general, ALT is defined as "the layer of ground subject to annual thawing and freezing in areas underlain by permafrost" (http://www.uspermafrost.org/glossary.php).Permafrost has also been defined as the region of the subsurface that remains at or below 0 • C for 2 or more years.The ALT defined that way would be the minimum of the maximum annual thaw depth over each 2-year moving window.We use a less arbitrary definition for the ALT here as the annual maximum thaw depth in accord with the general definition and similar to Koven et al. (2011).Given the discrete nature of our mesh, and the nonlinear nature of vertical soil temperature profiles near 0 • C, we determine ALT as the bottom of the deepest thawed mesh cell (temperature above 0 • C) for the year.

Annual thaw depth-duration (D)
ALT controls the amount of organic carbon experiencing thaw and thus microbially induced decomposition during a year.Because ALT is defined as the maximum thaw depth, it does not include information on duration of thaw.To quantify increasing duration of thaw in the future climate (i.e., the effects of earlier thaw and later freeze-up) as well as increasing depth, a new metric is introduced here: the annual thaw depth-duration D, defined as where H is the Heaviside function, 1 if T (z, t) is above 0 • C, 0 otherwise, z is depth in meters, and t is time in days.The fraction on the right side of Eq. ( 1) normalizes the metric by the 365 days in a year.We express D with units of m 3 m −2 to indicate that this metric defines the volume of thawed soil per unit area.D is a rough proxy for the potential for soil organic matter decomposition.It merges the amount of unfrozen soil and duration that soil is above freezing temperature for a given year.Therefore, the metric does not account for biological activity that occurs below 0 • C, which is generally considered to be greatly reduced (Mikan et al., 2002;Davidson and Janssens, 2006), but has been observed in permafrost soils (Sachs et al., 2008).It is noted that, while the annual amount of decomposition is likely correlated with D, the two quantities are not directly proportional because soil temperature and moisture will also change and affect the decomposition rates in the future climates.Nevertheless, uncertainty in D is of interest as it is indicative of uncertainty in future decomposition rates.

Annual mean liquid saturation (S l )
The annual mean liquid saturation S l is defined as where S l (z, t) is the liquid saturation as a function of depth and time.S l quantifies the spatially and temporally averaged liquid saturation in the unfrozen soil for a given year.Note that the denominator in Eq. ( 2) is the annual thaw depthduration metric D from above, except without dividing by 365.Liquid saturation within the active layer is of interest because of its control on decomposition rates, coupling hydrology to biogeochemical fluxes.

Stefan number (S T )
We propose an extension of the Stefan number from the form in Kurylyk et al. (2014) to one that incorporates intra-annual temporal changes and stratified soil properties.The Stefan number is the ratio of subsurface sensible to latent heat.In the current context, this refers to the amount of subsurface heat exchange that results in a change in temperature vs. the amount that is consumed in the isothermal conversion of ice to liquid water.The Stefan number provides information about the form of subsurface energy utilization in permafrost regions and is fundamental to a basic understanding of permafrost thaw mechanisms.
In its most basic form, the Stefan number is defined as where c b is the bulk specific heat of the material and L f is the latent heat of fusion of water (334 000 J kg −1 ).Kurylyk et al. (2014) define the Stefan number for the permafrost problem as where ρ b is the density of the thawed zone, T s is the surface temperature, T f is the temperature of freezing or thawing soil (taken as 0 • C), S wf is the liquid saturation in the thawed zone that was frozen, ρ w is the density of liquid water, and φ is porosity.Kurylyk et al. (2014) use this definition to evaluate the thermal regime of analytical solutions of soil thaw.
We expand this definition here to include the increased detail available in our numerical simulations as where S ice is ice saturation.The integrations are performed over the entire year (i.e., from 1 January through 31 December).Equation ( 5) expands on Eq. ( 4) to allow the consideration of details of transient heating and cooling throughout the year and stratified hydrothermal soil properties within the soil profile.

Comparison to climate uncertainty
To provide a reference point for the effect and magnitude of soil property uncertainty, we also perform ATS projections forcing the energy balance model with atmospheric projections from CESM, INM-CM4 (INM) (Volodin et al., 2010), BCC-CSM1-1 (BCC) (Ji, 1995), MIROC (Watanabe et al., 2010), CanESM2 (CAN) (Verseghy, 1991), and HadGEM2-CC (HAD) (Jones et al., 2011;The HadGEM2 Development Team, 2011;Collins et al., 2011) climate models based on RCP8.5 using the calibrated (fixed) soil parameters from Atchley et al. (2015).Using the calibrated soil parameters in these simulations isolates the effect of structural climate uncertainty.We compare permafrost projection uncertainty due to the NSMC ensemble of soil parameters (hydrothermal soil property uncertainty) and to the variability between climate models (structural climate uncertainty).
The soil property uncertainty in this analysis is parametric and can be considered more aleatoric/probabilistic in nature.The climate model uncertainty is epistemic in nature due to a lack of knowledge regarding modeling of atmospheric phenomena.These distinctions do limit comparisons that can be drawn between these two uncertainties.However, the comparison is relevant for our purposes to provide a frame of reference for soil property uncertainty to one of the other current, primary sources of permafrost thaw uncertainty.

Ensemble of calibration-constrained soil parameter combinations
In order to determine the effect that calibration-constrained soil property uncertainty can have on long term projections of permafrost conditions, we performed an uncertainty quantification around the calibrated soil parameters of Atchley et al.
www.the-cryosphere.net/10/341/2016/The Cryosphere, 10, 341-358, 2016 (2015).The strategy involved identifying a representative set of parameter combinations that all produce simulated temperatures that are consistent with observed temperatures.We use null-space Monte Carlo (NSMC) (Tonkin and Doherty, 2009), a form of calibration-constrained Monte Carlo, to accomplish this goal.NSMC was selected based on its sampling economy given the computational burden of the simulations involved.
A subset of the 16 soil parameters from the calibration of Atchley et al. (2015) are included here and presented in Table 1.The top pressures of the center and trough profiles from the calibration are not included here as these are internally calculated in the surface/subsurface ATS model.The van Genuchten water retention parameters are not included either as they were found to significantly exceed their physical boundaries during NSMC sampling.This is an indication that these are highly insensitive parameters and do not significantly effect simulated temperatures.This may be explained by the fact that these parameters control the shape of the water retention curve, but that this influences thermal properties of the soil only for a limited time near freeze up or thaw.
This leaves the 10 soil parameters listed in Table 1.The parameters r,peat and r,min are van Genuchten soil moisture characteristic residual saturations (Van Genuchten, 1980).K peat and K min are material thermal conductivities for peat organic matter and mineral grains within the soil layers.Bulk thermal conductivities are a function of material thermal conductivities and are sensitive to ice, liquid, and gas saturation, which is calculated within ATS as described in Appendix A of Atchley et al. (2015).A peat,fr , A peat,un , A min,fr , and A min,un are empirical exponents describing the dependence of frozen (fr) and unfrozen (un) Kersten numbers (i.e., ratios of partially to fully saturated thermal conductivities) to ice and liquid saturation states, respectively (Painter, 2011).The minimum and maximum parameter boundaries are modified from the calibration for the NSMC sampling (the parameter ranges are reduced in most cases) to physical limits identified through literature review and field observations from the BEO (Imnavait Creek and Kuparuk River, Alaska, Hinzman et al., 1991Hinzman et al., , 1998; large-scale pan-arctic modeling efforts, Beringer et al., 2001;Lawrence and Slater, 2008;Capricorn Fen, Northern Quebec, Letts et al., 2000; Gailbraith Lake, Northern Alaska, Overduin et al., 2006;Bonanza Creek, Delta Junction, and Washington Creek, Interior Alaska, O'Donnell et al., 2009;Siksik Creek, Northwest Territories, Quinton et al., 2000;Franklin Bluffs, West Dock, Imnavait Creek, Northern Alaska, Nicolsky et al., 2009;Fort Simpson, Scotty Creek, Northwest Territories and Wolf Creek, Yukon Territory, Zhang et al., 2010;Samoytov Island, Lena River delta, Siberia, Chadburn et al., 2015b).
To a lesser degree, other parameters were also found to exceed their physical boundaries during NSMC sampling.Therefore, we used the intersection of the null space and parameter boundaries as our criterion to accept samples.The generation of 20 000 samples within the null space resulted 0.7 0.8 0.9 Figure 1 presents histograms, while Fig. 2 presents paired plots of the NSMC ensemble soil parameters.In the matrix of plots in Fig. 2, parameter histograms are plotted along the diagonal (also presented in greater detail in Fig. 1), paired scatterplots in the lower triangle, and Pearson correlation coefficients are presented in the upper triangle.In Fig. 1, it is apparent that K peat followed by A peat,un are the most constrained parameter by the NSMC analysis.The rest of the parameters span significant portions of their range.This indicates that there are many combinations of parameters that result in calibrated temperatures.Many of the histograms are seen to butt up against their boundaries, indicating that these are parameters where the extent of the null space exceeds their range.
Applying NSMC to multiple calibration locations is often suggested (Tonkin and Doherty, 2009).In the calibration performed by Atchley et al. (2015), multiple local minima were identified.However, based on the broad range of parameter combinations with limited correlations and the fact that most parameters span most of their range, we conclude that the NSMC analysis from this single calibration point sufficiently captures the soil property uncertainty.
The correlations imposed by the NSMC sampling are evident by inspecting the Pearson correlation coefficients and scatterplots in Fig. 2. The strong correlations that are present are a result of a balancing act between parameters to achieve a least squares fit to measured temperatures.For example, the relatively strong negative correlation between K peat and K min (correlation of −0.81) is due to the fact that deeper temperatures in the soil profiles are controlled by the effective thermal conductivity.Therefore, there are numerous (negatively correlated) combinations of K peat and K min that produce similar effective thermal conductivities, resulting in good matches to measured temperatures.Many other correlated parameter pairs are also apparent, most with significantly lower correlations.There are also many uncorrelated parameter pairs (e.g., φ peat and K peat ) indicating a complete lack of interaction between the parameter pairs.The following analysis of permafrost projection uncertainty is conditional on the NSMC correlations presented here, and any conclusions take these correlations into account.References to Fig. 2 are made in the following sections explaining some of the impacts of these correlations.
The range in RMSE values is from around 0.55 to 0.65 • C. The accuracy of the temperature probes are ±0.1 • C. Therefore, the percentage of the RMSE that may be attributable to measurement imprecision is around 15-18 %.
Figure 3 presents the evaluation of the calibration against 2014 data and the 95 % confidence band of temperatures for the NSMC ensemble.The evaluation is presented as time series of temperatures where the fit between 2013 measured and calibrated temperatures can be compared to the 2014 measured and simulated temperatures.Since the 2014 measured temperatures are not included in the calibration, this comparison serves as an evaluation of the 2013 calibration.By inspection of the plots, it is apparent that the fit during the evaluation period is similar to the match during the cali-bration period (1st, 3rd, and 5th plots for the center, rim, and trough, respectively).This provides an initial indication that the calibration has extracted the hydrothermal relationships from the system and can be applied to years with different climate conditions than the calibration period.
The other plots in Figure 3 contain the corresponding 95 % confidence bands for 2013 temperatures.We performed a convergence analysis of the ensemble by calculating the ratio of measurements included in the 95 % confidence band as the number of ensemble members increased.We found that the ratio stabilized after the ensemble reached more than around 800 members.This indicates that the ensemble has converged and that more samples are not necessary.A plot of the convergence analysis is provided in the Supplement to this article, Fig. S2.
The measured temperatures are within the 95 % confidence band 79 % of the time for the center, 59 % for the rim, 46 % for the trough, and 61 % overall.The primary causes of these discrepancies are due to difficulties in capturing trends during the freeze up and thaw of the active layer.The low values are primarily due to the 95 % confidence band missing measured values at deep measurements apparent in Figs.S3, S4, and S5 in the Supplement.A lack of overlap is apparent during thawing (around day of year 150) and freeze up (around day of year 320), and is particularly evident in the rim profile in Fig. 3.These discrepancies are reduced in the decoupled calibrations (calibrations on individual profiles) (Atchley et al., 2015).We choose to use the coupled calibration parameters in order to obtain soil property values that provide a generalized characterization of the soil properties across the microtopography at the site.The expense of such a generalized characterization is a compromised fit across profiles.The discrepancies are also less pronounced in the center profile, which is the model used for the forward projections.Many physical processes may be leading to this result that become more pronounced in the coupled calibrations as parameter values are given less freedom to mask missing physical processes.For one, the exposed sides of the rim and subsequent lateral heat flow are not explicitly modeled and may at least partially explain the discrepancy.During the thaw, a lack of advective transport of heat by liquid water through the pore space created by sublimation during the winter (not included in the model) may result in warmer measured temperatures (Kane et al., 2001).
An initial ensemble created using Latin Hypercube Sampling with 1,000 samples postprocessed to include parameter combinations with RMSE's below various thresholds indicated that to achieve a convergent ensemble using Latin Hypercube Sampling would be computationally prohibitive.An additional NSMC analysis was performed with a more restrictive null space (only 2 eigenvectors out of 10 included in the null space).This ensemble did not require postprocessing based on RMSE, since all the RMSE values were deemed sufficiently small.This analysis resulted in over-correlated parameters.We therefore chose a loosely constrained NSMC (5 out of 10 eigenvectors included in null space) excluding samples with RMSE greater than 0.65 • C. We considered other RMSE cutoffs, but selected 0. sured temperatures within the 95 % confidence band, we consider all the unmodified NSMC samples to be calibrated and do not apply this step.These observations also led to the assumption that all NSMC samples are equally consistent with measured temperatures as opposed to using a weighting scheme.A datafile of the ensemble can be downloaded from http://dx.doi.org/10.5440/1236647.Boxplots of ALT are shown in Fig. 4a.The median ALT increased from approximately 30 cm in 2006 to nearly 0.9 m by the end of the century.The predictive uncertainty in ALT also increases significantly from the beginning to later years of the projections.The inter-annual variability of ALT projections is dependent on climate, as warmer years (e.g., 2094) have greater ALT and larger uncertainty than cooler years.This is apparent in Fig. 5 where the ensemble thaw depth statistics (median and 95 % confidence band), CESM8.5 air temperature, and ensemble snow depth statistics (95 % confidence band) time series are plotted together for comparison.

Permafrost thaw projection uncertainty
Boxplots of annual thaw depth-duration (D) are presented in Fig. 4b.The predictive uncertainty in D during the last decade of the projections is significantly greater than for the first year (2006).As expected, the inter-annual trends in D and ALT are similar.Also, the uncertainty of D is relatively larger during warmer years than cooler years, similar to ALT.
Boxplots of the annual mean liquid saturation (S l ) are presented in Fig. 4c.The predictive uncertainty in S l actually decreases slightly from the first year to the last decade.Also, in general, the last decade is slightly wetter than 2006, but only marginally so.Therefore, this hydrothermal analysis does not indicate that the soil moisture regime will change significantly as permafrost thaws.Soil moisture is one of the factors controlling the complex process of partitioning of carbon decomposition between CO 2 and CH 4 .However, other factors affecting carbon decomposition not considered here could affect the partitioning of carbon decomposition end products.
Boxplots of the Stefan number (S T ) are presented in Fig. 4d.In 2006 the soil profiles for the majority of the ensemble are latent-heat dominated.However, some Stefan numbers are greater than 1, with values ranging from around 0.3 to 1.4 (from around 3 times the latent heat as sensible heat to 1.4 times the sensible heat as latent heat).However, by the last decade, nearly all Stefan numbers are 0.2 or less (at least 5 times as much, and up to 20 times as much latent heat as sensible heat).This indicates a fundamental change in the way that the active layer processes energy between the beginning and later years of the projections.The thermal regime of the active layer becomes significantly more dominated by latent heat during the projections.The amount of energy that is utilized in creating a temperature gradient in the soil profile becomes proportionately smaller compared to the amount of energy consumed in the isothermal melting of ice.This is at least partially due to the approximately 3-times increase in the quantity of ice that is melted during later years of the projections.Predictive uncertainty appears to decrease from 2006 compared to the last decade, but this is likely due to the Stefan number approaching its lower limit.
To further illustrate predictive uncertainty of the ALT projections, temperature profiles at the time of ALT for year 2100 are presented in Fig. 6.Summary statistics (median and 5th and 95th percentiles) for 2006 are presented for reference.The discrete surface temperatures categorized by day of year (colors) reflect the fact that the surface temperature is highly dependent on the climate/air temperature for a given year, which is the same for all projections.Similarly, the day of ALT for 2006 do not all occur on the same day across realizations, occurring from day of the year 246 to 260.The increase in median ALT from around 30 cm to around 0.9 m from 2006 to 2100 is also apparent in this figure.The difference in the temperature regime within the profile is apparent in these figures as well by the curvature near the surface in most of the profiles in 2100 compared to 2006.This indicates that as the climate warms and the day of year when ALT occurs becomes later in the year, the surface temperature at that time will be cooler.This increase in lag time from the surface temperature to the active layer base is a result of the thermal wave traveling a greater distance to reach the permafrost.This may also be due to relative changes in the temperature gradient within the active layer and the permafrost as the ALT increases leading to delayed freeze from below.
Figure 7 shows similar plots to Fig. 6, but in this case, statistical measures of the ensemble are plotted.Statistical representation of the temperature profiles in Fig. 6 are plotted in Fig. 7a, along with bulk thermal conductivity (Fig. 7b) and ice (Fig. 7c), liquid (Fig. 7d), and gas (Fig. 7e) saturation profiles when ALT occurs in 2006 and 2100.The vari-   ation in thermal conductivity and saturation states further illustrates the predictive uncertainty due solely to soil properties.Substantial shifts in predictive uncertainty are also apparent from 2006 to 2100.In Fig. 7b, it is apparent that the thermal conductivity in the soil profile decreases from 2006 to 2100 due to the loss of the more thermally conductive ice from the profile, thereby inhibiting the propagation of the thermal wave.The deepening of the permafrost table is apparent in Fig. 7c as a deepening of the ice saturated region.Note that liquid saturations for mineral soil remain at its residual values below 0 • C and that residual liquid saturations ( r,peat and r,min ) are variable parameters within the uncertainty quantification (refer to Table 1).As a result, the ice saturation within the permafrost region is variable within the ensemble.In Fig. 7d and e, it is apparent that the liquid and gas saturations both increase as ice is converted to liquid and void space becomes available with the deepening of the permafrost table.This results in a potentially continuous gas phase, ranging up to at least 80 cm deep across the ensemble at the time of ALT, indicating the potential for aerobic conditions at these depths.Higher liquid saturations may result in lateral flow, a phenomenon not considered in our models.Given the polygonal micro-topography of the site, lateral flow may be less important than in hilly terrain.However, lateral flow may be important for the polygonal centers and rims.

Comparison to climate model structural uncertainty
In this section, we provide a frame of reference to the effect of soil property uncertainty on permafrost thaw projections by comparison to the uncertainty currently present in climate models.Without such a comparison, the relative contribution of soil property uncertainty would be difficult to gauge.Figure 8 presents histograms of projection metrics collected   This is expected since the calibration-constrained ensemble is forced using the CESM model.Similarly, the supports of calibration-constrained ensemble histograms for other climate models would be expected to encompass the calibrated soil parameter values (circles in Fig. 8) as well.This indicates that different climate models will result in different magnitudes of projection uncertainty due to soil property uncertainty.For example, if the calibration-constrained ensemble was simulated using MIROC, the magnitude of the projection uncertainty of D (Fig. 8b) could be as much as 4-5 times larger than for CESM.This indicates the interactive effect that soil property and structural climate model uncertainties have on projection uncertainty and that these forms of uncertainty are not easily decoupled.These plots present both the magnitude of projection uncertainty due to soil property uncertainty based on CESM atmospheric projections (histograms) and to structural climate model uncertainty (circles).By comparing the ensemble 95 % confidence bands for the metrics to the range of values across the climate models, it is apparent that structural climate model uncertainty has a greater impact on projection uncertainty than soil property uncertainty.The ratios of the ensemble 95 % confidence bandwidth and the range between the minimum and maximum values for climate models are 26 % for ALT, 9 % for D, 45 % for S l , and 80 % for S T .As explained above, if a different climate model had been used for the ensemble calculations, these percentages would be different.

Dependence of permafrost projections on soil parameters
Based on a correlation analysis, all the permafrost metrics are positively correlated, with lower correlations between annual mean liquid saturation and the other metrics.A paired plot of the permafrost metrics is provided in the Supplement to this article for additional detail (Fig. S6 in the Supplement).The correlation between ALT and D is expected given the definition of D as a metric defining the quantity and duration of unfrozen soil.The correlation of S l to ALT is a result of the deeper portions of the thicker ALT scenarios having slightly increased levels of saturation, which is apparent in the liquid saturation statistical profiles in Fig. 7d for year 2100.The correlation between D and S l can be explained by a similar argument.Increased levels of saturation lead to higher bulk thermal conductivity of the mineral soil layer, resulting in thicker ALT and larger D due to increased energy flux.
Correlations between S T and the other projection metrics indicate that as ALT increases, resulting in increased annual thaw depth-duration D and annual mean liquid saturation S l , the system becomes increasingly latent-heat dominated.This is due to the fact that more energy is required to thaw greater depths of frozen soil in later years.Figures 9,10,11,and 12 explore correlations between the calibration-constrained parameters and projected metrics.These figures contain scatterplots between hydro-thermal soil parameters and projection metrics for year 2100.The discrete nature of the samples with respect to ALT mentioned above due to the mesh discretization is also apparent in Fig. 9. Pearson correlation coefficients for each soil parameter/projection metric pair are presented on each scatterplot.The points are colored by D in Fig. 9 and by ALT in Figs. 10,11,and 12 to illustrate the correlations between metrics (see also Fig. S6 in the Supplement).Peat parameters are presented along the left column and mineral soil parameters along the right column of each figure.
Four strong correlations are apparent in Figs. 9, 10, 11, and 12 with coefficients greater than 0.9.Many of these correlations confirm our qualitative understanding of the model.It is apparent that in many cases, projection metrics have stronger dependencies on the mineral soil porosity (φ min ) and residual saturation ( r,min ) parameters compared to the corresponding peat parameters (φ peat and r,peat ).Dependence on the other parameters is less predictable.For example, decreasing mineral soil porosity (φ min ) increases the bulk thermal conductivity of the mineral soil due to the relatively large thermal conductivity of the mineral soil grains, leading to larger ALT (top right plot in Fig. 9).
We determine linear dependency coefficients of projection metrics to calibration-constrained parameters using ordinary least squares.We limit the analysis to soil parameter/projection metrics exhibiting moderate to strong correlation (|ρ| > 0.7).Table 2 presents the intercept and slope coefficients from the analysis, along with their 95 % confidence intervals.All coefficients in Table 2 are significant at the 1 % level.The coefficient of determination (R 2 ) is presented indicating the portion of the variance explained by the regression for each case.Note that since we use ordinary least squares including an intercept, the R 2 is simply the square of the correlation coefficients (ρ) presented in Figs. 9, 10, 11, and 12. Calibration-constrained parameters not included in Table 2 resulted in regressions with R 2 less than 0.5.
The slope coefficients are emphasized in bold in the table since these describe the first-order dependence of projection metrics on the calibration-constrained parameters.The slope coefficients describe the change in ALT given a unit change in the calibration-constrained parameter.For example, if φ min increases by 0.1, we would estimate that ALT will decrease by around 0.14 m.These coefficients can be useful in gaging the impact of soil parameter changes on projection metrics.

Discussion and conclusions
In summary, we extended previous calibration and model refinement work (Atchley et al., 2015) to quantify postcalibration uncertainty in soil properties and the impact of uncertainty on projections of permafrost thaw.Using a model with parameters calibrated against data from the BEO, driving the NSMC ensemble of models using the CESM climate model in the RCP8.5 scenario, and comparing against a set of other climate models in the RCP8.5 scenario, the following conclusions can be made.
-The median ALT and annual thaw depth-duration (D) of the calibration-constrained ensemble increase by around a factor of 3 by the end of the century.
-The effect of soil property uncertainty based on CESM atmospheric forcings is approximately 26 % of the uncertainty caused by climate model structural uncertainty for ALT, 9 % for D, 45 % for S l , and 80 % for Stefan number.
-Predictive uncertainty of ALT and D due to soil property uncertainty increase significantly from the first year to the last decade of the projections.
-Predictive uncertainty of soil moisture content due to soil property uncertainty is not significantly changed by the end of the century.
-Predictive uncertainty of the Stefan number due to soil property uncertainty decreases, but this is at least partially due to this metric approaching its lower boundary in the last decade.
-The manner in which the active layer processes incoming energy changes significantly.The active layer moves to an increasingly latent-heat-dominated system due to larger quantities of frozen ground thawed each year.
-ALT, D, and S T are highly dependent on φ min , while S l is highly dependent on r,min and moderately dependent on r,peat .
Efforts to quantify the relative roles of soil property vs. climate model uncertainty have only recently begun.We found that the effect of soil property uncertainties can be reduced to levels lower than the uncertainty generated by uncertainties in climate model structure through a process of calibration to field observations, model structural refinement (Atchley et al., 2015), and calibration-constrained uncertainty analysis.However, we had the advantage of high-resolution data from an unusually well-characterized site, which suggests that the residual uncertainty identified here using temperature data only is close to a practical limit.The quantitative results shown here are specific to the site, available data, RCP trajectory assumption, and climate model.Nevertheless, the approach presented here is anticipated to be useful for understanding the impact that additional data collection might have on reducing uncertainty associated with other high-latitude permafrost sites.Potential directions for future work include the investigation on the impact that longer data streams and other types of observation might have on reducing uncertainties.In particular, the calibration against borehole temperature data was uninformative concerning certain water retention properties of the soils (van Genuchten α and m parameters).Therefore, co-located measurements of soil moisture would be useful to help constrain those parameters, and may reduce the uncertainty associated with the other soil properties as well.Moreover, given the known spatial variability in soil properties across the pan-Arctic (Hinzman et al., 1998;Rawlins et al., 2013), calibration-constrained soil property uncertainty across larger spatial scales warrants further investigations.
The Supplement related to this article is available online at doi:10.5194/tc-10-341-2016-supplement.

Figure 1 .
Figure 1.Histograms of calibration-constrained hydrothermal soil parameter combinations obtained by NSMC sampling

Figure 2 .
Figure 2. Matrix of paired plots of calibration-constrained hydrothermal soil parameter combinations obtained by NSMC sampling.Parameter histograms are plotted along the diagonal, paired scatterplots in the lower triangle (2-D projections of the null space), and Pearson (linear) correlation coefficients in the upper triangle.The histogram counts for all histograms are indicated along the ordinate axis of the upper left plot.

Figure 3 .
Figure 3. Time-series of temperature at 40 cm depths plotted as a function of the day of the year for the polygonal center, rim, and trough profiles.Alternating plots include measured values from the BEO for 2013 (red line) and 2014 (grey line) and simulated temperatures from the 2013 calibration (blue line) and 2014 evaluation (black line).Every other plot contains the 2013 95 % confidence band for the NSMC ensemble as a shaded light blue region.

Figure 4
Figure 4 present boxplots of permafrost metrics for the first year (2006) and the last decade (2091-2100) of the projections.Individual boxplots for each year present the predictive uncertainty (due to parametric soil property uncertainty), while comparisons between boxplots for each metric indicate the inter-annual climate variability of the projections for the specified climate model.We present the year 2006 as an indication of the initial parametric uncertainty.Boxplots of ALT are shown in Fig. 4a.The median ALT increased from approximately 30 cm in 2006 to nearly 0.9 m by the end of the century.The predictive uncertainty in ALT also increases significantly from the beginning to later years

Figure 4 .
Figure 4. Boxplots of projected metrics including (a) ALT, (b) annual thaw depth-duration, (c) annual mean liquid saturation, and (d) Stefan number for year 2006 and from 2091 to 2100.The bottom and top of the boxes are the first and third quartiles, the red lines are medians, the whisker lengths are 1.5 times the interquartile range (50 %), and the plus symbols are outliers.

Figure 5 .
Figure 5. Thaw depth, air temperature, and snow depth time series for years 2006 and 2091 through 2100.The black line in the top plot is the median thaw depth of the ensemble and the blue shaded region is the 95 % thaw depth confidence band for the ensemble.The black region in the bottom plot is the 95 % snow depth confidence band for the ensemble.

Figure 6 .
Figure 6.Predictive uncertainty due to soil properties for depth profiles of temperature for the ensemble when ALT occurs for calendar year 2100.Color indicates the day of the year when ALT occurred for each realization.The 2006 median and 5th and 95th percentiles for the ensemble are plotted for reference.Day of year when ALT occurs for realizations in 2006 is from 246 to 260 (not indicated in the plot).

Figure 7 .
Figure 7. Predictive uncertainty due to soil property uncertainty for depth profiles of ensemble statistical quantities when ALT occurs for calendar years 2006 and 2100.The shaded regions are the 95 % confidence intervals for 2006 (red) and 2100 (blue).

Figure 8 .
Figure 8.Comparison of (a) ALT, (b) annual thaw depth-duration, (c) annual mean liquid saturation, and (d) Stefan number projection uncertainty due to soil property uncertainty (histograms) and structural climate model uncertainty (circles).Histograms include calibration-constrained ensemble values for years 2091 to 2100 (11 530 values) based on the CESM8.5 climate model.Open circles below the histograms are values for the various climate models for the same years using the calibrated soil parameters (10 values each, except for BCC which has 9).NSMC ensemble 95 % confidence band (CB) limits are indicated as vertical dashed lines.

Figure 9 .
Figure 9. Scatterplots between calibration-constrained parameters and projected ALT for year 2100.Soil parameters associated with peat are on the left and with mineral soil on the right (refer to column headings).Colors represent annual thaw depth-duration.The associated Pearson correlation coefficient ρ is indicated in each plot.The discrete nature of the ALT is due to the computational mesh discretization.

Figure 10 .
Figure 10.Scatterplots between calibration-constrained parameters and projected annual thaw depth-duration.Soil parameters associated with peat are on the left and with mineral soil on the right (refer to column headings).Colors represent ALT.The associated Pearson correlation coefficient ρ is indicated in each plot.

Figure 11 .
Figure 11.Scatterplots between calibration-constrained parameters and projected annual mean saturation.Soil parameters associated with peat are on the left and with mineral soil on the right (refer to column headings).Colors represent ALT.The associated Pearson correlation coefficient ρ is indicated in each plot.

Figure 12 .
Figure 12.Scatterplots between calibration-constrained parameters and projected Stefan number.Soil parameters associated with peat are on the left and with mineral soil on the right (refer to column headings).Colors represent ALT.The associated Pearson correlation coefficient ρ is indicated in each plot.

Table 1 .
NSMC parameter minimum and maximum bounds, units, and descriptions

Table 2 .
Linear regression intercept and slope coefficients for permafrost metrics as a function of calibration-constrained parameters