Development of perennial thaw zones in boreal hillslopes enhances potential mobilization of permafrost carbon

Permafrost thaw alters subsurface flow in boreal regions that in turn influences the magnitude, seasonality, and chemical composition of streamflow. Prediction of these changes is challenged by incomplete knowledge of timing, flowpath depth, and amount of groundwater discharge to streams in response to thaw. One important phenomenon that may affect flow and transport through boreal hillslopes is development of lateral perennial thaw zones (PTZs), the existence of which is here supported by geophysical observations and cryohydrogeologic modeling. Model results link thaw to enhanced and seasonally-extended baseflow, which have implications for mobilization of soluble constituents. Results demonstrate the sensitivity of PTZ development to organic layer thickness and near-surface factors that mediate heat exchange at the atmosphere/ground-surface interface. Study findings suggest that PTZs serve as a detectable precursor to accelerated permafrost degradation. This study provides important contextual insight on a fundamental thermo-hydrologic process that can enhance terrestrial-to-aquatic transfer of permafrost carbon, nitrogen, and mercury previously sequestered in thawing watersheds.


Introduction
In boreal watersheds underlain by permafrost, most hydrologic, biological, biogeochemical, microbial, and pedogenic activity occurs in the active layer that undergoes seasonal freeze/thaw (Hinzman et al 2005). Permafrost inhibits drainage, contributing to activelayer conditions characterized by cold, saturated soils with large organic carbon stores (Grosse et al 2011). State-of-the-art permafrost models predict considerable spatial and vertical degradation across much of the circumboreal region by the end of the 21st century (Lawrence and Slater 2005, Lawrence et al 2008, Jafarov et al 2012 that will impact the storage of soil carbon (McGuire et al 2018). In response to surface warming, depth of thaw eventually exceeds depth of seasonal frost, and a perennial thaw zone (PTZ), or supra-permafrost talik, develops. Though vertical talik development beneath surface-water bodies is well established through theory and observations (Plug and West 2009, Minsley et al 2012, Wellman et al 2013, Parsekian et al 2013, Roy-Leveillee and Burn 2017 and has been linked to enhanced stream discharge from sub-permafrost flow , Jafarov et al 2018, lateral taliks in terrestrial landscapes and their potential role as conduits for groundwater and solute transport are less studied. Early work focused on characterizing large lateral taliks associated with aufeis (icings) formation (Boikov et al 1984). More recent observational studies suggest the presence of thin perennial partial or fully thawed zones, unconnected with aufeis occurrence or other apparent surficial expression, which raises questions about the prevalence of such features in permafrost environments in response to warming. For example, interpolation of historical soil temperatures at sites in Russia suggest a concomitant permafrost-table deepening and seasonal frost shallowing supportive of PTZ development in recent decades (Streletskiy et al 2015). Likewise, comparison of permafrost and frost depth measurements suggests thin supra-permafrost taliks that have developed beneath forested plateaus in subarctic Canada (Connon et al 2018). Presently, these types of deep observational data are rare. Characterization of PTZs and determination of their prevalence remain important gaps in permafrost mapping (Walvoord and Kurylyk 2016) with implications for (1) fate and transport of vast stores of permafrost carbon (C) (Hugelius et al 2014) and mercury (Hg) (Schuster et al 2018), the former of which may be chemically labile and subject to decomposition even at temperatures 0°C (Waldrop et al 2010), and (2) infrastructure instability from changes in mechanical strength of thawing ground (Instanes and Anisimov 2016). Current state-of-theart remote sensing permafrost mapping technologies have depth limitations (typically <1 m) (National Research Council 2014, Pastick et al 2015 that impede PTZ identification. Ground-based geophysical field observations of PTZs, including novel results from borehole nuclear magnetic resonance (NMR) data, motivated the generalized modeling analysis and are presented here. This study addresses, from a theoretical perspective, (1) conditions conducive for PTZ development, (2) quantification of thaw rate, acceleration, and sensitivity to influential physical factors, and (3) effect of thaw development on baseflow magnitude, seasonality, and flowpaths.

Groundwater flow in permafrost systems
There is growing recognition of enhanced groundwater circulation, storage, and contribution to streamflow in thawing permafrost landscapes (Smith et al 2007, Walvoord and Striegl 2007, St Jacques and Sauchyn 2009, Velicogna et al 2012, Duan et al 2017. Shifts in water partitioning from surface to subsurface flow can substantially alter stream chemistry and chemical exports (Striegl et al 2005, Frey and McClelland 2009, O'Donnell et al 2012, O'Donnell et al 2014, Vonk et al 2015, Toohey et al 2016. Previous groundwater modeling studies indicate that thickening of the seasonally active (thawed) layer enhances groundwater discharge to streams  (2015)). The modeling analysis presented here extends prior work by addressing permafrost thaw induced changes, which include the development of PTZs, in baseflow seasonality and flowpath depth. Due to large uncertainties in the fate and transport of C (as gaseous or dissolved phase) from thawing permafrost (Lawrence et al 2015, Abbott et al 2016 and potential positive feedback to atmospheric warming (Schuur et al 2015), results are here viewed in the context of mobilization potential, as freeze/thaw state and export potential, for permafrost C. Recent work demonstrates the larger potential yield of dissolved organic carbon (DOC), total dissolved nitrogen (TDN), and mercury (Hg) from near-surface Holocene permafrost soils relative to active-layer soils (Wickland et al 2018, Schuster et al 2018, thereby elevating the importance of understanding subsurface hydrologic fluxes in thawing permafrost landscapes. Though abiotic and biotic processes affect the terrestrial to aquatic delivery of DOC, TDN, and Hg and require consideration for comprehensive study, the approach presented here offers insight on hydrologic processes and subsurface conditions that influence the seasonality of waterborne transport in permafrost-bounded hillslopes. Results highlight the need for coupled geochemical and cryohydrogeologic approaches in integrated observational and modeling studies to constrain the fate of globally important sources of C, N, and Hg, currently sequestered in permafrost. , acquired in September 2016 shows three distinct zones: low resistivity values characteristic of deeply thawed ground beneath both (1) a shallow tributary and (2) a small burned area from a wildfire that occurred in 2004 with a 2-4 cm thick organic layer, and high resistivity values indicative of frozen ground beneath (3) unburned black spruce forest with a 6-12 cm thick organic layer. Comparison of thaw characteristics in the latter two zones are consistent with the concept that thicker organic soil provides enhanced insulation and protection from permafrost thaw (Jafarov et al 2013, Minsley et al 2016. Although ERT measurements enable high-resolution mapping of permafrost, these data provide only a single snapshot in time near the point of maximum thaw and do not directly quantify the amount of liquid water.

Geophysical evidence for PTZs
Repeat NMR measurements collected in September 2016 and May 2017 were used to quantify seasonal changes in in situ liquid water content with depth in the burned and unburned areas. NMR data at the burned site (figure 1(c)) show seasonal changes in shallow water content in the upper ∼75 cm of the active-layer, beneath which a PTZ is evident from ∼87-200 cm depth, and a return to partly frozen conditions at 225 cm depth. In September 2016, measured temperature at this location ranged from +3.8°C at 111 cm depth to 0.0°C at 231 cm depth. At the unburned site (figure 1(d)), no significant thaw at depth is observed in either season, consistent with ERT results and characteristics typical of ecosystemprotected permafrost. While wildfire was the primary mechanism leading to conditions that favored PTZ development in this example, numerous other climatic and environmental factors can lead to PTZ formation, several of which are identified and discussed in this study.

Model methods
To examine active-layer dynamics and PTZ development broadly applicable to boreal hillslopes, coupled fluid-flow and energy transport is simulated through various hillslope conditions basally bounded by permafrost and subjected to seasonal air-temperature variation with a long-term warming trend superimposed on the seasonal cycle. Model results are analyzed to quantify seasonal shifts in groundwater flux attributed to warming and thaw. In contrast with a strict cryotic definition of permafrost (ground 0°C for at least two consecutive years) that translates to a sharp boundary, PTZs represent a continuum of thaw states and lack an established quantifiable definition. PTZ results are therefore presented in terms of liquid water (or ice) saturation in the porespace at seasonal maximum frozen conditions (i.e. at the onset of topdown active-layer thaw). Complexities in subsurface properties, freezing functions, and surface-energy balance are kept to a minimum in this analysis to highlight fundamental processes associated with permafrost dynamics and seasonal groundwater discharge relevant to permafrost DOC, TDN, and Hg potential mobilization.

Model description
An enhanced version of the US Geological Survey numerical code, SUTRA (Voss and Provost 2002), is employed here. The enhanced SUTRA code simulates coupled fluid-flow and heat transport incorporating the nonlinear phase change between ice and liquid water and ice-liquid phase transitions account for latent heat of fusion, accompanied by changes in permeability, thermal conductivity, and heat capacity depending on ice/liquid saturation. Details, including governing equations, have been described previously The model domain represents a generalized twodimensional 100 m long fully-saturated crosssectional hillside that slopes downward to the left from ERT transect from September 2016 shows clear differences in permafrost thaw in areas impacted by surface water and wildfire compared with undisturbed black spruce forest. Borehole NMR data indicate seasonal differences in in situ liquid water content at depth, with a PTZ evident at the burned site (c) that is not present at the unburned site (d).
the hill top (at elevation 5 m) to a small stream (at elevation 0 m) (figure S1 is available online at stacks.iop. org/ERL/14/015003/mmedia). A gauge pressure of zero Pa is specified along the top boundary to represent atmospheric pressure at the water table and at the intersection with a stream. The stream is assumed to be centrally located in a symmetric valley, implying no-flow fluid and fully insulating energy boundary conditions below the stream along the left side of the modeled domain. The hill crest represents the location of another symmetry boundary, with no-flow fluid and fully insulating energy conditions on the right side of the modeled domain. The model bottom is set at a large depth of 100 m to minimize the influence of basal boundary conditions on shallow permafrost dynamics, which are the focus of this study. A no-fluid flow condition and a constant upward geothermal energy flux with a typical value, 0.04 W m −2 , are applied along the basal boundary. Temperature is specified at the top boundary with time-varying values internally calculated via programming of SUTRA subroutine BCTIME by specifying a seasonal air-temperature function and applying 'n-factors' (described below) that represent buffering between air and groundsurface temperatures. Vertical discretization is fine (0.05 m) near the top to adequately resolve freeze/ thaw dynamics and coarsens with depth. Lateral discretization is constant (1 m; results for 0.1 m lateral discretization were comparable). Here, a simplified, yet flexible, approach is employed to capture variability in factors that affect surface-energy balance. N-factors, the ratio between cumulative ground-surface temperature and air-temperature degree days, are applied to compute ground-surface temperature from specified air-temperature (T air ) (Lunardini 1978). N-factors are site-specific. Values closest to 1 represent no difference between air and ground-surface temperatures; values approaching 0 indicate progressively more thermal buffering. Separate n-factors for freezing (n F when T air <0°C) and thawing (n UF when T air 0°C) periods are used due to the strong insulating effect of snow during the freezing period; thus typically, n F <n UF . To evaluate sensitivity of mean annual ground-surface temperature (MAGST) to n F and n UF and, by proxy, address the influence of snow, vegetation cover, shading, and other surface-energy factors, MAGST is calculated by applying T air seasonality for combinations of various values of n UF , n F , and mean annual air-temperature (MAAT). n UF is varied from 0.8 to 0.4, and n F is varied between 0.2 and 0.6, all at 0.1 increments, to reflect a range of literature-derived values Kreig 1988, Karunaratne andBurn 2003). MAAT is varied between −5°C and 0°C at 1°C increments, representative of conditions found in much of Alaska ( figure 1(a)).

Model simulations
The basecase simulation incorporates typical subsurface properties for boreal Alaska (table S1) and is used for comparison to other cases (table S2). The horizontal two-layer hydrogeologic configuration consists of organic soil (0.1 m thick for basecase) overlying mineral soil. Initial conditions for all cases reflect a dynamic steady-state condition (simulated as a repeated yearly cycle for 1000 years) for a sinusoidal yearly T air variation with −5°C MAAT and 20°C amplitude. Ground-surface temperatures are moderated from T air by n UF of 0.6 when T air >0°C and n F of 0.3 when T air 0°C, consistent with setting-appropriate values (e.g. Kreig 1988, Karunaratne andBurn 2003). Transient simulations are run from initial conditions for 100 years at 0.01 year timesteps. For times >0 years, a trend of +0.05°C yr −1 is superimposed on the seasonal air-temperature cycle. This warming trend is compatible with CMIP5 RCP8.5 trajectory for atmospheric warming projections that assumes no climate mitigation measures and is therefore a high-end constraint (Stocker et al 2015). In recognition that −5°C MAAT is not a universally appropriate present-day condition for hillslopes containing permafrost, the 1°C MAAT increase per 20 year specification additionally facilitates a convenient translation between simulated time and current MAAT for a given study area. For example, simulated results in the 40-60 year timeframe may best apply to current conditions near Fairbanks, Alaska, where MAAT=∼−2.5°C. Such translation is applicable for these simulations as surface-warming rates do not outpace subsurface equilibration times.
Although numerous variations on the basecase were conducted as part of a larger scoping study, only a small subset of critical simulations is presented here to concisely convey the most relevant and broadly applicable thermal and hydrologic responses observed. Table S2 and text S1 describe select cases and provide justification for the range of values evaluated.

N-factor matrix
MAGSTs are presented for combinations of n UF , n F , and MAAT ( figure 2(a)). Values in each of the six n-factor matrices relate to MAAT ranging from −5°C to 0°C in 1°C increments. Advancement from one matrix to the one below represents a warming progression through time or space, allowing for flexible interpretation. The large (∼5°C) range in calculated MAGST within each single-value MAAT matrix highlights the importance of factors that moderate air-toground temperature exchange. Boxes shaded in light red and light blue represent a zone of thermallysensitive conditions where MAGST falls between −1°C and 1°C. Results indicate conditions that favor positive MASGT even for MAAT<0°C. Results emphasize how subtle changes in variables that determine n-factors (for example, vegetation type, snow, shading, aspect) can shift the MAGST above and below 0°C across the full spectrum of MAAT shown here ranging from −5°C to 0°C. To demonstrate even greater subsurface thermal effects resulting from subtle n-factor shifts, the following section compares basecase results with 'reduced snow' (n F =0.4) and 'reduced shading' cases (n UF =0.7).

Permafrost thaw and annual groundwater discharge
To evaluate effects of air warming on subsurface thermal and hydrological conditions, the basecase and its variations of reduced snow, reduced shading, 'thick organic layer' (0.2 m thick organic soil), and 'no organic layer' (0.0 m thick organic soil) (table S2) were run to 100 years. In all cases, permafrost depth and groundwater discharge increase with time and MAAT (table S3); yet, these increases (figures 2(b), (c)) are nonlinear despite the linear increase in T air . For the basecase, permafrost-table deepening accelerates from 0.5 cm yr −1 at early times to 6.5 cm yr −1 at late times ( figure 2(b)). These simulated thaw rates are consistent with observations of 1-6 cm yr −1 in eastern Siberia (Streletskiy et al 2015). For the simulations considered, rates are similar (∼0.5 cm yr −1 ) at early times, but reach differing late-time rates ranging from 1.5 to 2 cm yr −1 for reduced snow and thick organic layer cases to 11-12.5 cm yr −1 for reduced shading and no organic layer cases. Whereas rate of permafrost-table deepening ( figure 2(b)) primarily reflects thermal vulnerability, the concomitant response in baseflow (figure 2(c)) additionally reflects sensitivity of the shallow-hydrologic system. Over the century-long model simulation period, the reduced snow case doubles baseflow, whereas the basecase nearly triples and the reduced shading case more than triples baseflow magnitude. The no organic layer case shows the greatest overall change in baseflow, increasing by almost an order of magnitude owing to the greatest loss of permafrost. The thick organic layer case shows the least change but supports the greatest total groundwater discharge among all cases throughout the century-long simulation due to greater ease of flow through the permeable organic layer.

PTZ development and impact on baseflow seasonality
Model results for ice saturation at the seasonal maximum frozen condition show the development of PTZs that progressively expand parallel to the surface over time for reduced snow, basecase, and reduced shading cases (figures 3(a)-(c)), as well as for thick and no organic layer cases ( figure S2(a)-(c)). Non-parallel patterns develop when mineral soil permeability is an order of magnitude greater than the basecase value due to advective heat transport that enhances thaw in areas influenced by warm groundwater recharge (text S1, figure S3). Conditions with <0.9 ice saturation indicate the presence of liquid water (in addition to specified 0.1 residual liquid water saturation). Compared to the basecase, PTZ timing and rate of development is later and slower, respectively, for the reduced snow case reflecting permafrost protection qualities of thin snowpack that allow cold winter temperatures to propagate deeply. In turn, the reduced shading case exhibits earlier and faster PTZ development relative to the basecase reflecting the influence of warmer ground-surface summer temperatures.
Groundwater discharge to the stream boundary varies seasonally, peaking in summer and diminishing during winter when connective discharge pathways are frozen (  these contours, liquid saturation remains above the contoured levels through the annual cycle. PTZs initially develop below the top of permafrost where temperatures remain just below 0°C throughout the year and serve as a precursor for accelerated permafrost thaw. Later in the thaw evolution, as the permafrost top declines and seasonal frost layer thins, PTZs extend above permafrost ( figure 4). The timing and extent of PTZ development varies considerably among cases that invoke only subtle differences in n-factor input (figures 4(a)-(c)). Results emphasize the importance of factors affecting near-surface-energy balance (including snow cover (figure 4) and summer shading (figure 4(c))) in influencing permafrost vulnerability.
Mid-slope spring (maximum frozen) and fall (maximum thawed) ice saturation conditions in the upper 10 m are shown to the right of figure 4 to illustrate the seasonal extremes for each case at t=60 years. The reduced snow case has not developed a PTZ according to spring ice saturation conditions near 0.9. In contrast, spring ice saturations at t=60 years for the other cases indicate PTZs by the presence of liquid water beneath seasonal frost.

Conclusions and implications
Previous work notes increased baseflow as a primary watershed response to permafrost thaw. This study provides physics-based support for such response and expands on existing cryohydrogeologic modeling studies by exploring the dynamics of seasonal groundwater discharge as PTZs develop. Geophysical measurements in boreal watersheds of interior Alaska provide support for PTZ occurrence and are an impetus for this modeling analyses. Model results indicate that relative increases in groundwater discharge associated with climate warming and PTZ evolution are greatest in spring and fall. The timing of increased baseflow and flowpath depths are significant for evaluating DOC, TDN and Hg mobilization as well as for the transport of mineral weathering products. The location of PTZs coincides with permafrost C, N and Hg sources. Furthermore, mineral weathering products may be concentrated in PTZs due to solute exclusion via ice formation. In early spring, flowpaths in still-frozen active layers will be quiescent except in the uppermost soil once seasonal thaw has begun. At this time, groundwater discharge to streams with hydrologically-connected PTZs along contributing hillslopes will therefore be dominated by activation of deeper flowpaths within PTZs that may also serve as conduits for transport of permafrost DOC, TDN, and Hg, and concentrated mineral weather products. In late fall, active-layer flowpaths progressively become quiescent due to seasonal frost development. Consequently, deeper groundwater flowpaths including those through PTZs, become progressively more important and sustain baseflow longer into the fall season. Thus, both spring and fall shoulder seasons provide cryohydrogeologic conditions that may be conducive for streamflow detection of permafrost thaw via aged DOC and permafrost-sourced mercury. Though previous studies provide compelling support to link observed changes in stream chemistry to permafrost thaw-related changes in hydrology, observations of aged DOC from ancient permafrost in aquatic systems remains somewhat elusive due in part to dilution effects, particularly for high-order streams . The strongest signal of aged DOC is thus expected in thawing headwater watersheds, in which a component of baseflow is derived from water that has moved through the zone of thawing permafrost. The greatest likelihood of detecting aged DOC in headwater streams is conventionally assumed to be during the early fall coincident with the seasonal thaw maximum. Based on insight gained from this modeling study's examination of supra-permafrost groundwater flow, aged DOC may also be likely to be observed in spring, post freshet, when flowpaths through PTZs are activated in hydrologically-connected systems well before seasonal frost thaws. Paired seasonal chemistry and age data from streams and their contributing permafrost sources are broadly needed to further interrogate this hypothesis. Geophysical measurements, as presented here, may aid in identifying watersheds with developed PTZs and source-to-sink hydrologic connectivity required to promote deep flowpath activation.
Model results demonstrate how subtle changes in factors affecting near-surface-energy exchange can exert great influence on permafrost thaw rate and PTZ development. Global climate projections suggest an overall reduction in snow extent and depth in northern high latitudes, but predict snow depth increases in Arctic tundra (Bring et al 2016). Current results show that locations with thin snowpack or thick organic layers would experience reduced thaw rate and vulnerability to PTZ development, even under continued warming scenarios; whereas regions with snowpack increases or thin organic layers would experience increased vulnerability to PTZ development. Increased wildfire activity in boreal regions (Kasischke and Turetsky 2006) could lead to reduced summer shading and thin/absent organic layer conditions, thereby increasing vulnerability to PTZs. Arctic shrub expansion (Tape et al 2006) may produce both a warming effect from localized snow entrapment (Sturm et al 2001), increasing vulnerability to PTZs, and a cooling effect from increased summer shading (Briggs et al 2014), decreasing vulnerability to PTZs. Understanding which countering effect may dominate is an ongoing research challenge.
Model results suggest that PTZ presence is a strong indication that an area is on the verge of accelerated permafrost thaw in coming decades, should warming trends continue. These findings are consistent with Lawrence et al (2008) that noted a link between suprapermafrost talik formation and thaw acceleration under four different MAAT warming scenarios, including one with a linear warming trend. Such sensitive thaw responses and corresponding hydrologic responses have important implications for hazard assessments of existing and planned infrastructure in arctic and boreal settings. Detection of PTZs is also important for assessing the global permafrost C feedback, because PTZ development increases the seasonal period (from summer only to perennial) during which newly thawed C can undergo decomposition and atmospheric release as greenhouse gases (Commane et al 2017). Improved methods for PTZ detection over large areas are needed for vulnerability mapping of shallow permafrost.
The model simulation analyses presented here consider complex, coupled thermo-physical processes in idealized permafrost-bounded hillslopes subject to idealized boundary conditions. Heterogeneity and other complexities in structure are purposefully minimized in order to examine fundamental physically-based seasonal flow changes in boreal permafrost hillslopes undergoing thaw. Current results underscore a need for advanced cryohydrogeologic models with coupled solute transport that constrain DOC, TDN, and Hg transit times to streams in the context of complex biogeochemical reactivity at ∼0°C. Continued modeling efforts to represent complexities in spatio-temporal boundary conditions that more explicitly account for the thermal and hydrologic roles of vegetation and snow as well as interannual variability in temperature and precipitation are crucial (Briggs et al 2014, Fisher et al 2016, Jafarov et al 2018. It is anticipated that subsequent site-based research will incorporate enhanced subsurface detail as well as surface-energy and water-balance information to best utilize the growing richness of data in boreal watersheds. Results from such future studies can be viewed in the light of the fundamental framework presented here.