Hydrologic modeling of a perennial firn aquifer in southeast Greenland

Abstract

A conceptual model, based on field observations and assumed physics of a perennial firn aquifer near Helheim Glacier (southeast Greenland), is evaluated via steady-state 2-D simulation of liquid water flow and energy transport with phase change. The simulation approach allows natural representation of flow and energy advection and conduction that occur in vertical meltwater recharge through the unsaturated zone and in lateral flow within the saturated aquifer. Agreement between measured and simulated aquifer geometry, temperature, and recharge and discharge rates confirms that the conceptual field-data-based description of the aquifer is consistent with the primary physical processes of groundwater flow, energy transport and phase change. Factors that are found to control simulated aquifer configuration include surface temperature, meltwater recharge rate, residual total-water saturation and capillary fringe thickness. Simulation analyses indicate that the size of perennial firn aquifers depends primarily on recharge rates from surface snowmelt. Results also imply that the recent aquifer expansion, likely due to a warming climate, may eventually produce lakes on the ice-sheet surface that would affect the surface energy balance.

1. Introduction

Ice-sheet mass loss represents a primary uncertainty in future projections of sea-level rise (Cazenave, 2006; Le Bars, 2017). Meltwater on an ice sheet can either run off or remain stored within the ice sheet as internal accumulation. Perennial firn aquifers temporarily store liquid water, much of which may eventually flow to the ocean (Forster and others, 2014; Koenig and others, 2014; Poinar and others, 2017; Miller and others, 2020). Firn aquifers form when meltwater at the snow surface infiltrates into the firn and fills available pore space without refreezing during winter. This requires that sufficient pore space exists to store meltwater and that high snow accumulation insulates meltwater and prevents refreezing (Forster and others, 2014; Kuipers Munneke and others, 2014). In ice-sheet accumulation zone areas unfavorable to aquifer formation (insufficient melt or accumulation), meltwater refreezes and does not contribute to mass loss to the ocean. Where firn aquifers form, they allow movement of accumulation zone meltwater along an elevation gradient, and permit meltwater discharge to the ocean should a connected hydrologic system develop. In Greenland, perennial firn aquifers cover ~20,000–90,000 km² within the ~1200–2000 m elevation range (Forster and others, 2014; Miège and others, 2016; Steger and others, 2017a). In a firn aquifer, ‘groundwater’ refers to liquid water below the surface of and within the ice sheet or englacial water.

Firn aquifers affect multiple processes related to ice-sheet mass balance, including meltwater transport and possibly ice-sheet dynamics, although the magnitude of their effects is poorly constrained. Firn aquifer total liquid water volume estimates would cause ~0.4 mm of sea-level rise (Koenig and others, 2014) if completely discharged to the ocean; however, sea-level rise resulting from the observed expansion (Miller and others, 2020) and continued flow of meltwater through firn aquifers would be greater. Where meltwater reaches the ice-sheet base, it may reduce basal friction, allowing ice to move downslope more quickly (Poinar and others, 2017, 2019). Thus, firn aquifer dynamics may also be of importance to glacial movement, which is an important contributor to sea-level rise. Hence, firn aquifers may play an important role in ice-sheet and glacier hydrology and processes of glacial mass transport and balance. Predicting the future of Greenland ice sheet (GrIS) mass balance thus requires improved understanding of the role that firn aquifers currently play in this balance, and how this may evolve under a warming climate.

Field work on the firn aquifer upslope from Helheim Glacier in southeastern Greenland (Fig. 1) as well as firn and regional climate modeling have provided basic initial data and analyses of general firn aquifer conditions including aquifer geometry and location (Forster and others, 2014; Koenig and others, 2014; Miège and others, 2016; Montgomery and others, 2017; Steger and others, 2017b), firn stratigraphy, density, porosity and temperature (Koenig and others, 2014; Killingbeck and others, 2020), hydraulic properties (Miller and others, 2017), liquid water flow rates and residence times (Miller and others, 2018, 2020), liquid water volume (Legchenko and others, 2018a, 2018b), impacts of aquifer discharge to crevasses (Poinar and others, 2017), subglacial hydrology (Poinar and others, 2019) and climatic conditions necessary for aquifer formation (Kuipers Munneke and others, 2014; Steger and others, 2017a, 2017b).

Fig. 1. (a) Firn aquifer site map in southeastern Greenland and conceptual model of hydrology adapted from Miller and others (2020). Landsat 8 composite image of Helheim Glacier (21 August 2014) showing profile simulated and conceptual model applied in simulations in 1-D and 2-D. Elevation contours from Cryosat-2 DEM (Helm and others, 2014). (b) Firn aquifer geometry along simulated profile determined from mean 1-D velocity profiles (red line) obtained from seismic (solid black line, inversion of travel times) and radar (water table, upper dashed black line) surveys over a 15 km profile. Aquifer base was dertermined by probability of seismic veloicty increase (lower dashed black line); details are presented in Montgomery and others (2017). Simulated profile extends upslope beyond the seismic profile.

Meltwater production on the ice sheet is projected to increase as global temperatures rise, most likely increasing recharge and perhaps causing firn aquifers to expand. Firn-aquifer expansion has indeed been observed (Montgomery and others, 2017; Miller and others, 2020) and liquid water volume in the Helheim aquifer has increased with time (Legchenko and others, 2018a). Increased meltwater production may also increase the supply of liquid water available to hydrofracture crevasses to the base of the ice sheet, pressurizing liquid water at the ice-sheet base and thereby possibly increasing ice-sheet velocity (Alley and others, 2005; Koenig and others, 2014; McNerney, 2016; Poinar and others, 2017, Poinar and others, 2019). Thus, understanding the fate of meltwater and the dynamics of liquid water flow through the firn aquifer is critical for predicting the ice sheet's response to climate change and, in turn, its impact on sea-level rise.

1.1. Prior firn hydrologic modeling

Typically, studies using numerical analyses to study liquid water flow in firn have used 1-D vertical models of firn thermodynamics and liquid water movement, often developed for modeling firn densification or other processes. However, these approaches (discussed below) are limited because they cannot consider lateral flow of liquid water, which is a key process that transports liquid water mass and energy in firn aquifers.

A range of snow and firn models have been applied to meltwater flow through firn. Regional climate models such as the Modèle Atmosphérique Regional (MAR, Gallée and Schayes, 1994) or the Regional Atmospheric Climate MOdel (RACMO, Van Meijgaard and others, 2008; Noël and others, 2015) can be coupled to subsurface models to resolve meltwater retention and refreezing within snow (e.g. Fettweis, 2007; Ettema and others, 2010; Reijmer and others, 2012; Steger and others, 2017a, 2017b). Many of these firn models are summarized in Stevens and others (2020) and Vandecrux and others (2020a, 2020b). For example, IMAU-FDM is a 1-D vertical model designed to couple with regional climate models to simulate firn density, temperature and liquid water content on polar ice sheets (Ligtenberg and others, 2011). Kuipers Munneke and others (2014) used IMAU-FDM to explain climatic and firn conditions that allow firn aquifers to form, namely high accumulation and high snowmelt rates.

To date, the application of firn models to firn aquifers has yielded useful, but rather limited success (Vandecrux and others, 2020b). Generally, model predictions of vertical firn temperature and liquid water content profiles do not match firn aquifer observations well. Firn model limitations stem from the fact that these models were primarily developed to represent processes such as densification in relatively dry firn, and not firn aquifers. The application of existing models to firn aquifers has several challenges including (1) Models do not fully or accurately represent the physical laws governing liquid water flow through porous media, instead using mathematical simplifications to represent such flow. (2) They are limited to 1-D-vertical-profile representation of what are really 3-D systems. Even where firn aquifers do not exist, lateral differences in meltwater recharge distribution and lateral heterogeneity in firn hydraulic properties would cause 1-D vertical firn models to provide poor representations of hydraulic and freeze/thaw processes in firn. (3) Existing models don't account for differences in liquid water flow physics between saturated and unsaturated zones (e.g. hydraulic conductivity variation with liquid saturation, and the impact of fluid pressure on water saturation). (4) They do not incorporate the physics, and impacts, of lateral liquid water flow on firn-aquifer configuration and dynamics. Some firn models (e.g. the DTU model) do not allow for the formation of firn aquifers at all (Sørensen and others, 2011; Simonsen and others, 2013; Vandecrux and others, 2020b).

To move liquid water through layers in a vertical column of firn, some firn models (e.g. Ligtenberg and others, 2011; Van Pelt and others, 2012; Marchenko and others, 2017) employ variations of ‘bucket-tipping’ schemes in which liquid water in a snow layer is moved to a deeper layer when the liquid water content exceeds a specific level (e.g. 2–6% of pore volume). This movement continues until the liquid water passes out of the model domain, or ‘runs off’, generally at depth (above an impermeable layer or at the bottom of the model domain). However, this ‘run off’ process does not include the physics of liquid water flow through porous media (either saturated or unsaturated), and therefore cannot represent flow through an aquifer or changes in the elevation of an aquifer's water table across any distance or time.

Some models do allow liquid water to move through firn according to more physics-based principles, although they are limited to vertical 1-D analysis. For example, SNOWPACK is a physically based snowpack model that simulates snowpack evolution, and can employ either a bucket-tipping or Richard's unsaturated-zone physics-based equation to simulate liquid water flow (including preferential flow paths) through snow (Wever and others, 2015, 2016). Meyer and Hewitt (2017) developed a continuum model for meltwater flow through snow following Darcy's law that allows for surface runoff if the snow is fully saturated. The subsurface models developed by Langen and others (2017) and Vandecrux and others (2018, 2020a, 2020b) both include Darcy flow. Specifically, the subsurface scheme of the regional climate model HIRHAM5 (Langen and others, 2017) was modified to account for evolution of firn density, albedo and hydraulic conductivity, and allows for liquid water retention in excess of the irreducible liquid saturation. It also incorporates a delayed runoff mechanism that allows liquid water in excess of the irreducible saturation to remain in a layer until it runs off or refreezes (Langen and others, 2015, 2017). For the application of the Community Firn Model to a firn aquifer, both slow matrix flow (Richard's equation) and preferential flow paths can be implemented, and at depth, no runoff is prescribed so liquid water accumulates to form an aquifer (Verjans and others, 2019; Stevens and others, 2020; Vandecrux and others, 2020b). However, vertical 1-D models that do allow the formation of a firn aquifer must include empirical terms that allow for a balance between liquid water accumulation and lateral flow, but such terms are not constrained by physics.

The groundwater flow model SEEP2D has been used to simulate 2-D liquid water flow in firn aquifers (Miège and others, 2016; Montgomery and others, 2020). Miège and others (2016) applied the model to the firn aquifer upslope from Helheim Glacier to simulate flow but used hydraulic conductivity values from mountain glaciers because this parameter had not yet been measured in the aquifer. Subsequent fieldwork showed that the measured hydraulic conductivity values in the aquifer upslope from Helheim Glacier are an order of magnitude higher than the values in mountain glaciers (Miller and others, 2017). Montgomery and others (2020) used SEEP2D to evaluate recharge scenarios to a newly discovered firn aquifer in the Wilkins Ice Shelf. However, SEEP2D does not incorporate equations for representing heat transfer nor does this code simulate liquid water flow in the unsaturated zone. Although these efforts were able to represent general patterns of meltwater flow through the saturated zone, they neglected thermodynamic (such as freezing and thawing) and unsaturated zone processes, which are critical to understanding firn aquifer dynamics and to developing reliable forecasts of future firn aquifer behavior.

1.2. Study goals

This study is intended to (1) provide a test of the overall 2-D cross-sectional conceptual model (that includes a quantitative description of firn-aquifer physics including groundwater flow with energy transport and phase change) of the firn aquifer located upslope from Helheim Glacier (presented in Miller and others, 2020) based on simulation analysis, (2) identify the primary processes that control aquifer configuration, and (3) evaluate how the aquifer configuration may change under continued warming of the Arctic that will likely cause increased recharge to the aquifer. The aquifer conceptual model describes the framework for understanding the system of liquid water movement from the snow surface through the firn. The framework includes the relationship between groundwater system components that affect the movement of liquid water, including recharge rates, firn stratigraphy and hydraulic properties, and water table elevations. Model results that differ substantially from field observations would indicate that the assumed physics and parameterizations or data are incorrect or incomplete. In contrast, model results that resemble field observations likely indicate that the model represents the appropriate physical processes that occur in a firn aquifer.

Given the relatively recent discovery of firn aquifers in Greenland, there were many uncertainties about their key physical processes at the outset of this study that warranted testing via numerical simulation. Does vertical infiltration of meltwater recharge the aquifer or does that liquid get stopped in the unsaturated zone by freezing? Using the conceptual model of the aquifer, can simulated lateral discharge rates match observations? Does the integration of observed vertical temperature profiles, water levels, recharge rates, firn physical properties such as porosity, permeability and stratigraphy, within a multi-dimensional numerical model that includes fluid and energy transport under both unsaturated and saturated conditions with phase change allow for simulation of a firn aquifer or do data and process representation deficiencies exist?

Although some of the modeling approach and results of this study evaluate basic processes and assumptions that are often assumed valid for mineral-based aquifers, it is not known that any of these apply to firn aquifers that are composed only of liquid water and ice. Therefore, given the lack of field and modeling studies of firn aquifers, and given that (in contrast with mineral-based aquifers) there is phase change between the frozen porous medium and the flowing liquid, a fundamental initial analysis is required to evaluate the physics and controlling parameters of firn aquifers. This analysis is the focus of the current study.

This is the first groundwater model of a firn aquifer in two spatial dimensions (2-D) that incorporates groundwater flow, heat transfer and phase change. This study employs a numerical simulation code called SUTRA-ICE, a modified version of SUTRA (Voss and Provost, 2002), that simulates variable-density groundwater flow through unsaturated and/or saturated porous media and transport of solutes or energy. The code accounts for water freezing and thawing, allowing evaluation of whether or not the conceptual model represents self-consistent physics of the fluid mechanics and thermodynamics of groundwater flow and phase change in a firn aquifer.

Transient and steady-state simulations are based on, and compared with, field observations of key features of the firn aquifer upslope of Helheim Glacier, that include snow surface temperature and melt rates, firn temperature as it varies with depth, firn hydraulic conductivity, water-table depths and lateral aquifer extent. ‘Steady state’ refers to simulations where total liquid water recharge equals total liquid water discharge and no changes occur with time. The intention here is not to try matching the details of all field observations, but rather to match the general characteristics of the firn aquifer system as interpreted from the field measurements. This approach allows development of understanding of broad physical controls on aquifer shape and dynamics, despite the fact that there is a lack of temporally and spatially detailed data that would be required to underpin a detailed simulation of complex transient aquifer conditions in a heterogeneous firn aquifer.

This test of the conceptual model through numerical simulation is the first quantitative assessment of a firn aquifer conceptual model and application of SUTRA-ICE code to full snow and ice conditions. This use of SUTRA-ICE to simulate meltwater flow through firn is hoped to be a substantial advance in ice-sheet hydrologic modeling, and this paper presents a new type of simulation with both vertical and horizontal (fully 2-D) liquid water flow and heat transport under both unsaturated and saturated conditions with phase change. The primary controls on the existence and configuration of firn aquifers are identified on the basis of both 1-D and 2-D simulations, and implications of the simulated steady-state 2-D firn aquifer behavior in warmer climates with greater meltwater recharge are elucidated.

2. Methods

2.1. Firn aquifer conceptual model

The conceptual understanding for the formation and function of the Helheim Glacier perennial firn aquifer was initially formulated on the basis of field measurements and observations (2010–2017) and prior firn/crevasse modeling work. A description of this conceptual understanding is presented in Miller and others (2020) and is summarized just below and in Figure 1. The conceptual model developed in this study numerically describes that conceptual understanding, including a hydrogeologic framework, groundwater recharge, flow and discharge to crevasses, hydraulic, physical and thermal properties, and boundary conditions of the aquifer as well as the physical processes of saturated and unsaturated groundwater flow and freeze–thaw defined in SUTRA-ICE. Developing the conceptual model using SUTRA-ICE enables testing of how consistent the conceptual understanding is when constrained by internally consistent representation of physically based hydrologic processes including water mass and total energy balance in multiple spatial dimensions. Inconsistencies between the conceptual model and field observations in any part of the system would indicate field data or process representation deficiencies. The primary question to be answered by testing is: Can the conceptual model predict firn aquifer configuration and dynamics that are consistent with field observations, when providing the model with boundary conditions based on field measurement (primarily firn surface temperature and snow-melt recharge rate)?

Perennial firn aquifers form when sufficient surface meltwater infiltrates to depth, accumulates and saturates open firn pore space (Kuipers Munneke and others, 2014; Miller and others, 2020). The heat contained in the infiltrating meltwater warms the shallow firn to the melting point through sensible and latent heat exchange, which keeps some water in its liquid state, permitting aquifer recharge to occur, and sustaining the liquid saturation within the saturated zone perennially. The upper ~10 m of the firn column, having the lowest liquid water content, experiences wide temperature variations as the firn cools through the fall, winter and spring months, and warms to ~0°C during the summer melt season (Koenig and others, 2014; Miller and others, 2020). Within the aquifer, temperatures remain near the melting point throughout the year as the existence of liquid water is sustained. Below the aquifer, temperatures decrease below 0°C and water is fully frozen (Miller and others, 2020). The liquid water in the firn aquifer overlies fully frozen sub-aquifer ice.

The Helheim firn aquifer is recharged (9–30 cm a⁻¹) during the summer melt season, as surface melt rapidly infiltrates through the firn unsaturated zone (0.4 m h⁻¹) and into the saturated zone (Miller and others, 2020). Mean specific discharge within the aquifer has been measured at 4.3 × 10⁻⁶ m s⁻¹ (std dev. = 2.5 × 10⁻⁶ m s⁻¹), and the aquifer has an average hydraulic conductivity of 2.7 × 10⁻⁴ m s⁻¹ (Miller and others, 2017, 2018). The average lateral hydraulic gradient of the aquifer is 0.01 m m⁻¹ (slope = 0.8°) (Miller and others, 2018). Liquid water flows through the firn aquifer and discharges downslope, likely to crevasses near the edge of the ice sheet (Poinar and others, 2017; Miller and others, 2020). The water table depth ranges between 5 and 50 m from the surface, and upslope from Helheim Glacier, the aquifer (i.e. the liquid-saturated zone above the low-permeability icy base of the aquifer) has an average thickness of 11 m (Forster and others, 2014; Montgomery and others, 2017). Porosity of the firn column decreases from ~60% at the snow surface to between 10 and 20% at the base of the aquifer (Koenig and others, 2014). The model assumes this is disconnected porosity and therefore the aquifer base is a no-flow boundary. These hydraulic and aquifer physical properties were often derived from short-term observations and therefore knowledge of temporal variation is poorly constrained.

2.2. SUTRA-ICE numerical simulations

In this study, the pre-existing field data, observations and understanding, and the physics of freeze/thaw and liquid water flow represented by the SUTRA-ICE code comprise the conceptual model (described above and depicted in Fig. 1) that is evaluated in 1-D and 2-D (Fig. 2). The goals of the simulation analyses are to (1) use measured field data to calibrate and constrain the model in order to elucidate major controls on the aquifer and to confirm the viability of the combined descriptive and physical process-based conceptual model, and (2) to use the calibrated model to understand the general effects on the aquifer of increasing recharge that is expected to accompany increasing melt as the Arctic warms. The design of the current study emphasizes parsimonious modeling that is focused on capturing overall behavior in the face of limited field data, rather than attempting to apply complex and highly parameterized approaches, following the method described by Voss (2011). This approach supports understanding of the primary controls on the firn aquifer and effects of varying recharge.

Fig. 2. Two-dimensional cross-sectional model setup and boundary conditions. Top is located at upper surface of snow (‘ground’ surface), bottom is located at approximate depth where firn temperature is −1°C and liquid water has frozen into solid ice except irreducible liquid saturation with no connected porosity. Downhill boundary is located at a crevasse in which the water table is set at 20 m depth, and uphill boundary is located near upper extent of the observed firn aquifer. Snow surface gradient is 0.01 m m⁻¹. Vertical exaggeration is 50 times.

The SUTRA-ICE code, created by the U.S. Geological Survey, was one of the first of a group of cryohydrogeologic codes, developed over the past two decades (Rühaak and others, 2015; Grenier and others, 2018), that couple groundwater flow equations to heat transfer equations with dynamic freeze–thaw capabilities. SUTRA-ICE represents the physics of the freeze–thaw process in groundwater where the liquid water can fully or partially saturate the porous medium. Details and governing equations are fully described in previous studies (McKenzie and others, 2007a; McKenzie and Voss, 2013; Wellman and others, 2013; Kurylyk and others, 2016; Evans and Ge, 2017; Evans and others, 2018; Evans and others, 2020) and are not repeated here. Readers may refer to these publications for complete information about the code and its functioning.

SUTRA-ICE was originally developed for studies of groundwater flow in permafrost regions and in regions with seasonal subsurface ice. The code has been applied to a variety of cases considering seasonal ground ice and permafrost, how the formation and loss of subsurface ice interacts with the flowing liquid groundwater, and how both are impacted by hydrologic and thermodynamic processes at the ground-surface (McKenzie and others, 2007b; Ge and others, 2011; McKenzie and Voss, 2013; Wellman and others, 2013; Briggs and others, 2014; Kurylyk and others, 2014, 2016; Evans and others, 2018, 2020; Lamontagne-Hallé and others, 2018; Zipper and others, 2018; Walvoord and others, 2019; Rey and others, 2020; Rushlow and others, 2020). The current study documents the first application of SUTRA-ICE to firn, in which only water in liquid and ice phases (no mineral grains) is present.

SUTRA-ICE has several advantages in simulating firn aquifers over other codes that have been applied to firn aquifers. SUTRA-ICE solves groundwater mass-balance equations for the unsaturated and saturated zones and allows for liquid water flow in one, two or three space dimensions (1-D, 2-D or 3-D). Representation of both vertical and horizontal flow in firn aquifer simulation is a significant advance beyond 1-D bucket-tipping, or even physics-based firn-water 1-D flow models because multiple dimensions allow for a dynamic, laterally extensive water table, representation of hydrologic and thermodynamic processes in both the unsaturated and saturated zones, and representation of liquid water flow through, and out of, the aquifer. Simulation of lateral liquid water flow is particularly important for representing a firn aquifer's saturated zone and aquifer discharge. This multidimensional capability can also be applied to simulate vertically and laterally heterogeneous firn, preferential flow paths and impermeable ice.

The SUTRA-ICE freeze–thaw capability accounts for the latent heat of fusion and varies thermal properties with changing total-water, liquid and ice saturations. Porosity is the fractional volume of void space in the porous medium (that can host liquid water and ice). ‘Saturation’ refers to the volume of total water (ice and/or liquid water), as a fraction of the volume of void space in the porous medium. ‘Total-water saturation’ is the sum of liquid water saturation and ice saturation. The code allows effective permeability to change as a result of phase change, with effective permeability decreasing as liquid water saturation decreases (due to drainage or freezing). The code does not simulate other possible processes that may occur in firn including cryosuction, air flow, vaporization and sublimation in unsaturated zones, and these represent sources of uncertainty in the model. Cryosuction has not been well studied in firn and is not expected to occur during the recharge period of the 1-D simulations or at all under the 2-D steady-state simulations, where no freezing occurs in the unsaturated zone. While cryosuction may be important during the transition to an isothermal unsaturated zone at 0°C, recharge during this transition is minimal and so cryosuction effects are expected to be minor. However, should cryosuction prove to be a significant process in firn aquifers, it could lead to some redistribution of liquid and ice contents in the unsaturated zone of a firn aquifer. Future experiments in firn might provide more certainty regarding the role of these other processes in firn aquifers. (Note that the version of SUTRA-ICE used was specifically modified for this project to output specific discharge, rather than liquid water velocity.)

In the case of a firn aquifer there is no rock-mineral phase and only liquid water and solid water (ice) exist, so special choices of SUTRA-ICE parameter values must be made. The mineral phase of the model domain is assigned as a permanent ‘ice backbone’ fraction that has the properties of ice. Seasonal snow cover porosity ranges between 0.4 for dense snow and 0.98 for freshly fallen, dry snow (Armstrong and Brun, 2008), and in the current study, a range of backbone values from 0.3 to 0.6 (as a fraction of the total porous medium volume) is considered. The permanent ‘backbone’ of ice can never melt or change its value in a simulation. Saturation values refer to only the non-backbone void space – and ice saturation (which can vary during a simulation) does not include the backbone fraction. Lower backbone values can be simulated for cases in which more firn can melt than considered in this study (i.e. a lower backbone value can be specified to allow more firn melt). This approach can simulate very low backbone fractions (e.g. 0.001) but cannot simulate complete melting of the ice (a zero-valued backbone fraction), because if this were to occur there would be no porous medium (which requires the presence of both liquid and ice grains) to simulate. SUTRA-ICE simulations that approach this all-liquid state could be done by setting the backbone fraction to lower values; however, such results are not reported in the current study.

Within the pore space (e.g. for a backbone fraction of 0.4, the total non-backbone initial porosity is 0.6, which is among the lower values observed in firn), the saturations of liquid and ice vary spatially depending on simulated steady-state fluid pressure and temperature, according to freeze/thaw unsaturated-zone physics and prescribed state functions. The initial pore space can fill with liquid water and/or ice as the simulation progresses through time, depending on temperature, fluid pressure and available water.

The assigned properties of liquid water and ice are shown in Table 1. For discussions related to SUTRA-ICE simulations, ‘total water’ refers to both liquid and ice, ‘liquid’ refers to the liquid-water phase and ‘ice’ refers to the solid-water phase. These all refer to water/liquid/ice mass in the non-backbone volume of the firn. The backbone volume of ice remains constant throughout a simulation.

Table 1. Table showing general model specifications, boundary conditions and physical parameters used in 1-D and 2-D simulations

Three user-defined state functions are required for the physics description contained in the firn-aquifer conceptual model to describe how: (1) liquid water saturation (for fully saturated conditions) varies as a function of only temperature (providing the saturation of liquid and ice under fully saturated conditions); (2) total-water saturation varies as a function of only fluid pressure (defining unsaturated conditions, the saturations of total water and air in the firn pores); and (3) relative permeability varies as a function of liquid saturation (defining the fraction of total permeability of firn that remains as the liquid saturation decreases due to either drainage or freezing). Currently, no direct data exist to accurately define the saturation functions for Helheim firn aquifer, so estimates are based on observations of the system as well as results of 1-D simulations (Table 2). A variety of function types have been applied in coupled water and heat transport models with freeze–thaw, including simplified linear functions when detailed laboratory-derived functions are not available (Kurylyk and Watanabe, 2013). The functions, and their shapes, remain a matter of ongoing research and debate.

Table 2. Table showing model inputs and parameters for 1-D sensitivity analysis and simulations

In addition to the firn aquifer base cases, parameters were uniquely varied for each base case with high, medium, low and very low values. Bold values indicate base case variations.

The low salinity of aquifer liquid water (<50 μS cm⁻¹) indicates that salinity that would allow liquid water to exist at temperatures much below 0°C is not present (Miller and others, 2018). Therefore, the liquid water saturation function and values were selected because they enable freezing within a narrow range of temperatures; liquid water freezes according to a steep exponential freezing function (see Evans and Ge, 2017, Eqn (1), where freezing temperature is 0°C).

Liquid water saturation for fully saturated firn conditions (S _Lsat) varies with temperature according to the model adapted from Kozlowski (2007) as:

(1)$$S_{{\rm Lsat}}( T ) = ( {1-\;S_{\rm r}} ) {\rm e}^{-( {T-T_{\rm f}/w} ) } + S_{\rm r}\;{\rm for}\;T < T_{\rm f}$$

where S _r is the residual liquid saturation, T is temperature, T _f is freezing temperature, and w is the exponential model parameter.

Total-water saturation (S _w, the sum of liquid and ice phase saturations) is calculated with a piecewise-linear function that varies with pressure:

(2)$$S_{\rm W}( p ) = \left\{{\matrix{ {1\;{\rm for}\;p > p_{\rm e}} \hfill \cr {( {1-S_{\rm t}} ) \left({\displaystyle{{\,p-p_{\rm w}} \over {\,p_{\rm e}-p_{\rm w}}}} \right) + S_t\;{\rm for}\;p_{\rm w} \le p \le p_{\rm e}} \hfill \cr {S_{\rm t}\;{\rm for}\;p < p_{\rm w}} \hfill \cr } } \right.$$

where S _t is the residual total-water saturation, p is pore pressure, p _w is the pressure at which the linear segment of the function gives S _w equal to S _t, and all pressures are negative. This equation follows Evans and others (2020), Eqn (2), where the air entry pressure (p _e) is 0. When p > 0, the firn is fully saturated (S _w = 1).

To determine the saturation of liquid, S _L, in the firn unsaturated zone, it must be noted that the saturation of liquid water, S _L, cannot exceed the total-water saturation, S _W, even if Eqn (1) indicates a higher value. Thus,

(3)$$S_{\rm L} = \min \{ {S_{{\rm Lsat}}^{} ( T ) , \;\;S_{\rm W}( p ) } \} $$

The ice saturation, S _I, is then given by

(4)$$S_{\rm I} = S_{\rm W}-S_{\rm L} = \max \{ {S_{\rm W}( p ) -S_{{\rm Lsat}}^{} ( T ) , \;\;0} \} $$

Relative permeability (k _r) is calculated as a function of liquid-water saturation as:

(5)$$k_{\rm r} = ( {1-k_{\rm s}} ) S_{\rm s}^{} + k_{\rm s}$$

where k _s is the relative permeability at residual liquid-water saturation and S _s is a scaled liquid-water saturation defined as:

(6)$$S_{\rm s}^{} \equiv \displaystyle{{S_{\rm L}-S_{\rm r}} \over {1-S_{\rm r}}}$$

Relative permeability decreases linearly from a value of 1 when liquid saturation has a value of 1, to a minimum relative permeability value when the liquid saturation drops to the selected liquid saturation value at which this minimum occurs. For lower liquid saturations, the relative permeability remains at its minimum value (values are given in Table 1). Decreases in permeability due to freezing are calculated using a linear function to simulate observed permeability, which is relatively homogeneous (Miller and others, 2017).

There is a current absence of field and laboratory data to constrain the functions for firn aquifers (Eqns (1–5)). However, this modeling study is not intended to represent the details of unsaturated zone water saturations and resulting relative permeabilities, as would a study that is focused on the details of an unsaturated zone in a typical hydrogeologic system. Rather, the current selections of simple linear relations for total-water saturation and relative permeability allow the current model to represent the primary saturation and drainage processes in a firn unsaturated zone. According to the functions, liquid water saturation decreases when temperature decreases, total-water saturation decreases when pressure drops and relative permeability decreases when liquid saturation drops. These simplifications represent the primary drainage and permeability changes and provide a basis for addressing the study goal of evaluating the major water and energy-balance processes that control the overall regional hydrologic behavior of a firn aquifer. However, such simplifications also represent a source of uncertainty in the modeling approach that future work may address.

The ice saturation, Eqn (4), in unsaturated conditions is determined by both the liquid water saturation, Eqn (3), and total saturation, Eqn (2), state functions (and is thus a function of both temperature and fluid pressure). The ice saturation state function appears as follows. For lower total water contents in the state function due to lower fluid pressure at constant temperature (for which the porous medium becomes more unsaturated), the liquid water content remains constant and ice content becomes lower. Kurylyk and Watanabe (2013) explain that in laboratory measurements (Jame, 1977; Watanabe and Wake, 2009; Zhou and others, 2014) the total unsaturated water content (which is a function of pressure) does not affect the liquid water content over a range of freezing temperatures. For example, this physical behavior can be observed in Figure 2 of Watanabe and Wake (2009), where, for a wide range of frozen soil types, the measured liquid water content depends only on temperature, and does not depend on the total-water saturation. Then, for lower total-water saturations at which there is no longer ice in the pores, the liquid saturation becomes equal to the total saturation and decreases with decreasing fluid pressure, as occurs in normal unfrozen mineral aquifer fabrics. This implies that the liquid-water saturation in partially frozen soils is independent of the total-water saturation, provided that the total-water saturation is greater than the maximum possible liquid saturation (the value that occurs for fully saturated conditions as a function of temperature). The current study assumes that similar state functions control unsaturated freeze–thaw behavior in firn, which contains no mineral grains, and this may represent a source of uncertainty in firn-aquifer modeling. Future studies could test this assumption in laboratory experiments on firn.

2.3. 1-D simulations

Calibration of model parameters employed in the 2-D cross-sectional steady-state model, described below, was achieved through analysis of 1-D time-variant simulations of firn temperature profiles. For the 1-D calibration, a variety of firn-aquifer situations were considered (Table 2). Three base cases were defined (drier firn, icier firn and wetter firn) and then input values were applied for high, medium, low and very low parameter-value variation cases (Table 2), resulting in a series of simulated vertical profiles of spring firn temperatures. The parameters and input values for total-water saturation, relative permeability and liquid water saturation equations from the combination that resulted in the closest match to measured firn temperature profiles were then employed in the 2-D steady-state simulations, described in the next section.

The 1-D vertical profile model analysis focuses on unsaturated zone freeze/thaw processes. Simulated temperature profiles were compared to measured temperature profiles for April 2016, based on a temperature string at FA15_1 (latitude: 66.3622, longitude: −39.3119, elevation: 1664 m) installed in April 2015 (Miller and others, 2020). By August 2016, the surface temperature sensor had become buried ~1.5 m below the surface by new snow accumulation. A burial rate of 0.1 m month⁻¹ was applied to the temperature sensor depths from 2015 to determine the sensor depth in April 2016. In reality, each sensor was deepened by ~1.2 m over this period. However, SUTRA-ICE observation nodes were set to occur at 1 m intervals, and so the measured depths are increased by only 1 m to coincide with SUTRA-ICE temperature observation nodes.

Transient April temperature profiles were simulated for a 1-D firn column consisting of a simple vertical stack of 2-D rectangular finite elements that provide the same 1-D vertical numerical results along the left and right edges of the stack of elements. The firn column simulation begins during the prior summer, when the firn column is isothermal at 0°C, and the water table is located at 20 m depth. Recharge of 0°C liquid water at the top of the column initially occurs in the simulation for 1.5 months, then stops. The liquid water is allowed to drain through the unsaturated zone for the first month, after which the snow surface temperature is lowered to a constant specified temperature value for 7 months, based on weather station data. The average air temperature for 2015 was −13°C and the average winter temperature for 2015 (January 1–April 30, and October 1–December 31) was −19°C (Reijmer and others, 2019). At the column base for all three periods, the temperature is held at 0°C (a specified-temperature boundary condition) and the fluid pressure is specified at a constant value of 294 300 Pa (equivalent to 30 m of hydraulic head, placing the simulated water table at 20 m depth, as observed in the field).

To characterize the effect of liquid water flow on the unsaturated zone temperatures under different parameterizations, dry, wet and icy firn base cases were established (Table 2). These base cases represent variations in parameters that control total and liquid water saturations. In the dry base case, the residual total-water saturation (ice + liquid) is low (0.01). In the icy base case, the residual total-water saturation is higher (0.20). In the wet base case, the minimum liquid saturation is high (0.20; this keeps more liquid water in the firn as the temperature decreases). For each of these base cases, the surface temperature during the winter, amount of summer recharge, fraction of backbone ice, unsaturated function, relative permeability function, and freezing function were varied over a range of values (see Table 2), providing a range of scenarios for each base case. Final simulated temperature profiles for each scenario were then compared to the measured firn temperatures.

2.4. 2-D simulations

2-D simulations were also developed and the aquifer response to changes in recharge was evaluated. The 2-D simulations represent the firn aquifer along a flowline of Helheim Glacier, where a series of sites provided surface-based measurements (Fig. 1) (Montgomery and others, 2017; Legchenko and others, 2018a). Each simulation was run in transient mode from arbitrary initial fluid pressure and temperature conditions until it reached a steady-state solution for liquid water pressure and temperature (same temperature value for liquid, ice and backbone), the two primary variables that are being solved for. The initial conditions and initial ice density vertical profile has no impact on the steady-state solution. The system of governing equations is inherently non-linear and running a transient simulation from arbitrary initial conditions to steady state resolves the non-linearities and provides the self-consistent numerical solution. The simulated 2-D cross-sectional domain consists of a 16 km long, 50 m thick vertical 2-D cross-section of the firn aquifer system (Fig. 2). The model configuration, the temperature and fluid pressure boundary conditions, as well as the freezing, unsaturated and relative permeability functions are explained in the following.

Regarding the spatial configuration and position of the model domain, there are two ways to accommodate glacial motion in the 2-D cross-sectional firn aquifer model. (1) It may be considered that the model domain moves downslope at the velocity of glacial movement with no new ice entering or leaving the cross-section; thus, the energy and water mass contribution of the moving glacier in which the firn aquifer exists would be zero. (2) If the model domain is set in a fixed location in space, the glacial ice moves through this fixed location, but in fact, the contribution of ice entering the upstream boundary and exiting the downstream boundary to the water mass and energy balances is negligible relative to that of the meltwater flux and vertical energy transport resulting from top-surface and bottom-surface temperature values. Therefore, the ice flow could be ignored in model simulations, and a fixed location of the model domain (2) was used for simplicity.

Boundary conditions and other details of firn-aquifer model setup are as follows. Meltwater generated at the surface of the snow was applied along the top boundary of the model as a specified fluid source boundary condition. Although non-uniform recharge occurs in the field, for simplicity in this initial firn-aquifer modeling study, recharge was applied uniformly across the top of the model, and the recharge rate is held constant during each simulation. Following confirmation of the conceptual model of a firn aquifer in the first 2-D simulations, the recharge rates, based on melt and recharge rates estimated from weather station data (Miller and others, 2020), were varied in another set of runs (5 and 30 cm a⁻¹) to determine their impact on the configuration of the firn aquifer.

The downslope model domain boundary was placed at the upslope edge of an open crevasse in the Helheim Glacier, where liquid water can flow out of (or in to) the aquifer. The hydrologic condition assumed at the crevasse was based on the measured water-table location at 20 m depth near the crevasse (Miège and others, 2016). The water table along the upslope side of the crevasse was simulated by setting a specified fluid pressure of 0 Pa at the water table location (20 m depth) and by specifying hydrostatic fluid pressure values everywhere along the crevasse boundary (with hydrostatic pressures becoming more negative along the crevasse boundary above the water table and becoming greater positive values with greater distance below the water table). The linear hydrostatic pressure distribution along the crevasse wall was calculated assuming a constant density of liquid water at 0°C, of 1000 kg m⁻³. (Note that liquid water density was held constant for all simulations, in order to preserve simplicity by avoiding natural convection of liquid water due to density differences.) Liquid water could theoretically flow in the upslope direction into the aquifer at the simulated crevasse boundary during any period in which the water level in the crevasse was higher than the water level in the aquifer, although this scenario was not evaluated and the simulations reported here have only aquifer outflow at this boundary.

The temperature specified at the base of the cross-section (top of assumed impermeable ice) was held constant at −1°C, based on vertical temperature profile measurements (Miller and others, 2020). The temperature specified at the snow surface was held at 0°C to match the temperature (also 0°C) of the liquid-water recharge added at the snow surface. Although measured surface temperatures decrease below 0°C when there is no meltwater recharge, the surface temperature remains at about 0°C during recharge periods. This choice was deemed appropriate for steady-state simulations.

Initial permeability values were set at 2.7 × 10⁻¹¹ m² based on field measurements of hydraulic conductivity (Miller and others, 2017). The longitudinal dispersivities were set to the lowest value (usually about one-fourth of the element size measured in the direction of flow) that avoids instability in the numerical energy-transport solution.

The user-defined function values for total-water saturation, relative permeability and liquid water saturation in the 2-D simulations were determined based on the 1-D simulation results (Table 3). The dry base case function values were selected for 2-D simulations because they produce firn temperature profiles similar to observed temperature profiles and they are more similar to field observations than the wet base case. In the assumed functionality, irrespective of how low the fluid pressure becomes, a 1% total-water saturation (liquid plus ice) was retained in the firn pore space (backbone fraction was set at 0.4) as the total residual water saturation. Liquid water freezes according to a steep exponential freezing function ranging from 100% liquid at 0°C to the residual liquid saturation at −0.2°C. For possible future simulation analyses, note that a steep linear function would likely give similar simulation results as were obtained using this function. As liquid freezes, the effective permeability, given by the product of total permeability for fully saturated firn and the relative permeability, decreases linearly from the total permeability value when fully saturated to a value 100 times lower when liquid saturation is at its minimum-possible value. The decrease by a factor of 100 is based on experience in mineral soils, and although this value is arbitrarily selected for this analysis, it provides a simple but significant decrease in the ability of liquid water to flow through firn when the liquid saturation deceases to very low values.

Table 3. Table showing backbone and user-defined function values used in 2-D simulations

3. Analyses and results

3.1. 1-D simulations

By producing a variety of firn vertical temperature profiles, the 1-D simulations provide insight into the controlling parameters of the system and the major processes occurring (Fig. 3). The simulated temperatures for the dry and wet base cases most closely match the measured April temperature profile (Fig. 3a). Inspection of differences and similarities among these results, as discussed below, indicates that the most important input parameters that control simulated firn temperature profiles are surface temperature, residual total-water saturation (minimum saturation of liquid water plus ice) and thickness of the capillary fringe (only if the capillary fringe reaches the snow surface so there is no zone at irreducible total-water saturation). The capillary fringe is the zone between the unsaturated and saturated zones where pores are saturated, but the pressure head is less than atmospheric. Further inspection indicates that less important parameters controlling firn temperature are the backbone fraction and freezing function parameters, but only for the icy cases. The simulated firn temperature profile in most cases is not sensitive to changes in recharge from the prior summer, backbone fraction or relative permeability function (Fig. 3b).

Fig. 3. Simulated and measured firn temperature. (a) Simulated firn temperature from the icy, wet and dry base cases and icy case with winter surface temperature of −5°C. Dry and wet base case results overlay one another. (b) Simulated dry base case compared to measured firn temperature profile for a range of parameters that simulated firn temperature is insensitive to. All results overlay one another. (c) Simulated firn temperature from the dry base case with a range of winter surface temperature compared to measured firn temperature, (d) and simulated firn temperature from the dry base case with a range of residual total-water saturation (unsaturated function) values compared to measured firn temperature. Simulated firn temperature is most sensitive to winter surface temperature and the residual total-water saturation in the total-water saturation function.

Simulated firn temperature is sensitive to the winter surface temperature (Fig. 3c). The simulated firn temperature becomes significantly warmer or colder than the measured firn temperature, depending on the snow surface temperature in the winter (Fig. 3c). The coolest winter surface temperature tested, −20°C, results in a substantially cooler firn than observed, although the differences are greatest at the snow surface where measurements are lacking, making comparison difficult (Fig. 3c). However, assuming the measured temperature line slope is constant to the snow surface, surface temperature differences are estimated to be >5°C. Simulated firn temperatures with a surface temperature of −13°C align closely with measured firn temperatures (Fig 3c). Large temperature differences between measured and simulated temperatures also occur with the icy base case, and a winter surface temperature of −5°C (Fig. 3a) in which simulated firn temperature is up to 4.9°C warmer than that measured, indicating that total water content (liquid and ice) is important as well.

Simulated firn temperature is also sensitive to parameters controlling total water content in the unsaturated firn including the residual total water content and the thickness of the capillary fringe (Fig. 3d). Higher residual total-water saturation values result in warmer simulated firn temperature compared with measured temperature. This is shown in Figure 3a, where simulated temperatures under the icy base case are warmer than dry base case simulated temperatures and measured temperatures. This is also shown in Figure 3d, which shows how variation in residual total-water saturation influences simulated temperatures. When the capillary fringe is set to include the entire unsaturated zone (i.e. the ‘high’ variation of each base case), the increased liquid water content results in warmer simulated firn temperature of up to 2°C for the dry case, 3°C for the icy case and 1°C for the wet case.

3.2. 2-D simulations

Two primary results are obtained from the 2-D simulations. First, agreement between modeled and observed aquifers demonstrates that the conceptual model captures the primary processes and parameters that control the firn aquifer. Second, as recharge rates increase, the aquifer expands.

3.2.1. Description and comparison of modeling results to field observations

Using the total-water saturation, relative permeability and liquid water saturation function parameters determined from the 1-D analysis, the 2-D cross-sectional model of the firn aquifer was run in a transient simulation from arbitrary initial conditions until it reached approximate steady-state conditions. A stable firn aquifer was formed in the simulation (see Figs 4–6 for a recharge rate of 15 cm a⁻¹). The system has an unsaturated zone of porous firn overlying a saturated zone where liquid fills pores above an ice zone where pores are filled with ice. The aquifer saturated zone (measured from its base, where ice saturation is at its maximum value and liquid water saturation is at its lowest possible value) to the water table (which occurs at some depth below the top of the cross-section) thins upslope. Temperatures within the simulated 2-D aquifer are at or just below 0°C, and are below 0°C below the aquifer base. The recharging snowmelt liquid water flows nearly vertically downwards (while undergoing freezing/thawing) through the firn's unsaturated zone overlying the aquifer, until the liquid reaches the saturated zone. At the firn-aquifer water table, the simulated meltwater naturally becomes part of the saturated zone of the aquifer, and the liquid water flows laterally (approximately parallel to the glacial surface slope) within the saturated zone, also with concomitant freezing or thawing. Simulated specific discharge rates (volume of liquid water in m³ flowing through a 1 m² cross-sectional area of firn, per second of time) range from 1.5 × 10⁻⁹ to 7.6 × 10⁻⁶ m s⁻¹. The liquid water discharges at the crevasse location along the downslope edge of the model domain.

Fig. 4. Simulated fluid pressure and firn temperature for three recharge rates (5, 15 and 30 cm a⁻¹). The water table occurs where the fluid pressure is zero (shown in blue-white boundary on the pressure plot). Snow surface gradient is 0.01 m m⁻¹. Vertical exaggeration is 50 times.

Fig. 5. Simulated liquid saturation, ice saturation and total-water saturation (liquid plus ice) for three recharge rates (5, 15 and 30 cm a⁻¹). Liquid and ice saturations are the volume of liquid and ice per volume of pore space, respectively. Snow surface gradient is 0.01 m m⁻¹. Vertical exaggeration is 50 times.

Fig. 6. Simulated firn aquifer showing liquid saturation, specific discharge and logarithm of specific discharge, to visualize near-vertical flow through the unsaturated zone and lateral flow through the saturated zone for three recharge rates (5, 15 and 30 cm a⁻¹). The region of transition from maximum to minimum liquid saturation around the aquifer is very narrow. Specific discharge is a product of fluid velocity, firn total porosity, and, liquid-water saturation, sometimes referred to as ‘Darcy velocity’. Snow surface gradient is 0.01 m m⁻¹. Vertical exaggeration is 50 times.

Agreement between measured and simulated aquifer geometry, temperatures and discharge confirms that the conceptual understanding of the Helheim Glacier firn aquifer is reasonable (see 4.2 for discussion of this comparison). Note that observations represent a system not at steady state and therefore may differ from results of the system simulated at steady state. Therefore, a general comparison of spatial patterns and rates is provided.

The 2-D steady-state simulation (recharge rate = 15 cm a^–1) produces a saturated zone that has generally the same geometry as the observed Helheim Glacier firn aquifer, where recharge has been estimated to range between 10 and 30 cm a⁻¹. In the observed system, an unsaturated zone overlies a saturated zone, and the base of the aquifer is formed by ice. The aquifer thins upslope. In the simulated system, an unsaturated zone overlies a saturated zone, and the base of the aquifer is formed by ice. The aquifer thins upslope. Observed water table depths range from 5 to 50 m, and simulated water table intersects the snow surface. Observed water table undulations were not simulated, potentially because simulations do not include spatial variation in recharge along the snow surface slope, snow surface slope variations or 2-D heterogeneity in firn physical properties. However, the overall patterns of liquid, ice and total saturation agree in that the unsaturated zone overlies the saturated zone, and ice forms the aquifer base. For more detail, see a comparison of Figures 4 and 5 (15 cm a⁻¹ recharge rate) with Figure 1b and also with Figure 10 in Miège and others (2016) and with Figure 4 in Montgomery and others (2017).

Vertical temperature profiles within the 2-D simulations agree with observations. Temperatures within the simulated 2-D aquifer are perennially at or just below 0°C, and cool below 0°C below the aquifer base. Observations of vertical firn temperature profiles indicate that within the aquifer, temperatures remain at or near 0°C throughout the year and sustain liquid water, and below the aquifer temperatures decrease below 0°C and water is fully frozen (Miller and others, 2020).

Simulated fluid flow rates and patterns agree with observations. Simulated specific discharge rate ranges (1.5 × 10⁻⁹ to 7.6 × 10⁻⁶ m s⁻¹) are similar to measured values, and show similar vertical profiles as are reported in Miller and others (2018). Measured specific discharge was highest at the top of the aquifer (~5 × 10⁻⁶ m s⁻¹) and generally decreases with depth toward the bottom (~1 × 10⁻⁷ m s⁻¹), with an average specific discharge of 4.3 × 10⁻⁶ m s⁻¹ (σ = 2.5 × 10⁻⁶ m s⁻¹) (Miller and others, 2018). The agreement between modeled and observed aquifer characteristics of geometry, temperature and discharge using a relatively simple representation of the physics and physical characteristics of the aquifer instead of mere fitting of many parameter values or an overparameterized model representation, suggests that this model of the firn aquifer includes and adequately represents the primary processes controlling the aquifer.

3.2.2. Aquifer response to changes in recharge rate

A series of simulations was run with a range of mean annual recharge rates to understand the effects of changing recharge on the general configuration of the steady-state saturated zone in firn aquifers. Results (Figs 4–6) show that the geometry of the firn aquifer's saturated zone varies with recharge rate. This is similar to what would be expected to occur in mineral aquifers, but for firn aquifers, where liquid and ice saturations depend not only on total porosity and recharge rates, but also on freezing and thawing and, thus, temperature, this relationship was untested. Should increased recharge cause an increased total porosity of the firn via melting, then the amount of water-table rise or aquifer expansion upslope from an increased recharge rate might not be as great (maybe none) as would be expected in a mineral aquifer, should some (or all) of the extra recharge fit into the newly created pore space. This is one of the unique behaviors that may occur only in firn aquifers and the relation of porosity dependence on flow rate and temperature is a main functionality of the numerical model.

Simulations show that the steady-state saturated zone reaches higher upslope and the water table is higher at increased recharge rates, whereas the saturated zone occurs lower downslope with a deeper water table at reduced recharge rates (compare 5–30 cm a⁻¹ in Figs 4–6). For cases with high-enough recharge (≥15 cm a⁻¹), the steady-state simulated water table intersects the snow surface. For high-enough recharge rates, the water table intersects the snow surface along the entire slope (see case of 30 cm a⁻¹).

4. Discussion

4.1. Insight into the controls on long-term meltwater retention and flow within firn

Scenarios representing variations in base model conditions and aquifer recharge rates were evaluated to learn the factors that control long-term meltwater retention and flow within firn. The cryohydrogeologic factors that control firn temperatures and aquifer configuration (water-table location and saturated-zone volume) are found to be as follows: (1) surface temperature, (2) meltwater recharge rate, (3) residual total-water saturation (minimum saturation of liquid water plus ice) and (4) thickness of the capillary fringe (only if it extends to the snow surface).

1. (1) The results of the 1-D transient simulations highlight the importance of surface temperature controls on firn temperature, and the sensitivity that firn may have to rising temperatures under a warmer climate. Surface temperature determines melt rates and influences firn column temperatures throughout the year. Warmer surface temperatures allow for more surface melting and reduce wintertime firn column cooling.

   The sensitivity of firn temperatures to surface temperatures in simulations can be explained by the established physics included in model code. Phase change is an important process that controls firn temperature profiles. Without phase change, cold winter surface temperatures would propagate downwards into the firn via heat conduction, dropping temperatures at depth to many degrees below 0°C. However, phase change from liquid water to ice requires thermal energy (latent heat) and occurs in a very small temperature range just below 0°C (here considered to −0.2°C), which keeps the local firn temperature in this range as long as winter freezing (and summer thawing) is occurring. Warmer air temperatures can therefore warm firn temperatures by both increasing liquid content (and associated thermal energy) from increased snowmelt rates, as well as by reducing firn cooling during winter months due to higher surface temperatures.

   Subsequent increased liquid water content, in the capillary fringe and/or in the saturated zone, retains more thermal energy in the firn column throughout the year via both sensible heat and latent heat, and buffers the winter cooling, resulting in warmer firn temperatures. Thus, the elevated liquid-water content due to increased snowmelt rates also impedes winter firn cooling, resulting in warmer simulated unsaturated firn temperatures. Dry firn has much larger yearly temperature variation than does firn that contains more liquid water. The firn in the saturated zone remains at or just below 0°C throughout the year due, in part, to the high liquid water content.

2. (2) Recharge rates correlate with surface temperatures. Aquifer response to potential climate change can be forecast because changes in snowmelt recharge rates corresponding to changes in surface temperature and solar radiation cause aquifer volume changes. Snowmelt rates are expected to rise under a warmer climate, potentially leading to increasing overall recharge rates, which would result in thicker firn aquifers. With enough additional recharge, the water table is projected to rise to the snow surface, forming lakes or rivers.

3. (3) Increased residual total-water saturation results in warmer simulated firn temperatures at the end of the winter cold season. The effect of increasing the residual total-water saturation is to add total water mass that stores more thermal energy (both sensible and latent). This causes slower seasonal temperature changes in the unsaturated firn column.

4. (4) Similarly, when the unsaturated zone properties are such that the capillary fringe reaches the snow surface (higher fluid pressure at residual liquid saturation), the increased liquid content throughout the firn column retains more thermal energy, resulting in greater thermal inertia. This causes firn temperatures at depth to lag behind surface temperature changes.

A range of approaches to formulating liquid-water saturation and relative permeability have been applied in unsaturated flow simulations and there is uncertainty in optimal function specification. However, sensitivity analysis demonstrates that firn temperature profiles are not sensitive to these functions. For example, no changes in the simulated vertical temperature profiles result from increasing the minimum liquid saturation by a factor of 10 or increasing the exponential model parameter by a factor of 100 in the liquid water saturation function or increasing the relative permeability at residual saturation in the relative permeability function by a factor of 500 (Fig. 3). This indicates that these choices of simplified state functions for liquid water saturation and relative permeability do exactly what was intended: to represent the most important freezing and permeability processes that control the aquifer configuration and primary dynamics. For this initial multidimensional analysis of firn aquifer physics, these selections have allowed a successful basic evaluation of firn aquifer configuration and response.

4.2. Surface water expression and implications of higher recharge rates

Although water levels vary within the aquifer with time, generally corresponding to melt rates (Miège and others, 2016; Miller and others, 2020), the difference between model and observed water table depths for the steady-state simulation with 15 cm a⁻¹ recharge suggest that the aquifer is not in equilibrium currently (i.e. not at steady state) for current climatic conditions. This is similar to a firn aquifer in the Wilkins Ice Shelf that does not appear to be in equilibrium based on simple fluid flow modeling that does not consider thermodynamics (Montgomery and others, 2020). As surface melt and recharge increases under a warming climate, the firn aquifer saturated zone is expected to grow in volume. The simulation-predicted inland aquifer expansion with greater recharge rates broadly agrees with SNOWPACK firn simulations (Steger and others, 2017a) as well as with observed inland expansion (Montgomery and others, 2017; Miller and others, 2020).

For recharge rates similar to those observed in the Helheim Glacier field area (15 cm a⁻¹), the modeled steady-state water table intersects the snow surface (Figs 4–6) and this would generate surface water features such as streams and lakes. In land-based hydrologic systems, lakes and streams occur at locations where water tables intersect the land surface. The snow surface topography, which features slope variation between 0 and 2° (Miège and others, 2016), would control the surface water feature type. Snow surface topography can be controlled by underlying subglacial topography or basal friction (Gudmundsson, 2003; Sergienko, 2013). Lakes could form in local topographic depressions and streams could form where the snow surface slopes more continuously downhill. The water table is generally closer to the snow surface at local surface topographic lows where lakes would form (Miège and others, 2016).

Currently, firn aquifers and supraglacial lakes rarely co-occur (Koenig and others, 2015; Miège and others, 2016). Southeast Greenland generally has few lakes compared to other parts of the ice sheet (Sundal and others, 2009; Selmes and others, 2011) and the absence of lakes in southeast Greenland has been attributed to steeper snow surface slopes (Sundal and others, 2009). At the Helheim field site, there is no evidence of streams or lakes. The ice sheet in southeast Greenland is characterized by thick permeable firn and relatively steep snow surface slopes that facilitate meltwater infiltration and aquifer formation and hinder surface water formation. These permeable firn conditions hinder the formation of streams and lakes by preventing meltwater from ponding and flowing along the low permeability ice surface, which, in contrast, occurs elsewhere on the GrIS (Echelmeyer and others, 1991; Selmes and others, 2011; Smith and others, 2015, 2017; Chen and others, 2017). Over time however, under current recharge conditions, the steady-state simulation results indicate that the water table in the Helheim firn aquifer may rise to the snow surface, forming lakes or streams. Thus, should the current climatic conditions persist and climate warming continue, streams or lakes may eventually appear in areas where firn aquifers intersect the snow surface.

Formation of such surface streams or lakes would dramatically alter the surface energy balance, possibly contributing to increased melt via reduced albedo. A surface hydrologic system of streams would presumably enhance flow of meltwater from the ice sheet to the ocean as well. Cold winter air temperatures could also freeze shallower liquid water in a firn aquifer, potentially forming thick ice slabs (MacFerrin and others, 2019), which may locally impede deep infiltration of snow melt.

This prediction of groundwater discharge at the snow surface highlights the critical link long-observed in land-based hydrology between surface water (supraglacial) and groundwater (englacial) hydrology, where they are considered a single ‘resource’ (Winter and others, 1998), that has not yet been emphasized in ice-sheet hydrology. Much emphasis has been placed on the connection between supraglacial and subglacial hydrology and the associated relevance for ice dynamics (Joughin and others, 1996; Zwally and others, 2002; Das and others, 2008), but as surface mass-balance processes have become the dominant control on ice-sheet mass loss and contribution to sea-level rise (Enderlin and others, 2014), the link between supraglacial and englacial hydrology also becomes more important.

Increased surface melt will likely increase discharge from firn aquifers, possibly resulting in a greater direct contribution of this fluid to sea-level rise as mentioned above. Increased aquifer discharge could also generate new and deeper crevasses via hydraulic fracturing due to elevated standing water levels in existing cracks and crevasses. Increased aquifer discharge to crevasses that reach the glacier base could also lead to increases in glacier movement via increases in sub-glacial fluid pressure and concomitant decreases in sub-glacial friction due to groundwater water discharge to crevasses.

4.3. Future work

As an approximation of general conditions and behaviors of the firn aquifer system that is based on full multidimensional hydraulic and thermodynamic physics, the steady-state simulations provide a robust test of the physical self-consistency of the conceptual firn-aquifer model and provide insight into the main factors that control the configuration of firn aquifers. However, the steady-state 2D simulations considered here do not incorporate time-varying or spatially varying surface temperature and recharge, snow accumulation, firn burial, preferential flow paths through heterogeneous firn or possible localized perching of meltwater upon discontinuous ice lenses above the water table that occur in reality. As a result, these results do not account for complexities of percolation zone processes including the heat and liquid that goes into warming the cold unsaturated zone in the spring or winter accumulation and cold content addition. Future work could implement time-variable surface temperatures and recharge and the impact of ice layers in firn on infiltration and firn aquifers. Field and lab studies to constrain the total-water saturation, liquid water saturation and relative permeability functions, as well as the possible effects of cryosuction on liquid and ice saturation distributions, would allow future modeling to study firn aquifers with more local detail, especially in the unsaturated zone. Such advances would improve estimates of aquifer volume changes and discharge, aid in study of crevasse hydrofracturing and improve estimates of pressurization of basal (sub-glacial) meltwater that may increase glacial motion, which in turn may increase the global rate of sea-level rise.

5. Conclusions

This study presents the first numerical simulation of multi-dimensional liquid water flow and thermal energy transport with a code (SUTRA-ICE) where groundwater can freeze and thaw depending on ambient temperature and also on fluid pressure when freezing occurs in an unsaturated zone. Firn aquifers should be represented by multi-dimensional models that allow for the substantial lateral fluid flow and energy transport that define them. This code provides the requisite self-consistent hydraulic and thermodynamic physics basis, and it was used in this initial study to test the conceptual model of the firn aquifer upslope from Helheim Glacier, identify primary processes and factors that control overall firn aquifer configuration, and evaluate long-term firn aquifer configuration as a function of meltwater recharge rate.

As in all mathematical descriptions of physical processes, a description can and should be tested for its ability to reproduce observations, and if it does so successfully, this particular description is an acceptable representation of the modeled system. This result never proves that the description is the only one that will reproduce observations – but it indicates that this particular model stands as a feasible and useful description of the system. Perhaps, in the future, this particular description will be unable to reproduce new data that may be collected, limiting the usefulness of the description. There may also be other different mathematical descriptions that reproduce the currently existing data, and as such, these too may be considered as acceptable representations of the modeled system. Thus, the conceptual model presented here that includes a mathematical description of firn-aquifer behavior may be considered as an acceptable description.

Another ambiguity in postulating a mathematical description of a physical system is that, given enough parameters in the description, the ability of the description to fit measurements may be due primarily to the large number of parameters, which may allow even a model that only poorly describes the main physical processes to fit the measured data. For the current firn-aquifer model, this (over-parameterization) was avoided by keeping the number of parameters as low as possible and by representing the physical system as simply as possible, forcing the dynamics of the specific physical processes included in the conceptual model to fit the data as well as possible. The result here is that the conceptual model, and the physics it represents, is a strong descriptor of firn-aquifer configuration (especially for steady-state conditions in 2-D cross-sectional representations of these systems).

This is the first use of a SUTRA-ICE multi-dimensional flow and transport type of code to simulate meltwater retention and unsaturated and saturated flow in a firn aquifer system. This type of code is expected to be widely applicable to firn aquifer simulations in the future, because it represents an improvement in modeling multi-dimensional meltwater flow through unsaturated and saturated firn. This successful use of SUTRA-ICE indicates that this code and other similar codes can be more broadly applied to future studies of firn aquifers or meltwater flow through snow.

The results of modeling analyses indicate that the following factors are important controls on the simulated aquifer: (1) surface temperature, (2) meltwater recharge rate, (3) residual total-water saturation (minimum saturation of liquid plus ice) and (4) thickness of the capillary fringe (only if it extends to the snow surface). Firn temperatures are not sensitive to the ‘backbone’ fraction of ice (the amount assumed to never melt in the simulation for the temperature range considered), nor to the liquid saturation range over which the relative permeability drops from its maximum to minimum value, nor to the minimum relative permeability value.

SUTRA-ICE simulation results indicate that the basic conceptual model of the physics of firn aquifer formation and function that was developed on the basis of field data is consistent with the fundamental physics of groundwater flow and with the thermodynamics of freezing and thawing as these occur within an aquifer consisting entirely of water in its liquid and solid states. Differences between simulated steady-state liquid water levels and observed liquid water levels may indicate that the Helheim Glacier firn aquifer is not at steady state for its current recharge rates. Should the aquifer attain steady state in the future, or if recharge increases under a future warmer climate, the water table will intersect the snow surface. This would generate lakes and/or streams co-located with the firn aquifer, representing a fundamental shift in ice-sheet hydrologic facies.