This study uses a three-dimensional groundwater flow model to investigate groundwater dynamics and groundwater—surface water (GW-SW) interactions considering the effects of permafrost distribution for the Tanana Flats Basin in interior Alaska. The Parameter ESTimation (PEST) code is used to calibrate the model with observed stream discharge data. A 36-year MODLFOW-USG regional simulation shows the following. (1) Permafrost impedes groundwater movement in all directions and through taliks provides a major pathway to connect the groundwater and surface water systems. More than 80% of the vertical groundwater flow occurs within the permafrost-free zones. (2) Permafrost holds a significant amount of water that cannot be easily released through groundwater movements; however, water above the permafrost table has much higher renewal rates than deep groundwater. (3) Groundwater upwelling supports the base flow for the Tanana River and its tributaries throughout the year and feeds water to the wetland ecosystems at the Tanana Flats through unfrozen zones. Stream leakage is also highly correlated with stream discharge. Our study suggests that cold regional hydrological cycle studies should consider the effects of permafrost distribution under future warming conditions. This study provides a robust three-dimensional hydrological modeling tool that can be applied for the regions underlain with either continuous or discontinuous permafrost.
Introduction
Groundwater and surface water, which are often considered as separate water sources, are essentially interdependent within the hydrological cycle (Winter, 1999; Kollet and Maxwell, 2006; Markstrom et al., 2008; Kraemer and Brabets, 2012). Surface water is the principal direct supply, while groundwater is the primary source of freshwater (Alley, 2006; Kundzewicz and Doell, 2009; Richey et al., 2015). Although groundwater in cold regions follows the same principles prevailing in temperate regions, freezing temperatures substantially modify water flow patterns through the frozen ground. When earth material remains at or below 0 °C for at least two consecutive years, it is called permafrost (Vaughan et al., 2013). Permafrost usually behaves like an aquiclude or an aquitard (Woo, 2012); it restricts groundwater recharge, discharge, and movement. It also affects subsurface water storage and movement above the permafrost table (Bense et al., 2009; Ge et al., 2011). In addition, groundwater and surface water (GW-SW) interactions are often obstructed by the presence of permafrost (Woo, 2012). Taken together, permafrost plays an important role in groundwater dynamics and GW-SW interactions. Groundwater in arctic or subarctic regions occurs in suprapermafrost, intrapermafrost, and subpermafrost zones. Permafrost distribution affects both groundwater distribution and its flow patterns. For example, when a talik is formed in continuous permafrost regions, subpermafrost groundwater can discharge directly to suprapermafrost groundwater (Walvoord et al., 2012; Kurylyk et al., 2014b).
To date, many studies have investigated the influences of permafrost distribution on groundwater dynamics and GW-SW interactions. Studies have further reported that river discharges in arctic and subarctic regions are changing under the warming climate, and permafrost degradation plays an important role (Walvoord and Striegl, 2007; Bense et al., 2009; Brabets and Walvoord, 2009; St Jacques and Sauchyn, 2009). For example, it has been shown that contribution to river discharge from groundwater upwelling is gradually increasing due to permafrost degradation (Walvoord and Striegl, 2007; St Jacques and Sauchyn, 2009). In addition, the upwelling of warming groundwater near streams has changed the river ice thickness and its environment (Burril et al., 2010; Herman-Mercer et al., 2011; Jones et al., 2015). Other studies have shown that wetland ecosystems are changing due to permafrost degradation (Viereck et al., 1993; Walters et al., 1998; Jorgenson et al., 2001; Yoshikawa and Hinzman, 2003). Due to the absence of groundwater measurements in remote regions, numerical models have been widely used to estimate ground-water dynamics (Frampton et al., 2011; Ge et al., 2011; van der Ploeg et al., 2012; Bosson et al., 2013; Kurylyk et al., 2014a; Johansson et al., 2015). For example, several studies have used the U.S. Geological Survey (USGS) standard groundwater flow model Modular Finite-Difference Ground-Water Flow Model (MODFLOW) to estimate the aquifer properties in permafrost regions (Nakanishi et al., 1998; Koch et al., 2011; Walvoord et al., 2012). However, most of these studies have not explicitly considered the effects of permafrost distribution, especially near stream channels. They have also omitted the dynamics of aquifer properties in different seasons. Further, they have used one- or two-dimensional approaches at relatively coarse spatial resolutions (Zhuang et al., 2001; Adam and Lettenmaier, 2007; Frampton et al., 2011; Ge et al., 2011). Consequently, great uncertainties in existing studies still remain.
This study uses the three-dimensional (3D) standard groundwater flow model MODFLOW-USG (Unstructured Grid) to estimate the influences of permafrost distribution on groundwater dynamics and GW-SW interactions in the Tanana Flats Basin, Alaska. We explicitly consider the stream effects on permafrost distribution and dynamics of aquifer properties. This study fills a gap in systematically modeling groundwater dynamics and GW-SW interactions considering the effects of permafrost distribution at regional scales. This study provides an effective 3D modeling tool to quantifying the role of permafrost in the arctic hydrological cycle.
Methods
Overview
First, we collected the data of surface elevation, hydrologic networks, permafrost distribution, and river stage/discharge prior to setting up the model. Second, we built the MODFLOW-USG model considering the effects of permafrost distribution Third, we calibrated and evaluated the model with observational data. Fourth, we applied the model to the study area. Finally, we analyzed the influences of permafrost distribution on groundwater dynamics and GW-SW interactions.
Description of Study Area and Data
The study area, the Tanana Flats Basin, is located southeast of the Yukon River Basin in interior Alaska (Fig. 1). Based on the latest geologic map of Alaska, approximately 68% of the study area is covered by unconsolidated sediments (e.g., soil and surficial deposits) of Quaternary age (Wilson et al., 2015). Most sediments are deposited along the alluvial fans of the Alaska Range (Fig. 1). At the foot of the Alaska Range, sedimentary rocks, mainly Nenana gravel and coalbearing rocks with ages from late Miocene to Pliocene, cover up to 11% of the study area. At the Alaska Range, land surface and bedrock are mainly composed of igneous rocks.
The northern part of the alluvial fans and the glacial outwash from the Alaska Range is commonly recognized as the Tanana Flats near Fairbanks (Fig. 1). Surface elevation decreases from approximately 4000 m at Mount Hayes at the Alaska Range to around 120 m at the Tanana Flats in less than 50 km. However, average elevation gradient in the lower portion of the Tanana Flats is only 1.0 m km-1. Poorly developed stream channels run from south to north through the flats. Clear Creek and Wood River are two well-developed continuous tributaries of the Tanana River (Racine and Walters, 1994). Tanana River, which is the largest tributary of the Yukon River, runs from east near Big Delta to west near Nenana along the northern part of the study area (Fig. 1). Groundwater and the Tanana River near Fairbanks have been studied intensively for decades (Williams, 1965; Glass et al., 1986; Racine and Walters, 1994; Walvoord et al., 2012). The groundwater levels in the alluvial plain between the Tanana River and Chena River near Fairbanks have been documented, and the measured depths to the water table were mostly within 4 m (Glass et al., 1986). However, few studies have focused on the Tanana Flats Basin and its hydrologic system. The Tanana Flats includes a variety of habitats, including white spruce forests, black spruce bogs, grasslands, meadows, and wetlands. Many studies suggest that the warming-induced permafrost degradation has gradually led the ecosystem shifting from boreal forests to fens and bogs in this region (Viereck et al., 1993; Jorgenson et al., 2001).
FIGURE 1.
Spatial location and hydrologic networks of the study area. The red-filled polygon in the upper right features the spatial location of the study area in interior Alaska. HUC819040507 is the Hydrologic Unit Code (HUC) from the U.S. Geological Survey (USGS) Watershed Boundary Dataset (WBD). Colored lines with numbered indices are stream segments. The Tanana River is represented by a number of connected stream segments, and the basin outlet is at the last reach of segment 37. Pink dots are the USGS gage stations (Table 1).
The study area lies entirely within discontinuous permafrost zones under the current climate. Approximately 60% of the study area is underlain by permafrost. The simulated active layer thickness (ALT) or seasonal frozen layer thickness (SFLT) data from Geophysical Institute Permafrost Laboratory (GIPL) at the University of Alaska Fairbanks varies between 0.5 and 1.5 m (Appendix, Section 3). The permafrost thicknesses, however, could reach 100 m (Glass et al., 1986; Walters et al., 1998; Romanovsky et al., 2002). Permafrost is generally missing under the floating mat fens and moss bog (Walters et al., 1998). These permafrost-free zones are the major pathway for groundwater discharge to the land surface.
Digital elevation model (DEM) data were obtained from the USGS National Elevation Dataset (NED) (Gesch et al., 2002). Climatic data were from the Global Historical Climatology Network (GHCN) through NOAA's National Climatic Data Center (NCDC) (Menne et al., 2012). Hydrologic data were from the National Hydrography Dataset (NHD)/Watershed Boundary Dataset (WBD) and National Water Information System (NWIS) (USGS, 2012). Among the four USGS gage stations, three were used as model inputs, and the remaining one at the watershed outlet was used for model parameterization (Table 1). Surface runoff and infiltration data were from a 36-year Precipitation-Runoff Modeling System (PRMS) simulation (Markstrom et al., 2015). The ALT&SFLT data were from the numerical simulation by GIPL (2011). Soil data were from the Natural Resources Conservation Service soil survey geographic database (SSURGO data sets) (Soil Survey Staff, 2015). Geology data were from the geologic map of Alaska (Wilson et al., 2015).
TABLE 1
U.S. Geological Survey gage stations used in this study.
Model Description
MODFLOW is a standard 3D finite-difference groundwater model (Harbaugh, 2005; Panday et al., 2013). For decades, it has been used extensively for studying groundwater dynamics and GW-SW interactions. Most studies have used MODLFOW in temperate regions and only a few in cold regions (Walvoord et al., 2012). In our study, the MODFLOW-USG, which supports time-variant aquifer properties, is used (Fig. 2) (Niswonger et al., 2011; Panday et al., 2013).
Our study region is divided into a 3D mesh at a horizontal resolution of 500 m × 500 m in the X and Y directions (Fig. 3, part a). In the Z direction, each column is divided into three layers (Fig. 3, part b). The layer thicknesses are chosen to accommodate the thicknesses of corresponding aquifer properties. The top layer is either active layer or seasonal frozen layer. Its thickness (1.2 m on average) was retrieved from the ALT&SFLT data. This layer is composed of saturated shallow ground-water and a vadose zone and is subject to a dry-rewetting problem due to its shallow thickness; therefore, the Newton-solver method is used (Niswonger et al.,2011). The second layer represents the permafrost layer with a thickness of 100 m on average, according to field observations (Glass et al., 1986; Walters et al., 1998; Jorgenson et al., 2001; Romanovsky et al., 2002). The bottom layer is the subpermafrost aquifer above the impermeable bedrock with a thickness of 400–1200 m (Wilson et al., 2015). Each layer is subdivided into several zones based on permafrost distribution and geologic setting. For example, the top layer includes unconsolidated deposits and fractured siltstone, and the second layer includes unconsolidated deposits and permafrost, whereas the bottom layer includes permeable sand and gravel (Fig. 4). Aquifer properties of these zones are obtained from other studies and model calibration. More details of the discretization of space and time are presented in the Appendix (Section 1).
Due to seasonal freezing and thawing, thickness and the aquifer property of the top layer are constantly changing throughout the year (Hinzman et al., 1991). Ideally, a time-variant spatial discretization and corresponding aquifer properties are required to consider these dynamics; however, these requirements currently cannot be met easily due to modeling limitations and data availability. Therefore, we used time-variant “effective” aquifer properties for the top layer, which enables us to change the aquifer properties while maintaining the spatial discretization. To implement this, we used the Time Variant Material (TVM) package to estimate these “effective” aquifer properties (Panday et al., 2013).
Groundwater divide is used as the lateral boundary. First, a groundwater system underlying a large river such as the Tan ana River usually forms a sub-basin along the entire river (Wilkins, 1997). Second, groundwater divide is a close approximation to the groundwater boundary under natural conditions (Franke et al., 1987; Reilly, 2001). Because areas near the groundwater divide are mostly underlain by permafrost, this assumption should not introduce significant uncertainty to our model.
The spatial extent of the groundwater divide was produced through watershed delineation because their spatial extents are the same. Along with the watershed delineation, the hydrologic networks were also retrieved using the System for Automated Geoscientific Analyses (SAGA GIS) and Arc Hydro tool (Olaya, 2004; Esri Water Resources Team, 2011). Overall, 37 stream segments containing 1570 stream reaches were identified (Fig. 1). Details of these operations are presented in the Appendix (Section 4).
To consider the stream effects on permafrost distribution, we corrected the original ALT&SFLT data based on the hydrologic networks. For example, ALT&SFLT are increased in stream channels even though their physical meanings do not stand anymore while they are required by numerical simulation. More details of these corrections are discussed in the Appendix (Sections 3 and 5). A complete scheme of these configurations is illustrated in Figure 5. In the end, a distribution map of ALT&SFLT was produced (Fig. 6).
Additional model packages include groundwater recharge (RCH) and stream flow routing (SFR). In the Tanana Flats Basin, surface infiltration is the only groundwater recharge source other than stream leakage. Therefore, we estimated the spatially explicit surface infiltration using a 36-year PRMS model simulation. This simulation explicitly considered the glacier effects and snowmelt because the Tanana Flats Basin is a typical glacier-fed catchment (Yang et al., 2009; Wada et al., 2011). Simulated surface infiltration and segment lateral inflow were then used in our MODFLOW simulation. Details of the PRMS model simulation are discussed in the Appendix (Section 6). For the SLR package, all 37 stream segments including the Tanana River were simulated. Because our study area watershed receives inflow from upstream basins, observed stream discharge from three USGS gage stations along the Tanana River were used as inputs (Fig. 1). Besides, simulated segment lateral inflow from the PRMS simulation was also used as segment inflow.
FIGURE 2.
Data flow in the MODFLOW-USG model simulation: DEM data were from the National Elevation Dataset (NED); hydrological data were from the National Hydrography Dataset (NHD), Watershed Boundary Dataset (WBD), and National Water Information System (NWIS); soil data were retrieved from the Soil Survey Geographic database (SSURGO datasets); climate data were from the Global Historic Climate Network (GHCN) through the National Climatic Data Center (NCDC); permafrost data were retrieved from the numerical simulation results by Geophysical Institute Permafrost Laboratory (GIPL). Packages used in MODFLOW-USG simulation include Basic (BAS), Stream Flow Routing (SFR) and Time-Variant Material (TVM) packages. Model outputs include hydraulic head, drawdown, groundwater storage, groundwater flow, and stream leakage.
Initial aquifer properties were obtained from soil survey data (SSURGO) and other studies (Yager, 1993; Batu, 1998; Nakanishi et al., 1998). Hydraulic conductivities of fine grain unconsolidated deposits are approximately 2.4 m d-1 and 7.8 m d-1near the Tanana River where the grain size is coarse (Soil Survey Staff, 2015). Hydraulic conductivity of the permafrost is set to 1.0 × 10-6 m d-1(Walvoord et al., 2012). When unconsolidated deposits are completely frozen in winter, their hydraulic conductivities are extremely low and they were set to the same order as that of permafrost. Given that hydraulic properties in the active layer vary with depth and they are intricately interrelated with thermal properties, the “effective” aquifer properties are believed to be highly nonlinear (Hinzman et al., 1991; Bolton et al., 2008; Quinton and Baltzer, 2013). Therefore, we used a nonlinear interpolation approach to estimate time series of “effective” aquifer properties. Specifically, the sine function and curve fitting techniques are used to describe the aquifer properties with the lowest and highest values in winter and summer, respectively.
A limited amount of groundwater data is available due to inaccessibility, especially for deep groundwater at subpermafrost (van der Ploeg et al., 2012). Therefore, we started the MODFLOW simulation with a steady state simulation, followed by a 36-year transient simulation from 1980 to 2015 (Reilly and Harbaugh, 2004; Walvoord et al., 2012).
FIGURE 3.
Spatial discretization of the MODFLOW-USG model (not drawn to scale), (a) The three-dimensional mesh of the spatial domain (XYZ-coordinate system is rotated around the Z-axis for better visualization of the Alaska Range. The thin red line represents the cross section of the mesh (Fig. 4). (b) A magnified region in (a).
Model Parameterization and Evaluation
The model independent Parameter ESTimation code (PEST) was used to automate the MODFLOW-USG model parameterization (Doherty et al., 1994; Dahlstrom and Carter, 2013). In PEST, the cost function Φ is defined:
where c is the system simulated vector; c0 is the “linearized” system simulated vector; J is the Jacobian matrix of the linearization function M, which maps n-dimensional parameter space into m-dimensional observation space; b is the system parameter vector; and b0 is the corresponding system parameter vector using M. To minimize the cost function, the “linearized” system parameter vector b0 is updated through iterations until Φ is reduced to certain criteria. To reduce the computational demands, we used a special version of PEST, BeoPEST, which supports Message Passing Interface (MPI) communications during parameter estimation (Hunt et al., 2010).
To construct the parameter vector b0, we selected all the model parameters. In order to reduce the number of parameters, the three-dimensional domain is subdivided into a number of zones, within which the parameters are assumed as constant (Fig. 4). A sensitivity analysis was conducted using BeoPEST to identify the most sensitive parameters. Initial values of these parameters were obtained from laboratory and other studies (Table 2) (Batu, 1998; Nakanishi et al., 1998; Walvoord et al., 2012). We also used the Marquardt parameter and singular value decomposition (SVD) method to address the “hemstitching” and invertible matrix problems in the parameter estimation process.
FIGURE 4.
The geological profile of the cross section in Figure 3 (not drawn to scale). Horizontally, there are four rock types, including unconsolidated deposits, sedimentary rock, igneous rock, and metamorphic rock. Vertically, each column is divided into three layers. The top layer is either active layer or seasonal frozen layer. The second layer contains both permafrost and permafrost-free zones. The bottom layer contains permeable sand and gravel. Near the Alaska Range, rocks are assumed to be continuous in depth (fractured on surface) and they have the same composition with bedrock.
We used 10-year daily observed stream discharge from 1980 to 1989 at the watershed outlet (Site ID: 15515500) to parameterize the model, which is long enough to calibrate the model. Model calibration results are listed in Table 2. Then we evaluated the model performance with the remaining data from year 1990 to 2015. Comparisons between observed and simulated stream discharge rates for both calibration and evaluation periods have shown that our estimates are accurate with a Nash-Sutcliffe efficiency (NSE) of 0.90. Time series and scatter plots have also shown that our estimates match observations with a Pearson's coefficient (r) of 0.97 (Fig. 7, parts a and b).
Model Simulations
The parameterized MODFLOW-USG model was used to estimate the influences of permafrost distribution on groundwater dynamics and GW-SW interactions. The simulation runs from 1 January 1980 to 31 December 2015 at a daily time step. We then analyzed the groundwater dynamics including hydraulic head, groundwater movements, and storage as well as the groundwater and stream water interactions.
Results and Discussion
Influences of Permafrost on Hydraulic Head and Drawdown
Simulated hydraulic heads generally follow surface topography. In the top layer, horizontal head gradients vary with surface elevations (Fig. 8, part a). In the southern part of the study area, heads change rapidly due to topography effects. In the central part, heads slowly decrease (1.0 m km-1 on average) from south to north. In the northern part, heads gradually decrease (0.9 m km-1) from southeast to northwest. This spatial pattern is consistent with previous studies (Anderson, 1970; Glass et al., 1986; Grinevskii, 2014). Most heads in the Tanana Flats are within 2.0 m below the surface. However, in areas where groundwater feeds the fens system, heads are higher than the surface elevation. In the second and bottom layers, spatial variability is small due to groundwater movements, which is also consistent with previous findings (Fig. 8, parts b and c) (Walvoord et al., 2012). Vertically, heads gradients vary with locations. The vertical gradients in the Tanana Flats are positive whereas they are negative near the Alaska Range. As a result, the groundwater system receives recharge (e.g., surface leakage) from areas near the Alaska Range whereas there is discharge to the surface (e.g., upwelling or springs) in the Tanana Flats.
FIGURE 5.
Interactions between a stream and a groundwater system underlain by permafrost and permafrost-free zones (not drawn to scale). The red (blue) line represents the high (low) gage height, and red (blue) arrows represent principal groundwater flow components when the stream is losing (gaining) water to (from) the groundwater system. The unsaturated zone may form beneath the streambed when a losing stream recharges the groundwater system. A large stream in the discontinuous-permafrost zone usually produces a temperature anomaly that forms a thaw bulb. The buffer zones were created to account for the lateral water flow between the stream channels and its banks.
Simulated drawdowns in the top layer are usually within 1.0 m. However, their magnitudes are higher near stream channels. In the second layer, simulated drawdowns are less than 0.1 m. Both heads and drawdowns are directly affected by surface topography and permafrost distribution. Because groundwater movements within permafrost zones are extremely slow, the magnitude of drawdowns is very small. However, in permafrost-free zones, groundwater movements are much faster, and drawdowns are also larger. In other layers, the influences of permafrost distribution on heads and drawdowns are not significant.
Influences of Permafrost on Groundwater Flow
Simulated cell-by-cell water flow rates within permafrost zones are very slow. In the second layer, simulated average horizontal flow rate within permafrost zones is 3.6 × 10-3 m3 d-1.1 In contrast, the rate in permafrost-free zones is 7.7 × 10-2 m3 d-1, which is about 20 times higher. The spatial pattern of horizontal flow rates is highly correlated with permafrost distribution. In the top layer, simulated horizontal flow rates are much higher (1.9 m3 d-1), and their spatial patterns are determined mainly by surface topography and stream channels instead of permafrost distribution (Fig. 9).
Because vertical groundwater movements connect both the top and bottom layers with the second layer, simulated vertical flow rates are always influenced by per-1 mafrost distribution. The average flow rates above permafrost zones and permafrost-free zones are 0.4 m3 d-1 and 2.2 m3 d-1, respectively. Their spatial patterns are also correlated with permafrost distribution (Fig. 10). Vertical flow rates between the top and second layers have also confirmed that groundwater upwells (more than 3.0 m3 d-1 in summer) to the surface in the Tanana Flats. This is consistent with previous studies (Racine and Walters, 1994; Walters et al., 1998).
FIGURE 6.
The spatial distributions of the Active Layer Thickness (ALT) and Seasonal Frozen Layer Thickness (SFLT) (unit: meter). Approximately 61% of the study area is underlain by permafrost. It is assumed that the Tanana River and its major tributaries have formed the unfrozen zones that perforate the permafrost. Both ALTs and SFLTs are corrected near the stream channels (see Appendix).
TABLE 2
Estimated values of model parameters from BeoPEST calibration.
Due to seasonal freezing and thawing in the top layer, simulated groundwater flow rates have shown a significant seasonality. The average horizontal flow rates in winter and summer are 5.8 × 10-2 m3 d-1 and 3.5 m3 d-1, respectively. The vertical flow rates between the top and second layers in winter and summer are 0.8 m3 d-1 and 2.2 m3 d-1, respectively.
FIGURE 7.
Comparison between observed and simulated stream discharge rates at the watershed outlet (site ID: 15515500). (a) and (b) The time series and scatter plots of the observed and simulated discharge rates, respectively.
Our simulation shows that, because of its extremely low hydraulic conductivity, permafrost impedes groundwater movements in all directions. Therefore, groundwater flow within or near permafrost zones is very low, and most flow is within permafrost-free zones (e.g., through taliks), which are the major pathway for GW-SW interactions. Simulation also shows that influences on horizontal flow in the top and bottom layers are very different. In the top layer, the average horizontal flow rate above permafrost zones (2.7 m3 d-1) is slightly higher than that above permafrost-free zones (1.0 m3 d-1). In the bottom layer, the average horizontal flow rate under permafrost zones (2.4 × 10-2 m3 d-1) is also higher than that under permafrost-free zones (1.3 × 10-2 m3 d-1). This implies that restriction on vertical flow because of permafrost may enhance horizontal flow. In addition, our simulation shows that, although the top layer is the thinnest layer comprising only 0.1% of the total volume of the entire spatial domain, it contributes more than 10% and 68% to the total horizontal and vertical groundwater flow, respectively. Combined together, 15% of the groundwater flow occurs within the top layer.
Influences of Permafrost on Groundwater Storage
Simulated groundwater storage change rates within permafrost zones are very low (3.1 × 10-3 m3 d-1 on average). Therefore, a large amount of water stored within permafrost zones cannot be easily released. In contrast, the change rates in permafrost-free zones are much higher (2.2 × 10-2 m3 d-1). The spatial pattern of simulated change rates in the second layer is also highly correlated with permafrost distribution. As a result, the warming-induced permafrost degradation in the near future may thus significantly impact the groundwater distribution in discontinuous-permafrost regions.
The top layer has the largest groundwater storage change rate (4.4 m3 d-1 on average); however, it is still affected by permafrost distribution. The average change rates above permafrost and permafrost-free zones are 2.2 m3 d-1 and 6.7 m3 d-1, respectively. Similar to groundwater flow rates, simulated storage change rates also have shown strong seasonality. The average change rates in winter and summer are 0.2 m3 d-1 and 25.0 m3 d-1, respectively.
FIGURE 8.
The spatial distributions of simulated hydraulic heads in (a) the top, (b) middle, and (c) bottom layers in summer 2015. The spatial patterns of heads are identical to that of the surface elevation. The spatial variabilities of simulated heads in the lower layers are smoothed due to groundwater movements.
Given that the top layer has the lowest volume (0.1% of the total volume) but contributes the highest storage change (62% of the total storage change), the water renewal rate in this layer is the fastest. Therefore, it is important to understand the special role of the active layer and its relationship with permafrost in regional hydrological dynamics.
Influences of Permafrost on Groundwater and Stream Interactions
Groundwater and stream water interactions are important components within the hydrological cycle, especially for arctic and subarctic hydrological systems under a warming climate. Recent studies also have suggested that changes in stream discharge (e.g., increasing winter base flow) in high latitudes are the results of changing groundwater and stream water interactions due to permafrost degradation (Walvoord and Striegl, 2007; Brabets and Walvoord, 2009; St Jacques and Sauchyn, 2009; Lyon and Destouni, 2010; Niu et al., 2011; Walvoord et al., 2012). Our simulated ground-water and stream water interactions, specifically stream leakage, have shown significant spatial-temporal variations over the decades. First, it has been shown that these interactions are active throughout the year. Second, while some streams (e.g., the Tanana River) are recharging the groundwater system, others (e.g., tributaries of the Tanana River) may be receiving ground-water upwelling. Even within the same stream, some reaches may be recharging while others are receiving groundwater at the same time (Fig. 11).
In winter, groundwater upwelling dominates the interactions, and it supports the base flow for the Tanana River and its tributaries. For the Tanana River in our study area, the total groundwater upwelling rate (negative stream leakage) is about 1.8 × 10)6 m3 d-1. In contrast, the total stream leakage rate is less than 1.0 × 103 m3 d-1. As the river discharge increases rapidly in spring and summer, the total stream leakage also increases significantly (1.0 × 107 m3 d-1). As a result, the total stream leakage often overruns the total groundwater upwelling, and the Tanana River recharges the groundwater system (Fig. 12).
FIGURE 9.
The spatial distribution of simulated horizontal flow rates in the top layer (Y direction) on 31 August 2015 (units: cubic meters per day). Horizontal flow rates in the top layer are directly affected by the hydrologic networks.
FIGURE 10.
The spatial distributions of simulated vertical flow rates (Z-direction) between the top and second layers on 31 August 2015 (units: cubic meters per day). Vertical flow rates are correlated with permafrost distribution.
FIGURE 11.
The spatial distribution of stream leakage rates on 31 August 2015 (units: cubic meters per day). In general, lowland streams recharge the groundwater system, and high-altitude streams receive the groundwater up welling.
TABLE 3
The average groundwater flow rates in different directions in permafrost and permafrost-free zones* (units: cubic meters per day).
FIGURE 12.
Time series of the total stream leakage rate in the Tanana River in our study area from 1 January 1980 to 31 December 2015. The pink-filled plot is a magnified portion (from 1 March 2015 to 31 October 2015) of the whole plot.
Because we assumed that the Tanana River produces a temperature anomaly that forms a thaw bulb, groundwater flow beneath the Tanana River is very fast. Our simulation has shown that stream leakage in the Tanana River is much larger than that in streams in remote regions. In addition, spatial-temporal variations of stream leakage in the Tanana River are significantly larger. Therefore, the warming-induced permafrost degradation will also change the groundwater and stream water interactions in discontinuous-permafrost regions.
Conclusions and Limitations
Our model simulation indicates that hydraulic heads and drawdowns as well as groundwater movements are influenced significantly by the permafrost distribution in our study area. Groundwater flow in permafrost-free zones contributes 44% and 83% to the horizontal and vertical flow, respectively, even though these zones cover only 39% of the study area. Most groundwater and surface water interactions (e.g., groundwater upwelling) are associated with permafrost-free zones (e.g., through taliks). The fast storage renewal rates in the top layer are indicative of the importance of the active layer. This layer contributes 10% and 80% to the total horizontal and vertical groundwater flow although it comprises only 0.1% of the total volume. In contrast, permafrost holds a large amount of water that cannot easily flow out the groundwater system due to its low hydraulic conductivity. The stream leakage rate above permafrost-free zones is much higher than that above permafrost zones. Stream leakage rates are correlated with stream discharge rates. Therefore, changes in surface hydrology including snow dynamics will influence groundwater and stream water interactions as well.
Our future research will address several limitations of this study. First, we will collect data of the cross section profiles of major stream channels to create stream geometry, since it plays an important role in groundwater and stream water interactions. Currently, we can only approximately infer stream geometry based on remote sensing images and hydrologic data sets. Second, we will consider more detailed surface hydrological processes in mountainous areas (e.g., the Alaska Range) where groundwater systems are very sensitive to incoming recharge. Snowmelt and glacier melting are the major recharge to groundwater systems in glacier-fed catchments. In our study, surface hydrology dynamics are simulated independently using a 36-year PRMS model simulation. Therefore, uncertainties in PRMS simulations could be propagated into our MODFLOW-USG model simulations.
Acknowledgments
This study is supported through projects funded to Zhuang by the NASA Land Use and Land Cover Change program (NASA-NNX09AI26G), Department of Energy (DE-FG02-08ER64599), the NSF Division of Information and Intelligent Systems (NSF-1028291), the NSF Carbon and Water in the Earth Program (NSF-0630319), and a project from the USGS focusing on quantification of carbon and methane dynamics in Alaska. We gratefully thank the technical support received from USGS PRMS, MODFLOW, and GSFLOW teams. We also appreciate assistance from Richard Winston, Richard Niswonger, and Damian Merrick on model development and Willem Schreuder on model calibration.
References Cited
Appendices
Appendix
Supporting Information
Discretization of Space and Time
Different combinations of cell size and layer configuration with consideration of computational demands were tested and the most appropriate settings were chosen for the numerical simulation.
Horizontally, the spatial resolution of the MODFLOW-USG model is 500.0 m in X (longitude) and Y (latitude) directions, which is fine enough to distinguish all major hydrologic features.
Vertically in Z direction, each column is divided into 3 layers. The layer thicknesses are defined on the basis of aquifer thicknesses (Fig. 3). The top layer is either active layer or seasonal frozen layer (see the definition in Appendix, Section 3). Its thickness was defined using the ALT&SFLT data. Because the average thickness of this layer is only ∼1.2 m, cells in the top layer are particularly subject to the drying and rewetting problem (Harbaugh, 2005). Therefore, the Newton Solver was used (Niswonger et al., 2011). The second layer contains both permafrost and permafrost-free zones. Several studies have shown that the permafrost thickness in this region varies between a few meters and more than 150 m (Jorgenson et al., 2001; Romanovsky et al., 2002). However, the spatial distribution data of the permafrost thickness is currently unavailable. Therefore, we assumed that the permafrost thickness is 100 m on average. The bottom layer is the deep groundwater aquifer above the impermeable bedrock. Its thickness was defined based on the geologic data (Wilson et al., 2015).The topography effects on the aquifer thickness were also considered. The three-dimensional (3D) domain discretization is illustrated in Figure 4, parts a and b.
MODFLOW supports various combinations of stress period and time step; however, a daily stress period and single time step combination is commonly used since it is sufficiendy small for MODFLOW to converge (Markstrom et al., 2008). In our study, most of the inputs including stream discharge and infiltration are at a daily time step. Therefore, the daily stress period and single time step were used.
Boundary and Initial Condition
The specification of the boundary and initial conditions is an essential part in numerical groundwater modeling. Direct groundwater measurements across the study area are generally missing, and boundary conditions were defined using a few MODFLOW packages.
First, the groundwater divide was used as the simulation boundary. It is the horizontal boundary of the 3D mesh (Fig. 4, part a). We assumed that no external groundwater flows into the basin at the boundary. Under natural conditions, groundwater divide is a close approximation to the groundwater flow boundary (Reilly and Harbaugh, 2004). Moreover, areas near the boundary are mostly underlain by permafrost, which has extremely low hydraulic conductivity. As a result, external groundwater flow can be omitted. The boundary was defined in the Discretization (DIS) file (Harbaugh, 2005).
Second, groundwater recharge was defined using the RCH package. Groundwater recharge in the study area mainly comes from surface infiltration other than stream leakage. In the lowlands, surface infiltration comes from snowmelt in spring and precipitation in summer. While in high altitudes, surface infiltration also comes from glacier melting. These processes were explicitly considered and simulated using a 36-year PRMS model simulation (See Appendix, Section 6).
Third, the Tanana River near Fairbanks receives upstream discharge. Therefore, the U.S. Geological Survey (USGS) gage station observational data were used as inflow data for corresponding stream segments in the SFR package. Moreover, stream segments also receive overland runoff and subsurface soil interflow. Therefore, simulated lateral segment inflow from the PRMS model simulation was also used.
Details of the PRMS model simulation were described in Appendix, Section 6.
Permafrost Distribution
Initial permafrost distribution was defined using the ALT&SFLT data. By definition, active layer is the area underlain by permafrost, whereas the seasonal frozen layer is the area underlain by permafrost-free zones but subject to seasonal freezing and thawing.
The ALT&SFLT data were obtained from numerical simulation results by GIPL (Jafarov et al., 2012). GIPL is a quasi-transitional, spatially distributed, analytical model used to estimate the active layer thickness considering air temperature, snow water equivalent, and other factors. These data were obtained from https://catalog.data.gov/dataset/gipl1-3-simulated-maximum-active-layer-thickness-alt-in-meters-averaged-for-particular-decade-f. Originally, the ALT&SFLT in GIPL are merged into one raster data set where positive and negative values represent the ALT and SFLT, respectively. In our study, their magnitudes were used to define the top layer thickness. In addition, we used an additional flag or color bar scale (Fig. 7) to illustrate the permafrost distribution. The original data were also resampled to 500.0 m spatial resolution.
FIGURE A1.
Spatial distribution of rock types based on geologic data. This map is the result of a reclassification of the raw geologic data.
Because GIPL does not explicitly consider the stream effects on permafrost distribution, corrections are made to the resampled ALT&SFLT data (Appendix, Section 5).
Stream Flow Routing
The hydrological networks were retrieved from watershed delineation using the 1:100,000-scale (medium resolution) NHD flowline and the downscaled DEM data. The System for Automated Geoscientific Analyses (SAGA GIS) and Arc Hydro tool (Esri Water Resources Team, 2011) were used to delineate the watershed using the following steps: (1) “burn” the NHD flowline through DEM reconditioning (Hellweger and Maidment, 1997); (2) fill the depressions in DEM with a minimum slope; (3) calculate the flow direction and accumulation; (4) define the drainage line and subbasin; and (5) correct potential errors (e.g. spikes on the edges).
The watershed delineation results were further validated with the raw NHD/WBD data. The corrected watershed boundary was then used as the groundwater divide. The topological relationships among stream segments, as well as the properties of each stream segment and its reaches, were retrieved with the aid of GIS analysis. Overall, 37 stream segments, which contain 1570 stream reaches, were identified and prepared for the model simulation.
Stream Effects
Permafrost distribution is directly influenced by stream effects. Not only do stream effects change the permafrost distribution around stream channels, but they also affect the 3D mesh layer configuration.
First, in the discontinuous-permafrost zone, unfrozen zones lie beneath most of the large rivers and streams. The unfrozen zones usually perforate permafrost throughout the year. The borings data beneath the Tanana River have indicated the existence of this type of unfrozen zone (Williams, 1965). Therefore, we assumed that major stream channels, including the Tanana River, Clear Creek, and Wood River, are underlain by unfrozen zones that perforate the permafrost. As a result, the initial ALT near-stream channels were converted to SFLT. However, for small streams near the Alaska Range, the permafrost may still prevail under certain depths. Although the physical meanings of ALT&SFLT do not stand anymore within the stream channels, they were retained for numerical simulation (Fig. 7).
Second, GW-SW interactions alter the layer configuration near stream channels. Streams can either gain or lose water from or to the groundwater system under different water table conditions. Studies also have observed the lateral flow between the Tanana River near Fairbanks and its banks during gaining or losing (Glass et al., 1986; Nakanishi et al., 1998). Therefore, the following additional corrections were made to address the stream effects on layer configuration: (1) ALT&SFLT in stream channels were modified on the basis of gage height observations, if available. Because stream segments were defined in the top layer, the ALT&SFLT data (less than 2.0 m) must be higher than the maximum gage height; (2) unsaturated zones beneath stream and associated vertical flow were considered to be in the top layer (Krause et al., 2013); and (3) buffer zones were created to account for the lateral flow between stream water and shallow groundwater. This is also required for numerical simulation convergence stability. A completed scheme of these settings is illustrated in Figure 6. After all, the updated spatial distribution of ALT&SFLT is illustrated in Figure 7.
PRMS Simulation
A 36-year PRMS model simulation was conducted ahead of our MODFLOW model simulation. PRMS is a deterministic, spatially distributed parameter—a physical process-based modeling system developed to evaluate the response of various combinations of climate and land use on streamflow and general watershed hydrology (Markstrom et al., 2015).
To provide spatially explicit inputs for our MODFLOW simulation, the PRMS model was developed in the same spatial-temporal resolution as the MODFLOW model. Input data for the PRMS simulation include grid-based climate data, vegetation type, and soil type. Grid-based climate data are produced using the thin plate smoothing spline surface-fitting technique. The ANUSPLINE package from Australian National University is used to implement this technique (Hutchinson and Xu, 2013).
Special attention has been paid to snow dynamics and glacier effects because the Tanana Flat Basin is a typical glacier-fed catchment. Elevation-dependent parameters were used to account for the snow dynamics and glacier effects (Molotch et al., 2005).
The PRMS model was calibrated with BeoPEST using observed stream discharge and snow data. Afterwards, simulated surface infiltration and surface/subsurface runoff were prepared for the MODFLOW simulation.
