Articles | Volume 22, issue 5
Research article
07 May 2018
Research article |  | 07 May 2018

A simple temperature-based method to estimate heterogeneous frozen ground within a distributed watershed model

Michael L. Follum, Jeffrey D. Niemann, Julie T. Parno, and Charles W. Downer

Frozen ground can be important to flood production and is often heterogeneous within a watershed due to spatial variations in the available energy, insulation by snowpack and ground cover, and the thermal and moisture properties of the soil. The widely used continuous frozen ground index (CFGI) model is a degree-day approach and identifies frozen ground using a simple frost index, which varies mainly with elevation through an elevation–temperature relationship. Similarly, snow depth and its insulating effect are also estimated based on elevation. The objective of this paper is to develop a model for frozen ground that (1) captures the spatial variations of frozen ground within a watershed, (2) allows the frozen ground model to be incorporated into a variety of watershed models, and (3) allows application in data sparse environments. To do this, we modify the existing CFGI method within the gridded surface subsurface hydrologic analysis watershed model. Among the modifications, the snowpack and frost indices are simulated by replacing air temperature (a surrogate for the available energy) with a radiation-derived temperature that aims to better represent spatial variations in available energy. Ground cover is also included as an additional insulator of the soil. Furthermore, the modified Berggren equation, which accounts for soil thermal conductivity and soil moisture, is used to convert the frost index into frost depth. The modified CFGI model is tested by application at six test sites within the Sleepers River experimental watershed in Vermont. Compared to the CFGI model, the modified CFGI model more accurately captures the variations in frozen ground between the sites, inter-annual variations in frozen ground depths at a given site, and the occurrence of frozen ground.

1 Introduction

Frozen ground (also known as frozen soil or soil frost) is important to predicting stormflows produced by certain watersheds (Shanley and Chalmers, 1999; McNamara et al., 1997; Prèvost et al., 1990; Woo, 1986). Several plot-scale studies have shown that frozen ground can impede infiltration and thus enhance runoff (Bayard et al., 2005; Dunne and Black, 1971; Stähli et al., 1999). Several of these studies have also shown that frozen ground is highly variable temporally and spatially (Campbell et al., 2010; Shanley and Chalmers, 1999; Stähli, 2017), which affects the amount and type of runoff (Wilcox et al., 1997). The presence, spatial pattern, and depth of frozen ground are driven by mass (water) and energy balances. The energy available from the atmosphere to thaw the soil is subject to the insulation of the snowpack (Pearson, 1920; Willis et al., 1961) and ground cover, including any vegetation, woody debris, and leaf litter (Brown, 1966; Diebold, 1938; Fahey and Lang, 1975; Sartz, 1973; Stähli, 2017). MacKinney (1929) found that ground cover reduced the depth of frost penetration by 40 % at a test site in Connecticut. Additionally, the presence and depth of frozen ground is affected by soil moisture (Fox, 1992; Willis et al., 1961) and the thermal conductivity of the soil (Farouki, 1981; Johansen, 1977).

Frozen ground has proven difficult to simulate within hydrologic models due to complex interactions of energy and water between the atmosphere, snowpack, and soil (Dun et al., 2010; Kennedy and Sharratt, 1998; Lin and McCool, 2006). Physically based models of frozen ground, such as the simultaneous heat and water (SHAW) model (Flerchinger and Saxton, 1989), the coupled heat and mass transfer model for soil–plant–atmosphere systems (COUP) (Jansson 2001; Jansson and Karlburg, 2010), and the distributed water–heat coupled (DWHC) model (Chen et al., 2007) have large parameter and forcing data requirements – such as wind speed, relative humidity, and short- and long-wave radiation – which restricts their applicability in many watershed. Additionally, these types of models either include, or are tightly coupled to soil moisture models, which can limit their applicability in models that do not explicitly simulate soil moisture content. To reduce data and parameter requirements and increase applicability, simple temperature-index or degree-day methods (Molnau and Bissell, 1983; Rekolainen and Posch, 1993) remain widely used within watershed models, including LISFLOOD (De Roo et al., 2001; Van Der Knijff et al., 2010), CREAMS (Rekolainen and Posch, 1993), and the gridded surface subsurface hydrologic analysis (GSSHA) model (Downer and Ogden, 2004). Degree-day approaches typically accumulate the daily average temperature as a frost index (C-days). When the frost index exceeds a threshold, the soil is considered frozen and impermeable to infiltration. The sudden restriction on infiltration can be an incorrect assumption, especially in forested environments where frozen soils often still experience infiltration (Lindstrom et al., 2002; Nyberg et al., 2001; Shanley and Chalmers, 1999). A limitation of degree-day approaches is that they are often untested against observed frost data because the frost index is not a physical property that can be compared to measurements. However, degree-day methods have been successful in capturing increased runoff from frozen ground events (Molnau and Bissell, 1983), and higher frost index values have been shown to correlate to deeper frost depths (Vermette and Christopher, 2008; Vermette and Kanack, 2012). Spatial variations of frozen ground within degree-day methods are typically based on variations in temperature (which are estimated from an elevation–temperature relationship) and variations in snowpack insulation (which are also typically inferred from an elevation–temperature relationship). Such reliance on elevation may lead to errors because Stähli (2017) found no clear connection between elevation and presence of frozen ground at test sites in the Swiss pre-Alpine zone.

The objective of this paper is to develop a model for frozen ground that (1) captures the spatial variations of frozen ground within a watershed, (2) allows the frozen ground model to be incorporated into a variety of watershed models, and (3) allows application in data sparse environments where limited forcing data may prohibit use of energy balance methods. In this paper, we use the GSSHA watershed model and develop the frozen ground model by modifying the commonly used conceptual frozen ground index (CFGI) (Molnau and Bissell, 1983) method in four ways. First, the CFGI method is coupled to an improved snowpack model that more accurately captures the spatial heterogeneity of the snowpack. In past applications of GSSHA, the CFGI method has been coupled with a temperature-index (TI) snowpack model based on SNOW-17 (Anderson, 1973, 2006). However, Follum et al. (2015) proposed a radiation-derived temperature index (RTI) snow model that uses a proxy temperature instead of air temperature to represent the energy available to the snowpack. Compared to the TI model, the RTI model more directly includes the effects of shortwave radiation and canopy cover and was shown to better represent the spatial variations of snow cover and snow water equivalent (SWE) in the Senator Beck Basin in Colorado. The RTI model is adopted to simulate the snowpack in the present study. Second, the effects of shortwave radiation and canopy cover are included in the CFGI model when calculating the energy available at the snow or ground surface. These effects are included by using a similar radiation-derived proxy temperature when calculating the frost index. Third, the insulation effects of ground cover are included by modifying the frost index equation. Fourth, an option is included to compute frost depth as a function of the frost index value. The modified Berggren equation and similar Stefan equation have been previously used to estimate frost depth from degree-days (Carey and Woo, 2005; DeWalle and Rango, 2008; Fox, 1992; Woo et al., 2004); a similar approach is used here to convert the frost index to frost depth.

The following sections first describe the existing TI and CFGI models within GSSHA. The combination of these two models serves as the baseline or control case for the experiments. Then, the RTI snow model and the modified CFGI frozen ground model (referred to as modCFGI) are described. Finally, the results of the TI/CFGI model and RTI/modCFGI models are compared to each other and to observations of snow depth, SWE, and frost depth at the Sleepers River experimental watershed (SREW) in Vermont.

2 Methodology

2.1 TI snowpack model

The TI snow model was implemented into GSSHA by Follum et al. (2014), who provides additional information about the model. Although GSSHA allows a variable time step for multiple processes, it always uses an hourly time step t) for snow calculations. GSSHA utilizes a structured grid in which each cell can have a different air temperature Ta (C) and precipitation P (m h−1). Air temperature is the primary driver of snowpack dynamics in the TI model and is estimated as follows:

(1) T a = T g + E g - E c ,

where Tg (C) is the air temperature at a gage, is a linear lapse rate (C km−1), and Eg and Ec (m) are the elevations of the temperature gage and the grid cell where Ta is being calculated, respectively. Precipitation accumulates as SWE (m) when TaTpx, where Tpx is the freezing point (0C by default). The precipitation P is multiplied by a uniform multiplication factor (Scf), which crudely represents snowpack sublimation and redistribution of snow due to wind (Anderson, 2006). The resultant effective precipitation (Peff) is added to the SWE.

Before the snowpack begins to melt, its heat deficit (or cold content) must be overcome. The change in heat deficit ΔDt (mm of SWE), due to a temperature difference between the snow surface and air, is calculated using the following equation:

(2) Δ D t = N mf , max Δ t / 6 M f / M f , max A TI - T sur ,

where Tsur is the snow surface temperature, and ATI is the antecedent temperature index (C), which is calculated using Ta and the antecedent snow temperature index parameter ATIPM (see Anderson, 2006, for details regarding Tsur and ATI). Nmf,max is the maximum negative melt factor (mm C−1 (6 h)−1), which is a parameter. Mf is the melt factor (mm C-1Δt-1), which is calculated as follows:

(3) M f = Δ t / 6 S v A v M f , max - M f , min + M f , min ,

where Sv and Av are seasonal melt adjustments that change by Julian day, and Mf, max and Mf, min are the maximum and minimum melt factors (mm C−1 (6 h)−1), which are parameters.

Once the heat deficit is overcome, SWE decreases as melt occurs. During normal conditions, the melt M (mm of SWE) is

(4) M = M f T a - T mbase + 0.0125 P eff f r T r Δ t ,

where Tmbase is the temperature at which melt begins (0 C by default), fr is the fraction of any precipitation that is rain (assumed equal to 1 when Ta>0C, otherwise set to 0), and Tr is the precipitation temperature (assumed equal to Ta or 0C, whichever is greater). During rain-on-snow events (more than 1.5 mm of rainfall in the previous 6 h), M is calculated from a simple energy balance:


where σ is the Stefan–Boltzmann constant, fu is the average wind function (mm mb−1 (6 h)−1) (see Anderson, 2006, for details), rh is the relative humidity (assumed to be 0.9 during rain-on-snow events) (Anderson, 1973, 2006), Pa is atmospheric pressure (mb) (either measured or calculated from elevation) (Anderson, 2006), and esat is the saturation vapor pressure (mb) (calculated based on Smith, 1993). The ripeness of the snowpack affects the amount of melt that is released and is controlled by the liquid holding capacity Lhc, which is a specified percentage of the ice in the snowpack (Anderson, 2006).

For frozen ground calculations, the snow depth is needed from the snow model. The snow depth Ds (cm) is found from the SWE and the snowpack density. GSSHA uses the single-layer snow density functions from SNOW-17 (Anderson, 1976, 2006). The density of newly fallen snow ρn (gm cm−3) varies between 0.05 (Ta-15C) and 0.15 (Ta=0C) according to the following equation:

(6) ρ n = 0.05 + 0.0017 T a + 15 1.5 .

Increases in snowpack density ρx from compaction, destructive metamorphism, and melt metamorphism due to the presence of liquid water are calculated as (Koren et al., 1999):

(7) ρ x , t = ρ x , t - 1 e B 2 B 2 e B 1



The variable t is an index for time, W is the ice portion of the snow pack (cm, W=100Sswe,t-1) where Sswe is the snow water equivalent on the ground in m, Ts is the average snow pack temperature (C, calculated based on Anderson, 2006), and ρd is the threshold density above which destructive metamorphism decreases (ρd is set to 0.15 gm cm−1 based on Anderson, 2006). Finally, β=0 if ρx,t-1ρd, and β=1 if ρx,t-1>ρd, c1=0.026 cm−1 h−1, c2=21 cm3 gm−1, c3=0.005 h−1, c4=0.10C−1, and c5=2.0 if there is liquid water in the snowpack and c5=1.0 if there is not (see Anderson, 1976, 2006, for details).

2.2 CFGI frozen ground model

The CFGI model was originally developed as a lumped model for flood forecasting in the Pacific Northwest, but it has been used in distributed models as well (De Roo et al., 2001; Van Der Knijff et al., 2010). The rationale of the CFGI method is that air temperature ultimately controls the ground temperature, but its impact is moderated by the insulating effects of any snowpack. The presence of frozen ground is determined by the frozen ground index F (C-days), which is calculated by

(10) F t = F t - 1 A - T a , d e - 0.4 K s D s ,

where Ta,d is the average daily air temperature (C), A is a daily decay coefficient, and Ks is the snow reduction coefficient (cm−1). The daily decay coefficient (A) controls the persistence of the F values, and Ks controls the insulation from the snowpack. Molnau and Bissell (1983) recommended changing Ks depending on whether Ta,d is above or below freezing (denoted as Ks,Ta>0C and Ks,Ta<0C, respectively).

Higher values of F indicate a higher likelihood that the ground is frozen. Once F exceeds a specified threshold (Fthreshold), the ground is considered frozen and infiltration is restricted. Molnau and Bissell (1983) found the ground to be frozen when F>83 C-days and thawed when F<56C-days. When F is between these values, the ground could be either frozen or thawed. It is worth noting that F does not depend on soil moisture, which is known to affect the initialization and depth of frozen ground (Kurganova et al., 2007; Willis et al., 1961).

2.3 RTI snowpack model

The RTI model makes two modifications to the TI model: (1) it uses a radiation-derived temperature Trad (C) to better describe the available energy, and (2) it estimates spatially varying snowpack sublimation based on solar radiation approximations.

The RTI model replaces Ta in Eqs. (4) and (5) with a radiation-derived proxy temperature Trad (C). In those equations, Ta is used to conceptually represent the energy available to the snowpack. Trad has a similar purpose but is intended to improve the estimation of available energy. Trad is calculated by assuming that the radiation terms dominate the energy balance at the snow surface so that outgoing longwave radiation balances the net incoming shortwave and longwave radiation (Follum et al., 2015). Thus

(11) R LW = R SW , net + R LW ,

where RLW↑ is outgoing longwave radiation, RSW,net is the net incoming shortwave radiation, and RLW↓ is the downwelling longwave radiation. The right side of Eq. (11) represents the energy that is supplied to the snowpack via the atmosphere. RLW↑ (W m−2) is the radiative response of the snowpack to that energy. Using the Stefan–Boltzmann law, RLW↑ can be written in terms of temperature Trad:

(12) T rad = R SW , net + R LW ε snow σ 1 / 4 - 273.15 ,

where εsnow is the emissivity of snow (assumed to be 0.97) and σ is the Stefan–Boltzmann constant.

RSW,net is calculated as follows:

(13) R SW , net = ( 1 - α s ) R SW ,

where αs is the albedo of the snowpack, which is calculated based on the time elapsed since the most recent snowfall and whether melt is occurring (Henneman and Stefan, 1999). RSW↓ is the incident shortwave radiation, which is calculated using the following:

(14) R SW = R SW , 0 φ r φ atm φ c φ v φ s φ t ,

where RSW,0 is the solar constant (Liou, 2002), φr accounts for distance from the Earth to the sun (based on Julian day, TVA, 1972), φatm accounts for atmospheric scattering (based on elevation, Allen et al., 2005), φc accounts for absorption by clouds (based on fractional cloud cover, TVA, 1972), φv accounts for vegetation (set equal to the vegetation transmission coefficient Kv (Bras, 1990), a vegetation-specific parameter ranging from 1.0 for no canopy coverage to 0.0 for complete canopy coverage), φs accounts for the slope/aspect of the terrain (based on latitude, slope, and azimuth angle, Duffie and Beckman, 1980), and φt accounts for topographic shading (based on elevation, azimuth angle, and solar elevation angle).

RLW↓ is calculated from the contributions of the atmosphere (including clouds) and the canopy:


where εa is the air emissivity (0.757 when snow is present based on Bras, 1990), N is the fractional cloud cover, Fc is the fractional canopy cover (estimated from leaf area index LAI following, Liston and Elder, 2006; Pomeroy et al., 2002), εc is the canopy emissivity (assumed equal to 1 following Sicart et al., 2004), and Tcanopy is the canopy temperature (C) which is assumed equal to Ta following DeWalle and Rango (2008).

Because the TI model uses Ta to drive snowpack dynamics, those dynamics are only directly associated with the downwelling longwave radiation from the air, which is a component of RLW↓. Furthermore, the spatial variations in the available energy only depend on the variations of Ta, which are inferred from elevation. Trad in the RTI model considers both RSW,net and RLW↓ and thus accounts for heterogeneity in topographic orientation and shading as well as canopy cover. The TI model partially accounts for seasonal variation in solar radiation and snow albedo by empirically adjusting Mf as shown in Eq. (3). In the RTI model, seasonal variations in solar radiation and snow albedo are included in Trad, so a constant melt factor Mf is used (Follum et al., 2015).

The TI model uses a uniform multiplication factor (Ssf) that is applied to the precipitation to account for sublimation, but sublimation is known to vary spatially (Musselman et al., 2008; Rinehart et al., 2008; Veatch et al., 2009). Most sublimation methods depend on relative humidity and wind speed (e.g., Pomeroy, 1988; Liston and Elder, 2006), which are often unavailable in data sparse environments. However, Gustafson et al. (2010) linked differences in sublimation rates to the amount of solar radiation a location receives. In the RTI model a simple approach is used to estimate hourly sublimation rates Ssub (cm h−1) as follows:

(16) S sub = S sub , d R SW R SW , flat ,

where Ssub, d (cm d−1) is the watershed-average daily maximum sublimation amount (a parameter), and RSW,flat is the daily shortwave radiation for a flat cell within the watershed on a cloud-free day. Thus, locations with higher RSW↓ (e.g., open areas and south-facing slopes in the Northern Hemisphere) will have higher values of Ssub. The method neglects wind speed and relative humidity, but does vary sublimation rates based on spatial patterns of solar radiation.

2.4 modCFGI frozen ground model

The CFGI model is modified in three ways to create the modCFGI model. First, the average daily proxy temperature Trad, d is used in place of Ta,d to represent available energy. Second, ground cover (leaf litter, woody debris, etc.) is included as an insulator in the frozen ground index. And third, an option is included to estimate frost depth based on the frozen ground index. The frost depth calculation is optional because it requires soil moisture estimates and may not be needed in many hydrologic models that only require the occurrence (not depth) of frozen ground.

The CFGI uses Ta,d in Eq. (10) to represent the energy that is available to heat the ground surface. In the modCFGI model, Ta,d is replaced with Trad,d. Trad,d is calculated using αs (see Eqs. 12 and 13) when snow is present, and the albedo of the land cover when snow is not present. By using Trad,d, the modCFGI model is expected to better represent the spatial heterogeneity of energy supply due to variations in the topography and canopy cover within a watershed.

The insulation by the ground cover is included by modifying Eq. (10) which becomes

(17) F t = F t - 1 A - T rad e - 0.4 K s D s + K gc D gc ,

where Kgc is the ground cover reduction coefficient (cm−1) and Dgc is the depth of ground cover (cm). This formulation retains the original form of the CFGI model but includes insulation from both snowpack and ground cover. F can still be used to identify the occurrence of frozen ground, which may be sufficient for many hydrologic models. However, because F is not a measurable quantity, an option to extend modCFGI to calculate frost depth is also needed.

Frost depth is calculated using F and the modified Berggren equation. As originally proposed (and described by DeWalle and Rango, 2008), the Berggren equation relates the number of degree days in the freezing/thawing period U (C-days) to the maximum frost depth Zmax (m) as follows:

(18) Z max = λ 48 U δ - 1 Ω m 1 / 2 ,

where λ is a dimensionless coefficient that accounts for changes in sensible heat of the soil, δ (J m−3) is the latent heat of fusion of the soil, and Ωm (J m−1 h−1 C−1) is the mean thermal conductivity of the frozen and unfrozen soil layers. The derivation and corresponding assumptions (i.e., linear soil temperature gradients, Aldrich Jr., 1956) do not reveal any major impediments to adapting this equation for a shorter time step. In addition, Fox (1992), Woo et al. (2004), and Carey and Woo (2005) have used a layered version of the Stefan equation, which is similar to Eq. (18) to simulate daily frost depths with daily input data. Thus, the modified Berggren equation is applied at a daily timescale and revised to become:

(19) Z d = λ [ 48 F - F threshold δ - 1 Ω m ] 1 / 2 ,

where Zd is the depth of frozen ground (m). By using the difference between F and Fthreshold, the degree-days of the current freezing/thawing period is utilized, which is similar to the use of U in the original equation. Zd is only calculated once the ground begins to freeze (when F>Fthreshold). Zd deepens as F becomes increasingly larger than Fthreshold. When F decreases (due to increasing Trad), so does the thickness of frost depth. No frost remains when F falls below Fthreshold.

For the original modified Berggren equation, λ can be estimated annually from Aldrich Jr. (1956) using U, the mean annual air temperature, and the soil water content ω (% of dry weight). Here, λ is calculated using daily differences between F and Fthreshold, the mean annual air temperature, and daily ω values. Thus, soil moisture is included in the calculation of Zd even though it is not included in the calculation of F. Furthermore, δ is estimated daily as

(20) δ = δ f ρ ω / 100 ,

where δf is the latent heat of fusion of water (0.334 MJ kg−1 at 0C) and ρ is the dry soil density. Ωm is estimated as follows (Farouki, 1981; Johansen, 1977):

(21) Ω m = Ω sat - Ω dry ω + Ω dry ,

where Ωdry and Ωsat are the thermal conductivity of dry and saturated soil, respectively. Ωsat is calculated as the geometric mean of the conductivities of the materials within the soil profile (Farouki, 1981; Johansen, 1977):

(22) Ω sat = Ω s ( 1 - n total ) Ω ice ( n ice ) Ω water ( n total - n ice ) ,

where Ωs, Ωice, and Ωwater are the thermal conductivity of solids, ice, and water, respectively (Farouki, 1981). ntotal is the porosity, and nice is

(23) n ice = n total Z d / H ,

where H (m) is the soil thickness.

Figure 1W-3 sub-basin in the Sleepers River experimental watershed. Sites used in this study are identified with red triangles and blue snowflakes. Basin delineation and elevation contours (m) are based on the 1/3 arcsec National Elevation Dataset, land cover classification is based on the 2006 National Land Cover Database, and sources of the background imagery include ESRI, DigitalGlobe, Earthstar Geographics, CNES/Airbus DS, GeoEye, USDA FSA, USGS, Getmapping, Aerogrid, IGN, IGP, and the GIS User Community.


3 Model application

3.1 Study area

The TI/CFGI and RTI/modCFGI models are tested at the W-3 sub-basin (Fig. 1) of the SREW. The study period is 1 October 2005 through 30 September 2010, which is water year (WY) 2006 through 2010. The SREW was founded in 1958 primarily for studies of snow accumulation, melt, and runoff (Anderson, 1973, 1976; Dunne and Black, 1970a, b; Dunne and Black, 1971; Shanley, 2000; Shanley and Chalmers, 1999). The W-3 sub-basin is located at 4429 N and 7209 W. Elevations range between 348 and 697 m, and the area is approximately 8.5 km2 (based on the National Elevation Dataset, Gesch et al., 2002). The basin is primarily forested with deciduous (57.7 %), evergreen (7.8 %), and mixed (15.3 %) trees (based on the 2006 National Land Cover Database (NLCD), Fry et al., 2011). Approximately 14.6 % of the land cover is pasture/hay and cultivated crops. These open areas are typically below an elevation of 525 m, which is the approximate limit for viable agriculture (Shanley and Chalmers, 1999). The W-3 sub-basin is extensively gaged for both hydrometeorology and hydrology by the US Geological Survey (USGS) and collaborators from federal agencies and universities. Additional basin information and data are provided by Shanley et al. (1995), Shanley and Chalmers (1999), and the USGS website (, last access: 7 November 2016).

Two snow sites and 35 frost sites within W-3 were monitored by the Vermont Field Office of the USGS. At the snow sites, SWE and snow depth were measured approximately weekly, and both sites are used in the present study. At the frost sites, snow depth and frost depth were measured periodically (between 0 and 14 measurements in a given winter). Frost depth was measured using CRREL-Gandahl frost tubes (Ricard et al., 1976), which are filled with a methylene blue solution. The frost depth is identified by a change in color within the tube (blue indicates thawed, clear indicates frozen). Vermette and Kanack (2012) provide images and descriptions of similar frost tubes, and Shanley and Chalmers (1999) provides detailed descriptions of the frost tubes at SREW. The frost sites (labeled FS in Fig. 1) are clustered in six parts of the watershed. For this paper, one site from each cluster (FS4, FS11, FS21, FS24, FS30, and FS40) was selected for analysis. The selected sites are far enough apart to be relatively independent but still capture the variations in elevation and land cover classification within the watershed.

3.2 Model inputs

The TI and CFGI models require hourly precipitation and temperature data, which were obtained from the USGS. Precipitation was measured at the W9 weir and R3 snow site (Fig. 1). The USGS then creates a single spatially averaged precipitation time series by weighting the measurements using the distribution of elevation (based on personal communication with James Shanley of the Vermont Field Office of the USGS on 14 November 2016). The W9 gage receives more weight because the watershed includes elevations both above and below this site. Hourly temperature was measured at the W9 site, which has an elevation of 520 m.

Table 1Allowable ranges and calibrated values for the TI and RTI model parameters using a model-independent parameter estimation and uncertainly analysis (PEST). Dashes indicate parameters that are not required in the associated model. The sensitivity ranking for each parameter is shown in parentheses.

Download Print Version | Download XLSX

The RTI and modCFGI models also require cloud cover data, which were obtained from the National Centers for Environmental Information (NCEI,, last access: 7 November 2016). The hourly cloud cover classification data (clear, few clouds, broken sky, etc.) were collected at the Edward F. Knapp State Airport (44 km southwest of the basin) and the Morrisville-Stowe State Airport (36 km west of the basin). The classification data were converted to cloud cover percentages using the method from Follum et al. (2015). Cloud cover data are routinely measured at most airports in the US (data archived at NCEI) as well as many meteorological stations. For simulation of frost depth (and comparison to frost depth observations), soil moisture and evapotranspiration were also simulated. These two components additionally require hourly relative humidity, wind speed, and atmospheric pressure data, all of which were obtained from a meteorological station at the Fairbanks Museum in Saint Johnsbury, VT (11 km southeast of the basin) with missing values replaced with hourly data from the two airports.

All the models require elevation data to determine the spatial patterns of snow and frozen ground. W-3 was delineated using the 1/3 arcsec ( 9 m) National Elevation Dataset (Gesch et al., 2002). The RTI and modCFGI models additionally require land cover classifications, which were obtained from the 2006 NLCD (Fry et al., 2011) and have a 30 m resolution. The classifications of some grid cells were changed to match the land covers observed in the field. In particular, the grid cell containing R3 was changed from deciduous forest to pasture/hay, FS11 was changed from mixed forest to evergreen forest, and FS21 was changed from developed to mixed forest. Both FS24 and FS30 are classified as pasture/hay, where FS24 is a managed pasture and FS30 is an unmanaged pasture (Ann Chalmers, Vermont Field Office of the USGS, personal communication, 15 November 2016). For example, during field observations in November 2016, FS24 had manure spread throughout the field, while FS30 was not fertilized.

Soil classification data are also required for calculating frost depth, and were obtained from the Digital General Soil Map of the United States (Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, Web Soil Survey, available online at, last access: 10 August 2016). Almost the entire W-3 basin is classified as fine sandy loam. The Watershed Modeling System (Aquaveo, 2013) was used to develop the GSSHA model with a 30 m structured grid. This resolution is adequate to capture the spatial heterogeneity of the basin while remaining computationally efficient.

3.3 Parameter estimation and calibration

The model-independent parameter estimation and uncertainly analysis (PEST) method (Doherty et al., 1994) was used to calibrate seven parameters in the TI model and eight parameters in the RTI model. PEST is a nonlinear local search parameter estimator that calibrates numerous parameters simultaneously to produce the best fit between simulated results and observations. WY 2006 and 2007 were used as the calibration period. The TI and RTI snow models were calibrated first to minimize the sum of the squared residuals between simulated and observed snow depths at the eight sites (six frost sites and two snow sites).

Table 1 displays the allowable range, calibrated value, and sensitivity ranking for the calibrated snow parameters. Goodness of fit statistics as well as description of the affects each parameter has on the snow simulations are described in the results and discussion section of this paper. The allowable ranges for ATIPM, fu, Lhc, Nmf, max, Mf , Mf, max, and Mf, min are based on physical limitations and typical ranges in the literature (Follum et al., 2015). LAI can be estimated from seasonal and annual relationships to remotely sensed normalized difference vegetation index (NDVI) values (Wang et al., 2005). However, snowpack affects the measurement of greenness in high latitude regions (Beck et al., 2006). Thus, LAI and Kv values were calibrated based on land cover classifications with forested land covers being categorized as deciduous forest (including deciduous forest, woody wetlands, and mixed forest) or evergreen forest. LAI and Kv values for non-forested land cover classifications were set to 0.0 and 1.0, respectively. Tpx and Tmbase were not calibrated (both set to 0 C) because the temperature data were post-processed by the Vermont USGS and are expected to be accurate. By comparing the temperature measurements at W9 and the Fairbanks Museum (elevation of  212.4 m), was estimated at 6.6 C km−1. All snow density parameters are set based on Anderson (1973, 2006).

The PEST results indicate that the TI model's snow depths are most sensitive to Scf, Mf, max, ATIPM, and Mf,min. For the RTI model, snow depths are most sensitive to Kv (deciduous), ATIPM, LAI(evergreen), and Mf (Table 1). The calibrated deciduous Kv is near the top of the allowable range (1.0) and LAI is near the bottom (0.103), indicating that the snow in the deciduous forest behaves similarly to the open pasture areas where Kv= 1 and LAI= 0.

The CFGI and modCFGI frozen ground models were calibrated to minimize the sum of squared residuals between the simulated and observed frost depths at the six frost sites. For purposes of comparison the modified Berggren equation was also added to the CFGI model to calculate frost depth. Table 2 displays the allowable range, calibrated value, and sensitivity ranking of each calibrated frozen ground parameter. Goodness of fit statistics as well as description of affects each parameter has on the frost depth simulations are described in the results and discussion section. Fthreshold was calibrated for both the CFGI and modCFGI models with the upper range based on Molnau and Bissell (1983). Three Kgc values were calibrated for the modCFGI frozen ground model: one for the managed pasture site FS24 (Kgc,FS24), one for the unmanaged pasture site FS30 (Kgc,FS30), and one for all other frozen ground sites (Kgc).

Table 2Allowable ranges and calibrated values for the CFGI and modCFGI model parameters using PEST. Dashes indicate parameters that are not required in the associated model. The sensitivity ranking for the modCFGI parameters are shown in parentheses.

Download Print Version | Download XLSX

Table 3Values of soil parameters used to calculate soil moisture in the single-layer Green and Ampt infiltration model.

Download Print Version | Download XLSX

Following Molnau and Bissell (1983), multiple combinations of A (0.8 and 0.97), and Ks,Ta<0C and Ks,Ta>0C (0.08, 0.2, and 0.5) values were tested with A= 0.97, Ks,Ta<0C= 0.08, and Ks,Ta>0C= 0.5 producing frost indices that best replicate the rise and fall of the frost depth as well as the timing of the peak frost depth. Depth of ground cover for each land cover type was obtained from field observations in November 2016. Specifically, Dgc=6 cm for deciduous forest (fallen leaves), Dgc=2 cm for evergreen forest (fallen leaves), Dgc=4 cm for pasture (grass), and Dgc=0 cm for all other land cover types.

The modified Berggren equation requires soil moisture, which can be simulated using several methods in GSSHA (Downer and Ogden, 2006). To facilitate extension of these results to other hydrologic models, the commonly used single-layer Green and Ampt infiltration model (Green and Ampt, 1911) with soil moisture redistribution between rainfall events (Ogden and Saghafian, 1997) is utilized to calculate infiltration. Soil moisture is tracked using a simple bucket approach, accounting for infiltration, evapotranspiration, and groundwater recharge as described in Downer (2007). The soil layer thickness (H) is set to 0.5 m for both the soil moisture calculations and frost depth equations. Soil infiltration parameters are set based on published values for the W-3 soil type (Downer and Ogden, 2006; Rawls et al., 1982, 1983; Rawls and Brakensiek, 1985) and are shown in Table 3. Evapotranspiration, which can reduce the soil moisture, is simulated using a Penman–Monteith approach (Monteith, 1965, 1981) with parameters estimated based on land cover (Downer and Ogden, 2006). The dry soil density (ρ= 1137 kg m−3) and dry soil thermal conductivity (Ωdry= 792 J m−1 h−1C−1) are set based on measurements of fine sandy loam by Nikolaev et al. (2013).

For the CFGI model, the calibrated Fthreshold value (Table 2) is relatively close to the lower bound value of 56 C-days found in Molnau and Bissell (1983). For the modCFGI model, the calibrated Fthreshold value is at the lower bound. The Fthreshold value is expected to be lower for the modCFGI model than the CFGI model. The modCFGI model incorporates the insulation by ground cover directly using Kgc and Dgc, whereas the CFGI model can only account for those effects by adjusting the Fthreshold value. It is also worth noting that Kgc,FS30 has a very low value (minimum of allowable range), which suggests that insulation from grass in an unmanaged pasture is very small. This could be the result of snow falling within the grass of the unmanaged pasture, thus making any insulating contribution from the grass very small.

4 Results and discussion

4.1 Snow depth and SWE (TI vs. RTI)

Figure 2 shows maps of simulated snow depth on 23 February 2007 from the TI and RTI snow models. The spatial variability in the TI snowpack is entirely based on elevation (due to the inference of local air temperature from elevation). Higher elevations have deeper snowpack due to lower air temperatures. The RTI snowpack also varies with elevation but shows variation due to land cover as well. In particular, pasture areas have slightly shallower snowpack than surrounding areas due to higher sublimation rates and higher absorbed shortwave radiation. North-facing slopes also have more snow than south-facing slopes due to lower absorbed shortwave radiation. Although no maps of observed snow depth are available for comparison, large-scale distributions of snowpack are known to be controlled by elevation, land cover, and slope/aspect (Fassnacht et al., 2017; Jost et al., 2007), which is more consistent with the RTI model.

Figure 2Simulated maps of snow depth (TI and RTI models) within the W-3 watershed for 23 February 2007. No observed maps of snow depth are available, but the map shows the differences between the temperature-based (TI) model and the modified (RTI) model.


Figure 3TI and RTI simulated snow depth at all eight test sites within the W-3 watershed.


Table 4Statistics for TI and RTI snow depth values at all eight test sites, and statistics for TI and RTI SWE values at the R3 and R25 snow test sites. Values are shown for calibration period (WY 2006–2007), validation period (WY 2008–2010), and overall (WY 2006–2010). RMSE values closer to zero and NSE values closer to one indicate better fit.

Download Print Version | Download XLSX

Figure 3 shows the snow depths from the TI and RTI models at all eight test locations and compares them to the observations. Root mean squared error (RMSE) and Nash–Sutcliffe efficiency (NSE) are shown in Table 4 for the calibration period (WY 2006–2007), validation period (WY 2008–2010), and complete period (WY 2006–2010). The TI and RTI models track closely together at the eight test locations despite differences in the snow depth shown in Fig. 2. Differences between the TI and RTI snowpack at the test sites are small (Fig. 3 and Table 4). The RTI model performs slightly better than the TI model in overall average RMSE (15.69 vs. 15.71 cm), while the TI model performs slightly better in overall average NSE (0.58 vs. 0.53). The observed snow depth is relatively low in WY2008 and 2009 at two of the pasture sites (FS24 and FS30) compared to the other sites. Specifically in WY2008 the small snow depth observations are not captured within either model. The R3 site is also classified as pasture yet has a higher snowpack in WY2008 and 2009. The higher snowpack at this pasture site may be explained by the proximity of R3 to forested areas, which may reduce the wind and help preserve the snowpack. Neither model considers wind effects.

Figure 4Ta and Trad values at all eight test sites within the W-3 watershed between 1 and 15 March 2005.


The snow depths from the two models are similar at each location (Fig. 3) because, on average, the available energy to melt snow (Ta in the TI model and Trad in the RTI model) is similar (Fig. 4). However, the diurnal variation of Trad is typically greater than that of Ta. Trad is derived from a simple radiation balance (i.e., neglecting other terms in the thermal energy balance). Thus, Trad is higher than Ta during the day due to high RSW↓ values, and it is typically lower than Ta at night because RSW↓ reduces to 0 and εa (set to 0.757) in Eq. (15) limits the affect Ta has on RLW↓ and therefore Trad. As shown in Fig. 4, the available energy is also similar between these locations. The elevation difference between the highest and lowest elevation site is approximately 300 m, corresponding to a maximum temperature difference of approximately 2 C between the sites. Also, the test sites are typically located on shallow slopes so topographic aspect has little influence on the energy available to melt the snowpack (i.e., Trad). All land cover classifications except evergreen forest (FS11) have Kv values at or near 1 and LAI values at or near 0, which reduces any variations due to land cover. Trad at FS11 (evergreen forest) is different from the other seven sites because its low Kv value (0.308) reduces RSW,net during the day, and a high LAI value (1.0) increases RLW↓ during day and night.

Figure 5 shows the simulated (both TI and RTI models) and observed SWE values, and Table 4 shows the associated performance metrics at the R3 and R25 snow sites. The TI and RTI models are only calibrated to snow depth, but SWE is calculated first and then combined with snow density to determine snow depth. Both models use the same method to calculate snow density. Both models also exhibit similar behavior and performance at the two sites, which is consistent with their similar snow depths discussed earlier (Fig. 3 and Table 4). Overall, this suggests that the snow density equations used within GSSHA are relatively accurate at the W-3 watershed. Thus, accurate estimates of snow depth typically correspond to accurate estimates of SWE as well.

Figure 5TI and RTI simulated SWE at R3 and R25 snow sites within the W-3 watershed.


4.2 Frost depth (CFGI vs. modCFGI)

Figure 6 shows simulated frost depth maps for 23 February 2007 using the CFGI and modCFGI models (no maps of observed frost depths are available for comparison). In the CFGI model, the frost depths mainly depend on elevation. Colder temperatures at higher elevations generally result in greater snowpack, which insulates the ground and produces smaller frost depths. However, at the beginning of the snow season when the snowpack is shallow, low temperatures at high elevations create deep frost in the higher elevations of the watershed. Later, deeper snowpack at high elevations insulate the ground, while the frost depth increases at lower elevations. This reversal in the elevation dependence can produce an inversion (localized minima in frost depth), as seen between the 500 and 650 m contour lines in Fig. 6. The modCFGI frost depth also has some elevation dependence, but the spatial variation mainly follows land cover classification, which is similar to observations of frozen ground in the Swiss pre-Alpine zone (Stähli, 2017). This variation is partly due to the use of Trad and the increased heterogeneity in the snow depth. The effect of snowpack can be seen by comparing hillslopes with the same land cover but different orientations, such as along the 500 m contour south of FS11. Lower Trad values on northeast-facing slopes result in deeper snowpack than the southwest-facing slopes (Fig. 2). This deeper snowpack produces shallower frost depths on the northeast-facing slopes due to insulation by the snow. However, the spatial pattern of frost depth is more heavily affected by the land cover. Land cover's impact largely occurs through the associated ground cover. This effect can be seen by comparing the deep frost at the unmanaged pasture (near FS30) with the shallower frost depth at the deciduous forest areas near FS4, FS21, and FS40. The low ground cover reduction coefficient at the unmanaged pasture (Kgc,FS30) reduces the insulation from the ground cover, creating deeper frost compared to the deciduous forest areas. The larger than expected role of ground cover in the modCFGI model may occur because ground cover is present during the initiation, deepening, and decrease of frost depth, while the snowpack is much more variable throughout the season.

Figure 6Simulated maps of frost depth (CFGI and modCFGI models) within the W-3 watershed for 23 February 2007. No observed maps of frost depth are available, but the map shows the differences between the temperature-based (CFGI) model and the modified (modCFGI) model.


Figure 7Observed frost depth compared against simulated (CFGI and modCFGI) frost depth at all six selected frozen ground test sites within the W-3 watershed.


Table 5Statistics for CFGI and modCFGI frost depth at all six frost sites. Values are shown for calibration period (WY 2006–2007), validation period (WY 2008–2010), and overall (WY 2006–2010). RMSE values closer to zero and NSE values closer to one indicate better fit. No frost was present at FS4 and FS21 during the validation period, resulting in an inability to calculate NSE. Statistics for a recalibrated modCFGI model without ground cover (labeled as ”modCFGI no gc”) are also shown.

Download Print Version | Download XLSX

Table 6Number of true positive (both simulated and observed data show frost depth), true negative (both simulated and observed data show no frost depth), false positive (simulated data shows frost depth but observed data does not), and false negative (simulated data shows no frost depth but observed data shows frost depth) occurrences during the entire test period. The Accuracy is the sum of the true positive and true negative divided by the total number of observations.

Download Print Version | Download XLSX

Figure 7 shows the frost depths from the CFGI and modCFGI models along with the frost depth observations. The RMSE and NSE values during the calibration, validation, and overall periods are shown in Table 5. The simulated frost depth remains more constant amongst the sites when using the CFGI model, which produces similar maximum frost depths for a given year independent of the land cover. The modCFGI results deviate considerably from the CFGI results, producing greater frost depths at the unmanaged pasture (FS30) and evergreen (FS11) sites and smaller frost depths at the deciduous (FS4, FS21, and FS40) and managed pasture (FS24) sites. These simulated differences between the sites are consistent with the observations. The decreased frost depth in the deciduous forest and managed pasture result from their high measured litter depth (Dgc=6 cm) and high reduction coefficient (Kgc,FS24=1.887 cm−1), respectively. The two pasture sites (FS24 and FS30) differ considerably in the observed frost depth with FS30 consistently having deeper frost. This difference likely occurs because FS24 is managed and FS30 is not. With the exception of the validation period at FS30, the modCFGI model performs better (lower RMSE and higher NSE values) than the CFGI model. The difference in performance is most pronounced at the deciduous sites (FS4, FS21, and FS40) where the average overall NSE value is 11.9 for the CFGI model and 0.20 for the modCFGI model.

In hydrologic models, capturing the presence of frozen ground is important because even shallow frost with high moisture content (concrete frost) has the potential to impede infiltration (Dunne and Black, 1971). Therefore, the ability of the CFGI and modCFGI models to accurately capture the presence of frozen ground is evaluated. Whenever frost observations are available, the simulated frost depths are categorized as follows: true positive (both simulated and observed data show frost), true negative (both simulated and observed data show no frost), false positive (simulated data shows frost but observed data shows no frost), or false negative (simulated data shows no frost but observed data shows frost). Table 6 shows the number of observations in each category for each test site. The table also shows the model accuracy, which is calculated as the percent of the observations that are correctly classified (true positive or true negative). The CFGI and modCFGI models perform similarly in capturing true positives at FS4, FS21, FS24, and FS40, while modCFGI has more true positives at FS11 and FS30. The lower true positives and higher false negatives indicate that the CFGI model tends to underestimate the presence of frozen ground at FS11 and FS30. Overall, both the CFGI and modCFGI models capture most of the frozen ground events, with the modCFGI model performing better than the CFGI model at five sites and worse at one site (FS21). The average accuracy of the modCFGI model is 15.2 % higher than the CFGI model, with the largest increase in accuracy at FS11 (29.8 %).

Figure 8Observed frost depth compared against simulated (modCFGI with and without ground cover included) frost depth at all six selected frozen ground test sites within the W-3 watershed. The modCFGI model without ground cover is labeled as ”modCFGI no gc”.


A simple test is employed to explore the modification that contributes most to the increased accuracy of the modCFGI model. This test removes ground cover from the modCFGI model, recalibrates, and then compares the results to observations. When ground cover is removed, the calibrated Fthreshold value is 83 C-days, which is at the top of the calibration range. This change indicates that ground cover has a large impact on the appropriate value of this threshold. Figure 8 shows the simulated frost depths using the modCFGI model with and without ground cover for each test site. Performance metrics for the modCFGI model with and without ground cover are shown in Table 5. Variability in frost depth between the sites is diminished when ground cover is removed, leading to large errors between simulated and observed frost depth. When ground cover is removed, the frost depth results decrease in accuracy (higher RMSE values and lower NSE values) compared to the complete modCFGI model. The only exception is the overall period at FS30, which is also the only site where the CFGI model outperforms the full modCFGI model. These results suggest that inclusion of ground cover is an important reason why the modCFGI model outperforms the CFGI model.

Figure 9Simulated frost depths from the modCFGI model using simulated soil moisture (θsim), a constant high soil moisture (θhigh), and a constant low soil moisture (θlow) at all six selected frozen ground test sites within the W-3 watershed.


The sensitivity of the modCFGI results to soil moisture is also examined. Soil moisture does not affect the calculation of F, but it is included within the modified Berggren equation (Eqs. 18 and 19) in the calculation of δ (Eq. 20) and Ωm (Eq. 21). Soil moisture was simulated using a single-layer Green and Ampt approach. However, no soil moisture measurements are available at any of the test sites to evaluate the accuracy of the simulated values. Sensitivity of the modCFGI model to volumetric soil moisture is tested by artificially setting the soil moisture to either the residual water content (θlow) or the effective porosity (θhigh), which are the lower and upper bounds for soil moisture values within the model. Figure 9 shows the modeled frost depths from the modCFGI model using θlow, θhigh, and the soil moisture from the Green and Ampt approach (θsim, which is identical to modCFGI in Figs. 7 and 8). Also shown are the observed frost depths for reference only. The frost depth from the θsim case is similar to the frost depth from the θhigh because the simulated soil moisture is usually close to the effective porosity. Frost depth increases when θlow is used, which coincides with other studies (Fox, 1992; Willis et al., 1961). The timing of the frozen ground (when it begins and ends) is identical in all three of the simulations. The consistent timing occurs because soil moisture is not used to calculate F and the same Fthreshold (which controls when frozen ground begins) was used for all three simulations. This result highlights a deficiency in the modeling framework. Specifically, soil moisture should be considered for determining the initiation of frozen ground because wet soils have a higher specific heat capacity and require more energy loss to cool and freeze the soil (Kurganova et al., 2007).

5 Conclusions

The main purpose of this paper was to better estimate the spatial pattern of frozen ground for distributed watershed modeling by modifying an existing degree-day frozen ground model (CFGI), which uses a frost index value to determine whether the ground is frozen or not. The modifications to the CFGI model include (1) use of a radiation-derived temperature index (RTI) snow model instead of a standard temperature-index (TI) snow model, (2) use of a radiation-derived proxy temperature (Trad) instead of air temperature (Ta) in the calculation of the frost index, (3) inclusion of ground cover (litter, debris, grass, etc.) as an insulator of the ground from air temperatures, and (4) an option to use a version of the modified Berggren equation to calculate frost depths based on the frost index values. The CFGI and modCFGI models were tested using the GSSHA hydrologic model over a five-year period within the W-3 watershed, which is part of Sleepers River experimental watershed in Vermont. The model results were compared against snow depth at eight sites, snow water equivalent at two sites, and frost depth at six sites. The primary conclusions of the paper are as follows:

  1. The RTI snow model produces much more complex spatial patterns of snow depth than the TI snow model for the W-3 watershed. The TI model, which is based on SNOW-17 (Anderson, 2006), only produces spatial variation using elevation. The RTI model accounts for elevation, hillslope orientation, canopy shading, and longwave radiation from the canopy through the use of the radiation-derived proxy temperature. It also includes a simple sublimation method based on solar radiation. Thus, its snow depths exhibit spatial heterogeneity based on elevation, slope/aspect, and land cover, all of which are known to affect the large-scale distribution of observed snow depths (Fassnacht et al., 2017; Jost et al., 2007).

  2. Both the RTI model and TI model produce accurate results for the eight snow depth sites at W-3. Two of the eight sites also measure snow water equivalent, where the RTI and TI model also show similarly accurate results. The eight test sites have similar topographic attributes and primarily differ in their land covers, which include pasture, deciduous forest, and evergreen forest. Because the leaves have typically fallen prior to snow accumulation, all but the evergreen site behave similarly in snow accumulation and ablation.

  3. The modCFGI frost model produces more complex spatial patterns of frost depth than the CFGI frost model for the W-3 watershed. The CFGI model uses elevation to infer the spatial variation of air temperature. It additionally uses the TI model for snow depth, which also depends on elevation. Thus, the simulated frost depths at W-3 primarily reflect the watershed elevations. In contrast, the modCFGI model uses the radiation-derived proxy temperature to infer the energy available to heat the ground and the RTI model to simulate snow depth. Furthermore, it accounts for the insulating effects of ground cover (in addition to snowpack), which also depends on the land cover. Thus, the frost depths simulated by the modCFGI model at W-3 depend on the local elevation, hillslope orientation, and land cover, all of which are known to affect the distribution of frozen ground (Fox, 1992; MacKinney, 1929; Wilcox et al., 1997; Willis et al., 1961).

  4. The modCFGI model produces more accurate frost depths than the CFGI for all but one of the six test sites in the W-3 watershed. Overall, the modCFGI model more accurately captures the inter-annual variability in frost depth at a given site and variability of frost depth between sites. Although both the CFGI and modCFGI capture the majority of frozen ground events observed, the modCFGI model has 15.2 % better accuracy in capturing the presence of frozen ground, which is expected to be important for capturing runoff that is produced by frozen ground.

  5. A key reason for the difference in performance between the two frost models is that the modCFGI model includes the insulation of the ground by ground cover while the CFGI model does not. When ground cover is removed from the modCFGI model its results for W-3 are less accurate and the variability in simulated frost depth between the sites is limited. Ground cover is likely important in this watershed because it is relatively thick and is also present at all stages of the winter while snowpack is not.

Overall, the modCFGI model provides improved spatial representation of frozen ground while requiring only cloud cover estimates as additional forcing data (more forcing data may be required if soil moisture is simulated to obtain frost depth). Limited data requirements should make modCFGI well suited for data sparse environments. Hydrologic models often need to account for the presence of frozen ground, which in data-sparse environments often means using simple degree-day approaches that typically vary frozen ground with elevation only (as was shown with the CFGI model). To calculate Trad the modCFGI model does require cloud cover data, which are collected operationally at most airports within the US. If soil moisture is explicitly simulated within the hydrologic model the modCFGI model can also be used with the modified Berggren equation to simulate frost depth, which requires information on soil type and an estimate of the thermal conductivity of the soil.

Five main avenues are available for future research. First, the modCFGI model should be generalized to include the effects of wind (as it relates to the snowpack) and more completely consider the role of soil moisture. Soil moisture is not considered when calculating the frost index, so it does not impact the initiation or duration of frozen ground. This limitation results from using a degree-day approach and may be important in some cases (Kurganova et al., 2007; Willis et al., 1961). Second, the modCFGI model should be tested further. Additional testing should consider other areas where snow and frozen ground are known to affect runoff, such as the upper Midwest region of the US. Additional testing should also better characterize the insulation properties of ground cover under different management scenarios. Third, the calculation of Trad is simple and applicable in data-sparse environments, but other approaches for adjusting a temperature value based on topography and land cover are available (Fox, 1992; Kang, 2005; Webster et al., 2017) and could be further tested. Fourth, future research should also determine the effects of spatial heterogeneity of snow and frost depth on runoff and streamflow at both the local and watershed scales. Similar to Campbell et al. (2010), the RTI and modCFGI models could be used in data-sparse watersheds to investigate how changes in historic and future climate affect snow, frozen ground, and runoff. Finally, although this paper focuses on the simulation of frost depth in the context of watershed modeling, the methods described could also be used for agriculture, overland mobility modeling, and infrastructure where snow and frost depth are major concerns.

Data availability

Hourly precipitation and temperature data used in this study are available by contacting the Vermont Field Office of the USGS (; USGS, 2018a). Cloud cover, relative humidity, wind speed, and atmospheric pressure data used in this study are available from the National Centers for Environmental Information (NCEI, 2018;; last access: 7 November 2016). Land cover data (2006 National Land Cover Database; Fry et al., 2011) and elevation data (National Elevation Dataset; Gesch et al., 2002) used in this study are available from the USGS (; USGS, 2018b). Soil data used in this study were obtained from the Digital General Soil Map of the United States (Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, Web Soil Survey, available online at; last access: 10 August 2016). Snow and frozen-ground data used in this study are available by contacting the Vermont Field Office of the USGS (; USGS, 2018a).

Information regarding the GSSHA model is available at (GSSHA, 2018). Additional information and output data sets are available by contacting the corresponding author.

Competing interests

The authors declare that they have no conflict of interest.


This project was funded by the Flood and Coastal Research Program at the U.S. Army Corps of Engineers, Engineering Research and Development Center, Coastal and Hydraulics Laboratory in Vicksburg, Mississippi. We thank the Vermont Field Office of the USGS (specifically James Shanley and Ann Chalmers) for their diligent efforts to collect, process, and analyze numerous hydrologic data for research. We also thank three anonymous reviewers and the editor who provided helpful comments that substantially improved this article.

Edited by: Markus Weiler
Reviewed by: three anonymous referees


Aldrich Jr., H. P.: Frost penetration below highway and airfield pavements, Highway Research Board Bulletin, 135, 124–149, 1956. 

Allen, R. G., Walter, I. A., Elliott, R., Howell, T., Itenfisu, D., and Jensen, M.: The ASCE Standardized Reference Evapotranspiration Equation, ASCE, New York, 2005. 

Anderson, E. A.: National Weather Service River Forecast System – Snow Accumulation and Ablation Model, Technical Memorandum NWS Hydro-17, November 1973, 217 pp., Silver Spring, Maryland, 1973. 

Anderson, E. A.: A Point Energy and Mass Balance Model of a Snow Cover, U.S. National Oceanic and Atmospheric Administration NOAA Technical Report NWS 9, Silver Spring, Maryland, 1976. 

Anderson, E. A.: Snow Accumulation and Ablation Model – SNOW-17, NWSRFS User Documentation, U.S. National Weather Service, Silver Springs, MD, 2006. 

Aquaveo: Watershed Modeling System Version 9.1.4, Provo, Utah, 2013. 

Bayard, D., Stähli, M., Parriaux, A., and Flühler, H.: The influence of seasonally frozen soil on the snowmelt runoff at two Alpine sites in southern Switzerland, J. Hydrol., 309, 66–84, 2005. 

Beck, P. S., Atzberger, C., Høgda, K. A., Johansen, B., and Skidmore, A. K.: Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI, Remote Sens. Environ., 100, 321–334, 2006. 

Bras, R. L.: Hydrology: An Introduction to Hydrologic Science, Addison-Wesley, Reading, Massachusetts, 643 pp., 1990. 

Brown, J.: Influence of vegetation on permafrost, in: Permafrost International Conference, edited by: Woods, K. B., National Academy of Sciences-National Research Council, NRC Div. Bldg., Ottawa, Canada, 6 pp., 1966. 

Campbell, J. L., Ollinger, S. V., Flerchinger, G. N., Wicklein, H., Hayhoe, K., and Bailey, A. S.: Past and projected future changes in snowpack and soil frost at the Hubbard Brook Experiment Forest, New Hampshire, USA, Hydrol. Process., 24, 2465–2480, 2010. 

Carey, S. K. and Woo, M.-K.: Freezing of subarctic hillslopes, Wolf Creek Basin, Yukon, Canada, Arct. Antarct. Alp. Res., 37, 1–10, 2005. 

Chen, R. S., Kang, E. S., Ji, X. B., Yang, Y., Zhang, Z. H., Qing, W. W., Bai, S. Y., Wang, L. D., Kong, Q. Z, Lei, Y. H., and Pei, Z. X.: Preliminary study of the hydrological processes in the alpine meadow and permafrost regions at the headwaters of Heihe River, Journal of Glaciology and Geocryology, 29, 387–396, 2007. 

De Roo, A., Odijk, M., Schmuck, G., Koster, E., and Lucieer, A.: Assessing the effects of land use changes on floods in the Meuse and Oder catchment, Phys. Chem. Earth Pt. B, 26, 593–599, 2001. 

DeWalle, D. R. and Rango, A.: Principles of Snow Hydrology, Cambridge University Press, Cambridge, England, 2008. 

Diebold, C. H.: The effects of vegetation upon snow cover and frost penetration during the March 1936 floods, J. Forest., 36, 1131–1137, 1938. 

Doherty, J., Brebber, L., and Whyte, P.: PEST: Model-Independent Parameter Estimation, Watermark Computing, Corinda, Australia, 122, 1994. 

Downer, C. W.: Development of a Simple Soil Moisture Model in the Hydrologic Simulator GSSHA, ERDC-TN-SWWRP-07-8, U.S. Army Engineer Research and Development Center, Vicksburg, Mississippi, 2007. 

Downer, C. W. and Ogden, F. L.: GSSHA: Model to simulate diverse stream flow producing processes, J. Hydrol. Eng., 9, 161–174, 2004. 

Downer, C. W. and Ogden, F. L.: GSSHA Users' Manual. in: SR-06-1, U.S. Army Engineer Research and Development Center, Vicksburg, Mississippi, 2006. 

Duffie, J. A. and Beckman, W. A.: Solar Engineering of Thermal Processes, 3, Wiley, New York, 1980. 

Dun, S., Wu, J., McCool, D., Frankenberger, J., and Flanagan, D.: Improving frost-simulation subroutines of the Water Erosion Prediction Project (WEPP) model, T. ASABE, 53, 1399–1411, 2010. 

Dunne, T. and Black, R. D.: Runoff processes during snowmelt, Water Resour. Res., 7, 1160–1172, 1971. 

Dunne, T. and Black, R. D.: An experimental investigation of runoff production in permeable soils, Water Resour. Res., 6, 478–490, 1970a. 

Dunne, T. and Black, R. D.: Partial area contributions to storm runoff in a small New England watershed, Water Resour. Res., 6, 1296–1311, 1970b. 

Fahey, T. J. and Lang, G. E.: Concrete frost along an elevational gradient in New Hampshire, Can. J. Forest Res., 5, 700–705, 1975. 

Farouki, O. T.: The thermal properties of soils in cold regions, Cold Reg. Sci. Technol., 5, 67–75, 1981. 

Fassnacht, S. R., López-Moreno, J. I., Ma, C., Weber, A. N., Pfohl, A. K. D., Kampf, S. K., and Kappas, M.: Spatio-temporal snowmelt variability across the headwaters of the Southern Rocky Mountains, Front. Earth Sci., 11, 1–10, 2017. 

Flerchinger, G. and Saxton, K. E.: Simultaneous heat and water model of a freezing snow-residue-soil system I, Theory and development, T. ASAE, 32, 565–0571, 1989. 

Follum, M. L., Downer, C. W., and Niemann, J. D.: Simulating the spatial distribution of snow pack and snow melt runoff with different snow melt algorithms in a physics based watershed model, AGU Fall Meeting, American Geophysical Union, San Francisco, CA, 2014. 

Follum, M. L., Downer, C. W., Niemann, J. D., Roylance, S. M., and Vuyovich, C. M.: A radiation-derived temperature-index snow routine for the GSSHA hydrologic model, J. Hydrol., 529, 723–736, 2015. 

Fox, J. D.: Incorporating freeze-thaw calculations into a water balance model, Water Resour. Res., 28, 2229–2244, 1992. 

Fry, J. A., Xian, G., Jin, S., Dewitz, J. A., Homer, C. G., Limin, Y., Barnews, C. A., Helold, N. D., and Wickham, J. D.: Completion of the 2006 national land cover database for the conterminous United States, Photogramm. Eng. Rem. S., 77, 858–864, 2011. 

Gesch, D., Oimoen, M., Greenlee, S., Nelson, C., Steuck, M., and Tyler, D.: The national elevation dataset, Photogramm. Eng. Rem. S., 68, 5–32, 2002. 

Green, W. and Ampt, G.: Studies on Soil Physics, J. Agr. Sci., 4, 1–24, 1911. 

GSSHA (Gridded Surface Subsurface Hydrologic Analysis): available at:, (last access: 25 April 2018), 2018. 

Gustafson, J. R., Brooks, P. D., Molotch, N. P., and Veatch, W. C.: Estimating snow sublimation using natural chemical and isotropic tracers across a gradient of solar radiation, Water Resour. Res., 46, 1–14, 2010. 

Henneman, H. E. and Stefan, H. G.: Albedo models for snow and ice on a freshwater lake, Cold Reg. Sci. Technol., 29, 31–48, 1999. 

Jansson, P. E.: Coupled heat and mass transfer model for soil-plant-atmosphere systems, Royal Institute of Technology, Stockholm Sweden, 2001. 

Jansson, P. E. and Karlberg, L.: Coupled heat and mass transfer model for soil-plant-atmosphere systems, Royal Institute of Technology, Stockholm, 454 pp., 2010. 

Johansen, O.: Thermal Conductivity of Soils, CRREL-TL-637, Cold Regions Research and Engineering Laboratory, Hanover, NH, 1977. 

Jost, G., Weiler, M., Gluns, D. R., and Alila, Y.: The influence of forest and topography on snow accumulation and melt at the watershed-scale, J. Hydrol., 347, 101–115, 2007. 

Kang, D. H.: Distributed Snowmelt Modeling with GIS and CASC2D at California Gulch, Colorado, Colorado State University, Fort Collins, CO, 195 pp., 2005. 

Kennedy, I. and Sharratt, B.: Model Comparisons to Simulate Soil Frost Depthwi, Soil Sci., 163, 636–645, 1998. 

Koren, V., Schaake, J., Mitchell, K., Duan, Q. Y., Chen, F., and Baker, J. M.: A parameterization of snowpack and frozen ground intended for NCEP weather and climate models, J. Geophys. Res.-Atmos., 104, 19569–19585, 1999. 

Kurganova, I., Teepe, R., and Loftfield, N.: Influence of freeze-thaw events on carbon dioxide emission from soils at different moisture and land use, Carbon Balance and Management, 2, 2–2, 2007. 

Lin, C. and McCool, D.: Simulating snowmelt and soil frost depth by an energy budget approach, T. ASABE, 49, 1383–1394, 2006. 

Lindstrom, G., Bishop, K., and Lofvenius, M. O.: Soil frost and runoff at Svartberget, northern Sweden – Measurements and model analysis, Hydrol. Process., 16, 3379–3392, 2002. 

Liou, K. N.: An Introduction to Atmospheric Radiation, second edition, Academic Press, New York, 2002. 

Liston, G. E. and Elder, K.: A distributed snow-evolution modeling system (SnowModel), J. Hydrometeorol., 7, 1259–1276, 2006. 

MacKinney, A.: Effects of forest litter on soil temperature and soil freezing in autumn and winter, Ecology, 10, 312–321, 1929. 

McNamara, J. P., Kane, D. P., and Hinzman, L. D.: Hydrograph separation in an Arctic watershed using mizing model and graphical techniques, Water Resour. Res., 33, 1707–1719, 1997. 

Molnau, M. and Bissell, V. C.: A continuous frozen ground index for flood forecasting, Proceedings 51st Annual Meeting Western Snow Conference, Canadian Water Resources Association, Cambridge, Ontario, 109–119, 1983. 

Monteith, J.: Evaporation and environment, Symposia of the Society for Experimental Biology, 4, Cambridge, England, 1965. 

Monteith, J.: Evaporation and surface temperature, Q. J. Roy. Meteor. Soc., 107, 1–27, 1981. 

Musselman, K., Molotch, N. P., and Brooks, P. D.: Effects of vegetation on snow accumulation and ablation in a mid-latitude sub-alpine forest, Hydrol. Process., 22, 2767–2776, 2008. 

NCEI (National Centers for Environmental Information): Cloud cover, relative humidity, wind speed, and atmospheric pressure data, available at:, (last access: 7 November 2016), 2018. 

Nikolaev, I. V., Leong, W. H., and Rosen, M. A.: Experimental investigation of soil thermal conductivity over a wide temperature range, Int. J. Thermophys., 34, 1110–1129, 2013. 

Nyberg, L., Stähli, M., Mellander, P. E., and Bishop, K.: Soil frost effects on soil water and runoff dynamics along a boreal forest transact: 1. Field investigations, Hydrol. Process., 15, 909–926, 2001. 

Ogden, F. L. and Saghafian, B.: Green and Ampt infiltration with redistribution, J. Irrig. Drain. E., 123, 386–393, 1997. 

Pearson, G. A.: Factors controlling the distribution of forest types, Part I, Ecology, 1, 139–159, 1920. 

Pomeroy, J.: Wind Transport of Snow, PhD Thesis, University of Saskatchewan, 1988. 

Pomeroy, J., Gray, D., Hedstrom, N., and Janowicz, J.: Prediction of seasonal snow accumulation in cold climate forests, Hydrol. Process., 16, 3543–3558, 2002. 

Prèvost, M., Barry, R., Stein, J., and Plamondon, A. P.: Snowmelt runoff modelling in a balsam fir forest with variable source are simulator (VSAS2), Water Resour. Res., 25, 1067–1077, 1990. 

Rawls, W. J. and Brakensiek, D.: Prediction of soil water properties for hydrologic modelling, in: Watershed Management in the Eighties, edited by: Ward, T. J., ASCE, 293–299, 1985. 

Rawls, W., Brakensiek, D., and Saxton, K.: Estimation of soil water properties, T. ASAE, 25, 1316–1320, 1982. 

Rawls, W. J., Brakensiek, D. L., and Miller, N.: Green-Ampt infiltration parameters from soils data, J. Hydraul. Eng., 109, 62–70, 1983. 

Rekolainen, S. and Posch, M.: Adapting the CREAMS model for Finnish conditions, Hydrol. Res., 24, 309–322, 1993. 

Ricard, J. A., Tobiasson, W., and Greatorex, A.: The field assembled frost gage, United States Army Cold Regions Research and Engineering Laboratory, Hanover, New Hampshire, 1976. 

Rinehart, A. J., Vivoni, E. R., and Brooks, P. D.: Effects of vegetation, albedo, and solar radiation sheltering on the distribution of snow in the Valles Caldera, New Mexico, Ecohydrology, 1, 253–270, 2008. 

Sartz, R. S.: Effect of Forest Cover Removal on Depth of Soil Freezing and Overland Flow, Soil Sci. Soc Am. J., 37, 774–777, 1973. 

Shanley, J. B.: Sleepers River, Vermont: A water, energy, and biogeochemical budgets program site, Fact Sheet-166-99, U.S. Department of the Interior, U.S. Geological Survey, 2000. 

Shanley, J. B. and Chalmers, A.: The effect of frozen soil on snowmelt runoff at Sleepers River, Vermont, Hydrol. Process., 13, 1843–1857, 1999. 

Shanley, J. B., Sundquist, E. T., and Kendall, C.: Water, energy, and biogeochemical budget research at Sleepers River Research Watershed, Vermont, 2331–1258, U.S. Geological Survey, Open-File Reports Section [distributor], 1995. 

Sicart, J. E., Pomeroy, J. W., Essery, R. L. H, Hardy, J., Link, T., and Marks, D.: A sensitivity study of daytime net radiation during snowmelt to forest canopy and atmospheric conditions, J. Hydrometeorol., 5, 774–784, 2004. 

Smith, J. A.: Precipitation, in: Handbook of Hydrology, edited by: Maidment, D. R., McGraw-Hill Inc., New York, 3.4, 1993.  

Stähli, M.: Hydrological significance of soil frost for pre-alpine areas, J. Hydrol., 546, 90–102, 2017. 

Stähli, M., Jansson, P.-E., and Lundin, L.-C.: Soil moisture redistribution and infiltration in frozen sandy soils, Water Resour. Res., 35, 95–103, 1999. 

USGS: Sleepers River Research Watershed, Danville, Vermont, Vermont Field Office, available at:, (last access: 23 April 2018), 2018a. 

USGS: The National Map, available at:, (last access: 23 April 2018), 2018b. 

TVA: Heat and Mass Transfer Between a Water Surface and the Atmosphere, Tennessee Valley Authority, Norris, TN, 1972. 

Van Der Knijff, J., Younis, J., and De Roo, A.: LISFLOOD: a GIS-based distributed model for river basin scale water balance and flood simulation, Int. J. Geogr. Inf. Sci., 24, 189–212, 2010. 

Veatch, W., Brooks, P. D., Gustafson, J. R., and Molotch, N. P.: Quantifying the effects of forest canopy cover on net snow accumulation at a continental, mid-latitude site, Ecohydrology, 2, 115–128, 2009. 

Vermette, S. and Christopher, S.: Using the rate of accumulated freezing and thawing degree days as a surrogate for determining freezing depth in a temperate forest soil, Middle States Geographer, 41, 68–73, 2008. 

Vermette, S. and Kanack, J.: Modeling frost line soil penetration using freezing degree-day rates, day length, and sun angle, 69th Eastern Snow Conference, Frost Valley YMCA, Claryville, New York, USA, 2012. 

Wang, Q., Adiku, S., Tenhunen, J., and Granier, A.: On the relationship of NDVI with leaf area index in a deciduous forest site, Remote Sens. Environ., 94, 244–255, 2005. 

Webster, C., Rutter, N., and Jonas, T.: Improving representation of canopy temperatures for modeling subcanopy incoming longwave radiation to the snow surface, J. Geophys. Res.-Atmos., 122, 9154–9172, 2017. 

Wilcox, B. P., Newman, B. D., Brandes, D., Davenport, D. W., and Reid, K.: Runoff from a semiarid ponderosa pine hillslope in New Mexico, Water Resour. Res., 33, 2301–2314, 1997. 

Willis, W., Carlson, C., Alessi, J., and Haas, H.: Depth of freezing and spring run-off as related to fall soil-moisture level, Can. J. Soil Sci. Sci., 41, 115–123, 1961. 

Woo, M.: Permafrost hydrology in North America, Atmos. Ocean., 24, 201–234, 1986. 

Woo, M., Arain, M., Mollinga, M., and Yi, S.: A two-directional freeze and thaw algorithm for hydrologic and land surface modelling, Geophys. Res. Lett., 31, L12501,, 2004. 

Short summary
Spatial patterns of snow and frozen ground within watersheds can impact the volume and timing of runoff. Commonly used snow and frozen ground simulation methods were modified to better account for the effects of topography and land cover on the spatial patterns of snow and frozen ground. When tested using a watershed in Vermont the modifications resulted in more accurate temporal and spatial simulation of both snow and frozen ground.