Hydrologic Impacts of Thawing Permafrost—A Review

Permafrost hydrology is a rapidly progressing research field, and a number of new discoveries and questions have emerged in recent years. Research interest in cold regions has been spurred in part by surface temperature warming rates in high latitudes (McBean et al., 2005) and high altitudes (Pepin et al., 2015) that are greater than the global average. This warming has produced changes to the cryosphere, including permafrost (ground that is ≤0°C year round), that impact hydrologic processes and conditions (ACIA, 2005; Hinzman et al., 2013). Despite increased attention, there are still critical limitations in hydrologic data coverage, subsurface characterization, process-level understanding, and integrated modeling approaches. Due to its low hydraulic conductivity (K), permafrost strongly affects the movement, storage, and exchange of surface and subsurface water. In turn, subsurface flow can influence permafrost distribution by enhancing the transfer of thermal energy via heat advection (de Grandpré et al., 2012; Sjöberg et al., 2016). This interplay, together with additional feedbacks among physical, chemical, and biogeochemical processes, creates complex, and often nonintuitive, dynamics in permafrost regions. Understanding hydrologic changes in response to climate-induced permafrost thaw is critical for anticipating changes in ecosystem services and dynamics at local to regional scales (Jorgenson et al., 2013) and for constraining the strength of the permafrost–carbon feedback at the global scale (Schuur et al., 2015; Lawrence et al., 2015).

In temperate climates, near-surface hydrologic analyses must account for the fact that the porous medium contains variable amounts of liquid water, with the balance of the pore space occupied by gas. In cold regions, much of the “vadose” or shallow zone may be fully occupied by water, but that water is variably partitioned between liquid and solid (ice). Some of the concepts from gas–liquid unsaturated-zone hydrology, such as soil–water retention curves and unsaturated K relations, have been applied to liquid–water–ice systems as reviewed by Kurylyk and Watanabe (2013). Many vadose zone hydrologic problems in temperate regions can be idealized as one-dimensional vertical flow problems. In contrast, in permafrost regions hydrologic flow problems are inherently three-dimensional and in some cases tied to subtle variations in topography (Painter et al., 2013). Also in permafrost hydrology applications, the transfer and transformations of heat are of paramount importance.

Previous papers in Vadose Zone Journal, including those in a 2013 special section (see Toride et al., 2013 and papers highlighted therein), have focused on the generalizations of unsaturated flow theory and advances in characterization and modeling of frozen soil processes. Hayashi (2013) provides an overview of hydrologic, ecologic, and mechanical processes occurring in seasonally frozen soil and offers perspective on pressing science needs in cold vadose zone research. The objective of this review is to expand these lines of discussion to pan-Arctic permafrost and address the hydrologic consequences of climate warming with focus on the discontinuous permafrost zone due to its enhanced susceptibility to near-term thaw. The section that immediately follows introduces fundamental concepts and terminology related to hydrological processes and phenomena that occur in permafrost regions. The next section discusses how the distribution of permafrost affects surface and subsurface routing of water through landscapes and the consequent hydrologic impacts of permafrost thaw. Then, the state of permafrost mapping and methods for multiple scale characterization are reviewed. Lastly, the modeling sections highlight recent improvements in simulating hydrological and hydrogeological processes in permafrost environments and conclude by identifying several challenges and opportunities for future advancements in this field.

Current estimates indicate that permafrost regimes cover approximately 24% of the exposed land surface of the Northern Hemisphere (Brown et al., 1997, 2002; Fig. 1). Permafrost can be present in soil, sediment, or rock and is defined by ground that is cryotic (≤0°C) for at least two consecutive years (e.g., Dobinski, 2011). The permafrost zone extends vertically from the permafrost table, encountered on the order of tens of centimeters to several meters below the ground surface, to the permafrost base several meters to 10³ m deep (Fig. 2) depending on conditions. Above the permafrost table, ground temperatures exceed 0°C for some duration during the summer, and below the permafrost base, perennially noncryotic conditions exist due to the influence of the geothermal heat flux. Strictly speaking, the active layer thickness is defined as the lesser of the maximum seasonal frost depth and the maximum seasonal thaw depth. In discontinuous permafrost zones, where depth of unfrozen soil at the end of summer can exceed the frost depth, a perennially unfrozen zone can exist above the permafrost. Freezing-point depression, induced by adsorption and capillary forces, high overburden pressure, and pore water solute concentration, causes liquid water to persist in porous media and fractured rock at temperatures <0°C (Watanabe and Mizoguchi, 2002; Rempel, 2012). As such, the permafrost table can exist slightly above the maximum depth of the active layer, and the permafrost base extends below the maximum depth of the perennially frozen zone where most of the pore water exists as ice (Woo, 2012) (Fig. 2). The zones where a significant fraction of liquid water persists in the pore space at ground temperatures slightly below 0°C are called the cryopeg (base of permafrost) and the transition zone (top of permafrost). The permafrost transition zone is of interest from the perspective of thaw vulnerability and potential carbon decomposition and mobilization even at temperatures below 0°C (Waldrop et al., 2010).

The transition from unfrozen to frozen ground coincides with a reduction in K of several orders of magnitude between 0 and −0.5°C for saturated porous media (Burt and Williams, 1976; McCauley et al., 2002). Thus, permafrost distribution can substantially influence subsurface water pathways and fluxes (Walvoord et al., 2012; Frampton and Destouni, 2015). Frozen ground is often conceptualized as an impermeable barrier that inhibits infiltration and promotes surface and near-surface runoff. However, unsaturated frozen and partially frozen ground may allow for considerable flow through macropores (e.g., Mackay, 1983; Boike et al., 1998; Scherler et al., 2010) over short durations until the infiltrating water freezes and further reduces the soil K. Infiltration into frozen soils may be further restricted due to cryosuction-induced upward migration of water that eventually freezes and blocks the downward flow of water (Stähli et al., 1999; Scherler et al., 2010).

The active layer exerts control on surface and near-surface water storage, drainage, and routing. The thawing front in the active layer moves progressively downward into late summer, separating thawed and frozen soil. If the frozen soil beneath the thawing front is ice-saturated (thus relatively impermeable), the active layer can function as a very shallow perched aquifer that controls runoff and streamflow response to snowmelt and summer precipitation (Carey and Woo, 2005; Yamazaki et al., 2006; Wright et al., 2009; Koch et al., 2014), as well as isotope and nutrient transport and cycling (Koch et al., 2013; Tetzlaff et al., 2015). Thus, the depth to the thawing front is an important variable in permafrost hydrology.

Taliks are unfrozen zones within permafrost and may be characterized as “closed” or “open,” depending on whether the talik is bounded by permafrost (closed) or allows for a connection between unfrozen zones (open). Taliks can form beneath surface water bodies due to the high heat capacity of water, the influence of the water surface on the surface energy balance, and the reduced heat transfer (e.g., wind shear mixing) during winter ice-covered conditions. The propensity for open talik development increases for surface water bodies that do not freeze to their beds in the winter (Burn, 2005). Upwelling from deep groundwater, including thermal and saline springs, can give rise to open taliks (Woo, 2012). Taliks can also form in response to land disturbances, such as wildfires (Jafarov et al., 2013; Minsley et al., 2016) and infrastructure development (Nelson et al., 2001). Talik evolution initiates as heat is conducted into permafrost from a heat source at a temperature above 0°C. Once subsurface water flowpaths become connected through permafrost (talik breakthrough), advection can potentially add to the transfer of heat and accelerate thaw given sufficient subsurface K (Rowland et al., 2011; Wellman et al., 2013; McKenzie and Voss, 2013; Fig. 3b).

Subsurface water flow in permafrost environments can occur above permafrost (seasonally and perhaps perennially depending on conditions), below permafrost, and within taliks (e.g., Kane et al., 2013). These flow zones are known as suprapermafrost aquifers, subpermafrost aquifers, and intrapermafrost groundwater, respectively (Woo, 2012; Fig. 3). Fluxes and exchanges in these zones will depend on typical hydrogeologic considerations, including the unfrozen K of the substrate and hydraulic gradients in the system. The generalized framework results in (i) vertical partitioning between shallow and deep flow systems (i.e., supra- and subpermafrost aquifers) supported by notable seasonal variability in water river chemistry as the dominant source of water changes (i.e., shallow organic-rich soil water in summer–autumn vs. deep mineral-rich groundwater in winter; O’Donnell et al., 2012a), and (ii) lateral partitioning that limits subsurface hydrologic connectivity among inland waters (Connon et al., 2014; Jepsen et al., 2015). Because of its hydrogeologic function as an aquitard, permafrost influences water and solute subsurface flowpaths, residence times through organic and mineral soils (Frampton and Destouni, 2015), and the magnitude of groundwater discharge to streams (Walvoord et al., 2012), thereby impacting aquatic chemical exports (Suzuki et al., 2006; Striegl et al., 2007; Walvoord and Striegl, 2007; Vonk et al., 2015).

There is ample evidence that permafrost is warming and thawing in the pan-Arctic basin (Osterkamp, 2005; Harris et al., 2009; Romanovsky et al., 2010). It is projected that these trends will continue and lead to large-scale losses of near-surface permafrost (Slater and Lawrence, 2013; Koven et al., 2013; Pastick et al., 2015). However, nonuniform rates of permafrost degradation and irregular spatial distribution of thaw are anticipated, thereby imposing important sources of uncertainty in estimates of future conditions.

The most readily observed mode of permafrost thaw is increasing active layer thickness (ALT). Many studies report large interannual ALT variability and note a surprising lack of consistent observed ALT increases (see references in Table 1), despite ubiquitous positive air temperature trends. This apparent discrepancy underscores the fact that the thermal regime in permafrost soils is additionally mediated by interactions among soil properties, moisture, snow, and vegetation that positively or negatively influence permafrost stability (Jorgenson et al., 2010). For example, pore ice has a higher thermal conductivity than pore water, and thus soil heat transfer is typically more efficient in winter than summer. This leads to a “thermal offset” between the mean annual temperatures at the permafrost table and the ground surface (Smith and Riseborough, 2002). A pan-Arctic assessment by Park et al. (2013), using ALT observational data and a land surface model, highlighted the importance of hydrologic variables, including snow depth and summer soil moisture, which counterbalance or amplify ALT increases expected from atmospheric warming alone. As an example, a thinner snowpack may impede soil warming by offering weaker winter insulation to cold atmospheric air, whereas a thicker snowpack may enhance soil warming through the opposing effect. Variability of snow depth and summer soil moisture may help explain positive ALT trends in Eurasia and the contrasting lack of consistent trends observed in North America.

Although large-scale permafrost models generally agree that permafrost spatial extent has decreased in the past several decades with southern boundaries moving northward and altitudinal boundaries moving upward (Zhang et al., 2008a; Harris et al., 2009), observational support is limited by the lack of baseline information and established methods for effectively mapping large-scale permafrost at sufficient resolution to detect change. As an alternative, landscape and vegetation changes strongly linked to permafrost degradation and detectable via aerial or satellite imagery have been used to infer reductions in permafrost spatial extent in study areas of interior Alaska, USA (Jorgenson et al., 2001) and Northwest Territories, Canada (Quinton et al., 2011).

An increase in the abundance of open vertical taliks, particularly beneath water bodies and disturbed zones, may be expected in some areas and has been inferred from observed drainage of thermokarst lakes (Yoshikawa and Hinzman, 2003; Smith et al., 2005). However, lack of adequate historical context and limited current data presents a challenge in distinguishing between the classic evolution of thermokarst lakes and an actual trajectory of climate-induced widespread change. An increase in suprapermafrost and intrapermafrost taliks is expected and modeled (e.g., Zhang et al., 2008a), yet observational support for this contention is also limited. A study by Jepsen et al. (2013) provides indirect support for the growth of supra- and intrapermafrost taliks in the Yukon Flats of Alaska over the past several decades based on a large-scale examination of lake area dynamics and subsurface permafrost characterization.

Permafrost degradation generated through the modes described above will likely produce large changes in surface and subsurface hydrology in some areas. The input of permafrost meltwater contributes a relatively small and transitory pulse of subsurface water that is unlikely to be a dominant, long-term signal in streamflow records in the watershed or basin undergoing permafrost thaw. The more impactful hydrologic modification from permafrost thaw is in the concomitant change of the hydrogeologic framework, particularly in permeable settings that allow for the opening of previously blocked (permafrost-limited) vertical and lateral flowpaths that can transmit large groundwater fluxes when thawed. Permafrost thaw can also generate rapid landscape changes in certain settings (e.g., thermokarsting and plateau subsidence; Quinton et al., 2011) that in turn influence surface water storage, routing, and runoff (Connon et al., 2014).

Changes in the three-dimensional distribution of permafrost have the potential to influence surface and subsurface water fluxes and flowpaths (Fig. 4) at local, regional, and circum-Arctic scales. Large-scale Arctic assessments project a transition away from a surface-water dominated system to a more groundwater-dominated system. Hydrologic and geochemical support for this contention is derived primarily from sub-Arctic and low Arctic basins (Table 1) with warm, discontinuous permafrost. Permafrost thaw may enhance subsurface fluxes, including soil drainage and recharge, suprapermafrost flow, groundwater–surface water exchange, subpermafrost flow, and baseflow. Direct and indirect evidence for such changes in water fluxes induced by permafrost thaw has been observed throughout northern regimes (Table 1). However, due to inherent geologic heterogeneity and the hydrodynamics (transient water table elevations and hydraulic gradients) associated with permafrost thaw, as well as additional hydroclimatic variables (i.e., changes in precipitation, air temperature, evapotranspiration, and snow), actual spatial and temporal changes in subsurface water fluxes present a challenge for prediction. For example, water tables in the active layer may decline with permafrost degradation, and the reduction in hydraulic gradient may or may not be countered by increases in transmissivity affected by permafrost thaw.

Increased ALT resulting in enhanced flow through the suprapermafrost aquifer is often postulated to be the cause for the increased baseflow detected in the historical record in many permafrost basins (Table 1). However, the saturated K of a vertical soil profile typical of boreal forests with a thick organic layer overlying mineral soil is known to decrease exponentially with depth (Quinton and Baltzer, 2013; Streletskiy et al., 2015). Thus, the thickening of the active layer can drive more water to be transmitted through a lower K zone (Koch et al., 2014). If water table elevation is maintained in the thicker active layer, then suprapermafrost flux will increase, but if water table elevation declines with an increase in soil drainage, then the enhancements in suprapermafrost transmissivity due to thickening may be countered by reductions in bulk K (Fig. 5). An alternative explanation for commonly observed positive trends in baseflow is wholesale permafrost loss spanning decades that leads to the enhancement of regional groundwater circulation and discharge to major rivers (Walvoord et al., 2012). However, the attendant increase in baseflow derived from the latter mechanism is also strongly influenced by potential and nonstraightforward changes in the water table elevation affected by thaw. Wholesale permafrost loss is more likely to contribute to baseflow increases in discontinuous permafrost than in continuous permafrost, where permafrost tends to be cold and thick.

Concurrent changes in water fluxes and flowpaths resulting from permafrost thaw influence the distribution and volume of water as soil moisture, lakes, wetlands, groundwater storage, and, in the winter, as aufeis and river ice. The projected trajectories of change in water distribution are less uniform than those of water flux (Fig. 4). Observational support for permafrost thaw–related changes in soil moisture, the distribution of lakes and wetlands, and groundwater storage is often inconsistent in terms of the magnitude and even direction of changes (trajectory information and supporting references in Table 1). This observed variability underscores the importance of location and scale of observation, as well as complexities in feedbacks and interactions between climatologic, thermal, ecologic, and hydrologic processes. A relevant research question currently posed at regional to pan-Arctic scales relates to identifying the areas that will become drier or wetter.

In addition to potential large-scale redistribution of near-surface water, changes in water stored in the form of ice during winter months are also expected. Aufeis, or icing, is layered ice that accumulates along streams and rivers as a result of groundwater discharge during freezing temperatures (Woo, 2012). Aufeis volume and spatial distribution may be expected to increase concomitantly with observed increases in baseflow, yet warming winter temperatures and the thermal energy from increased groundwater discharge may counteract, to some degree, expected aufeis build up. Observational studies of changes in aufeis distribution are limited, but show no substantial change (Yoshikawa et al., 2007). However, winter river ice thickness has been observed to be decreasing (Table 1), fueling concerns for winter transportation hazards. Also, river ice breakup, which normally occurs in the spring, can be triggered by mid-winter warm periods. These mid-winter warm periods have caused destructive ice jams and river flooding in temperate regions of Canada (e.g., Beltaos, 2002), and such events may become increasingly common at higher latitudes in the coming decades. Evaluating how much river ice thinning is due to warming from below via enhanced groundwater temperature and flux vs. increased heat exchange from above due to warming winter air temperature remains untested. Less is known about historical trends in river ice thickness than in river ice phenology because the latter is more readily derived from remote sensing techniques.

Surface–subsurface partitioning of water fluxes and distribution of surface–subsurface water storage are fundamentally linked to factors that affect ecosystem dynamics. Ecosystem variables sensitive to permafrost thaw-induced changes in hydrology include, but are not limited to, vegetation, hydrologic connectivity between inland waters, and seasonal variability in the stream hydrograph and thermograph (Fig. 4).

Observed changes in vegetation composition ascribed to changes in soil moisture and/or surface water conditions derived from permafrost thaw include large ecosystem shifts from birch forests to fens and bogs (Jorgenson et al., 2001), black spruce to bogs (Baltzer et al., 2014), shrub to graminoid dominance (Christensen et al., 2004), and other more subtle compositional shifts (Jorgenson et al., 2013). Also, impacts of climate change on vegetation density and health, such as the browning of boreal forests (Beck and Goetz, 2011; Verbyla, 2011), have been observed in high-latitude regions. However, variations in vegetation structure and composition are the result of complex interactions between climatic, hydrologic, thermal, and microbial conditions and disturbance factors (e.g., fire and insects), which make it extremely difficult to project how vegetation may respond in the future or even to assess trends that are primarily attributable to permafrost thaw effects.

Changes in hydrologic connectivity in permafrost lakes and/or wetland rich lowlands have been linked to ecosystem changes caused by permafrost thaw. Through changes to the landscape initiated by permafrost degradation and ground subsidence, the static surface hydrologic connectivity of the landscape can be notably augmented (Connon et al., 2014). Likewise, increases in subsurface hydrologic connectivity may be expected with shallow permafrost thaw, although minimal research has been conducted to date on this topic. River floods offer a means of enhancing functional hydrologic connectivity through increased water levels and the lateral transmission of water between inland waters. Ice-jam floods, occurring periodically during river ice breakup, are of particular importance in sustaining some cold region lowland lakes because water levels produced by ice jams are typically much greater than those caused by open-water floods (Beltaos and Prowse, 2009). Though hydroclimatic conditions play a primary role in controlling ice-jam flooding frequency and intensity (Goulding et al., 2009), we speculate that substantial increases in baseflow from permafrost thaw and the resulting thermal degradation of river ice (reducing pre-breakup ice cover thickness and strength) may exert some influence toward reducing the severity of ice-jam floods. However, the effects of potential reductions in ice jams severity in sustaining lowland lakes may be countered to some degree by the enhanced subsurface hydrologic connectivity that would result from permafrost thaw in permeable channels. Jepsen et al. (2015) described and modeled an expected mechanistic shift in floodwater propagation to lowland lakes from classic “fill-and-spill” typical of permafrost systems (Woo, 2012) to “fill-and-seep” associated with such enhanced subsurface hydrologic connectivity. Lowering of subsurface thresholds via permafrost table decline to form lateral suprapermafrost taliks between inland water bodies would allow for the subsurface propagation of lower intensity floods that may be expected due to climate warming effects described above. Such a shift toward increased subsurface hydrologic connectivity could bring about more synchronous water level changes between nearby lakes. Likewise, the enhanced water solute and nutrient exchange among inland waters may serve to reduce variability in lake chemistry and isotopic composition.

Reduction in the seasonal variability of streamflow magnitude and temperature can be expected with permafrost thaw (e.g., Ye et al., 2009; Liu et al., 2005) as proportionally more water is routed through the subsurface, which serves as a buffer, before entering stream networks. With this shift comes the potential for aquatic ecosystem changes. For example, enhanced and warmer groundwater discharge in winter months in southern permafrost regions may create larger and warmer in-stream thermal anomalies, such as those that are utilized by aquatic ectotherms for providing thermal refuge at lower latitudes during winter months (Cunjak and Power, 1986). These groundwater-generated thermal refugia could facilitate the winter residence of fish species currently constrained by the colder water temperatures in rivers and streams in permafrost environments.

Pan-Arctic permafrost (Fig. 1) spans a wide range in characteristics, thus resulting in considerable spatial variability in the hydrologic influence of permafrost, vulnerability to thaw, and potential hydrologic impacts of thaw. Permafrost characteristics relevant to these aspects include spatial coverage, ice content, temperature, depth to permafrost table and base, and potential magnitude of change in K on thaw.

Spatial coverage of permafrost is very broadly mapped for the pan-Arctic drainage basin (Brown et al., 2002), with higher resolution mapping available for smaller regions (e.g., Hegginbottom and Radburn, 1992). Coverage is typically mapped in general categories as: (i) continuous (90–100%), (ii) discontinuous (50–90%), (iii) sporadic (10–50%), and (iv) isolated (<10%). Although large-scale mapping indicates broad patterns of permafrost that fall into these categories, the actual distribution of permafrost, especially beyond the continuous permafrost zone, may be quite heterogeneous and complex, creating considerable scaling challenges. Some promising empirical approaches have been developed to generate baseline maps of permafrost probability at the local to regional scales using basal temperature of snow (Bonnaventure et al., 2012; Bonnaventure and Lewkowicz, 2013), though these approaches can be limited by methodology assumptions and availability of ground-truth information. Both the general classification and the typically low resolution of widely available permafrost spatial distribution maps are insufficient for robust permafrost hydrology modeling applications (as baseline input or calibration/validation) and thereby limit realistic model representation.

Volumetric ice content of permafrost in the upper 20 m at the pan-Arctic scale is also broadly mapped; general categories include high (>20%), medium (10–20%), and low (0–10%) (Brown et al., 2002). Ice content of permafrost is relevant for constraining rates of thaw and evaluating land deformation impacts of thaw. Ice-rich permafrost thaws more slowly than ice-poor permafrost due to the thermal inertia of the latent heat of fusion. However, upon thaw, ice-rich permafrost supports a greatly enhanced potential for ground subsidence and landscape disturbance compared with ice-poor permafrost. Ice content is primarily dependent on geologic properties of the permafrost and its depositional history. Syngenetic permafrost, formed synchronously with sediment deposition, tends to have high ice content relative to epigenetic permafrost, formed subsequent to sediment deposition (Jorgenson et al., 2010). As such, efforts to broadly map ice content have relied heavily on surficial geology mapping products combined with limited borehole reference data.

Permafrost temperature is also relevant to vulnerability assessments because permafrost just below 0°C is most susceptible to near-term thaw as minimal sensible heat is required to first raise the temperature to 0°C (Kurylyk and Hayashi, 2016). Although correlations between permafrost temperatures and air temperatures can be drawn, direct measurements of permafrost temperature are fairly limited, and long-term monitoring sites across the pan-Arctic basin are even more sparse (Romanovsky et al., 2010).

Characterization of the permafrost table is critical for many cold region studies due to seasonally dynamic thermal and hydrologic processes that occur in the active layer where subsurface biological, biogeochemical, and pedogenic activity is concentrated. As a consequence, mapping ALT has become a high research priority. The Circumpolar Active Layer Monitoring (CALM) network currently includes more than 125 sites globally (www.gwu.edu/∼calm, accessed 23 May 2016; Brown et al., 2000). These sites and other point ALT measurements commonly obtained via manual techniques, such as frost probing, have been used for extrapolation via empirical and statistical methods (e.g., Pastick et al., 2013). They also provide essential ground-truth information for calibration of equilibrium and transient thermal models to estimate ALT at large scales (Riseborough et al., 2008, and references therein).

Compared with ALT, much less information is known on the distribution of the permafrost base, and detection is restricted to areas of focused study where deep borehole and/or deep-seeing geophysical information is available (e.g., European Union PACE permafrost monitoring network; Harris et al., 2009). Similarly, geologic mapping of permafrost regimes typically consists of surficial geologic information and very deep geologic investigations (energy and economic explorations), with substantially less information relating to the depths that are of interest for hydrogeologic characterization. This intermediate depth information is critical for projecting changes to the subsurface flow system that would accompany warming. Subsurface characterization, including baseline permafrost distribution and the thawed hydrogeologic framework of high latitude systems, represents a major gap in efforts to project thaw rates and evaluate impacts. Permafrost characterization efforts span the spectrum from point-scale observations obtained via both invasive and noninvasive techniques to pan-Arctic-scale assessments inferred from satellite remote sensing approaches, the latter of which requires the former for calibration and validation. Along the spectrum of increasing spatial coverage comes a compromise in resolution and loss in ability to detect increasingly deep frozen ground features.

Geophysical methods, including ground-based and airborne approaches, fill an important niche in bridging the large gap in permafrost characterization between point-scale data and remote sensing analyses. A variety of geophysical methods (radar, electromagnetics, electrical resistivity, and seismic) have been developed to map permafrost and distinguish frozen–unfrozen ground transitions relying on the differences in the electrical properties and compressibility of ice and liquid water. The presence of ice vs. liquid water in the pore space of sediments affects contrasts in electrical resistivity, seismic and electromagnetic wave velocity, and dielectric permittivity. Due to inherent tradeoffs between data resolution and spatial coverage (Fig. 6a) and between data resolution and depth of investigation (Fig. 6b), many current permafrost studies using geophysics implement a multimethod strategy accounting for scale and features of interest.

Despite early, seminal work that first applied geophysical methods to map permafrost (e.g., Daniels et al., 1976; Hoekstra et al., 1975; Hoekstra, 1978), cold region geophysics remains an active area of research in method development, field application, data interpretation, and implementation of results into hydrologic modeling. A comprehensive review of geophysical approaches in permafrost terrain is beyond the scope of this paper, and readers are referred to broad reviews by French et al. (2006), Kneisel et al. (2008), Harris et al. (2009, section 5), and Hauck (2013). Yet, because permafrost characterization is such an integral component of permafrost hydrology studies and inadequate subsurface characterization contributes a great source of uncertainty in such analyses, we highlight select recent advances in geophysical applications to permafrost regimes.

Active layer characterization is perhaps the most common target for cold regions geophysical approaches, and several shallow-seeing methods are applicable (Fig. 6). Westermann et al. (2010) describe an approach using multi-channel ground-penetrating radar (GPR) to determine ALT and monitor the seasonal evolution of thaw depth, as well as the average water content in the soil profile in Svalbard. Several recent studies have combined GPR with other geophysical tools, including electromagnetic induction (EMI), electrical resistivity tomography (ERT), and seismic refraction (SR) to optimize the characterization strengths of each method and enhance subsurface interpretation (e.g., Langston et al., 2011; Hubbard et al., 2013; Sjöberg et al., 2015). In general, GPR provides a useful approach to capture ALT with good spatial coverage, and EMI also shows promise in this direction (Briggs et al., 2016). Permafrost investigations that require information at intermediate depths benefit from the use of ERT surveys. Though setup is labor intensive, ERT transects yield resistivity cross-sections that can be used to infer detailed two-dimensional information on permafrost morphology and lithology. Close electrode spacing along ERT transects allows resolution of ALT for comparison with GPR and frost probe estimates, whereas more distant electrode spacing allows for deeper imaging of the subsurface; thus, a nested approach for electrode spacing in ERT surveys may be employed to optimize resolution of ALT and support deeper permafrost characterization (Briggs et al., 2016). A range of depths for permafrost characterization is relevant for vadose zone and groundwater hydrology studies. Though less popular than other methods mentioned here, SR is another approach shown to be valuable for detecting changes in ground ice over time (Hilbich, 2010).

New research shows that enhanced information on water mobility near the permafrost table may be obtained in the field through the combined use of ERT and in situ downhole nuclear magnetic resonance (NMR) (Minsley et al., 2016). Nuclear magnetic resonance is directly sensitive to liquid water content, as mobile, capillary, and bound water. Nuclear magnetic resonance results indicate notable liquid water content persisting below the permafrost interface (presumably the 0°C isotherm) detected via ERT and frost probing. This characterization is relevant in light of recent results that demonstrate the potential for rapid mineralization of dissolved organic carbon and release of carbon dioxide following permafrost thaw (Drake et al., 2015).

Appropriate geophysical methods for detecting the bottom of permafrost depend on the target depth, which is site specific. In relatively thin permafrost (<10–20 m), ERT may be able to capture the interface between permafrost and unfrozen sediment beneath. For thicker permafrost, deeper seeing ground-based methods such as time-domain electromagnetics (TDEM), which offers one-dimensional vertical information, are required.

In general, ground-based geophysical techniques are time intensive and spatially limited. Geophysical data collection via airborne platforms offers a means of increasing spatial coverage and accessing remote and challenging terrain. For example, airborne electromagnetic (AEM) surveys are useful for characterizing permafrost distribution on the order of 10³ to 10⁴ km² in sedimentary basins to depths of 50 to 100 m (Minsley et al., 2012). Recent forward and inverse modeling analysis by Minsley et al. (2015) highlighted the utility and limitations of AEM in the interpretation of taliks beneath lakes in discontinuous permafrost landscapes. Determining whether open or closed taliks exist below surface water bodies is a critical characterization component required for assessing deep subsurface hydrologic connectivity and the potential for groundwater–surface water exchange. Direct measurements of liquid water content in the subsurface (i.e., NMR) can substantially reduce the uncertainty associated with interpreting resistivity images below surface water bodies. A recent novel study conducted near Fairbanks, Alaska, USA demonstrated the value of using surface NMR, together with TDEM resistivity data, to determine talik dimensions beneath thermokarst lakes (Parsekian et al., 2013).

Unlike other components of the terrestrial cryosphere, permafrost is inherently a subsurface feature that is unable to be directly detected via current remote sensing capabilities used to infer near-surface properties. Therefore, permafrost mapping efforts utilizing large-scale remotely sensed data typically employ empirical and statistical approaches to extrapolate point-scale measurements or geophysical survey estimates of ALT (e.g., Shiklomanov and Nelson, 1999; Pastick et al., 2013, 2014). Extrapolation of permafrost characteristics is particularly challenging due to limited field measurements in the remote and vast regions of the pan-Arctic basin. In addition, the resolution of satellite remote sensing data may be incompatible with small-scale heterogeneity in permafrost characteristics that are common, particularly in discontinuous permafrost landscapes. Techniques to improve permafrost models derived from remotely sensed data sets are evolving, yet much research is needed to advance toward the ultimate goal of establishing permafrost baseline information against which future change can be assessed (National Research Council, 2014). A recent review of permafrost remote sensing techniques to (i) identify surface cryologic features linked to permafrost and (ii) map physical variables directly or indirectly related to subsurface thermal conditions can be found in Westermann et al. (2015).

Some emerging remote sensing methods for permafrost mapping include airborne P-band synthetic aperture radar to estimate ALT (Tabatabaeenejad et al., 2015) and Interferometric Synthetic Aperture Radar (InSAR) to measure seasonal subsidence from thawing of pore ice in the active layer and to infer ALT (Schaefer et al., 2015). Current methods of mapping permafrost using satellite remotely sensed data apply best to near-surface permafrost (in the upper 1 m), and methods to estimate ALT using airborne sensed data are most applicable for continuous permafrost with shallow ALTs (i.e., High Arctic). Considerable demand persists for high-spatial resolution permafrost mapping capabilities that are applicable for sub-Arctic regions where discontinuous to sporadic permafrost coverage is typical and where the permafrost table may be greater than 1 m below the ground surface. In addition to ALT and spatial distribution, mapping permafrost ice content is of great value, because ice content influences thaw vulnerability, rate of thaw, and geomorphic and hydrologic impacts. To date, most large-scale permafrost ice-content mapping has relied on identification of surface features, such as ice wedge polygons, to infer information of ground ice content in near-surface permafrost. Continued research efforts to develop and apply cross-scale permafrost characterization techniques that combine ALT, spatial, vertical, and ice-content information from geophysical surveys and airborne and/or satellite remote sensing data with thermal permafrost models may be the most viable path forward to address scaling concerns.

Hydrological modeling in permafrost environments presents challenges that are seldom addressed in permafrost texts or review papers. Because of the influence of the ground thermal regime on water storage and transmission in permafrost basins, many model developers have incorporated heat transfer equations into physically based hydrological models. Herein, we focus on these heat transfer approaches rather than the water routing algorithms per se because explicitly accounting for the ground thermal regime is the distinctive challenge of permafrost hydrological modeling.

This approach invokes many assumptions, including no horizontal heat transfer, constant moisture content (latent heat) with depth, uniform thermal conductivity, negligible soil heat capacity, and no heat advection. However, the simplicity of the Stefan equation facilitates modifications to relax the limiting assumptions, and recent studies have demonstrated how the equation can be modified to accommodate temporally variable soil moisture conditions (Hayashi et al., 2007), spatially variable moisture content and thermal properties (Kurylyk, 2015), two-directional freezing and thawing (Woo et al., 2004), heat advection (Kurylyk et al., 2014b), and soil heat capacity (Kurylyk and Hayashi, 2016). Due to the simplicity and flexibility of the Stefan equation, variations have been incorporated in permafrost hydrology models (Fox, 1992; Carey and Woo, 2005) and land surface schemes (Li and Koike, 2003; Yi et al., 2006). Semi-empirical methods based on the Stefan equation (e.g., Nelson et al., 1997) have also been used to predict soil freeze–thaw, but these approaches require data for calibrating equation coefficients and lack spatiotemporal transferability.

One difficulty with the Stefan equation and other approaches for determining soil freeze–thaw is that the ground surface temperature is a required input, but often only meteorological data are available at high latitudes or altitudes. Air and ground surface temperature can be strongly decoupled, especially under an insulating snowpack (Zhang, 2005). Cold regions engineers and scientists have long employed the empirical n-factor, which is essentially the ratio of the temporal integral of the surface temperature (e.g., the degree-days) to the temporal integral of the air temperature during either a thawing or freezing period (Lunardini, 1981). This n-factor can exhibit interannual variability (Juliussen and Humlum, 2007) and is expected to change with climate warming and concomitant snowpack changes. Thus, others have proposed quasiempirical means to determine ground surface temperatures in permafrost regions from meteorological data and an understanding of the surface energy balance (Hwang, 1976; Williams et al., 2015).

The Kudryavtsev et al. (1977) equation is an alternative analytical solution to determine the depth of seasonal freezing or thawing (Romanovsky and Osterkamp, 1997), which performs better than the standard, unmodified Stefan equation because it accounts for soil thaw retardation caused by the heat capacity of the soil. This approach also requires the determination of the temperature at the top surface of the permafrost (T_TOP), which can be calculated using the T_TOP model approach proposed by Smith and Riseborough (1996). The Kudryavtsev et al. (1977) equation is not presented here due to the complexity of the approach and the number of equations. Further details were given by Bonnaventure and Lamoureux (2013). Alternatively, the modified Stefan equation proposed by Kurylyk and Hayashi (2016) can be used to account for the influence of soil heat capacity on soil freeze–thaw. Results from this method, which does not require the prior calculation of T_TOP, have been shown to be very similar to those obtained using the more complex Kudryavtsev et al. (1977) approach (Yin et al., 2016).

Finally, Semenova et al. (2014) described an alternative analytical algorithm to calculate vertical soil freeze–thaw. This approach allows for phase-dependent snow and soil thermal conductivity. To date, this algorithm has only been incorporated into the hydrology model Hydrograph as described by Lebedeva et al. (2014).

In general, the computational savings achieved by incorporating analytical, rather than numerical, approaches for soil freeze–thaw into a hydrology model can be substantial, particularly for global-scale simulations.

Other hydrological models have incorporated numerical solution schemes (e.g., finite difference, finite element, or finite volume) to model ground freezing and thawing, particularly when simulating one-dimensional infiltration intro frozen soils (e.g., Zhao et al., 1997). Recent examples of distributed hydrological models employing numerical soil freeze–thaw schemes include Gouttevin et al. (2012), Rawlins et al. (2013), Zhang et al. (2013), and Clark et al. (2015). These powerful numerical approaches can accommodate complex processes or conditions, including soil heterogeneities, coupled heat and water transfer, complex temperature boundary conditions, intermittent freeze–thaw, and temporally variable thermal properties. Properly parameterized numerical approaches typically perform better than analytical approaches for predicting ground freeze–thaw (Zhang et al., 2008b), but there is a tradeoff between model fidelity and computational expense. Table 2 presents results from several modeling studies that have considered the influence of climate change on surface and near-surface hydrological processes in permafrost basins.

One primary advantage of numerical models is that freezing and thawing processes can occur over a range of temperatures rather than the sharp phase change assumed by analytical solutions. This freezing range allows for some subsurface water flow at sub-zero temperatures. However, it can be difficult to properly partition between pore water and pore ice based on the temperature and water (or ice) potential, and this partitioning is even more complex in unsaturated soils (Kurylyk and Watanabe, 2013). Often the “freezing = drying” approach is employed, which assumes that the remaining unfrozen water content in partially frozen soils can be determined in an analogous manner to how water content in drying soils is obtained via the soil water characteristic curve (Koopmans and Miller, 1966; Spaans and Baker, 1996). The soil freezing curve (relationship between temperature and unfrozen water content) and soil characteristic curve are related through a form of the Clapeyron equation (Kurylyk and Watanabe, 2013).

Despite the recent advances in permafrost hydrological modeling, there are still a number of challenges and unresolved questions. Five examples are provided herein, although this list is not comprehensive. First, all of the analytical soil freeze–thaw solutions and most of the numerical soil freezing models that have been coupled to hydrological models assume that heat transfer is restricted to the vertical dimension. However, lateral heat transfer to isolated permafrost bodies can be important in lowland discontinuous permafrost environments (McClymont et al., 2013; Kurylyk et al., 2016; Sjöberg et al., 2016), as well as steep, alpine environments (Noetzli et al., 2007; Noetzli and Gruber, 2009). Thus, the one-dimensional approaches described above are better suited for simulating seasonal freeze–thaw rather than multidecadal thawing of permafrost bodies and resultant hydrologic changes. For example, Connon et al. (2014) demonstrated that pronounced increases in streamflow in the Northwest Territories, Canada likely arose from changes to the hydrologic connectivity and contributing area of a watershed due to a reduction in the lateral extent of permafrost rather than ALT increases. A hydrological model only considering vertical heat transfer would not likely reproduce these multidecadal and multidimensional dynamics.

Second, many of the land surface schemes and hydrological models employing a version of the Stefan equation have not incorporated the equation modifications noted above that accommodate soil layering, heat capacity, and other factors that influence the rate of soil freeze–thaw. Neglecting these complicating factors can yield considerable errors in the calculation of the frost or thaw depths and the associated impact on subsurface water storage and routing.

Third, in most numerical models of soil freeze–thaw, thermodynamic equilibrium is tacitly employed via the utility of the Clapeyron equation to represent the interplay between water potential, temperature, and phase change. However, disequilibrium freezing and thawing can occur when temperatures change too quickly to allow for equilibrium ice formation or melt. Also, disequilibrium pressure can occur during snowmelt infiltration into partially frozen soils because ice content can increase despite the constant temperature. To our knowledge, no hydrologic model has attempted to represent disequilibrium phase change processes, and this represents an ongoing challenge.

Fourthly, predicting the ground surface temperature (the drivers of seasonal soil freeze–thaw) under snowpack and in snow-free conditions remains a persistent challenge in permafrost hydrology modeling. The common approach of employing the simplistic n-factor should be replaced by more physical or quasi-physical approaches (e.g., Williams et al., 2015), but these can be restricted to snow-free conditions. Accurately modeling ground temperatures can still be challenging even when applying more complex, physically based models as simulated ground thermal regimes can be extremely sensitive to model parameters (Gubler et al., 2013; Wu et al., 2016; Harp et al., 2016).

Finally, further research is warranted to mathematically represent the relationship between the soil K and air, ice, and liquid water contents in partially frozen, unsaturated soils. The K is normally represented with a relative K function similar to that used for unfrozen, unsaturated soils. However, considerable uncertainty remains as to whether K of partially frozen soils should also be further decreased through the means of an empirical impedance factor and whether dual porosity approaches should be alternatively employed (Kurylyk and Watanabe, 2013). These unknowns could be addressed through carefully designed laboratory experiments to better constrain soil freezing curves and the K of partially frozen soils for different soil types.

These five noted ongoing hydrological modeling challenges relate to subsurface conditions. Other aspects of hydrologic and hydraulic modeling that are unique to cold regions also merit consideration. For example, most currently available water routing algorithms for river channels are not well suited to permafrost regions. Given the high uncertainty in river ice-jam frequency, intensity, and impact on floodplain hydrology in permafrost settings, improvements in these routing algorithms are warranted. Applying flexible, modular hydrology simulators such as RAVEN (RAVEN Development Team, 2016) or the Cold Regions Hydrological Model (Pomeroy et al., 2007) may facilitate the incorporation of new routines designed specifically for permafrost environments.

In the early 1970s, a number of one-dimensional coupled groundwater flow and energy transport models that considered freezing and thawing began to emerge. Until recently, such simulations were restricted to the vertical direction; however, in the past 10 yr, a large number of multi-dimensional models have been developed to investigate groundwater flow in permafrost regions (e.g., McKenzie et al., 2007; Bense et al., 2012; Endrizzi et al., 2014; Karra et al., 2014; Frederick and Buffett, 2015). These models are typically formed by coupling a three-dimensional Richards-type equation for water flow to a three-dimensional heat transfer equation considering heat conduction, heat advection, thermal dispersion, and pore water phase change (Kurylyk et al., 2014a). These powerful simulators, known as cryohydrogeology models, consider the influence of the latent heat of pore water phase change on the subsurface effective heat capacity and also simulate the reduction in K due to pore ice formation.

Cryohydrogeology models have been applied in a range of scenarios to simulate the influence of climate warming on groundwater flow in permafrost basins (Table 3). In general, these studies have demonstrated that increased ALT may enhance groundwater storage during the recharge season and thus increase baseflow and decrease streamflow seasonality. Several of these studies have particularly focused on the role of groundwater flow in taliks and have demonstrated that open talik breakthrough can substantially enhance groundwater exchange between supra- and subpermafrost aquifers (Fig. 3).

In the past, these models have often been developed by independent teams of researchers. However, many of these codes have become recently integrated in the international cold regions groundwater modeling network known as InterFrost (Grenier et al., 2015), which is predominantly focused on model testing and developing inter-code comparisons (Kurylyk et al., 2014b; Rühaak et al., 2015).

A number of challenges or limitations associated with the development or application of cryohydrogeology models remain in addition to the already noted limitation of accurately representing the K of partially frozen soils. First, studies employing these models (Table 3) have typically been restricted to idealized environments. To our knowledge, few studies have been produced that attempt to reproduce field conditions (e.g., multidecadal permafrost thaw or groundwater flow rates) (Atchley et al., 2015; Kurylyk et al., 2016; Sjöberg et al., 2016). This is at least partially due to model parameterization difficulties arising from a scarcity of hydrogeologic data at high altitudes or latitudes. We expect that these models will begin to be tested in well-instrumented catchments to ground simulations in reality.

Second, most of the models listed in Table 3 do not include a land surface scheme; the Arctic Terrestrial Simulator (Atchley et al., 2015) and GeoTop2.0 (Endrizzi et al., 2014) are notable exceptions. Due to the lack of consideration of surface hydrologic and thermal processes, the boundary conditions to drive most of these models must be specified subsurface conditions (e.g., groundwater table location or recharge rates and soil temperature). Such data are seldom available, and this has led modelers to employ simplified boundary conditions, such as assuming ground surface temperature is equal to air temperature or that the groundwater table follows the ground surface topography. In the absence of a coupled land surface scheme, an alternative approach is to first run simulations in a surface model, and then employ the output from the surface model (e.g., groundwater recharge and near-surface soil temperature) as the boundary conditions for the hydrogeology model (Kurylyk et al., 2016).

Third, a remaining limitation of cryohydrogeology models is that simulations are typically restricted to two dimensions, simple structure and geometry, and relatively small spatial scales due to the computational resources required to obtain the numerical solution. The equations are highly nonlinear and generally demand small time steps and a fine mesh or grid. This challenge may be overcome by implementing codes with massively parallel solvers, although this presents further challenges (Painter et al., 2013). Computing capabilities are continuously improving, and we expect that these models will be increasingly applied in three dimensions and for larger and more structurally complex systems.

A final limitation of available permafrost hydrogeology models relates to the difficulty of their utility. As noted by Fritz et al. (2015), many permafrost heat transfer modeling tools are “too complex to be used by anyone other than modeling experts.” This limitation may help explain the lack of simulations conducted for site-specific application. Concerted efforts to provide more frequent and detailed model training workshops for interested researchers and develop more user-friendly model interfaces may help overcome this limitation.

Water resources and the ecosystems they support are particularly vulnerable to climate change in permafrost environments because (i) observed and projected climate warming rates are high in northern latitudes and high altitudes, (ii) air temperatures and hydrological processes are linked through the process of ground freezing and thawing, and (iii) substantial alteration in the hydrogeologic framework of regions may result from permafrost degradation. Recent climate warming has altered hydrologic conditions and processes in the pan-Arctic basin, including streamflow magnitude and seasonality, aquifer activation, and the spatial extent and distribution of wetlands and lakes. Although phenomenological explanations of the link between permafrost thaw and altered hydrological processes abound, we currently have limited capability to predict how future warming may influence local, regional, and pan-Arctic scale hydrology. This review identified hydrologic components that are expected to change or have changed with permafrost thaw but lack observational data to support such inferences (“data gaps” noted in Table 1). These components possess major data gaps in part due to logistical challenges in being able to measure and monitor subsurface water fluxes and storage in permafrost environments. Continued development toward subsurface hydrologic and physical characterization coupled with efforts to integrate field data into cryohydrogeologic models will improve our understanding of expected trajectories of change for various types of permafrost systems.

High-resolution characterization of baseline permafrost conditions and monitoring data represent critical components for assessing projected thaw rates and potential hydrologic consequences. Patterns of permafrost distribution in discontinuous permafrost are complicated by ecosystem-protected processes, thereby requiring advancements in refining cross-scale permafrost characterization and gaining knowledge of the complex factors controlling near-surface ground temperatures. Multi-method geophysical approaches show promise in cross-scale efforts to bridge point data and satellite sensor data. Recent advances have led to unprecedented knowledge of the spatial and vertical extent of permafrost as well as liquid–ice-content information. Most of the uncertainty in geophysical approaches to distinguish frozen–unfrozen ground interfaces is derived from the underdetermined inverse problem and from data and image interpretation. Multiple factors, in addition to ice and water presence, influence geophysical properties of the subsurface. Because lithology also influences electrical and electromagnetic properties, information about the region’s geologic framework is imperative for accurate interpretation of permafrost distribution using geophysical electrical methods. Research opportunities aimed at reducing uncertainties and ambiguity in the interpretation of geophysical surveys in permafrost include: (i) experimental laboratory studies to better quantify the geophysical properties of freezing soil, (ii) joint inversion and interpretation methods, and (iii) modeling approaches, such as error models and appraisal methods, to constrain inversions and assess reliability of interpretations (Hauck, 2013). In addition, efforts to characterize permafrost systems will require continued integration among various geophysical methods and with remote sensing approaches to assimilate information across scales with variable resolution. Adequate ground-truth information for calibration of remote sensing data remains a critical and typically underprovided component in cold regions.

This review identified several studies employing physically based models to consider how climate-induced changes in the thickness and lateral extent of permafrost may influence surface and near-surface hydrological processes (Table 2) and the dynamic interaction between permafrost and groundwater flow (Table 3). Very few of the cryohydrogeology modeling studies were conducted for alpine permafrost environments, which represent the headwaters for rivers that supply fresh water for much of the world’s population (Viviroli et al., 2007). The current gaps in model applications represent significant opportunities for cold regions hydrologic modelers. Potential model developments and applications include linking simulations with long-term field observations of permafrost and hydrologic changes, developing improved representation of hydraulic properties of frozen soils, and incorporating land surface schemes in cryohydrogeology models to simulate surficial water and energy balances. These research directions have the potential to increase confidence in simulation results and thus improve our capability for predicting how future warming will affect water fluxes and distribution, and thereby mediate the permafrost–carbon feedback and influence ecosystems.

We thank VZJ co-editor David Stonestrom (USGS) for the invitation to provide this review. We also thank Fred Day-Lewis (USGS) for the conception of Fig. 6 and for offering insightful comments on the manuscript. Three anonymous reviewers and the associate editor, Masaki Hayashi, provided useful suggestions for improving the manuscript. B. Kurylyk was funded by postdoctoral fellowships provided by the Natural Sciences and Engineering Research Council of Canada, the Killam Trust, and the University of Calgary Eyes High program. M. Walvoord gratefully acknowledges support from the USGS National Research Program and the NASA ABoVE campaign.