Modeling groundwater upwelling as a control on river ice thickness

Historically, many Alaskan villages and towns were established along waterways, which served as transportation corridors throughout most of the year. Owing to the limited road network, frozen river systems are still important transportation networks in rural Alaska. Collaborative partnerships were established with several local (non-academic) Alaskans who have extensive experience traveling on Alaskan rivers in all seasons. Their observations on the Tanana River (Jones et al. 2013; Schneider et al. 2013a, b) mirror those found by Herman-Mercer et al. (2011) who reported that rural Alaskans have observed that thin ice was becoming more common on the Yukon River in recent winters. Such ice thinning represents a significant danger to winter travelers.

One potential mechanism for thinning ice relates to changes in permafrost distribution. Ice-rich permafrost is an effective barrier to water transport and recharge (Burt & Williams 1976: Horiguchi & Miller 1980; Kane & Stein 1983), but warming temperatures have caused significant permafrost degradation across the Arctic (Jorgenson et al. 2001; Hinzman et al. 2005; Romanovsky et al. 2010). Permafrost degradation in discontinuous permafrost regions has been described as a potential mechanism for increased hydrologic connectivity between surface and ground water systems. The impermeability of permafrost to water infiltration prevents groundwater (GW) recharge and forces near-surface water to flow laterally and quickly into streams. With permafrost thawing, there is more infiltration to sub-permafrost GW, which delays the rapid runoff from a watershed and allows for more sustained subsurface discharge throughout the winter. Permafrost degradation has been associated with: lake drainage (Yoshikawa & Hinzman 2003); increased daily summer and winter minimum river flows (Smith et al. 2007); and increases in winter river baseflow (Walvoord & Striegl 2007; St Jacques & Sauchyn 2009; Brabets & Walvoord 2009; Lyon & Destouni 2010). During winter, Alaskan ice-covered rivers are fed entirely by GW. Therefore, increased winter baseflow is likely due to increased GW input as a consequence of permafrost degradation associated with the increased GW input corresponding with permafrost degradation (Walvoord & Striegl 2007; St Jacques & Sauchyn 2009, Brabets & Walvoord 2009, Lyon & Destouni 2010) in discontinuous permafrost regions.

Heat from GW upwelling can degrade river ice from below. Field observations show that GW heat flux exceeds atmospheric heat losses and that dangerously thin ice conditions can be maintained for extended periods despite very cold wintertime air temperatures (Figure 1). In recent decades, there have been several reports that thin ice and open water may be more prevalent on interior Alaskan rivers in winter (Herman-Mercer et al. 2011; Jones et al. 2013; Schneider et al. 2013b), which are consistent with independent data on changes in regional hydrology. These conditions appear to be caused by a combination of factors that include fast-moving turbulent water, warmer air temperatures, and GW upwelling in shallow areas.

The Tanana Flats, which are found to the south of the Tanana River, have a high GW table and have exhibited extensive permafrost degradation since the 1700s (Jorgenson et al. 2001). Much of the Tanana River and its neighboring regions appear to be fed by GW upwelling. The observations of increased thin ice and open water on the Tanana during winter are hypothesized to be caused by increased GW flow caused by permafrost degradation, which has been intensified by a warming climate.

In this paper, fluctuations in GW hydrology were modeled to assess the potential impacts on river ice dynamics. In collaboration with rural Alaskans, changes in hydrology were explored to examine how residents of interior Alaska could be affected by recurring dangerous ice conditions in the Tanana River. Field studies were initiated in locations identified by local residents as having recurring thin ice. Based upon field observations, a numerical model was developed to explore the relationship between seasonal GW flow and ice thickness under changing environmental conditions. The model was designed to address two primary research questions:

1. What physical factors have the greatest influence on seasonal dynamics between river ice thickness and GW upwelling on the Tanana River?

2. How do variations in environmental conditions change the capacity of GW to melt river ice?

Study sites were established in areas known to have dangerously thin ice. For modeling, we use data from Hot Cake Slough which is located within the Bonanza Creek Long-Term Ecological Research (LTER) area on the Tanana River, approximately 15 km southwest of Fairbanks, Alaska, USA (longitude: 148 ° 7′ 22″; latitude: 64 ° 43′ 26″) (Figure 2). This region is characterized by discontinuous permafrost and is bordered by the Tanana Flats region, which is part of the Tanana–Kuskokwim Lowlands (Anderson 1970). The Tanana Flats is an area with a high water table fed by GW originating from glaciers in the Alaska Range found to the south. The Alaska Range glaciers (typically found above 1,525 m) are considered as the primary source of water for the Tanana River (120 m elevation at Hot Cake Slough) (Anderson 1970).

Hot Cake Slough has a persistent open water feature that typically remains ice-free during much of the winter, but often attains a thin ice cover during extended cold (−25 to −40 °C) periods. A time-lapse camera taking daily photographs of the study area clearly showed the dynamic nature of the ice conditions. The water depth (measured from channel floor to the bottom of the ice), channel velocity, and ice-free winter discharge were measured using a wading rod and a Marsh McBirney Flow-mate 2000 electromagnetic flow meter. Channel slope was measured using a level over a distance of 100 m. The local hydraulic conductivity (K, m/day) of the sediment within the open water feature of Hot Cake Slough was measured at a streambed sediment depth of 20 cm using the falling head permeameter method (Todd & Mays 2005). Nested piezometers (installed from 20 to 120 cm below channel thalweg) were installed in the open water feature of Hot Cake Slough and were used to measure the hydraulic gradient (dh/dz, dimensionless). From the measured hydraulic conductivity and hydraulic gradient, the vertical GW flow velocity (Q_GWz, m/day) was calculated using Darcy's Law (Todd & Mays 2005).

Salinity tracers were used to measure the horizontal GW flow velocity (V_GWx, m/day) between shallow (1.6 m) GW wells (Todd & Mays 2005). Using horizontal and vertical flow velocities, the total GW flow (Q_GW, m³/d) into the channel at Hot Cake Slough was calculated per unit area (m²).

The vertical temperature profile from the GW through the water column, ice, snow, and air within the river channel was measured at two different points in Hot Cake Slough (referenced as the upstream and downstream study sites). The temperature profile was monitored using thermistor strings calibrated to 0.0 °C in a continuously flowing ice bath.

The physical system is represented by Figure 3, which illustrates how GW affects river ice thickness under observed field conditions. From Figure 3, a physically based numerical model was developed using MATLAB^® (version 2011b) to explore the heat balance between GW discharge and thermally degraded ice-covered areas. The model estimates the potential ice melt rate associated with GW under uniform and steady-state conditions (wind velocity, snow density, upwelling rate, GW temperature, air temperature, and snow depth). Ice-free circumstances are not considered in this model. Symbol definitions and units are summarized in the Appendix (available online at http://www.iwaponline.com/nh/046/026.pdf).

Unsteady flow was modeled similarly by examining how continuous GW inputs would affect ice thickness in a downstream direction in a channel under field conditions. Simulating movement in a downstream direction, water depth of the volumetric cell was estimated, assuming subcritical flow based upon the incoming water depth and the volumetric flow of GW into the cell. Flow velocity was calculated using Manning's Equation for SI units (U_w = R^2/3S^1/2n⁻¹, m⋅s⁻¹) (Sturm 2001) using the hydraulic radius (R), the calculated water surface slope (S), and an estimated composite Manning's coefficient (n = 0.03, dimensionless), which accounts for ice cover and channel roughness.

The initial air temperature, ice thickness, snow depth, and wind velocity were specified. Freezing air temperatures and snow depth could change over time, but wind velocity was assumed constant. The heat transfer coefficient was calculated from the average measured water column temperature (T_w = 1.4 °C), while GW temperature was assumed 3.2 °C based on field measurements.

To examine the steady-state model at equilibrium conditions, theoretical scenarios were used to examine the system when the total energy gained from GW equaled the energy transferred to the ice, which equaled the energy lost to the atmosphere (H_GW = H_wi = H_a). In the first assessment, H_GW input was based upon field measurements (T_GW = 3.2 °C and Q_GW = 0.85 m⋅day⁻¹). H_wi was adjusted by modifying T_w until H_wi was equal to H_GW. H_a was adjusted by modifying T_a until H_a was equal to H_GW. In the second assessment, T_a was assumed (−23 °C) based upon field observations when open water was just beginning to freeze (Δη_ice=0.001 m) and H_a was calculated. H_wi was adjusted by modifying T_w until H_wi was equal to H_a. H_GW was adjusted by modifying Q_GW until H_GW was equal to H_a. In a third assessment, the GW and surface water temperatures were assumed (T_w = 1.4 °C and T_GW = 3.2 °C), while T_a and U_GW were adjusted until H_GW = H_wi = H_a.

Analyses were performed to estimate the sensitivity of the model to changes in the dominant input variables. Specific initial conditions were used as model input for the standard analyses (Table 1). Each input variable was increased and decreased by 30% to determine the resulting sensitivity of the ice melt calculation. The 30% change in each input variable was selected as a reasonable, yet consistent proportional modification to assess the relative response of the model.

A confidence analysis was performed to estimate model performance relative to the estimated measurement accuracy of each input variable (Table 2). Modified input values were varied by increasing and decreasing the initial value by its estimated accuracy (%). All of the modified input values that would maximize (or minimize) the ice melt rate were used simultaneously to determine how uncertainty of the input variables affected model output.

The water depth showed little variation from December through early April each year (varying less than 1 cm until mid-April). The average water velocity of the upstream channel cross-section was 0.02 m⋅s⁻¹. The channel slope was measured to be less than 0.001 and when calculated from the Manning's Equation for the observed conditions (with Manning's n =0.03), the slope was calculated to be less than 0.00001. The vertical hydraulic conductivity (K) was measured to be 4.0 m⋅day⁻¹. The vertical hydraulic gradient varied with the depth of each piezometer. The 1.90 m deep piezometer had a consistently lower flow rate than the 1.10 m deep piezometer. For April 2013, the average piezometric head at 1.10 m depth was 0.085 m (s.d. = 0.022) across a vertical distance of 0.44 m (Figure 4) and the calculated vertical GW flow rate averaged 0.77 m⋅day⁻¹ (s.d. = 0.20) (Figure 4). One outlier of piezometric head was omitted from the analysis. On two different days, horizontal GW velocities were measured to be 0.5 and 0.3 m⋅day⁻¹ at a depth of 1.6 m, which corroborates measurements of 0.3 m⋅day⁻¹ reported by Clilverd et al. (2008) in similar substrates from nearby LTER sites. With these GW flow rates and a channel depth of 20 cm, the total average GW velocity into the channel was calculated to be 0.86 m⋅day⁻¹.

In March 2013, the vertical temperature profile measured nearly constant GW temperatures of 2.7 °C at a depth 0.5 m below the thalweg at the upstream site in Hot Cake Slough (Figure 5), while the downstream site in Hot Cake Slough had a GW temperature of 3.2 °C (Figure 6). At the upstream site, the water column temperature averaged 1.4 °C (Figure 7(a)), while the downstream site water column averaged 1.7 °C (Figure 7(b)).

Time-lapse photos were used to observe site conditions on a daily basis (Figure 1). Using the photo series, changes in site conditions were observed over three different winter seasons. Notably, despite much colder air temperatures in the winter 2011–2012, there was more open water in Hot Cake Slough than in the winter of 2012–2013. Data were insufficient from the 2010–2011 season to compare the period of ice cover to 2012–2013.

According to Equation (8), the potential ice melt rises with increasing air temperature, water column temperature, initial ice thickness, initial snow depth, or water velocity, but falls with rising water depths. The water column temperature increases with increased GW upwelling rates and GW temperatures, both of which increase the GW heat flux (Equation (3)). At a GW upwelling temperature of 3.2 °C, the measured GW heat flux (133 W⋅m⁻²) would melt up to 3.8 cm of ice per day under ideal conditions with no heat lost to the atmosphere or gained by the water column.

Snow-free ice that was 10 cm thick was modeled to melt at the measured GW upwelling rate (0.86 m⋅day⁻¹) when the average water temperature was 1.4 °C and air temperatures were −14 °C or greater (Figure 8(a)). However, the presence of snow radically changes these relationships (Figure 8(b)). At the temperatures down to −40 °C (while including conductive and convective heat losses), the ice would melt with only 5 cm of snow. Moreover, the variability in ice melt rates associated with colder air temperatures is tempered by increasing snow cover (Figure 8(b)) and ice melt is largely independent of air temperatures when snow depths are greater than 10 cm (measured snow density = 395 kg⋅m⁻³). This depth represents the ‘hiemal threshold’ (Pruitt Jr 1957), which varies according to snow density, but is defined as the snow depth at which the physical and biological systems under the snow are thermally isolated from oscillations in air temperature.

At a snow depth of 1 cm, the air temperature is directly proportional to the potential ice melt, while increasing GW flow results in a non-linear increase in the ice melt rate (Figure 8(a)). Similarly, the complexity of the response surface of the potential ice melt rate over a range of snow depths and air temperatures (Figure 8(b)) show that the importance of air temperature is attenuated with increasing snow depth. It also illustrates the linear response of ice melt relative to changes in air temperature and the steep response gradient at low snow depths.

In the assessment where GW measurements were assumed accurate (U_GW = 0.86 m⋅day⁻¹; T_GW = 3.2 °C) in a system at thermal equilibrium, H_GW = 133 W⋅m⁻². To be in equilibrium under these conditions (η_ice =2 cm and η_snow = 1 cm), T_w would be 1.4 °C at T_a of −19 °C. In an assessment assuming water column temperatures were accurate (mean T_w = 1.4 °C, mean T_GW = 3.2 °C, η_ice = 0.2 cm; η_snow = 1 cm), H_ice =136 W⋅m⁻². T_a would be −19 °C and U_GW would be 0.87 m⋅day⁻¹ at thermal equilibrium. In the assessment where field observations of T_a were used while water was starting to freeze (T_a = −23.0 °C; η_ice = 0.1 cm; η_snow = 0 cm), H_a =230 W⋅m⁻². At thermal equilibrium under these conditions, U_GW would be 1.48 m⋅day⁻¹ with T_GW = 3.2 °C and T_w would be 2.4 °C.

Table 1 shows that the calculated ice melt rate was most sensitive to increases in the water column temperature, water velocity, air temperature, water depth, and snow depth, respectively. The sensitivity of each variable is dependent on the initial conditions. At the conditions specified (Table 1), a 30% change in the water temperature, water velocity, air temperature, water depth, or snow depth would affect the ice melt calculations by more than 10% (Table 1).

A confidence analysis was performed to gain an understanding of model performance based upon the measurement confidence of the input variables. Under the assumption that each input variable has errors associated with its measurement, confidence in model performance would decrease with increasing measurement errors. Table 2 suggests that if the estimated minimum and maximum errors of all input variables were realized simultaneously, the modeled ice melt rate could range from ±16% of the initial output value of 18 mm⋅day⁻¹ (e.g., 15–20 mm⋅day⁻¹).

Degrading permafrost, warming air temperatures, forest fires, or changes in precipitation are factors that could affect GW flow rates. Walvoord & Striegel (2007) hypothesized that permafrost degradation has caused the 20% increase in observed winter GW discharge in the Tanana River between 1963 and 2005. Increased air temperatures associated with climate have increased permafrost temperatures and can transform discontinuous permafrost zones into areas without permafrost (Rowland et al. 2010; Romanovsky et al. 2010). These factors increase water exchange rates between surface and ground water (Hinzman et al. 2013) and in areas with a positive hydraulic gradient (GW pressure), upwelling rates will increase.

Field observations suggested that GW upwelling decreased river ice thickness under winter conditions in interior Alaska (Jones et al. 2013; Jones 2014). Increased snow accumulation would increase the thermal insulation of the ice, decrease the heat losses from the water column to the air, and leave more heat energy available for degrading ice. While GW flow into the channel maintained a 2.7 to 3.2 °C temperature gradient between the channel bottom and the bottom of the ice, the thermal gradient between the water column and the air temperature was tempered by the constant temperature at the bottom of the ice. Modeling results suggest that GW temperatures played a minor role in melting the ice relative to other parameters (Table 1), contrary to field observations. Water column temperature and flow velocity were modeled to be the most important drivers of ice degradation. Water column temperature and velocity are both affected by GW input, but the GW measurements suggest that only 0.25% of channel flow is attributable to GW for each meter of the flow path. The cumulative downstream effect of GW flow on an ice-covered channel was not modeled, but could be a future extension of the research.

Groundwater velocities are inherently difficult to measure during winter field conditions. The modeling results suggest that GW flows were not measured with sufficient spatial or temporal resolution to capture their variability. Whereas GW flow was measured, preferential GW flow paths may exist that were not sampled during field measurements. The locations of maximum vertical GW flux were sought, but could only be identified by the presence of open water. In addition, hydraulic conductivity was measured at a shallower depth than the maximum vertical flux was measured. Therefore, it is probable that the hydraulic conductivity of preferential flow paths was higher than measured and winter GW flow rates were higher than calculated. If the flow rate was higher, the GW and surface water may have mixed less and GW temperatures may have been higher than measured, which would increase the GW heat flux (Kurylyk et al. 2013, 2014) and the ice melt rate (Jones et al. 2013).

Of the dominant environmental factors (Table 1), air temperature and snow depth change most rapidly. Modeling results suggest that changes in the timing or magnitude of snow accumulation influenced the potential ice melt rates affiliated with GW. The effect of air temperature on ice melt rate was attenuated as snow accumulated on the ice (Figure 8(b)). During late winter in interior Alaska, snow depths typically exceed the hiemal threshold, thereby thermally isolating the water column from air temperature oscillations. Assuming stable channel morphology, water velocity and depth were not expected to vary significantly through the winter months after stabilizing in November or December, although the proportion of GW to surface water would increase slowly through the winter.

Rural Alaskans have observed that winter travel conditions have become more dangerous due to thin or open ice in recent memory (Herman-Mercer et al. 2011; Schneider et al. 2013a, b). Researchers observed slow-moving open water at temperatures as low as −30 °C (Figure 1), but modeling results suggest that field measurements of GW flow (hydraulic conductivity, vertical gradient, and horizontal velocity) were not sufficient to maintain open water conditions. This research suggests that more extensive characterization of GW flow paths would be necessary to model the influence of GW even in areas of limited spatial extent. It was hypothesized that degrading permafrost in distant recharge areas could allow for enhanced infiltration and greater fluxes in GW upwelling areas to melt ice from below to generate hazardous ice conditions; however, the results of this study were inconclusive and suggest that additional study would be warranted.

In sub-arctic environments, winter travel on rivers is common, especially in the context of subsistence activities. Predictable and reliable environmental conditions help maintain the quality of life for subsistence peoples. While subsistence peoples are adaptable, environmental variability increases their vulnerability to cultural disruptions that include individual or household decisions to migrate from subsistence lifestyles in remote communities to jobs in cash-based economies in cities. This study illustrates the importance of using applied science to investigate the complex interactions between climate, hydrology, and society to understand how increasing variability in climate and natural processes may affect the lives of subsistence users in the north.

This research was made possible by the generous partnership with the Village of Tanana, Alaska and the individuals that assisted with this research. We especially would like to express our gratitude for the support and efforts of our collaborators, including Ruth Althoff, Karen Brewster, Charlie Campbell, Wallie Carlo, Sam Dementieff, Ronnie Evans, Kristy Johnson, Dave Norton, Karl Olsen, Neil Scannell, Marla Statscewich, Matthew Sturm, and Charlie Wright. This research was supported by the National Science Foundation's (NSF) Division of Polar Programs (grants OPP-0422068, ARC-0517762, and ARRA ARC-0909517), NSF's ‘Resilience and Adaptation of Social-Ecological Systems in a Rapidly Changing North’ IGERT program (grant #0654441), and the Bonanza Creek Long-Term Ecological Research program (DEB-10226415). Additional support was provided by the International Arctic Research Center, Alaska EPSCoR NSF (OIA-1208927) with the state of Alaska, Alaska Climate Science Center, UAF Water and Environmental Research Center, UAF Center for Global Change Student Research Grant, and the USGS through the National Institutes for Water Research program. The views expressed in the paper are those of the authors and do not represent the funding agencies.