2.1. Soil Model

[6] Figure 1 depicts a schematic representation of the updated PWBM described in this paper. The updated PWBM soil model discretizes a 60 meter soil column which contains 23 layers, with layer thickness increasing with depth. The model simulates snow/ground temperature dynamics in a more physically based manner, using the 1-D heat equation with phase change

C∂T∂t+Lζ∂θ∂t=∂∂z(λ∂T∂z), z∈[zs,zb](1)

and diffusive and gravitational movement of water in thawed ground by solving Richard's equation

∂ζ∂t=∂∂z(k[ ∂ψ∂z+1 ]), z∈[0,zb].(2)

[7] Here T = T(z, t) is the temperature, ζ=ζ(z,t) is the volumetric water content, ψ=ψ(z,t) is the soil matrix potential. The quantities C=C(T,z)[Jm−3K−1] and λ=λ(T,z)[Wm−1K−1] represents the volumetric heat capacity and thermal conductivity of soil, respectively; L [Jm−3] is the volumetric latent heat of fusion of water, and θ=θ(T,z) is the so-called “unfrozen liquid pore water” fraction, and k = k(T, z) is the hydraulic conductivity. Details of the numerical solution of the heat equation can be found in Appendix Modeling of Soil Temperature. The interested reader may consult the CLM 4.0 technical description for details of the numerical solution of the Richard's equation. The equations are solved implicitly with a daily time step. At the upper boundary condition air temperature and precipitation are prescribed as described below. We emphasize that equation (2) is applicable only for the thawed ground material. In order to extend to simulations of water motion in frozen ground, some models (e.g., CLM 4.0 [Oleson et al., 2010]) propose that ψ=ψ(T,θ,x) if T < T_p, where T_p is the so-called freezing-point temperature depression. In this work, we assume that the water migration in the frozen ground is negligibly small. The latter can be modeled by assuming that the coefficient of hydraulic conductivity k(T, z) = 0, if T < T_p. Thus, when a layer of the ground material becomes frozen the total water content ζ in it stays constant until the moment when the layer becomes thawed again. Note that for frozen soil layers the matric potential can be arbitrarily defined, since it does not enter into calculations, and the water flux boundary condition is imposed at the bottom of the thawed region. Another distinction of the proposed model from the CLM is the way soil thermal properties are parameterized. Additional information is provided in Appendix Modeling of Soil Temperature.

2.2. Snow Model

[8] Seasonal snow cover has a significant influence on the ground thermal regime [Zhang, 2005]. For example, snowcover provides an insulating layer and limits the degree of soil cooling during cold winter months. Classification maps derived from wind, precipitation, and air temperature data [Sturm et al., 1995] can be used to infer snow thermal and physical properties for simulations with numerical models. Land surface models (LSMs) capable of simulating soil freeze/thaw dynamics have recently begun to focus on the effects of snow on underlying soil temperatures.

[9] Formulations in the PWBM for snow accumulation, sublimation, melt, and associated processes are described in Rawlins et al. [2003]. The simulated snowpack contains both a solid and a liquid portion, providing a total model estimate for snow water equivalent (SWE). New routines to account for seasonal changes in snow density follow on recent improvements to the CLM and to the Simple Biosphere/Carnegie-Ames-Stanford Approach (SiBCASA). The CLM version 4 includes a snow model which simulates processes such as accumulation, melt, compaction, snow aging [Lawrence et al., 2012]. The Simple Biosphere/Carnegie-Ames-Stanford Approach (SiBCASA) uses the snow classification system of Sturm et al. [1995], and a snow model derived from the CLM. We take here a similar approach in modeling the temporal evolution in snowpack density. Details of our new snow density model are described in Appendix Modeling of Soil Temperature.

2.3. Organic Content and Parameterizations

[10] In many parts of the Arctic where soil carbon is high, the first 40–50 cm of soil is nearly 100% organic, with a transition from organic to mineral soil and 100% mineral soil below. In the middle transition zone, thermal and hydraulic properties of mixed mineral and organic soil material can be approximated as a weighted combination of the mineral soil and organic soil properties. As described in Rawlins et al. [2003], the previous version of PWBM contained an upper organic layer and a lower mineral layer imposed in the two layer profile.

[11] In order to better account for the thermal and hydrologic properties of soils, we parameterized carbon density in each soil layer. We take a similar but not identical approach to the one described in Lawrence and Slater [2008]. The Global Soil Data Task (GSDT, 2000) data set contains soil-carbon density (C, kg m⁻³) across the depth interval of 0–1 m. To obtain C across the pan-Arctic basin, we averaged the five arc-second GSDT data for each EASE-Grid cell. We applied the soil profile for polar and boreal soils from Zinke et al. [1986] to obtain carbon storage over the top 11 model soil layers (1.4 m depth). Soil carbon or organic fraction for each layer was then determined as

fsc,i=ρsc,i/ρsc,max(3)

where fsc,i is the carbon fraction of each layer i, ρsc,i is the soil carbon density, and ρsc,max is the maximum possible value (peat density of 130 kg m⁻³, Farouki [1981]). Soil properties for each layer are specified as a weighted combination of organic and mineral soil properties

P=(1−f)Pm+fPo(4)

where f is the fraction of organic material in the soil layer, Pm is the value for mineral soil, Po is the value for organic soil, and P is the weighted average quantity. Thermal and hydrologic parameters as a function of soil class are listed in Table 1.

Table 1. Soil Parameters Used in the PWBM Simulationsa

Soil Type	λ (W m−1 K−1)	C (J m−3 K−1 × 106)	Θsat	ksat (m s−1 × 10−3)	Ψsat (mm)	β
a Each grid cell in the model is characterized by one of the five mineral soil types. Model parameters are defined through a weighted combination of organic and mineral soil properties. See Rawlins et al. [2003] for more detail on the PWBM soil routine.
Sand	λm = 3.6	Cm = 3.0	0.25	0.023	−47	3.4
Loam	λm = 3.0	Cm = 3.0	0.35	0.042	−207	6.1
Clay	λm = 2.3	Cm = 3.0	0.45	0.020	−390	12.1
Sandy loam	λm = 3.3	Cm = 3.0	0.40	0.071	−132	4.5
Clay loam	λm = 2.6	Cm = 3.0	0.39	0.028	−289	8.2
Organic soil	λo = 1.5	C0 = 1.9	0.9	0.02	−120	2.7

3.1. Forcing Data

[12] Many of the static input data fields (e.g., soil properties, landcover type, and snow class) and meteorological forcings used in the prior and present version of the PWBM are available within the ArcticRIMS project archive (http://rims.unh.edu/). For the two pan-Arctic simulations, we draw daily 2 m air temperature, precipitation, and wind speed from two atmospheric reanalyses data sets. Atmospheric reanalyses are retrospective forms of numerical weather prediction using a fixed model and data assimilation system. The first (hereafter referred to as NNR) was derived from the NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) effort [Kalnay et al., 1996] (in ArcticRIMS see: http://rims.unh.edu/data/read_me.cgi?category=7&subject=0 and http://rims.unh.edu/data/read_me.cgi?category=2&subject=0). The second set of reanalysis data (hereafter ERA-40) was drawn from the European Center for Medium Range Weather Forecasts (ECMWF) ERA-40 reanalysis [Uppala et al., 2005] (in ArcticRIMS see: http://rims.unh.edu/data/read_me.cgi?category=7&subject=8 and http://rims.unh.edu/data/read_me.cgi?category=2&subject=5). No ERA-40 wind speed data is available in Arctic-RIMS and so NNR wind speed are used in their place for the simulations described below.

[13] Across the terrestrial Arctic, total precipitation over cold season months can be used as a proxy for snowfall. As a first step, we assessed biases in NNR and ERA-40 total precipitation for November–March using precipitation data available from the University of Delaware (UDel) [Matsuura and Willmott, 2009; Willmott and Matsuura, 2000]. The UDel data sets were developed through interpolations of meteorological station data records. For this assessment, average November–March precipitation was taken over the period 1980–1999. Bias maps are shown in Figures 2a and 2b. The distribution of biases is largely positive across the pan-Arctic. Integrated area average biases (excluding Greenland) are 0.17 and 0.24 mm day⁻¹ for NNR and ERA-40, respectively. Across Alaska, biases are also mostly positive. For the region west of 135°W, averages are 0.21 mm day⁻¹ for NNR and 0.31 mm day⁻¹ for ERA-40.

3.2. Ground-Based Validation Data

[14] Evaluations of model simulated soil temperatures are made using multiple within-site observations from a subset of the Alaska Ecological Transect (ALECTRA) sites (Figure 3). We chose ALECTRA sites that have sufficient multiyear data over the period 2000–2005. The evaluation sites are Bonanza Creek, Coldfoot, and Dietrich (Table 2). Missing data due to logger or sensor issues limit the number of sites available for examination. The ALECTRA sites were conceived and implemented to capture spatial heterogeneity in soil and vegetation stem/branch temperature conditions across the landscape that largely reflect microclimatic variability. Many of the ALECTRA sites are situated near research stations where meteorological observations are available. The sites are distributed along a north-south latitudinal transect extending from Arctic coastal tundra through the boreal and into maritime forest. Soil temperature and snow data for Council were drawn from records archived at the University of Alaska Water and Environment Research Center (http://ine.uaf.edu/werc/). For the simulations at each of the four research sites (Table 2), we used the average 2 m air temperatures across available sites. Precipitation and wind speed were drawn from the NNR.

[15] Observations of ALT from the Circumpolar Active Layer Monitoring (CALM) data set [Brown et al., 2000] are used to evaluate PWBM simulated estimates. Field sampling in the CALM program typically involves measurements across a grid spanning 0.1 km² or in some locations 1 km² area. Each recorded ALT value represents maximum seasonal depth of thaw of the 121 samples. For our comparisons with PWBM-simulated values, we take the average of all nonmissing CALM ALTs over the period 1980–2011.

4.1. Site Comparisons

[16] Figure 4 shows the simulated snowpack behavior over several years for each of the four study sites. At each site, snow density increases through the cold season, starting at around 100 kg m⁻³. The profiles for Bonanza Creek and Coldfoot are similar, with snow densities spanning a range between 100 and 250 kg m⁻³. A larger range is evident for Dietrich and Council, where end-of-season snow density is over 300 kg m⁻³. Comparing simulated snow depth with observations reveals no consistent biases. In some years, there is good agreement, in some years the model overestimates snow depth, in other years it underestimates. This is worth noting given the positive bias (overestimate) in cold season precipitation (Figure 3).

[17] The thermal conductivity of organic soils is low, ranging from 0.50 W m⁻¹ K⁻¹ at saturation to 0.06 W m⁻¹ K⁻¹ under dry conditions [Farouki, 1981]. Figure 5 shows model simulated thermal conductivity profiles for May, June, September, and October at the four sites. For Bonanza Creek, Coldfoot, and Council the simulated organic density as parameterized is nearly 100% in upper 20 cm and transition to fully mineral soil at around 40 cm. Note the much thinner modeled organic horizon for the Dietrich site. Simulated thermal conductivity there increases sharply down to around 30 cm depth where it approaches 2.4 W m⁻¹ K⁻¹. For Dietrich, the reduced thermal conductivity in the soil layer centered at 0.55 m is largely attributable to lesser water content. This soil layer is the lowest that thaws over the spin-up period during which time it looses some water. During the simulation period, permafrost again develops at this layer, however, the water content and thermal conductivity then is lower. At each site, conductivities tend to be higher in spring versus autumn. This result reflects higher soil moisture amounts, and thus thermal conductivity, following spring snowmelt.

[18] Soil temperatures are influenced by a number of climate and landscape factors [Callaghan et al., 2011]. The seasonal cycle in air temperature and seasonal accumulation and melt of snowpack are perhaps the two most important processes that influence annual variations in soil temperatures. Examining simulated soil temperatures against observed data provides insights into model efficacy. Figures 6a and 6b show simulated and observed daily soil temperatures at 25 cm depth for Bonanza Creek and Coldfoot. Colored traces represent temperatures from measurements in each landscape units within the research sites. For Bonanza Creek, model simulated summer temperatures are well captured, falling between the warm south facing slope and the cooler spruce sites. Model temperatures during winter are clearly too warm. At Coldfoot, summer temperatures again are well simulated. In contrast with Bonanza Creek, 25 cm temperatures during the first 5 months of the year fall in the middle of the range of observations.

[19] Figure 7 depicts monthly average soil temperature at both 25 and 50 cm depth for Bonanza Creek, Coldfoot, Council, and Dietrich. For Bonanza Creek, simulated temperatures in summer fall well within the range of the measured data from the three landscape units (Table 2). In winter, the model is warmer than observations, with simulated temperatures at 25 and 50 cm remaining near freezing through the winter. Model simulated snow depth compares well with observations over the winters of 2000–2001 and 2001–2002 (Figure 4). The warm soil temperature bias in winter 2002–2003 is consistent with an overestimation in simulated snow depths. Over the period 1980–2008, precipitation averages 374 mm yr⁻¹ in the daily NNR data versus 286 mm yr⁻¹ for National Weather Service observations from Fairbanks International airport over the climate normal period 1971–2000.

[20] For Coldfoot, model simulated soil temperatures at 25 cm depth fall well within the measured range across the landscape units during most months (Figure 7). Simulated snow depths are lower than the observed values in two of the three years shown. The exception is winter 2002–2003 when the model overestimates temperatures. Simulated temperatures at 50 cm depth remain near 0°C through summer, whereas the observations are 1–4°C warmer.

[21] At Council, simulated temperatures are colder than the observations by some 5°C during the coldest months (Figure 7). Simulated midwinter snow depths in 2000–2001 were approximately half (40 versus 80 cm) of the observed values. This bias could explain some of the discrepancy in simulated soil temperatures. With no south facing slope (Table 2), observed temperatures at Council are notably cooler than those at Bonanza Creek and Coldfoot sites. Both the annual temperature cycle and the range in temperatures across the four landscape units at Council are small. Good agreement with the observations occurs in summer each year, where the model and observations at 25 and 50 cm are in the range 0–5°C. The small variations around 0°C in the winter 2000–2001 observations at Council may be a result of a deep snowpack and/or a reflection of a thick organic layer, both of which would tend to limit the influence of the overlying cold winter air.

[22] For Dietrich, simulated temperatures fall within the observed range during most months, excluding spring of 2004 and 2005 when simulated temperatures are approximately 1–2°C colder than the observations (Figure 7). No observed snow depth data are available for Dietrich. Like with Bonanza Creek, Coldfoot, and Council sites, simulated temperatures in summer agree well with observations. When compared with Bonanza Creek, the larger degree of cooling in the Dietrich simulation is largely a reflection of the thinner organic layer and, in turn, higher thermal conductivities (Figure 5).

[23] We performed a series of sensitivity experiments to help elucidate potential sources of bias in simulated soil temperatures at depth during both the cold and warm season. In particular, we examined the simulated soil temperatures at 25 cm at Bonanza Creek and Council sites. For Bonanza Creek, reducing snowfall over the cold season by 50% cools temperatures at the 25 cm depth over the cold season (November–March, 2001–2000) by 2.3°C. These results suggest that any errors in snow depth can explain some but not all of the bias in simulated winter temperatures at depth. We then analyzed the effect of errors in organic content (by carbon density in the model) of soil layers by decreasing the organic content by 50%. This change resulted in a cooling of temperatures by approximately 1°C. In view of these results, we speculate that a combination of errors in snow depth and organic layer thickness (carbon density) could well explain most of the bias in simulated soil temperatures at depth. For Council, doubling of cold season snowfall warms soil temperatures at 25 cm (1999–2000) by 3.6°C. This result suggests that errors in the approximation of snow depth can explain most but not all of the bias in the simulation results shown in Figure 7.

[24] We also performed a pair of experiments to better understand model sensitivity during the warm season. For these simulations we simultaneously perturbed soil layer carbon density and air temperature and focused on Bonanza Creek sites. First, in order to mimic typical conditions across the south-facing slope, we approximated an organic layer thickness of 20 cm and as described above reparameterized the resultant carbon density of each layer. At the same time, we scaled the input 2 m air temperatures upward by 2°C. In the second sensitivity experiment, we approximate condition in either the white spruce or black spruce stand by parameterizing a 60 cm organic layer and resultant carbon densities at depth, along with 2°C cooler temperatures. A difference of 4.2°C was found between the two simulations. Our results here demonstrate that the numerical soil freezing and thawing scheme, parameterizations and input data for organic content, and input air temperature forcings are able to capture much of the variations in observed summer soil temperatures across the Bonanza Creek sites.

4.2. Simulated Active Layer Thickness and Permafrost Extent Across the Pan-Arctic

[25] To further evaluate model performance, we ran additional simulations across the entire pan-Arctic drainage basin. Model spin-up was performed through 50 iterations over the year 1980, the first year of the transient simulation. We use the daily precipitation, air temperature, and wind data derived from the reanalysis products as described in section 'Forcing Data'. A grid cell is labeled as containing permafrost if at least one layer within the upper 15 soil layers (3.25 m depth) remains frozen throughout the year. Simulated ALT is estimated as the depth where soil temperature crosses the 0°C threshold. In these evaluations, observed CALM ALT represents all available nonmissing values for a given site over the period 1980–2011.

[26] Figure 8 shows comparisons between simulated and observed ALT over the period 1980–1999 from the prior (a,b) and updated (c,d) versions of PWBM. Simulated mean ALT is taken as the 30 year average for the grid cell encompassing each respective CALM location. Mean biases (simulated minus observed ALT) from the prior version are −36 cm and −27 cm for the NNR and ERA-40 simulations, respectively. With the updated version, biases are −4.0 cm and 5.9 cm for NNR and ERA-40, respectively. A systematic underestimation of ALT is evident in simulations with the prior version. It should be noted that a scale mismatch exists in any comparison involving grid-to-point observations. While the updated model simulates well the mean among the observed ALTs, it struggles with predictability. This shortcoming is not uncommon among numerical models which simulate active layer dynamics [Su et al., 2005; Lawrence et al., 2012]. In a study of ALT estimates across northern Alaska produced by three different models, Shiklomanov et al. [2007], found that large differences in ALT were attributable to the different approaches used for characterization of surface and subsurface conditions. Measured ALT can vary substantially over small distances. At any given location, ALT is influenced by several factors including thermal properties of soil, snow characteristics, vegetation type, and soil moisture amount. It has been shown that variations in ALT of up to 50% can occur over short distances [Nelson et al., 1997].

[27] At the pan-Arctic scale, the PWBM captures well the north-south ALT gradient (Figure 9). Estimates range between 30 and 40 cm along the Arctic coast and approach 1.5 m in parts of southern and western Siberia and central Canada. The simulation with NNR forcings (Figure 9a) produces slightly greater ALTs compared to the ERA-40 simulation (Figure 9b). Simulated permafrost areal extent is 14.6 (14.4–14.8) × 10⁶ km² and 13.3 (13.0–13.6) × 10⁶ km² from the NNR and ERA-40 simulations, respectively. This compares with an observed area of continuous plus discontinuous permafrost of 12.5 × (11.8–14.6) 10⁶ km² [Zhang et al., 2000]. Thus, our simulated areas are 7 and 17% above observed extent for the ERA-40 and NNR simulations, respectively. However, as Lawrence et al. [2012] noted, permafrost area obtained from a coarse model is likely biased high. If true, the 14.6 × 10⁶ km² (continuous plus discontinuous permafrost) extent estimated by Zhang et al. [2000] would represent a better goal for model simulations.

[28] Figure 10 shows the results of a comparison between simulated runoff in the present study (labeled PWBM₂₀₁₃) and those described in Rawlins et al. [2003] (PWBM₂₀₀₃). Marked differences are evident. For the Yenisei, Lena, and Yukon basins, annual runoff in the present simulation is closer to observed values from the R-ArcticNet v4.0 archive (http://www.r-arcticnet.sr.unh.edu/v4.0/) [Lammers et al., 2001]. For the Ob and Nelson basins, observed runoff is relatively low and simulated runoff is clearly overestimated. Annual observed runoff averaged over the pan-Arctic basin (1979–2001) is approximately 230 mm yr⁻¹ [Rawlins et al., 2010]. Comparing this recent estimate to the 180 mm yr⁻¹ (1980–2001) from Rawlins et al. [2003] gives a bias (predicted minus observed) of approximately −22%. Using the same NNR data as Rawlins et al. [2003] in the updated PWBM gives an annual runoff of 250 mm yr⁻¹, a bias of 9%. We note that the observed runoff/precipitation (R/P) ratio for the Nelson basins is low relative to other high latitude river basins.

[29] Figure 11 illustrates the differences in soil moisture distribution with depth between simulations with the updated and prior model versions. It shows soil moisture variations for the grid cell encompassing the Bonanza Creek site from simulations forced with the measured (site) daily air temperatures along with NNR precipitation and wind data. The soil profile in the new version exhibits expected patterns, with nearly saturated mineral soils above the permafrost (Figure 11a). Upper organic rich layers become wet following snowmelt and are considerably drier than the underlying mineral soils in summer. ALT averages 68 cm, close to the mean observed value of 55 cm from the CALM data set. In the prior version of the model (Figure 11b), the ALT reaches approximately 30 cm each year. Soil moisture is lower as well. In the prior version, soil ice thaw was calculated based on the fraction of the layer which thaws on a given day. This results in relatively higher soil ice and, hence, lower soil water amounts through summer.

4.3. Potential Future Changes in Near Surface Conditions

[30] Lastly, we explored potential future changes in soil thermal and moisture regimes using data from the North American Regional Climate Change Assessment Program (NARCCAP, https://www.narccap.ucar.edu/). The NARCCAP [Mearns et al., 2007] is archiving outputs from a set of regional climate model (RCM) simulations over a domain spanning North America. The NARCCAP includes six RCMs, each forced with boundary conditions from two atmosphere-ocean general circulation models (GCMs). For each model pairing, air temperature and precipitation are available for two 30 year periods: 1971–2000 and 2041–2070 (hereafter present and future, respectively). From the native RCM grids, we selected the cell which encompasses the Bonanza Creek area and averaged daily air temperature and precipitation across the available GCM-RCM pairings. Wind speed was taken as the daily climatology from the NRR data. Model spin-up for each 30 year simulation was performed by iterating over the first year of the present and future periods (1971 and 2041, respectively). For the grid chosen, a moderate cold and wet bias is present. To ameliorate the influence this bias would place on results we scaled each daily time series using monthly observations for Fairbanks International airport. Table 3 lists original and bias-adjusted values. Temperature change between the future and present periods is 2.9°C and precipitation change is just over 30%. As described below, our analysis of potential changes in ALT, thawed season length, soil temperature and moisture focuses on averages over the periods 1996–1999 and 2066–2069.

[31] Figure 12a shows simulated soil temperature as a function of depth for the period 1996–1999. Forcing the model with the bias-adjusted NARRCAP data results in a seasonally frozen (nonpermafrost) soil column. We note that average grid cell temperature over the 1971–1999 period with this simulation is just over 1.5°C warmer than the NNR and ERA-40 data. The absence of permafrost here is not entirely unexpected as the area around Bonanza Creek is in the discontinuous permafrost zone. For the present day simulation, the 30 cm soil layer is above 0°C for an average of 156 days (DOY 118–273). In the future simulation (Figure 12b), the thawed period increases to 166 days (DOY 115–280) and June–August temperature are warmer by 1.9°C. In the present simulation, soil temperatures at 45 cm are below 0°C for 158 days (late December to early June). In the future simulation frost does not reach 45 cm during winter. By the late 2060s, increased precipitation contributes to a deeper simulated snowpack, with depths in February and March greater by 26 and 28%, respectively (Figure 13b). The end of snowmelt advances by approximately 10 days (May 3 to April 23). For the present simulation soils between 30 and 40 cm thaw near DOY 160, experience recharge from above, and average 44% water content during the unfrozen period. In the future simulation, soil moisture averages 39%, with the 30 and 40 cm layer unfrozen all year. Soil moisture between 50 and 100 cm lowers from 58% to 54% in the future simulation.