A simple model is developed to evolve daily ground temperatures from surface air temperatures (SATs) in snow-dominated areas. Ground surface temperatures (GSTs) are calculated by propagating the daily SAT through the snowpack, and attenuating the signal amplitude. Subsequent subsurface heat transfer is then modeled using the analytical solution of the one-dimensional heat conduction equation. The thermal impacts of nonconductive heat transfer processes and seasonal freeze thaw are implicitly represented by the time-dependent apparent thermal diffusivity of the subsurface. The model is tested in four snow-dominated regions: Barrow, Council, Ivotuk (all in Alaska) and Reynolds Creek Experimental Watershed (in Idaho). The model captures the seasonal evolution of the ground temperature at all sites. The model demonstrates the feasibility of simulating subsurface temperatures using only air temperature and snow depth.
In contrast to other land features, snow has a number of unique properties: it is a winter-only feature, it has a substantially higher albedo than the surrounding vegetation, and has characteristically low thermal conductivity. Specifically, the surface albedo of snow ranges from 0.6 to 0.85, while that of vegetation ranges from 0.1 to 0.3, and the thermal conductivity of snow ranges from 0.1 to 0.5 Wm−1 K−1, while that of soil ranges from 0.8 to 2.2 Wm−1 K−1. In winter, therefore, snow strongly alters the surface energy budget (Yeh et al. 1983; Namias 1985; Barnett et al. 1989) and prevents effective heat exchange between the ground and the atmosphere. From an ecological standpoint, variations in timing and depth of snow cover impact soil processes such as winter net nitrogen mineralization (Schimel et al. 2004) and net ecosystem CO2 efflux (Welker et al. 2000). Evidence from long-term snow manipulation experiments suggest that an increase in soil temperature due to an increase in snow depth (SD), as projected during this century, can lead to degradation of the permafrost (Hinkel and Hurd 2006).
Models of varying degrees of complexity have been developed to describe snow dynamics and the coupling of land, snow, and boundary layer processes. Some models realistically capture the exchange of energy, mass, and momentum across the atmosphere–snow–ground system, and explicitly include a suite of snowpack processes such as the dynamics of gravitational settling, metamorphism, phase changes, and heat transfer through the percolation of water (Anderson 1976; Jordan 1991; Loth et al. 1993; Lynch-Stieglitz 1994; Stieglitz et al. 2001; Bartelt and Lehning 2002; Liston and Elder 2006). While these complex models can simulate ground temperatures with high fidelity (Lynch-Stieglitz 1994; Koster et al. 2000; Stieglitz et al. 1999, 2001), they are computationally expensive. Other models have been developed using simple empirical schemes (Bartlett et al. 2005; Pollack et al. 2005; Stieglitz and Smerdon 2007). They employ frameworks that implicitly represent snow insulation impacts on subsurface temperatures by use of simple governing equations that depend on few parameters. For example, Stieglitz and Smerdon (2007) employed a one-dimensional diffusion equation coupled to the surface air temperature (SAT) through a time varying flux boundary condition at the land surface. The time varying flux boundary condition is a function of the SAT, ground surface temperature (GST), and a coupling function. The temporal character of the coupling function implicitly represents the cumulative thermal effects of the processes operating at the land–atmosphere interface. These processes include snow insulation, vegetative insulation, freeze–thaw processes, vapor transport in soils, and evapotranspiration.
We develop a simple semiempirical model to evolve daily ground temperatures from daily SAT and SD in snow-dominated areas. We generate the daily GST by propagating the daily SAT through the snowpack and attenuating the SAT signal amplitude. Subsequent subsurface heat transfer is then modeled using the analytical solution of the one-dimensional heat conduction equation. The thermal impacts of nonconductive heat transfer processes and seasonal freeze–thaw are implicitly represented through a time-dependent apparent thermal diffusivity (ATD). The model is tested in four snow-dominated regions: Barrow, Council and Ivotuk, all in Alaska, and Reynolds Creek Experimental Watershed (RCEW) in Idaho.
a. Model description
The analytical solution to a sinusoidal signal of mean temperature, T, and amplitude Ao applied at the surface of a homogenous infinite half-space is (Carslaw and Jaeger 1959)
where z is distance from the surface of the half space. Here ω is the radial frequency, which is 2π times the actual frequency of the signal. Here A is the signal attenuation of the following form:
and ϕ is the phase lag:
where k is the wave vector:
and λ, D, τ, are the damping depth for the snow/ground, defined as the characteristic depth at which the temperature signal is attenuated to 1/e of the surface, the thermal diffusivity of the snow/ground, and the period of the forcing, respectively.
We modify Eq. (1) to evolve GST from SAT in the presence of snow cover. Based on observed daily SAT and GST data for Council, Ivotuk (Fig. 1) and at Fargo, North Dakota (Smerdon et al. 2003), the SAT signal is not phase lagged with depth (see section 4). Equation (1) then reduces to
where λsnow is the seasonal damping depth for snow, which is approximately 67 cm (Hillel 1998) for a snowpack of density 300 kg m−3. In our model, snow is an insulative material that only attenuates the daily SAT signals. The deeper the snowpack, the higher the attenuation of the daily SAT signals. During snow-free periods, SD(t) = 0, GST = SAT (see section 4).
To generate subsurface ground temperatures from GST, Eq. (1) is again modified such that T = Tg, where Tg is the observed annual mean ground temperature, A0 = GST(t) − Tg, and ω = 2π/365. At depth z, A0 is attenuated by a factor of exp[−z/λsoil(t)] and phase shifted by the calibration parameter, ϕd(z). Here λsoil(t) is the damping depth for the soil such that
where Dh is the time-dependent ATD (McGaw et al. 1978). The ATD represents the combined effects of conductive and nonconductive heat transport processes (Chen and Kling 1996; Hinkel 1997; Hinkel et al. 2001) resulting from freeze–thaw, soil water evaporation, and movement, etc. Latent heat effects associated with seasonal freeze–thaw increases the apparent volumetric heat capacity of the soil and decreases the ATD. Therefore, Dh will vary over the course of a year. For example, a typical value of specific heat capacity for water is 4200 kJ m−3 K−1 (Hillel 1998), while that of ice is 1900 kJ m−3 K−1 (Hillel 1998). As the soil thaws, the apparent heat capacity can become higher and reduce the ATD by at least an order of magnitude (Ochsner and Baker 2008). To capture this seasonal variation of Dh, we choose a simple sinusoidal function (see section 4):
where Dh is the mean ATD and B is the amplitude of the sinusoid. Optimal values of Dh and B are determined through calibration; Dh is constrained to be below the thermal diffusivity of pure ice (11.6 × 10−7 m2 s−1 or 0.1 m2 day−1) and above zero. This constraint is necessary since ATDs with values less or equal to zero will not yield physically meaningful solutions when used in the analytical solution [e.g., Eq. (6)]. ATD calculated using this methodology therefore represents the average bulk thermal diffusivity over the entire soil column. For high Dh, the damping depth is high, resulting in lower attenuation of the temperature signal at a given depth. For low Dh, the damping depth is low, resulting in higher attenuation of the temperature signal at the same depth.
b. Site and data description
The model behavior is evaluated using daily SAT(t) and SD(t) provided at four snow-dominated sites: Barrow, Alaska (71.3°N, 156.8°W); Ivotuk, Alaska (68°29′N, 155°44′W); Council, Alaska (64°53′N, 163°40′W); and RCEW, Idaho (43°05′N, 116°43′W). Barrow is located in northwestern Alaska on the coast of the Arctic Ocean. The climate is cold and dry with a mean annual air temperature of −12.2°C (1949–2003) and annual solid precipitation of 745 mm (water equivalent; more information is available online at http://www.wrcc.dri.edu/summary/climsmak.html). Barrow is snow covered for ∼270 days yr−1 (Zhang et al. 1996). Ivotuk is located 350 km (220 miles) south of Barrow, at the foothills of the North Slope of the Brooks Range. The site experiences mean annual air temperature of −10.9°C and mean annual precipitation of 202 mm (Riedel et al. 2005). Council is located in the central Seward Peninsula and represents a transitional area between the boreal forest and tundra. Mean annual air temperature at Council ranges from −4.06° to −4.62°C (Chapin et al. 2006). The RCEW is located in the Owyhee Mountains in southwestern Idaho. Mean annual air temperature ranges from 4.7° to 8.9°C (Hanson et al. 2001) while mean annual precipitation ranges from 230 mm to greater than 1100 mm (Slaughter et al. 2001). Details of the datasets for each site can be found in Table 1.
Data were inspected for missing values. Periods of less than one month with missing values were filled in through linear interpolation, while periods of missing data greater than one month were excluded from calculations. Daily SAT and SD time series records used in simulations for the four sites are shown in Fig. 2. Model performance is evaluated using root-mean-square errors (RMSEs), σ:
where Toi is the observed temperature, Tsi is the simulated temperature, and N is the total number of data points.
Figure 3 illustrates the best-fit ATDs for the four sites: for each site, ATDs peak toward the end of spring and decrease toward zero until fall. Values range from 0.0025 to 0.0775 m2 day−1 in Barrow, 10−5 to 1.99 × 10−3 m2 day−1 in Ivotuk, 3.5 × 10−4 to 1.05 × 10−3 m2 day−1 in Council, and 0.0015 to 0.0585 m2 day−1 in RCEW (Fig. 3).
Figure 4 illustrates observed and simulated ground temperature at 50 cm from 1 January 1990 to 31 December 1997 in Barrow. At Barrow, observed ground temperature at 50 cm ranges from −25° to 15°C. For the entire 8-yr period, simulated ground temperature is in good agreement with the observed ground temperature. RMSE for Barrow is 3.87°C. Figure 5 illustrates the observed and simulated ground temperature at 5 and 10 cm from 17 September 2003 to 31 December 2006 in Ivotuk. Although the model captures much of the seasonal trend of the ground thermal regime for the entire period, it consistently underestimates the fall ground temperatures by ∼2.0°C for all years. RMS error for Ivotuk are 3.49°C at 5 cm, and 3.09°C at 10 cm. Figure 6 illustrates the observed and simulated ground temperature at 5, 10, 15, and 25 cm from 1 January 1992 to 11 December 1992 in Council. While the model captures the seasonal trend, it overestimates spring ground temperatures at all depths. RMSEs for Council range from 2.16° to 2.84°C. Figure 7 illustrates the observed and simulated ground temperature at 10, 20, 30, 40, 50, 60, 90, and 120 cm from 12 June 1992 to 30 September 1996 for RCEW. At RCEW, simulated ground temperatures at all depths were remarkably similar to the observed ground temperatures. RMS errors for Reynolds Creek range from 2.39° to 3.49°C. RMSEs for all sites can be found in Table 2.
Simulation results depict clear seasonal trend in ATDs. The simulated ATDs are the lowest in midfall and the beginning of winter (Fig. 3). Ground temperatures during this period are falling and the soil water is freezing, which increases the apparent heat capacity of the soil and lowers the ATD to a minimum value. Ground temperatures continue to drop in winter and the freezing front propagates deeper into the soil. This in turn increases soil ice and the bulk ATD of the soil column. As temperatures rise above 0°C in spring, the frozen soil begins to thaw. As the soil thaws, latent heat is absorbed, which increases the apparent heat capacity, and lowers the ATD of the soil from a maximum value.
To justify our sinusoidal ATDs in Eq. (7), we directly calculate ATDs for Ivotuk, RCEW, and Council from observed ground temperature time series using a finite-difference scheme (McGaw et al. 1978; Outcalt and Hinkel 1989; Hinkel et al. 1990; Hinkel et al. 2001). We then compare them with our sinusoidal ATDs (Fig. 8). Results show that our sinusoidal ATDs in Eq. (7) capture the main seasonal variations in Ivotuk and RCEW: lowest from the end of summer to midwinter and the highest toward the end of spring. At Council, there are no obvious trends in the calculated ATDs. The seasonal trends depicted by our sinusoids are also consistent with the results of earlier works that calculated thermal diffusivity from temporal ground temperature and meteorological data. For example, Hinkel et al. (2001) calculated thermal diffusivity from time series ground temperature records at Barrow and observed that the ATD decreases toward zero around September/October. Likewise, Pollack et al. (2005) quantified 10 yr of daily thermal diffusivity of the shallow subsurface empirically from daily meteorological records at Fargo, North Dakota. For the period of their study, thermal diffusivities are lowest from the end of summer to mid winter and the highest toward the end of spring.
As expected, our sinusoidal ATDs in Eq. (7) do not capture the daily variability (Fig. 8). Daily variability in soil ATDs are dependent on daily fluctuations in soil moisture and nonconductive processes such as evapotranspiration and vapor transport in soils. Moreover, our sinusoidal ATDs do not capture the negative values seen in the calculated ATDs (Fig. 8). These negative calculated ATDs have been interpreted to reflect the dominance of the nonconductive heat transfer over conductive heat transfer in the soil column (Hinkel et al. 2001). In our modeling framework, however, ATDs with values less or equal to zero do not yield physically meaningful solutions [e.g., Eq. (6)], and therefore have been constrained to values greater than zero. Nevertheless, the validation strategies discussed above lend support to the notion that our sinusoidal ATDs capture the seasonal behavior of ATDs at our study sites.
Simulated mean ATDs differ by an order of magnitude across the sites (Fig. 3). At Barrow and RCEW, simulated mean ATDs are 0.04 and 0.03 m2 day−1, respectively, and are consistent with typical values of 0.0173 − 0.0432 m2 day−1 (Hillel 1998; Hinkel et al. 2001). At Ivotuk and Council, simulated mean ATDs are 10−3 and 7 × 10−4 m2 day−1, respectively, and are an order of magnitude lower than typical values. However, at Ivotuk and Council, simulated ATD are of the same order of magnitude as the calculated ATDs (Fig. 8). To ascribe why the differences in ATDs exist between the sites, a more complex model that allows for nonconductive heat transfer processes is needed.
Although our model captures the seasonal evolution of ground temperatures, there are caveats to consider: 1) the depth dependence of the ATD has not been considered. Mineral soil fraction and soil moisture content vary with depth and affect the magnitude of the ATD. ATD typically decreases as volumetric soil moisture increases (Hinkel et al. 2001) but increases with mineral soil fraction. Future work is needed to incorporate a depth-dependent ATD. 2) During periods of snow cover we assume a constant snowpack density of 300 kg m−3, which translates into a constant λsnow of 67 cm for the snowpack. This assumption of a fixed snow density therefore results in early and late season biases, which contribute to notable differences between the simulated and observed ground temperatures during fall in Ivotuk and during spring in Council. Studies have shown that snowpack density varies geographically and seasonally (Lynch-Stieglitz 1994; Sturm and Benson 1997; Marks et al. 2001): between 150 kg m−3 for fresh snow to 500–700 kg m−3 for the end-of-season snowpack (Zhang 2005). 3) To maintain the simplicity of the model, we assume that ground–atmosphere decoupling due to vegetation is zero (SAT = GST) during snow-free periods. Analysis of the SATs and GSTs for a number of sites where snow cover is significant (more information is available online at http://public.ornl.gov/ameriflux/), has shown that the ground–atmosphere decoupling due to snow is significantly greater than the ground–atmosphere decoupling due to vegetation (Stieglitz and Smerdon 2007). However, Stieglitz and Smerdon (2007) have also shown that ground–atmosphere decoupling due to vegetation can be significant in some sites. For example, at Campbell River, British Columbia, Canada, summer ground–atmosphere decoupling due to vegetation can attenuate the GST by as much as 25%. 4) We observe no phase lag between SAT and GST in any of our datasets and model the GST as described by Eq. (5). This seems to contradict the common understanding of the propagation of a sinusoidal signal through any conductive material as described by Carslaw and Jaeger (1959). For example, a hypothetical and permanent snowpack that is 1 m thick with λsnow of 67 cm, will result in ∼87 days of phase lag between SAT and GST. It seems that there are nonconductive processes such as vapor transport and water transport within the snowpack, which are responsible for significant reduction in phase lag between the SAT and GST. To isolate the thermal effects of these nonconductive processes, a more complex model is needed and is beyond the objectives of this work. 5) Both λsnow and λsoil(t) are modeled as a function of a single radial frequency with an annual period. This assumption implies that all temperature fluctuations with higher frequencies will not be damped out by the snow and soil. This in turn results in higher variability in the simulated ground temperatures than the observed ground temperatures. Despite all the caveats, we have constructed a useful modeling framework that captures the first-order controls of the ground–atmosphere heat transfer in snow-dominated regions. This simple scheme permits rapid prediction of daily ground temperatures over large areas using only daily SAT and SD.
A simple modeling framework is developed to evolve daily GST from SAT in snow-dominated areas and subsequently simulate ground temperatures using the analytical solution to the one-dimensional thermal diffusion equation. Complex latent heat effects of seasonal freeze–thaw were incorporated into the framework through a time-dependent ATD. Advantage of this modeling framework is that it permits rapid prediction of ground temperatures using only SAT and SD data. Coupling of this simple ground temperature model with a spatially distributed snow model capable of generating SD will provide a powerful tool for determining the magnitude of change in ground temperatures in high-latitude regions undergoing changes in winter precipitation/snow at various spatial scales.
This research was supported in part by the following NSF Grants 0439620, 0436118, and 0922100, as well as by a grant from the U.S. Dept. of the Interior. Yiwei Cheng was additionally supported in part by an Everglades Foundation Fellowship from the Everglades Foundation.
* Current affiliation: Department of Geography, University of North Texas, Denton, Texas
Corresponding author address: Yiwei Cheng, School of Civil and Environmental Engineering, Georgia Institute of Technology, 790 Atlantic Dr. NW, Atlanta, GA 30332-0355. Email: email@example.com