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

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 T_a (^∘C) and precipitation P (m h⁻¹). Air temperature is the primary driver of snowpack dynamics in the TI model and is estimated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{T}_{\mathrm{a}}=T_{\mathrm{g}}+\mathrm{\varnothing }\left(E_{\mathrm{g}}-E_{\mathrm{c}}\right),\end{array}$$

where T_g (^∘C) is the air temperature at a gage, ∅ is a linear lapse rate (^∘C km⁻¹), and E_g and E_c (m) are the elevations of the temperature gage and the grid cell where T_a is being calculated, respectively. Precipitation accumulates as SWE (m) when T_a≤T_px, where T_px is the freezing point (0^∘C by default). The precipitation P is multiplied by a uniform multiplication factor (S_cf), which crudely represents snowpack sublimation and redistribution of snow due to wind (Anderson, 2006). The resultant effective precipitation (P_eff) 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 ΔD_t (mm of SWE), due to a temperature difference between the snow surface and air, is calculated using the following equation:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }{D}_{\mathrm{t}}=N_{\mathrm{mf},\max }\left(\mathrm{\Delta }t/\mathrm{6}\right)\left(M_{\mathrm{f}}/M_{\mathrm{f},\max }\right)\left(A_{\mathrm{TI}}-T_{\mathrm{sur}}\right),\end{array}$$

where T_sur is the snow surface temperature, and A_TI is the antecedent temperature index (^∘C), which is calculated using T_a and the antecedent snow temperature index parameter A_TIPM (see Anderson, 2006, for details regarding T_sur and A_TI). N_mf,max is the maximum negative melt factor (mm ^∘C⁻¹ (6 h)⁻¹), which is a parameter. M_f is the melt factor (mm ^∘C${}^{-\mathrm{1}}\mathrm{\Delta }{t}^{-\mathrm{1}})$, which is calculated as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{M}_{\mathrm{f}}=\left(\mathrm{\Delta }t/\mathrm{6}\right)\left[S_{\mathrm{v}}{A}_{\mathrm{v}}\left(M_{\mathrm{f},\max }-M_{\mathrm{f},\min }\right)+M_{\mathrm{f},\min }\right],\end{array}$$

where S_v and A_v are seasonal melt adjustments that change by Julian day, and M_f, max and M_f, min are the maximum and minimum melt factors (mm ^∘C⁻¹ (6 h)⁻¹), which are parameters.

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

$$\begin{array}{}\text{(4)}& {\displaystyle }M=\left[M_{\mathrm{f}}\left(T_{\mathrm{a}}-T_{\mathrm{mbase}}\right)+\mathrm{0.0125}{P}_{\mathrm{eff}}{f}_{\mathrm{r}}{T}_{\mathrm{r}}\right]\mathrm{\Delta }t,\end{array}$$

where T_mbase is the temperature at which melt begins (0 ^∘C by default), f_r is the fraction of any precipitation that is rain (assumed equal to 1 when $T_{\mathrm{a}}>\mathrm{0}{}^{\circ }$C, otherwise set to 0), and T_r is the precipitation temperature (assumed equal to T_a or 0^∘C, 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:

$$\begin{array}{ll}{\displaystyle }M& {\displaystyle }=\mathit{\sigma }\left[{\left(T_{\mathrm{a}}+\mathrm{273}\right)}^{\mathrm{4}}-{\mathrm{273}}^{\mathrm{4}}\right]\mathrm{\Delta }t+\mathrm{0.0125}{P}_{\mathrm{eff}}{f}_{\mathrm{r}}{T}_{\mathrm{r}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }+\mathrm{8.5}{f}_u\left(\mathrm{\Delta }t/\mathrm{6}\right)\left[\left(r_{\mathrm{h}}{e}_{\mathrm{sat}}-\mathrm{6.11}\right)+\mathrm{0.00057}{P}_{\mathrm{a}}{T}_{\mathrm{a}}\right],\end{array}$$

where σ is the Stefan–Boltzmann constant, f_u is the average wind function (mm mb⁻¹ (6 h)⁻¹) (see Anderson, 2006, for details), r_h is the relative humidity (assumed to be 0.9 during rain-on-snow events) (Anderson, 1973, 2006), P_a is atmospheric pressure (mb) (either measured or calculated from elevation) (Anderson, 2006), and e_sat 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 L_hc, 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 D_s (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⁻³) varies between 0.05 ($T_{\mathrm{a}}\le -\mathrm{15}{}^{\circ }$C) and 0.15 ($T_{\mathrm{a}}=\mathrm{0}{}^{\circ }$C) according to the following equation:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{n}}=\mathrm{0.05}+\mathrm{0.0017}{\left(T_{\mathrm{a}}+\mathrm{15}\right)}^{\mathrm{1.5}}.\end{array}$$

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):

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\rho }}_{x,t}={\mathit{\rho }}_{x,t-\mathrm{1}}\left({\displaystyle \frac{e^{B_{\mathrm{2}}}}{B_{\mathrm{2}}}}\right)e^{B_{\mathrm{1}}}\end{array}$$

where,

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{1}}=c_{\mathrm{3}}{c}_{\mathrm{5}}\mathrm{d}t{e}^{c_{\mathrm{4}}{T}_{\mathrm{s}}-c_x\mathit{\beta }\left({\mathit{\rho }}_{x,t-\mathrm{1}}-{\mathit{\rho }}_{\mathrm{d}}\right)},\mathrm{and}\text{(9)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{2}}=W_{t-\mathrm{1}}{c}_{\mathrm{1}}\mathrm{d}t{e}^{\mathrm{0.08}{T}_{\mathrm{s}}-c_{\mathrm{2}}{\mathit{\rho }}_{x,t-\mathrm{1}}}.\end{array}$$

The variable t is an index for time, W is the ice portion of the snow pack (cm, $W=\mathrm{100}{S}_{\mathrm{swe},t-\mathrm{1}})$ where S_swe is the snow water equivalent on the ground in m, T_s 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⁻¹ based on Anderson, 2006). Finally, β=0 if ${\mathit{\rho }}_{x,t-\mathrm{1}}\le {\mathit{\rho }}_{\mathrm{d}}$, and β=1 if ${\mathit{\rho }}_{x,t-\mathrm{1}}>{\mathit{\rho }}_{\mathrm{d}}$, c₁=0.026 cm⁻¹ h⁻¹, c₂=21 cm³ gm⁻¹, c₃=0.005 h⁻¹, c₄=0.10 ^∘C⁻¹, and c₅=2.0 if there is liquid water in the snowpack and c₅=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

$$\begin{array}{}\text{(10)}& {\displaystyle }{F}_{\mathrm{t}}=F_{t-\mathrm{1}}A-T_{\mathrm{a},\mathrm{d}}{e}^{-\mathrm{0.4}{K}_{\mathrm{s}}{D}_{\mathrm{s}}},\end{array}$$

where T_a,d is the average daily air temperature (^∘C), A is a daily decay coefficient, and K_s is the snow reduction coefficient (cm⁻¹). The daily decay coefficient (A) controls the persistence of the F values, and K_s controls the insulation from the snowpack. Molnau and Bissell (1983) recommended changing K_s depending on whether T_a,d is above or below freezing (denoted as $K_{\mathrm{s},{\mathrm{T}}_{\mathrm{a}}>\mathrm{0}{}^{\circ }\mathrm{C}}$ and $K_{\mathrm{s},{\mathrm{T}}_{\mathrm{a}}<\mathrm{0}{}^{\circ }\mathrm{C}}$, respectively).

Higher values of F indicate a higher likelihood that the ground is frozen. Once F exceeds a specified threshold (F_threshold), 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<56 ^∘C-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 T_rad (^∘C) to better describe the available energy, and (2) it estimates spatially varying snowpack sublimation based on solar radiation approximations.

The RTI model replaces T_a in Eqs. (4) and (5) with a radiation-derived proxy temperature T_rad (^∘C). In those equations, T_a is used to conceptually represent the energy available to the snowpack. T_rad has a similar purpose but is intended to improve the estimation of available energy. T_rad 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

$$\begin{array}{}\text{(11)}& {\displaystyle }{R}_{\mathrm{LW}\uparrow }=R_{\mathrm{SW},\mathrm{net}}+R_{\mathrm{LW}\downarrow },\end{array}$$

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

$$\begin{array}{}\text{(12)}& {\displaystyle }{T}_{\mathrm{rad}}={\left({\displaystyle \frac{R_{\mathrm{SW},\mathrm{net}}+R_{\mathrm{LW}\downarrow }}{{\mathit{\epsilon }}_{\mathrm{snow}}\mathit{\sigma }}}\right)}^{\mathrm{1}/\mathrm{4}}-\mathrm{273.15},\end{array}$$

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

R_SW,net is calculated as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{R}_{\mathrm{SW},\mathrm{net}}=(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}})R_{\mathrm{SW}\downarrow },\end{array}$$

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). R_SW↓ is the incident shortwave radiation, which is calculated using the following:

$$\begin{array}{}\text{(14)}& {\displaystyle }{R}_{\mathrm{SW}\downarrow }=R_{\mathrm{SW},\mathrm{0}}{\mathit{\phi }}_{\mathrm{r}}{\mathit{\phi }}_{\mathrm{atm}}{\mathit{\phi }}_{\mathrm{c}}{\mathit{\phi }}_{\mathrm{v}}{\mathit{\phi }}_{\mathrm{s}}{\mathit{\phi }}_{\mathrm{t}},\end{array}$$

where R_SW,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 K_v (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).

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

$$\begin{array}{ll}{\displaystyle }{R}_{\mathrm{LW}\downarrow }& {\displaystyle }=\mathit{\sigma }{\mathit{\epsilon }}_{\mathrm{a}}{\left(T_{\mathrm{a}}+\mathrm{273.15}\right)}^{\mathrm{4}}\left(\mathrm{1.0}+\mathrm{0.17}{N}^{\mathrm{2}}\right)\left(\mathrm{1}-F_{\mathrm{c}}\right)\\ \text{(15)}& {\displaystyle }& {\displaystyle }+F_{\mathrm{c}}\mathit{\sigma }{\mathit{\epsilon }}_{\mathrm{c}}{\left(T_{\mathrm{canopy}}+\mathrm{273.15}\right)}^{\mathrm{4}},\end{array}$$

where ε_a is the air emissivity (0.757 when snow is present based on Bras, 1990), N is the fractional cloud cover, F_c is the fractional canopy cover (estimated from leaf area index L_AI following, Liston and Elder, 2006; Pomeroy et al., 2002), ε_c is the canopy emissivity (assumed equal to 1 following Sicart et al., 2004), and T_canopy is the canopy temperature (^∘C) which is assumed equal to T_a following DeWalle and Rango (2008).

Because the TI model uses T_a to drive snowpack dynamics, those dynamics are only directly associated with the downwelling longwave radiation from the air, which is a component of R_LW↓. Furthermore, the spatial variations in the available energy only depend on the variations of T_a, which are inferred from elevation. T_rad in the RTI model considers both R_SW,net and R_LW↓ 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 M_f as shown in Eq. (3). In the RTI model, seasonal variations in solar radiation and snow albedo are included in T_rad, so a constant melt factor M_f is used (Follum et al., 2015).

The TI model uses a uniform multiplication factor (S_sf) 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 S_sub (cm h⁻¹) as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }{S}_{\mathrm{sub}}=S_{\mathrm{sub},\mathrm{d}}\left({\displaystyle \frac{R_{\mathrm{SW}\downarrow }}{R_{\mathrm{SW}\downarrow ,\mathrm{flat}}}}\right),\end{array}$$

where S_sub, d (cm d⁻¹) is the watershed-average daily maximum sublimation amount (a parameter), and $R_{\mathrm{SW}\downarrow ,\mathrm{flat}}$ is the daily shortwave radiation for a flat cell within the watershed on a cloud-free day. Thus, locations with higher R_SW↓ (e.g., open areas and south-facing slopes in the Northern Hemisphere) will have higher values of S_sub. 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 T_rad, d is used in place of T_a,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 T_a,d in Eq. (10) to represent the energy that is available to heat the ground surface. In the modCFGI model, T_a,d is replaced with T_rad,d. T_rad,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 T_rad,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

$$\begin{array}{}\text{(17)}& {\displaystyle }{F}_{\mathrm{t}}=F_{t-\mathrm{1}}A-T_{\mathrm{rad}}{e}^{-\mathrm{0.4}\left(K_{\mathrm{s}}{D}_{\mathrm{s}}+K_{\mathrm{gc}}{D}_{\mathrm{gc}}\right)},\end{array}$$

where K_gc is the ground cover reduction coefficient (cm⁻¹) and D_gc 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 Z_max (m) as follows:

$$\begin{array}{}\text{(18)}& {\displaystyle }{Z}_{\max }=\mathit{\lambda }{\left(\mathrm{48}U{\mathit{\delta }}^{-\mathrm{1}}{\mathrm{\Omega }}_{\mathrm{m}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where λ is a dimensionless coefficient that accounts for changes in sensible heat of the soil, δ (J m⁻³) is the latent heat of fusion of the soil, and Ω_m (J m⁻¹ h^{−1 ∘}C⁻¹) 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:

$$\begin{array}{}\text{(19)}& {\displaystyle }{Z}_{\mathrm{d}}=\mathit{\lambda }[\mathrm{48}\left(F-F_{\mathrm{threshold}}\right){\mathit{\delta }}^{-\mathrm{1}}{\mathrm{\Omega }}_{\mathrm{m}}{]}^{\mathrm{1}/\mathrm{2}},\end{array}$$

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

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 F_threshold, the mean annual air temperature, and daily ω values. Thus, soil moisture is included in the calculation of Z_d even though it is not included in the calculation of F. Furthermore, δ is estimated daily as

$$\begin{array}{}\text{(20)}& {\displaystyle }\mathit{\delta }={\mathit{\delta }}_{\mathrm{f}}\mathit{\rho }\mathit{\omega }/\mathrm{100},\end{array}$$

where δ_f is the latent heat of fusion of water (0.334 MJ kg⁻¹ at 0^∘C) and ρ is the dry soil density. Ω_m is estimated as follows (Farouki, 1981; Johansen, 1977):

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathrm{\Omega }}_{\mathrm{m}}=\left({\mathrm{\Omega }}_{\mathrm{sat}}-{\mathrm{\Omega }}_{\mathrm{dry}}\right)\mathit{\omega }+{\mathrm{\Omega }}_{\mathrm{dry}},\end{array}$$

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):

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathrm{\Omega }}_{\mathrm{sat}}={\mathrm{\Omega }}_{\mathrm{s}}^{(\mathrm{1}-n_{\mathrm{total}})}{\mathrm{\Omega }}_{\mathrm{ice}}^{\left(n_{\mathrm{ice}}\right)}{\mathrm{\Omega }}_{\mathrm{water}}^{(n_{\mathrm{total}}-n_{\mathrm{ice}})},\end{array}$$

where Ω_s, Ω_ice, and Ω_water are the thermal conductivity of solids, ice, and water, respectively (Farouki, 1981). n_total is the porosity, and n_ice is

$$\begin{array}{}\text{(23)}& {\displaystyle }{n}_{\mathrm{ice}}=n_{\mathrm{total}}{Z}_{\mathrm{d}}/H,\end{array}$$

where H (m) is the soil thickness.

Figure 1 W-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 44^∘29^′ N and 72^∘09^′ W. Elevations range between 348 and 697 m, and the area is approximately 8.5 km² (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 (https://nh.water.usgs.gov/project/sleepers/index.htm, 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 1 Allowable 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.

The RTI and modCFGI models also require cloud cover data, which were obtained from the National Centers for Environmental Information (NCEI, https://www.ncdc.noaa.gov/, 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 http://websoilsurvey.nrcs.usda.gov/, 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 A_TIPM, f_u, L_hc, N_mf, max, M_f , M_f, max, and M_f, min are based on physical limitations and typical ranges in the literature (Follum et al., 2015). L_AI 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, L_AI and K_v 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. L_AI and K_v values for non-forested land cover classifications were set to 0.0 and 1.0, respectively. T_px and T_mbase 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⁻¹. 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 S_cf, M_f, max, A_TIPM, and M_f,min. For the RTI model, snow depths are most sensitive to K_v (deciduous), A_TIPM, L_AI(evergreen), and M_f (Table 1). The calibrated deciduous K_v is near the top of the allowable range (1.0) and L_AI is near the bottom (0.103), indicating that the snow in the deciduous forest behaves similarly to the open pasture areas where K_v= 1 and L_AI= 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. F_threshold was calibrated for both the CFGI and modCFGI models with the upper range based on Molnau and Bissell (1983). Three K_gc values were calibrated for the modCFGI frozen ground model: one for the managed pasture site FS24 (K_gc,FS24), one for the unmanaged pasture site FS30 (K_gc,FS30), and one for all other frozen ground sites (K_gc).

Table 2 Allowable 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.

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

Following Molnau and Bissell (1983), multiple combinations of A (0.8 and 0.97), and $K_{\mathrm{s},T_{\mathrm{a}}<\mathrm{0}{}^{\circ }\mathrm{C}}$ and $K_{\mathrm{s},T_{\mathrm{a}}>\mathrm{0}{}^{\circ }\mathrm{C}}$ (0.08, 0.2, and 0.5) values were tested with A= 0.97, $K_{\mathrm{s},T_{\mathrm{a}}<\mathrm{0}{}^{\circ }\mathrm{C}}=$ 0.08, and $K_{\mathrm{s},T_{\mathrm{a}}>\mathrm{0}{}^{\circ }\mathrm{C}}=$ 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, D_gc=6 cm for deciduous forest (fallen leaves), D_gc=2 cm for evergreen forest (fallen leaves), D_gc=4 cm for pasture (grass), and D_gc=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⁻³) and dry soil thermal conductivity (Ω_dry= 792 J m⁻¹ h^−1 ∘C⁻¹) are set based on measurements of fine sandy loam by Nikolaev et al. (2013).

For the CFGI model, the calibrated F_threshold 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 F_threshold value is at the lower bound. The F_threshold 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 K_gc and D_gc, whereas the CFGI model can only account for those effects by adjusting the F_threshold value. It is also worth noting that K_gc,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 2 Simulated 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 3 TI and RTI simulated snow depth at all eight test sites within the W-3 watershed.

Table 4 Statistics 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.

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 4 T_a and T_rad 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 (T_a in the TI model and T_rad in the RTI model) is similar (Fig. 4). However, the diurnal variation of T_rad is typically greater than that of T_a. T_rad is derived from a simple radiation balance (i.e., neglecting other terms in the thermal energy balance). Thus, T_rad is higher than T_a during the day due to high R_SW↓ values, and it is typically lower than T_a at night because R_SW↓ reduces to 0 and ε_a (set to 0.757) in Eq. (15) limits the affect T_a has on R_LW↓ and therefore T_rad. 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., T_rad). All land cover classifications except evergreen forest (FS11) have K_v values at or near 1 and L_AI values at or near 0, which reduces any variations due to land cover. T_rad at FS11 (evergreen forest) is different from the other seven sites because its low K_v value (0.308) reduces R_SW,net during the day, and a high L_AI value (1.0) increases R_LW↓ 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 5 TI 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 T_rad 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 T_rad 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 (K_gc,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 6 Simulated 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 7 Observed frost depth compared against simulated (CFGI and modCFGI) frost depth at all six selected frozen ground test sites within the W-3 watershed.

Table 5 Statistics 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.

Table 6 Number 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.

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 (D_gc=6 cm) and high reduction coefficient ($K_{\mathrm{gc},\mathrm{FS}\mathrm{24}}=\mathrm{1.887}$ cm⁻¹), 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 8 Observed 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 F_threshold 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 9 Simulated 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 F_threshold (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 (T_rad) instead of air temperature (T_a) 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 T_rad 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 T_rad 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.