Time-Domain Reflectometry Measurements and Modeling of Firn Meltwater Infiltration at DYE-2, Greenland.

Surface meltwater can be retained in an ice sheet if it infiltrates the firn and refreezes. This is an important mass balance process for the Greenland Ice Sheet, reducing meltwater runoff and associated sea-level rise. The processes of meltwater infiltration and refreezing are not fully understood, however, and remain difficult to monitor remotely. We deployed vertical arrays of thermistors and time-domain reflectometry (TDR) probes to 4-m depth in the firn to continuously monitor meltwater infiltration and refreezing processes at DYE-2, Greenland. The observations provide a detailed picture of the coupled thermal and hydrological evolution of the firn through the 2016 melt season, including estimates of firn water content. The thaw and wetting fronts reached a maximum depth of 1.8 m, with meltwater infiltration concentrated in four main pulses of melting and subsurface warming that reached progressively deeper into the firn. The observations were used to constrain a coupled model of firn thermodynamics and hydrology, which was then run over the period 1950-2020, driven by meteorological forcing from GC-Net station data and ERA5 climate reanalyses. Model results suggest that decadal-scale firn evolution at DYE-2 is strongly influenced by extreme melt seasons such as those of 1968, 2012, and 2019, when meltwater infiltration reached depths of 6-7 m. Extreme melt years drive increases in firn temperature, ice content, and density, reducing firn meltwater retention capacity. Such processes are likely to govern future meltwater retention as the percolation zone extends to higher elevations in Greenland in the coming decades.

Introduction

In the accumulation area of the Greenland Ice Sheet, meltwater percolates into the underlying seasonal snow and firn, where it can refreeze (Benson, 1962; Braithwaite et al., 1994; Pfeffer et al., 1991). Meltwater retention from this process is well documented in the Greenland percolation zone, and it is a critical mass balance process with the potential to buffer sea level rise associated with ice sheet melt (Harper et al., 2012; Rennermalm et al., 2013). Meltwater infiltration and refreezing are expected to increase in importance as the Greenland Ice Sheet percolation zone expands to higher elevations in a warming world (MacFerrin et al., 2019; Noël et al., 2017; van Angelen et al., 2013).

Meltwater that percolates and refreezes is difficult to account for in altimetric measurements of surface mass balance, as it is not possible to observe how much meltwater is retained within the system (Kuipers Munneke et al., 2015; Sørensen et al., 2011). Meltwater retention is also a source of uncertainty in mass balance models in Greenland, due to the complexities of meltwater infiltration and ice layer formation processes (e.g., Langen et al., 2017; Reijmer et al., 2012; Steger, Reijmer, van den Broeke, et al., 2017; Vernon et al., 2013). Processes of meltwater percolation and refreezing occur at fine scales, are spatially heterogeneous, and in situ observations are scarce (Humphrey et al., 2012; van As et al., 2016). This makes it difficult to calibrate and validate snow and firn models, particularly on the scale of polar ice caps and ice sheets (Steger, Reijmer, & van den Broeke, 2017; Vandecrux, Mottram, et al., 2020; Verjans et al., 2019).

Meltwater percolation in snow can take the form of relatively homogeneous infiltration (e.g., Marsh & Woo, 1984) or vertical ‘piping’ of water, via preferential flow pathways (Humphrey et al., 2012; Marsh & Woo, 1984; Pfeffer & Humphrey, 1998; Williams et al., 2010). Meltwater may percolate to a certain depth before it reaches an impermeable ice layer or refreezes (e.g., Colbeck, 1979). Subsurface grain size, density, or permeability contrasts can cause meltwater to spread horizontally, leading to the development of ice layers (Avanzi et al., 2016; Katsushima et al., 2013). With sufficient melt, individual ice lenses can merge into thick, horizontally extensive ice layers, which inhibit meltwater percolation (de la Peña et al., 2015; Gascon et al., 2013; MacFerrin et al., 2019; Machguth et al., 2016).

At high elevations in the Greenland Ice Sheet percolation zone, cooler temperatures and lower melt rates result in greater near‐surface pore space, cold content, and hydraulic permeability, promoting meltwater retention (Vandecrux, Fausto, et al., 2020; Verjans et al., 2019). Where melting is moderate, ice lenses and layers are typically discontinuous over horizontal length scales of less than 1 m (e.g., Dunse et al., 2008; Parry et al., 2007). There is evidence that meltwater can effectively infiltrate ice‐rich firn in this environment (Machguth et al., 2016), including ice layers that are up to ∼0.1 m thick (Samimi et al., 2020). However, climate warming is expected to cause increases in firn temperature, density, and near‐surface ice layers in the upper percolation zone of the ice sheet in the coming decades (de la Peña et al., 2015; MacFerrin et al., 2019; Machguth et al., 2016; Noël et al., 2017; van Angelen et al., 2013), reducing meltwater storage capacity and potentially driving a transition from meltwater retention to runoff. Future projections of Greenland Ice Sheet response to climate warming and associated sea‐level rise require a quantitative understanding of these processes (e.g., Langen et al., 2017; van As et al., 2016).

It is uncertain whether firn evolution is more sensitive to long‐term warming (e.g., Hanna et al., 2012, 2021) or the impacts of extreme melt seasons such as those witnessed in 2012 and 2019 (Hanna et al., 2014; Nghiem et al., 2012; Tedesco & Fettweis, 2020). These extreme melt seasons were associated with strong anticyclonic conditions that persisted for days to weeks (Fettweis et al., 2013; Hanna et al., 2014; Rajewicz & Marshall, 2014), causing anomalous melt which extended to high elevations of the ice sheet in both 2012 (Nghiem et al., 2012) and 2019 (Tedesco & Fettweis, 2020). Large volumes of meltwater during such events can create thick near‐surface ice slabs (Machguth et al., 2016; Noël et al., 2017), which are likely to be impermeable to meltwater infiltration in subsequent years. Extreme melt seasons also drive deep meltwater infiltration, firn warming, and densification, reducing meltwater retention capacity.

To advance understanding of meltwater infiltration processes in the Greenland percolation zone, we designed field experiments at DYE‐2 in southwest Greenland to continuously monitor the coupled thermal and hydrological evolution in the near‐surface firn. Vertical arrays of thermistors and time‐domain‐reflectometry (TDR) sensors were deployed to 4.3‐m depth to track meltwater infiltration through the 2016 melt season. The thermistors monitor firn temperature and the latent heat signature of meltwater refreezing (e.g., Humphrey et al., 2012; Charalampidis et al., 2016). TDR probes provide additional information because the dielectric permittivity of snow and firn is strongly sensitive to small amounts of liquid meltwater (Schneebeli et al., 1997; Techel & Pielmeier, 2011). Dielectric permittivity therefore provides a quantitative measure of liquid water content and the TDR arrays directly track the seasonal progression of the wetting front, in conjunction with its thermal signature.

Samimi et al. (2020) presented the temperature and dielectric permittivity data from the thermistor and TDR arrays at DYE‐2. Initial results in Samimi et al. (2020) demonstrated that the experimental setup provided a high‐quality continuous record of the coupled thermal and hydrological evolution of the near‐surface firn through the 2016 melt season, including evidence of meltwater percolation through continuous ice layers as much as 0.12‐m thick. This study builds on Samimi et al. (2020) through estimates of liquid water content from the dielectric permittivity data and assessments of the meltwater infiltration velocity and irreducible water content in the firn. The observations are used to constrain and tune a multilayer model which simulates the coupled subsurface hydrological and thermal evolution in the upper 20 m of snow and firn. The calibrated model is then run through historical reconstructions from 1950 to 2020 to place the 2016 melt season in context and examine decadal‐scale firn evolution at the site. Model simulations from 1997 to 2018 are forced by GC‐Net Automatic Weather Station (AWS) data from DYE‐2 (Steffen & Box, 2001) and are extended to the period 1950–2020 using bias‐adjusted meteorological forcing from the ERA5 climate reanalysis (Hersbach et al., 2020).

The field experiments and modeling studies were designed to improve understanding of meltwater infiltration and retention processes in the Greenland percolation zone. Our specific objectives are to evaluate the depth of seasonal meltwater infiltration in firn, to characterize the liquid water content and meltwater infiltration velocities in firn, to constrain firn hydrological modeling efforts, and to assess the interannual variability and decadal‐scale evolution of firn temperature, density, and ice content at DYE‐2. We use the calibrated firn model to examine the historical frequency of extreme melt events since 1950, their impacts on the firn, and the significance of such events in comparison with the long‐term climate and firn evolution.

Study Site and Research Methods

DYE‐2, Greenland

Our field study was part of the FirnCover project, which aimed to increase understanding of firn densification and meltwater retention processes in the Greenland Ice Sheet percolation zone (MacFerrin et al., 2019; Machguth et al., 2016). In April 2016, we established observation sites in the near‐surface firn at 66˚28′39″N, 46˚17′5″W, near DYE‐2 station on the southwestern flank of the ice sheet (Figure 1). The study sites were at an elevation of 2120 m, in the upper part of the percolation zone of southwestern Greenland. Summer melt in this zone is retained through meltwater refreezing, but commonly infiltrates beyond the annual snow layer (Benson, 1962). DYE‐2 has a slope of less than 1° and there is no evidence of surface water ponding or horizontal surface water flow.

Figure 1

(a and b) Study area at DYE‐2, Greenland, indicating the firn pit and Automatic Weather Station (AWS) locations. (c) AWS at site A. (d) The upper part of firn pit A showing the TDR sensors and examples of ice layers at this site, prior to filling in the pit. (e) Firn pit B after filling it in (foreground, with the datalogger box, ultrasonic depth gauge [UDG], and temperature sensor). Background map in Figure (a) courtesy of Google Earth.

A 100.2‐m ice core collected at DYE‐2 in 1977 indicated a mean annual accumulation rate of 0.34 m w.e.  for the period ∼1740–1976 (Clausen et al., 1988; Shoji et al., 1991). This has increased by ∼10% in recent decades, as a mean annual accumulation of 0.365 m w.e.  at DYE‐2 was inferred from firn cores acquired by Ohio State University (OSU) in 1998, representing the period 1969–1998 (Mosley‐Thompson et al., 2001). From 2013 to 2016, shallow firn cores were drilled each year by FirnCover project participants at DYE‐2. Machguth et al. (2016) analyzed the 2013 and 2015 firn density and stratigraphy and compared their results with the 1998 OSU ice cores. From 1998 to 2013, density increased by ∼140 kg m−3 in the top 7 m of firn and ∼100 kg m−3 from 8 to 16 m depth. The increases in density are associated with higher ice content (i.e., increases in refrozen meltwater), and meltwater penetration heavily modified the pre‐1998 firn stratigraphy recorded by Mosley‐Thompson et al. (2001). The Machguth et al. (2016) study indicated that DYE‐2 may be undergoing a transition to denser and more ice‐rich firn, characteristic of the lower percolation zone.

Our study sites were within 1 km of the Greenland Climate network (GC‐Net) automatic weather station established at DYE‐2 in 1996 (66°28′50″N, 46°16′59″W; Steffen & Box, 2001). Mean annual and mean summer (June through August, or JJA) 2‐m air temperatures at the GC‐Net station were T ann = −17.4°C and T JJA = −5.4°C from 1997 to 2018. Based on the GC‐Net data, the melt season at DYE‐2 typically runs from mid‐June to mid‐August. DYE‐2 experiences an average of 28 melt days each summer (JJA), as calculated from days with maximum hourly air temperatures exceeding 0°C. Annual positive degree days (PDD) at the site averaged 8.6°C d from 1997 to 2018.

On average, firn temperatures at 10‐m depth have been ∼2.5°C warmer than the mean annual air temperature at DYE‐2, due to latent heat release from refreezing meltwater (Steffen & Box, 2001). However, there is evidence of firn warming at DYE‐2 over the last four decades. Clausen et al. (1988) reported a 10‐m firn temperature of −17.2°C at DYE‐2, based on measurements made during the 1977 ice‐coring campaign. This warmed to −16.5°C in the late 1990s (Steffen & Box, 2001), and GC‐Net thermistor measurements at DYE‐2 indicate an additional ∼1.5°C of warming from 1998 to 2009, to a 10‐m temperature of −15°C. Consistent with this, Vandecrux, Fausto, et al. (2020) reported a warming of 1.1°C per decade at 10‐m depth at DYE‐2, based on firn modeling over the period 1998–2018. Firn warming exceeded atmospheric warming over this period, indicating an increase in meltwater production and refreezing at this site.

DYE‐2 appears to have remained within the accumulation area of the ice sheet for the entire period of GC‐Net observations, with seasonal snow persisting each year, although the end‐of‐summer snowline (ELA) advanced near to this elevation in 2012, the warmest summer on record. DYE‐2 experienced 66 melt days in 2012, with a mean summer temperature T JJA = −2.9°C and PDD totals of 37.2°C d. The mean JJA surface albedo in 2012 was 0.73, lower than normal (0.79) but diagnostic of persistent snow cover. The heavy melt season of summer 2012 may account for much of the deep meltwater infiltration, increased ice content, and densification recorded by Machguth et al. (2016) in the 2013 firn core, relative to the 1998 OSU cores.

Specific to the period of study discussed in this manuscript, the summer of 2016 was 0.8°C warmer than mean conditions from 1997 to 2018, with T JJA = −4.6°C and a total of 37 melt days. The winter snowpack was ∼0.9 m deep in mid‐April of 2016, representing an accumulation of 325 mm w.e. In addition to our TDR experiments to monitor meltwater infiltration through the summer 2016 melt season, Heilig et al. (2018) installed an upward‐looking ground‐penetrating radar system at 4.3‐m depth to monitor the evolution of the wetting front and the snow/firn liquid water content. The radar measurements indicate a total melt of ∼0.15 m w.e. at DYE‐2 in the summer of 2016 (Heilig et al., 2018). Of this, 0.09 m w.e. (60%) refroze within the seasonal snowpack and 0.06 m w.e. (40%) infiltrated the firn, with a maximum meltwater infiltration depth of 2.3 m below the surface.

Firn Measurements and TDR Methods

Thermistor and TDR Measurements

We established TDR and thermistor arrays at two locations near DYE‐2 in April, 2016, excavating firn pits to depths of 2.2 and 5.3 m (Figure 1). The two sites, referred to as A and B, are 400 m apart (Figure 1b). The 2015–2016 seasonal snowpack was at subzero temperature and free of ice layers at the time of installation of the sensors. We define firn as the multiyear snow that underlies the seasonal snowpack. Snow and firn density were measured at 10‐cm intervals in the firn pits, using a 100‐cm3 box cutter. The seasonal snow had an average density of 290 kg m−3. The underlying summer 2015 melt surface took the form of an ice layer (i.e., a refrozen crust), which was difficult to penetrate with an avalanche shovel. Below this, the firn was made up of a mixture of dense snow, ice lenses and ice layers. There were 32 discrete ice layers in the upper 5.3 m of snow and firn, with a total ice thickness of 1.08 m. This is equivalent to an average ice‐layer thickness of 3.4 cm. Roughly half of the ice layers (15) were less than 1‐cm thick and three layers exceeded 10 cm in thickness, with the maximum being a 33‐cm thick layer extending from 3.79 to 4.12 m depth. Several ice layers were horizontally continuous over the ∼2 m width of the firn pits (cf. Figure 1c).

The north‐facing vertical face of each firn pit was instrumented with eight thermistors and eight time‐domain reflectometry (TDR) probes to monitor snow water content (Table 1, Figure 1c). Sensors were inserted horizontally, with the thermistors and TDR probes extending 10 and 30 cm into the wall of the firn pit, respectively. Sensor spacing was irregular in order to concentrate observations near the surface as well as immediately above and below ice layers, to test whether these acted as impermeable barriers to water flow. The firn pits were filled in with snow after sensor installation. Thermistors and TDR probes were wired to Campbell Scientific CR1000 dataloggers and values were recorded each 30 min from May 11 to September 30, 2016, capturing the subsurface hydrological and thermal evolution over the complete melt season.

Table 1

Installation Depths (d) and Initial Snow/Firn Temperature (T ), Bulk Dielectric Permittivity (ε ), and Density (ρ 0) at Each Measurement Level for Sites A and B

Site A (installed May 11, 2016)					Site B (installed May 8, 2016)
66.47775°N, 46.28510°W					66.47505°N, 46.29115°W
Level	d (m)	T s0 (°C)	ε b0	ρ 0 (kg m−3)	d (m)	T s0 (°C)	ε b0	ρ0 (kg m−3)
1	0.3	−7.3	2.37	380	0.1	−3.7	2.31	280
2	0.6	−10.1	1.88	380	0.2	−4.2	2.49	310
3	0.9	−12.7	2.21	230	0.4	−6.6	2.09	280
4	1.4	−14.5	2.32	510	0.6	−9.2	2.06	320
5	1.8	−15.1	2.37	410	0.9	−11.9	2.47	640
6	2.1	−15.6	2.44	460	1.2	−13.1	2.54	620
7	2.8	−15.7	2.14	360	1.4	−14.4	2.65	560
8	3.7	−15.8	2.67	520	1.6	−14.9	2.36	480

Firn Cores

Firn cores reaching depths of 5–21 m were drilled with a Kovacs ice coring system in 2016 and 2017, including cores adjacent to each firn pit. Core densities were measured at 5‐cm vertical resolution and visual stratigraphies were carried out at 1‐mm vertical resolution. The ice‐layer stratigraphy in the firn cores supports the interpretation of the TDR data and also informs the total ice content in the firn, a measure of average melt/refreezing rates at the site.

TDR Data Analysis

Snow is made up of a mixture of ice, air, and liquid water, with differing dielectric properties for each constituent. We adopt fixed values for the relative dielectric permittivity of air, ice, and water:  = 1, , and , respectively (Evans, 1965; Lundberg, 1997). Because of these differences, the dielectric permittivity of snow increases strongly with liquid water content, making the bulk value, , a sensitive indicator of snow moisture. For porosity θ and volumetric liquid water fraction , the bulk density of snow or firn is where is the density of ice crystals in the snow matrix (917 kg m−3) and and are the densities of air and water (1.1 kg m−3 and 1000 kg m−3, respectively). The dry density, , is calculated from Equation 1 when : for dry snow/firn at temperatures below 0°C. This was the case on installation of our sensors, when we measured the initial snow and firn densities (Table 1). Snow density and water content both influence dielectric permittivity because solid and liquid water molecules are polarized in the presence of an electric field; a greater mass of these molecules increases polarization and relative permittivity (Schneebeli et al., 1997; Stein et al., 1997).

Measured dielectric permittivity needs to be converted to water content based on a dielectric mixing model. Empirical equations that relate TDR‐derived dielectric permittivity to liquid water content have been developed for soils (e.g., Ledieu et al., 1986; Topp et al., 1980), but these are not appropriate for snow or firn. Tiuri et al. (1984) and Denoth (1994) introduce empirical equations relating the bulk dielectric permittivity to snow density and water content. Stein et al. (1997), Schneebeli et al. (1997), Lundberg (1997), and Techel and Pielmeier (2011) discuss the specific application of TDR methods to snow, including empirical equations to separate the effects of snow density and liquid water content on the bulk dielectric permittivity.

Following work in soil water studies (Birchak et al., 1974; Roth et al., 1990), three‐component mixing models have also been applied to analyses of bulk dielectric permittivity in snow, :

Subscripts i, w, and a refer to ice, water, and air, for dielectric permittivty and volume fraction θ. Exponent β is equal to 0.5 for the Birchak et al. (1974) model, based on the relation between bulk dielectric permittivity and the volume‐averaged wave velocity in a dielectric medium. Several recent snow radar studies also apply an exponent of 0.5 in dielectric mixing models (e.g., Heilig et al., 2015, 2018; Schmid et al., 2014). This choice of exponent was found to work well for TDR monitoring of meltwater content in supraglacial snow (Samimi & Marshall, 2017), giving:

The porosity θ can be calculated from Equation 1 with dry density measurements (i.e., prior to the onset of the melt season), θ ≈ 1 − ρ /ρ . To update θ through the summer melt season, we include a model of snow and firn densification following Vionnet et al. (2012), for overburden pressure σ and snow/firn viscosity η. Vionnet et al. (2012) developed this parameterization of viscous deformation for the seasonal snowpack, but it has been applied to firn densification in several studies (e.g., Gascon et al., 2013; Verjans et al., 2019). The formulation includes a parameterization of the effects of temperature and liquid water content on snow viscosity. This model for densification applies to the “background” firn density (i.e., the firn matrix), and liquid water and refrozen ice content act in addition to this to increase bulk density, following Equation 1.

Meteorological Data

Automatic Weather Station Data, 2016

An automatic weather station (AWS) configured for surface energy balance monitoring was installed adjacent to site A in April 2016 (Figures 1b and 1c). The AWS measured air temperature, humidity, incoming and outgoing longwave and shortwave radiation, wind speed and direction, air pressure, and snow surface height every 10 s, with 30‐min averages saved to the datalogger. Snow surface height was measured by an ultrasonic depth gauge (UDG) mounted on an independent stake that was drilled 2 m into the firn. Separate meteorological sensors were installed at site B to measure air temperature, humidity, and snow surface height (Figure 1e). All sensors were left in place for one year, from April 2016 to April 2017.

The UDGs track both surface drawdown and fresh snow accumulation. This provides information on the evolving burial depth for the TDR probes and thermistors and the potential melt‐out of the sensors. Surface drawdown can occur through wind scour, sublimation, melting, and snow/firn densification, so interpretation is required to infer surface melt from the UDG records. Increases in surface height due to snowfall or windblown snow accumulation are also recorded by the UDG sensors, and can usually be corroborated through the measured increases in albedo that accompany fresh snow. While we did not see evidence of it, gravitational settling (sinking) of the UDG pole can also occur, masquerading as snow accumulation or causing an underestimate of the ablation. We drilled the UDG poles into dense firn at 2‐m depth to minimize this effect.

We combine different lines of evidence to infer summer melt totals (in m w.e.) from these records. Initial measurements of snow density provide data for the snow‐water equivalent of the snowpack and the density stratigraphy. We use the densification model (Equation 4) to estimate vertical motion due to settling of the initial snow and firn stratigraphy and assume that the remaining surface drawdown signal is due to melting. Where snow accumulation is recorded by the UDGs, we assume that this represents fresh snowfall with an initial snow density of 120 kg m−3, after Heilig et al. (2018). We assign an uncertainty 20% to the UDG ablation estimates, given the various sources of uncertainty in interpreting these records. As an additional constraint, thermistors and TDR sensors that melt out have distinctive signals (e.g., daytime warming well above 0°C, which cannot occur in snow). This provides an independent estimate of the amount of ablation that occurred at a site, based on the snow‐water equivalent that was removed above any exposed sensors. We assume that any such ablation occurred due to melt. This can be cross‐checked against the modeled melt from the surface energy balance, which provides an additional independent estimate of the melt.

AWS data are used as input to a surface energy balance model to estimate surface melt rates (Ebrahimi & Marshall, 2016) and inform the interpretation of the thermistor and TDR data. We also use the 30‐min AWS data to drive a model of firn evolution through the 2016 melt season, as described in Section 2.4. Thermistor and TDR data are used to constrain and tune the firn model, and we then apply the calibrated model to analyses of multidecadal firn evolution at DYE‐2, using meteorological forcing from the GC‐Net AWS (Steffen & Box, 2001) and ERA5 climate reanalyses (Hersbach et al., 2020).

GC‐Net Weather Station Data, 1997‐2018

For applications requiring more remote meteorological forcing (i.e., nearby weather stations; output from climate reanalyses or climate models) we use a different configuration of the surface energy balance model. Specifically, outgoing shortwave and longwave radiation are calculated internally in the surface energy balance and firn model, based on the modeled surface albedo, α , and the surface‐layer temperature. The snow albedo parameterization is described in Marshall and Miller (2020) and includes a linear melt‐season decay as a function of cumulative positive degree days, PDD, where  = 0.80 is the initial (spring) albedo and  = 0.002°C−1 d−1 is the decay constant used in the modeling, based on tuning to the observed albedo data at DYE‐2. This parameterization roughly accounts for the reduction in snow albedo that occurs through the melt season due to the combined effects of increasing snow grain size, grain rounding, liquid water content, and progressive concentration of impurities in melting snow (Brock et al., 2000). We also include a stochastic representation of summer snow events, which transiently refresh the surface albedo to the dry‐snow value through the summer melt season (Marshall & Miller, 2020).

Required input fields to drive the decadal simulations include air temperature, relative humidity, wind speed, air pressure, and incoming shortwave and longwave radiation. With the exception of incoming longwave radiation, these fields are available from the GC‐Net DYE‐2 AWS data from 1997 to 2018. These data contain numerous gaps (see Vandecrux et al., 2018); the data were quality‐controlled and gap‐filled to provide mean daily values of all missing or compromised meteorological variables. Gap filling was based on random sampling of the statistical distribution of all available data for a given day. For incoming longwave radiation, we use hourly ERA5 reanalysis data from the nearest 0.25° grid cell to DYE‐2, at 66.5°N and 46.25°W. ERA5 incoming longwave radiation is biased‐adjusted using daily mean bias corrections, based on comparison with our AWS data for the period April–October 2016. The average bias adjustments for April and October are used for the months of November to March. Precipitation data (monthly snow accumulation) are also taken from ERA5 for the period 1997–2018.

The energy balance and firn models are run on 30‐min or hourly time steps. Where only daily mean data are available (i.e., for the gap‐filled data), we parameterize the diurnal cycles of temperature and incoming shortwave radiation to construct 30‐min values from mean daily fields (Ebrahimi & Marshall, 2016). This parameterization requires minimum, maximum, and mean daily temperature and the mean and maximum daily incoming shortwave radiation.

ERA5 Climate Reanalyses, 1950–2020

Hourly ERA5 climate reanalyses are available at a resolution of 0.25° for the period 1950–2020 (Hersbach et al., 2020). We compiled this into daily mean data for the grid cell at 66.5°N and 46.25°W, which contains DYE‐2. Incoming longwave radiation was bias‐adjusted as noted above. For the remaining ERA5 data, daily bias adjustments are applied based on mean deviations in the ERA5 data from the gap‐filled GC‐Net data over the common period 1997–2018. The firn model is forced by the mean daily bias‐adjusted ERA5 meteorological fields, with parameterizations of the diurnal temperature and incoming shortwave radiation cycles and surface albedo and outgoing longwave radiation modeled internally, as noted above.

Daily mean meteorological forcing is used for the ERA5 simulations for pragmatic reasons. It is simpler to bias‐adjust and gap‐fill than 30‐min or hourly data and it is also readily available from a wide range of climate reanalyses and future climate projections. As this is of interest to follow‐up work with ice‐sheet scale firn modeling, we focus on tuning and calibration of the surface energy balance and firn model to these inputs. Ebrahimi and Marshall (2016) demonstrated that modeled surface energy balance and melt using daily mean forcing is in good accord with hourly forcing, as long as diurnal cycles of air temperature and incoming shortwave radiation are accounted for.

Surface Energy Balance and Firn Model

Model Physics and Parameterizations

Surface energy balance and melt rates are calculated following the model of Ebrahimi and Marshall (2016), including a subsurface model for heat conduction in the upper 20 m of firn and snow. The model has 43 layers, with resolution concentrated near the surface; layers are 0.1‐m thick from the surface to a depth of 0.6 m, 0.2‐m thick from 0.6 to 2 m, 0.4‐m thick from 2 to 10 m, and 1‐m thick from 10 to 20 m. Net energy, Q , is calculated at 30‐min intervals as the sum of the energy fluxes towards the surface layer, where the terms on the right represent (in order) the incoming and outgoing shortwave radiation, incoming and outgoing longwave radiation, and the sensible, latent, and conductive heat fluxes. The model requires inputs from either AWS data or climate reanalyses/models for the incoming shortwave and longwave radiation, wind speed, air temperature, humidity, and air pressure. Sensible and latent heat flux are calculated from these meteorological fields and conductive heat fluxes are calculated within the subsurface model, based on the temperature gradient in the upper three layers. Reflected shortwave radiation Q ↑ = Q ↓ (1 − α ) and outgoing longwave radiation, Q ↑, follows Stefan‐Boltzmann's equation, Q ↑ = E σT 4, for surface emissivity E  = 0.98, Stefan‐Boltzmann's constant σ, and absolute temperature of the surface layer, T . Albedo and surface temperature are calculated internally (Equation 5, Ebrahimi & Marshall, 2016), or AWS data can be used directly for these where data are available (i.e., summer 2016).

When the surface layer is at 0°C and , net energy goes to melting, following where is the melt rate (m s−1), and L is the latent heat of fusion. If net energy is negative and the surface layer is 0°C, any liquid water that is present will refreeze, releasing latent heat, and the surface layer will cool once all liquid water is refrozen. When surface layer temperatures are below the melting point, and surplus or deficit of energy drives warming or cooling, based on the energy balance solution within a one‐dimensional model of subsurface temperature evolution:

The right‐hand terms represent heat conduction, latent heat release from meltwater refreezing, and heat transport from advection of meltwater or rainwater, respectively. Here , c , and k are the bulk density, specific heat capacity, and thermal conductivity of the snow or firn, and c are the density and specific heat capacity of water, and is the vertical rate of water percolation, with units m s−1. We calculate the thermal conductivity of the snow and firn after Calonne et al. (2019). The specific heat capacity of ice is taken from Cuffey and Paterson (2010), c  = 152.2 + 7.122 T , for absolute temperature T . Bulk heat capacity is then calculated from .

The refreezing term in Equation 8 has units W m−3 and is calculated from where is the refreezing rate (m s−1) and this heat is spread across the layer thickness, Δz.

The subsurface temperature model is coupled with a simple treatment of meltwater percolation. For a snow or firn layer with thickness Δz and the volume fraction of liquid water the amount of water in the layer (m) is equal to . Conservation of mass in each subsurface layer gives the expression for local water balance,

Any water that refreezes is assumed to be distributed over the layer Δz. We assume that all meltwater flow is vertical (gravitational drainage with no horizontal advection), such that the flux divergence in Equation 10 is calculated from the vertical derivative and where q and q refer to the meltwater flux into (upper boundary) and out of (lower boundary) the layer. At the upper boundary, this is equal to the melt rate. Refreezing in Equation 11 is calculated within the subsurface thermal model, as a function of cold content in the layer.

Once a layer is temperate, liquid water can be retained within the pore space or it can percolate deeper into the snow or firn. Each layer can retain a certain amount of water within its pores due to capillary forces, the irreducible water content. We follow the formulation from Coléou and Lesaffre (1998), where the mass fraction of irreducible water, is a function of porosity, θ, following:

The porosity is calculated from the dry snow and ice densities in Equation 1, based on initial field measurements and the snow compaction model for melt‐season densification (Equation 4). For bulk snow/firn density ρ and water density ρ , the volume fraction of irreducible water, , is related to the mass fraction through:

We define the volumetric water content θ and the irreducible water θ based on the total snow or firn volume, rather than volume fraction within the pore space (e.g., Colbeck, 1974).

Water in excess of the irreducible water content is allowed to move downward. We calculate the water flux from Darcy's law, where ∇ϕ is the hydraulic gradient and the hydraulic conductivity k is a function of the porosity, θ, liquid water content, θ , and snow/firn grain size d . For unsaturated, gravitational percolation, the hydraulic gradient includes the combined effects of capillary and gravitational forces. Our calculation of hydraulic conductivity follows that of Meyer and Hewitt (2017), defining an effective permeability, k  = , where μ is the viscosity of the water and g is the gravitational acceleration. For water at 0°C, the viscosity is equal to 1.793 × 10−3 kg m−1 s−1. For unsaturated flow, the effective permeability is equal to the product of the intrinsic permeability, k (θ, d ), and the relative permeability, k (S): where S is the liquid water saturation (i.e., the fraction of the pore space filled with water, calculated from S = θ /θ) and P is the water pressure. Following Meyer and Hewitt (2017), we adopt a Carman‐Kozeny relation for the intrinsic permeability, k  = d 2/180·θ 3, with the relative permeability calculated from k  = S β, for exponent β = 2. We do not simulate the grain size or its evolution with temperature and depth, so treat this simplistically in this study by assuming a constant value d  = 0.001 m, based on typical grain diameters in Greenland firn (e.g., Linow et al., 2012). As an illustration with these values of β and d , a porosity of θ = 0.5 and a liquid water content of θ  = 0.02 give k  = 1.1 × 10−12 m2 and  = 6.1 × 10−6 m s−1.

In the unsaturated infiltration zone, the water pressure in Equation 15 is equal to the capillary tension, which accounts for the irreducible water content θ as parameterized in Equation 12. As this is implicitly accounted for, we directly model just the gravitational forces in the hydraulic potential, which allows the approximation q  ∼ k () in Equations 14 and 15. The TDR sensors provide a measure of the wetting front propagation speed of free water in the pore space. Following the Green and Ampt (1911) equation for water infiltration in soil, as recast by Mein and Larson (1973), the wetting‐front propagation speed v is related to water flux q through where θ is the initial liquid water content below the wetting front, which we assume to be 0 in the snow and firn. For wetting front depth z , v is equivalent to dz /dt and θ is the liquid water content of the advancing wetting front (i.e., the overlying layer). This allows us to relate measurements of the wetting front velocity with the Darcy flux and hydraulic conductivity.

Ice layers in the snow and firn can also impede meltwater infiltration. We parameterize this effect through an additional permeability factor , giving  =  .. We assign  = 1 where ice layers are less than 0.1 m thick, based on evidence of meltwater penetration through layers up to 0.12 m thick at DYE‐2 (Samimi et al., 2020). Thicker ice slabs act as impermeable barriers to meltwater infiltration (Gascon et al., 2013; MacFerrin et al., 2019), but their effective permeability depends on continuity and meltwater pathways through flaws or fractures (Humphrey et al., 2021). We assume that ice layers more than 0.5 m thick act as impermeable barriers (  = 0). For ice layers with thickness t between 0.1 and 0.5 m, we prescribe a nonlinear decrease in permeability,  = 10−γ ( − 0.1), with a reference value γ = 10 which reduces by four orders of magnitude as ice‐layer thickness increases from 0.1 to 0.5 m.

We solve the subsurface temperature and hydrological evolution in Equations 8 and 11 using 30‐min time steps through the summer melt season, permitting a direct comparison with the TDR probe and thermistor data at the different measurement depths. The simulations for summer 2016 use the observed spring snowpack as an initial condition and model densification following Equation 4, along with the effects of modeled ice layers that form from meltwater refreezing. For the multiyear simulations forced by GC‐Net and ERA5 meteorological fields, monthly snow accumulation at DYE‐2 is taken from the ERA5 precipitation. Layers in the firn model are added or removed as snow accumulates and ablates, effectively advecting accumulated snow and ice‐layer stratigraphy through the model grid.

Firn Model Initialization

Initial conditions are required for the model temperature, density, and ice‐layer stratigraphy. For the ERA5 simulations, we run the model through a 40‐year spin‐up over the period 1950–1990, then restart the simulation at 1950 with the spun‐up temperature, density, and ice‐layer stratigraphy as initial conditions. A 40‐year spin‐up is adopted because this is the approximate number of years represented in the upper 20 m of firn (Mosley‐Thompson et al., 2001). For numerical experiments forced by the GC‐Net AWS record, 1997–2018, we take the spun‐up firn conditions from the ERA5 simulations from 1950 to December 31, 1996, then apply the AWS forcing from the gap‐filled GC‐Net station record as of January 1, 1997.

For simulations of just the 2016 melt season to compare with the thermistor and TDR records, the model is initialized with the observed snow/firn thermistor temperatures, densities, and ice layers in the upper 3.5 m of the snow and firn (April 2016). Below this, densities and ice‐layer stratigraphy are based on measurements from the 21‐m firn core. Initial temperatures below 3.5 m are taken from the GC‐Net simulations for May 1, 2016.

Model Sensitivity Experiments and Evaluation

Model results are compared with the 2016 observational data set to examine one of the most uncertain process parameterizations in the firn model, the treatment of capillary water retention (Samimi et al., 2020; van As et al., 2016). Most previous firn modeling studies use either the parameterization of Coléou and Lesaffre (1998) for the irreducible water content or prescribe a fixed value such as θ  = 0.02 (e.g., Reijmer et al., 2012; Vandecrux, Mottram, et al., 2020; Verjans et al., 2019). Our reference model uses Equation 12 from Coléou and Lesaffre (1998), but we also carry out experiments with fixed values of the irreducible water content from 0.002 to 0.04 and with a simple scaling of Equation 12, from 50% to 150%. The latter approach captures the variation of capillary water retention with depth, in proportion with the empirical relation of Coléou and Lesaffre (1998).

Model results are evaluated against three different observations from the summer 2016 field studies: firn temperatures, T , maximum meltwater infiltration depth, z , and total summer melt, m. Firn temperatures are compared directly at all eight levels, and we evaluate model performance based on mean absolute error (MAE), mean error, and root‐mean‐square (rms) error of modeled versus observed temperatures. We also construct a composite model performance index, I, to evaluate the model results based on all three observed fields of interest, derived from the standardized model deviations (model minus observation) of firn temperature, ΔT , meltwater infiltration depth, Δz , and summer melt, Δm:

Mean absolute error is used as a measure of ΔT and σ is the mean standard deviation of observed firn temperature at all levels. An index of I = 1 indicates a perfect match with observations. I ∈ (0, 1) for all model experiments that we carried out, though this index can be negative if model deviations are large. Variables in Equation 17 can be weighted according to the relative importance of each field, but we assign them equal weights in our evaluation.

Results

AWS Observations and Surface Energy Balance

Figure 2 plots meteorological and surface energy balance conditions through summer 2016, spanning the melt season at DYE‐2. Snow surface temperature in Figure 2a is calculated from the measured outgoing longwave radiation, T  = (Q ↑/σE )1/4. There was a minor surface thaw while we were in the field in early May, with temperatures reaching 0°C, but the subsurface remained below the melting point through this period and there was no evidence of subsurface meltwater. A thin (∼1 mm) ice layer was observed at the snow surface. Two colder weeks followed, during which time the seasonal snow cooled to below −13°C. The snowpack remained dry throughout the month of May, while we were at the site. The melt season began the first week of June and ran through to late August (Figure 2a). Diurnal cycles of air and snow‐surface temperature indicate a daily freeze‐thaw cycle through most of the melt season, with surface refreezing at low solar radiation angles (“overnight”). Air and snow‐surface temperatures reach 0°C for a few hours on most days through the summer (JJA) melt period, with two occasions when melt conditions persisted through the night (July 18–19 and August 9–11).

Figure 2

Meteorological data and modeled surface energy balance at DYE‐2 over the summer melt season, May 24 through September 18, 2016. (a) Average 30‐min values of air and snow‐surface temperature, T and T , daily average air temperature, T ad, and daily maximum snow‐surface temperature, T . (b) Ultrasonic depth gauge (UDG) record of snow‐surface height. (c) Mean daily net radiation, Q*, turbulent flues, Q (sensible plus latent heat fluxes), and net energy, Q . (d) Daily total melt, net melt (ablation), and refreezing. Net melt refers to the surface melt that infiltrates below the surface layer and total melt includes refrozen near‐surface meltwater (i.e., “recycled” mass).

Net longwave radiation was consistently negative save for two overcast periods on the dates noted above, when snow and air temperatures remained near 0°C. This drove high values of incoming longwave radiation and net energy (Figure 2c), giving the highest melt rates of the summer (Figure 2d). Net radiation was positive over the summer (JJA), averaging 14 W m−2, while sensible, latent, and conductive heat fluxes averaged +3, −11, and +1 W m−2, respectively (Table 2). Net energy averaged 5 W m−2 for the summer and reached a maximum of 10 W m−2 in July. It declined to 3 W m−2 in August, due to the reduced shortwave radiation and atmospheric cooling, but early August meltwater production was nonetheless significant, exceeding that of June and representing ∼32% of the summer total. Conductive heat flux to the surface was more important in August and September, averaging ∼3 W m−2.

Table 2

Mean Monthly and Summer June Through August (JJA) Surface Energy Fluxes and Melt Totals, Summer 2016

Period	Q S ↓	Q S ↑	Α	Q L ↓	Q L ↑	Q H	Q E	Q C	Q N	T a	PDD	Melt
May	303	252	0.83	211	266	9	−9	0.1	0	−11.6	0.6	0.02
June	330	272	0.82	255	298	2	−14	−0.3	4	−4.0	1.8	0.09
July	341	256	0.75	227	296	6	−11	0.2	9	−4.1	2.7	0.17
Aug	252	189	0.75	237	289	2	−9	1.4	5	−5.7	2.4	0.12
Sept	143	117	0.82	222	250	2	−6	1.4	−1	−14.5	0.0	0.00
JJA	308	239	0.78	239	295	4	−12	0.4	6	−4.6	6.9	0.38

Note. All energy fluxes (Q) have units W m−2, α is the surface albedo, T is the mean air temperature (°C), PDD denotes positive degree days (°C d), and melt reports monthly and summer totals (m w.e).

UDG snow surface height measurements at the two firn pits indicate a total summer melting of ∼0.5 m of snow or 0.18 ± 0.04 m w.e. (Figure 2b). The error estimate includes uncertainties about the snow density and potential sinking of the UDG pole, which would cause the ablation to be underestimated. Summer snow accumulation of ∼0.25 m partially offset the melting, resulting in ∼0.25 m of ablation relative to the April snow surface (Figure 2b). The upper two sensors at site A, at depths of 0.1 and 0.2 m, melted out on July 25 and July 30, respectively. Data from these sensors is omitted from the analysis after these dates. Sensors at and below 0.3 m depth remained below the surface through the full melt season. The surface energy balance model gives an estimated total melt of 0.40 m w.e. from May through September (Figure 2d). About half of the melt occurred in July, when net radiation was the strongest. The total includes melting of refrozen near‐surface ice lenses and pore water, which we call “recycled” melt. The ablation or “net melt” is much less than the total melt, equal to 0.16 m w.e. with the reference model. This value depends on the parameter settings for capillary water retention (see Section 3.4).

All of the meltwater refroze in the snow and firn. Many days had strong diurnal melt‐freeze cycles, with identical daily melt and refreezing totals (Figure 2d). At other times, refreezing lagged melting by a few days, reflecting meltwater that infiltrated isothermal snow or firn and remained as liquid water for a short while. Table 3 lists the monthly and summer partitioning between melt, ablation, refreezing, and the different energy fluxes within the subsurface model. Of 234 MJ m−2 in total positive net energy over the summer, 100 MJ m−2 (43%) went to warming the near‐surface snowpack and 134 MJ m−2 (57%) was directed to melting. The latent heat released in the near‐surface snow and firn during refreezing was identical to the melt energy, 134 MJ m−2, as 100% of meltwater refroze.

Table 3

Monthly, Summer, and Melt‐Season (MJJAS) Mass and Energy Fluxes for the Reference Model Parameter Settings

Period	Mass balance (m w.e.)			Energy fluxes (MJ m−2)
	m	m s	r	E warm	E melt	E adv	E ref
May	0.02	0.01	0.02	30	7	0.0	7
June	0.09	0.03	0.09	18	31	0.1	31
July	0.17	0.08	0.14	23	55	0.3	47
Aug	0.12	0.05	0.14	23	41	0.2	47
Sept	0.00	0.00	0.01	24	2	0.1	7
JJA	0.38	0.15	0.37	64	127	0.6	125
MJJAS	0.40	0.16	0.40	117	135	0.6	135

Note. m, m , and r refer to the total monthly or summer melt, the net melt (surface ablation associated with melting, including the effects of “Recycled” meltwater), and refreezing. E warm is the total surface energy used to warm the surface layer of the snow (when T  < 0 and Q  > 0), E melt is the energy used for melting (when T  = 0 and Q  > 0), E adv is the sensible heat advection associated with meltwater infiltration, and Eref is the latent heat energy released by meltwater refreezing.

TDR and Thermistor Observations

Figure 3 plots measured subsurface temperatures and liquid water content in the two firn pits through summer 2016. The melt‐season progression was similar at the two sites, with nearly identical timing and depth for the meltwater infiltration and subsurface warming events. There were three main periods of melting in summer 2016, indicated by the grey shading. The first major melt event was in the second week of June, with meltwater infiltration to a depth of 0.2 m, as recorded by the two upper two sensors at site B (Figures 3e and 3f). The snow at 0.4 m depth remained below 0°C and dry during this initial event. This initial meltwater infiltration episode was also undetected by the upper sensors at site A, at 0.3 m (Figure 3c), so appears to have been confined to the upper 0.2 m of the snowpack. Sensors at 0.3 and 0.4 m depth warmed to near 0°C during this event, but an atmospheric cooling cycle followed and the entire snowpack refroze and cooled.

Figure 3

Observed (a) air and snow surface temperatures, (b, e) snow temperatures, and (c, f) liquid water content, along with (d) modeled melt and refreezing, summer 2016. (b and c) are for site A and (e and f) are for site B. Shaded areas indicate periods of meltwater infiltration, as indicated by the subsurface temperature and liquid water content measured at site B. Liquid water content data are smoothed using an 11‐point (5‐h) moving average.

Melting resumed on June 22, accompanied by an abrupt warming by up to 5°C in the upper ∼1.4 m of snow and firn and wetting to ∼0.5 m depth at both sites (Figures 3c–3f). The warming was recorded over less than 24 h, evidence of latent heat release from refreezing meltwater, followed by slower conductive warming at firn depths below 2 m. A detailed analysis of the meltwater infiltration event on June 22–23 provides a quantitative illustration of the coupled thermal and hydrological evolution during this abrupt warming and meltwater infiltration event (Figures 4a and 4b). Average air temperature was −1.1°C over this 48‐h period, with a maximum temperature of +1.2°C. Total modeled melt on these two days equals 13.6 mm  w.e., with the onset of melt at 11:00 on June 22. By that evening, wet, temperate conditions developed to 0.6 m depth at site A and 0.4 m depth at site B. Temperatures at 11:00 on June 22 at 0.1, 0.2, 0.4, and 0.6‐m depth were −2.6°C, −2.7°C, −3.5°C, and −5.2°C, respectively. Snow rapidly warmed at these four depths, reaching 0°C between 18:00 and 19:00 in the upper 0.6 m at site A and between 18:00 and 20:00 in the upper 0.2 m at site B. Sensors at 0.4 m depth at site B first recorded temperate, wet conditions at 08:00 on June 23. All sensors in the upper 0.5 m of snow then maintained temperate conditions for the next several days at the two sites.

Figure 4

Observed snow/firn liquid water content (top of each panel) and temperature (bottom of each panel) for select periods of interest in the summer 2016 melt season. (a and b) Melt onset at sites A and B for the period June 21–July 4, for all sensors at depths less than 1 m. (c) Records from 0.6 to 1.4‐m depth at site A from August 6–28, indicating the arrival of the thaw front (0°C isotherm) at 1.8‐m depth on August 11 and the subsequent refreezing and cooling as the melt season shut down. (d) Records from 0.6 to 1.6‐m depth at site B during the melt season shutdown, August 13–25.

Dielectric permittivity (liquid water content) increased from background values at the same time as snow became temperate at each site during this initial melt event, with a simultaneous arrival of the thaw and wetting fronts at each depth. This melt cycle ended on July 2 at each site, as atmospheric cooling caused meltwater refreezing followed by cooling of snow temperatures. The thermistor and TDR measurements were highly consistent through this melt event, with θ  = 1–2% under melting conditions and θ  ∼ 0, when temperatures were below 0°C.

Following a two‐week cool hiatus with limited surface melting, the melt season resumed in earnest on July 18 (Figure 3). The snow and firn had returned to dry, subzero conditions, so the next melt episode was again accompanied by abrupt warming. Melting conditions persisted for ∼one month, with temperate conditions developing to 1.8 m depth by August 10. All of the sensors at site B registered liquid water and melting temperatures (Figures 3e and 4d). At site A, the upper three sensors recorded wet and temperate conditions through most of this record, through to August 23, while the firn at 1.8 m depth reached the melting point for only a short interval (four hours on August 10), with a small, simultaneous increase in dielectric permittivity. Sensors below 1.8 m depth at site A recorded sub‐zero temperature with no evidence of meltwater. The thaw and wetting fronts therefore briefly reached a depth of 1.8 m. Volumetric water contents of 3%–4% were common in the upper 0.6 m during the main melt July‐August episode, with short‐lived values of up to 5% during the early stages of the main meltwater infiltration events. Water content in wetted firn was typically 1%–2% at depths greater than 0.6 m.

Declining shortwave radiation and cooler air temperatures in the second week of August triggered the end of the melt season, although it took several days for the firn to refreeze and cool. Figures 4c and 4d plot the end of the melt season at mid‐depths (0.6–1.8 m) at each of the sites. During the last stages of refreezing, an isothermal, wet layer of snow and firn from 0.6 to 1.1 m depth was sandwiched between colder snow above and below. Cooling fronts moved down through the snowpack from the surface and also upward from the cold, underlying firn. The step‐wise drying and refreezing at different depths is particularly clear in pit B in Figure 4d, with the deeper firn being the first to refreeze (at 1.6 m and then 1.4 m), followed by refreezing at 0.6, 0.9, and 1.1 m. It took ∼one week for the upper 1.8 m of snow and firn to refreeze and start cooling below 0°C, with the firn at ∼1 m depth being the last portion of the seasonal thaw layer to refreeze. Small melt events after August 21 (Figure 3) appear to have refrozen near the surface, with no indication of subsurface liquid meltwater at either site after this time.

Overall, four abrupt, step‐like warming events were registered at each site, associated with latent heat release from meltwater infiltration to progressively greater depths. The seasonal firn temperature evolution was driven by these abrupt events, along with a more gradual, diffusive response to the atmospheric forcing and the meltwater infiltration/latent heat episodes (Figures 3b and 3e). Firn warming from 2.1 to 3.7 m depth on ∼August 12 (Figure 3b) is interesting in that we see no evidence of meltwater infiltration below 1.8 m depth, but there are relatively sharp temperature responses below this (e.g., a warming of 4.2°C at 2.7 m over ∼4 days from August 8 to 12). This likely reflects thermal diffusion, given the timescale, but it is also possible that deep meltwater infiltration and refreezing occurred through preferential pathways in the vicinity of our sensors, driving local warming that was recorded by the thermistors.

Firn Modeling, 2016 Melt Season

Figure 5 plots the modeled subsurface temperature evolution in the upper 4 m of the snow and firn for the reference model parameter settings. Observed air and subsurface temperatures at site A are also plotted in Figure 5, for comparison. Diurnal freeze‐thaw cycles at the surface are well represented in the model, and are generally restricted to the upper 0.2 m of the snow. The abrupt subsurface warming events during the main pulses of melt are captured by the model (Figures 5b and 5c), and are part of a stepwise penetration of the thawing and wetting front through the summer melt season. The model also reproduces the isolated warm layer that was observed at 0.6–1.2 m depth during the August freeze‐up. However, temperate conditions reached a maximum depth of 1.2 m in the second week of August in the model, which is an underestimate relative to the observed thaw front penetration to 1.8 m. Modeled temperatures are highly correlated with the observations, with an average melt‐season linear correlation coefficient of 0.96 over the eight observations depths from site A (Table 4).

Figure 5

Results for the reference model parameter settings, with meltwater retention after Coléou and Lesaffre (1998). (a) Modeled melt and refreezing, (b) observed subsurface temperatures at site B, for comparison, and (c) modeled subsurface temperature evolution in the upper 4 m of snow and firn.

Table 4

Model Comparison With Observations for Different Treatments of Irreducible Water Content, θ

θ wi model	m (m w.e.)	z thaw (m)	r T	ME (°C)	MAE (°C)	RMSE (°C)	I
Reference	0.16	1.2	0.96	−1.07	1.44	2.01	0.84
CLp8	0.19	1.4	0.96	−0.82	1.27	1.78	0.88
CLp7	0.21	1.6	0.96	−0.67	1.16	1.63	0.89
CLp6	0.23	1.8	0.97	−0.38	0.95	1.32	0.89
CLp5	0.25	2.0	0.98	−0.05	0.77	1.06	0.85
0.04	0.16	1.0	0.96	−1.28	1.60	2.17	0.79
0.03	0.20	1.2	0.96	−1.00	1.51	1.91	0.85
0.02	0.25	1.8	0.97	−0.40	0.97	1.34	0.85
0.01	0.32	2.4	0.95	0.92	1.46	1.90	0.69

Note. The reference model uses the θwi parameterization of Coléou and Lesaffre (1998). Experiments CLp8, CLp7, etc. scale the reference values by 0.8, 0.7, etc., and numbers in column 1 refer to constant values prescribed for θwi. The observed summer melt and thaw depth were m = 0.18 m w.e. and z thaw = 1.8 m. ME, MAE, and RMSE are the mean, mean absolute, and mean rms errors in the modeled temperatures, averaged over the summer for all observation depths, rT is the linear correlation coefficient between the observed and modeled temperatures, and I is the composite skill index for the model. Bolded values indicate the optimal model performance for a given measurement.

The mean bias in the modeled temperature is −1.1°C, with a mean absolute error of 1.4°C. Modeled net summer melt equals 0.16 m w.e. with the reference parameter settings. This is within the uncertainties but ∼10% less than the estimated summer melt based on the ultrasonic depth sensors, 0.18 ± 0.04 m w.e. Both subsurface temperatures and meltwater infiltration are underestimated compared with the observations. We do not simulate preferential flow paths that can lead to enhanced infiltration (e.g., Wever et al., 2014), so meltwater may be retained or refrozen too near the surface, contributing to the cold bias in the firn temperatures. Modeled meltwater infiltration is also highly sensitive to the parameterization of capillary water retention, as explored in the next section.

Sensitivity of the Modeled Meltwater Infiltration to Capillary Water Retention

Figure 6 plots modeled summer melt and the maximum depths of the thaw and wetting fronts as a function of θ . Greater near‐surface meltwater retention (higher θ ) limits meltwater infiltration and also reduces net melting, since a higher proportion of meltwater is recycled in near‐surface freeze‐thaw cycles. The blue shading in Figure 6 indicates the observed values and uncertainty range for meltwater infiltration depth and summer melt. The thaw front reached 1.8 m in our TDR data, but Heilig et al. (2018) estimate a maximum meltwater infiltration depth of 2.2 m, so we accommodate that possibility in the uncertainties. Averaged over the snow and firn column, the reference model parameterization from Coléou and Lesaffre (1998) gives a value θ  = 0.0273 (2.73%). Results are plotted against the mean value in Figure 6, but θ varies with depth, following Equation 11, with higher values near the surface. For a given mean value of θ , melt totals with irreducible water content based on Equation 11 are systematically lower than for uniform θ (Figure 6), and are in better accord with the observed melt.

Figure 6

Sensitivity of modeled (a) maximum summer thaw (z ) and meltwater infiltration (z ) depths and (b) net surface melt to the parameterization of capillary water retention. The observed meltwater infiltration depth, 1.8–2.2 m, and net summer melt, 0.18 ± 0.04 mm w.e., are indicated with the blue lines, with shading indicating the estimated confidence range. Purple symbols are for the mean (depth‐averaged) values of θ from the parameterization of Coléou and Lesaffre (1998) and blue symbols are for uniform values of irreducible water content θ  ∈ [0.01, 0.04]. From left to right, the Coléou and Lesaffre (1998) cases that are plotted are CLp5, CLp6, CLp7, CLp8, CLp9, and the reference model parameterizaton (CL) (see Table 4).

Table 4 summarizes model results for a selection of the experiments in Figure 6. No treatment of the capillary water retention simultaneously optimizes the fit to the different observations. For instance, a scaling factor of 0.5 applied to the parameterization of Coléou and Lesaffre (1998), which we denote CLp5, gives the best fit to the observed subsurface temperatures, but this comes at the expense of a large overestimation of the net melt (0.25 m w.e.) and meltwater infiltration depth (2.4 m). We evaluate the model based on the composite index I, which includes the agreement with observed firn temperatures, maximum summer meltwater infiltration depth, and net summer melt. The optimal value of capillary water retention is lower than that of our reference model, through a scaling factor of 0.7 for the parameterization of Coléou and Lesaffre (1998), which we denote CLp7. This gives the best combined fit to the three observational criteria, I = 0.89. We adopt this as the optimized model for decadal‐scale simulations. The average irreducible water content with this setting is 0.019. For uniform values of the irreducible water content, the optimal results were for θ  = 0.025, but this setting gives both an over‐estimate of melt and an underestimate of meltwater infiltration.

Figure 7 plots modeled versus observed mean daily subsurface temperatures at five different depths for the CLp7 model configuration. This has a higher correlation of temperatures, lower temperature error, and greater meltwater infiltration depth than the reference model (Table 4), in better accord with the observations, but net summer melt is 0.21 m w.e., 0.03 m greater than observed. The model is slightly too cold below 0.6 m depth; it is possible that the model thermal conductivity is too low, such that near‐surface warming from latent heat release does not adequately conduct to depth. Consistent with this possibility, observed abrupt warming events from 2 to 4‐m depth are more diffuse in the model (Figures 7c and 7d). These sudden warmings, with a magnitude of ∼5°C, are well reproduced in the near‐surface firn (Figures 7a and 7b).

Figure 7

Observed (thick lines, site A) versus modeled (thin lines) mean daily snow/firn temperature through the summer melt season at depths of (a) 0.3 m, (b) 0.6 m, (c) 1.4 m, and (d) 2.1 and 3.7 m. Results are from a simulation with the optimized model parameters (see the text).

Multidecadal Firn Evolution

Multi‐year simulations with the CLp7 parameter settings provide insight into the interannual variability of meltwater infiltration depths and longer‐term trends in firn temperature, density, and ice content at DYE‐2. Over the period of available AWS data at the DYE‐2 GC‐Net station, 1997–2018, mean summer (JJA) temperature (±1σ) was −5.4 ± 1.1°C, with no statistically significant trend. This result is echoed in the modeled summer melt, with no significant trend and a mean value of 0.14 ± 0.13 mm w.e. from 1997 to 2018 (Table 5). The average modeled depths of the summer meltwater infiltration and thaw fronts are 1.7 ± 1.4 m and 1.4 ± 1.4 m, respectively, with temperate conditions and meltwater infiltration confined to the upper 2 m of snow and firn in most years. Interannual variability in summer melt was dominated by the exceptional summer of 2012 over this period, with a modeled melt of 0.64 m w.e., ∼4σ above normal. The summer of 2012 was also exceptional with respect to meltwater infiltration, with the modeled 2012 wetting and thaw fronts reaching depths of 7.2 and 6.8 m, respectively.

Table 5

Summer (June Through August [JJA]) Climate, Surface Energy Balance, and Firn Conditions at DYE‐2 From 1950 to 2020, Based on Firn Modeling Forced by the ERA5 Meteorological Reanalyses

Period	Q S ↓	Q*	Q H	Q E	Q N	T a	T JJA	PDD	Melt	Net melt	T 10	z t	z w
1950–2020E	307	15	4	−13	6	−18.0	−5.4	18.0	0.26	0.16	−16.7	1.6	1.9
1997–2018E	306	16	3	−14	5	−17.5	−5.4	14.3	0.24	0.14	−16.6	1.5	1.7
1997–2018G	306	16	2	−14	5	−17.6	−5.4	11.0	0.23	0.14	−16.3	1.4	1.7
2012E	304	26	5	−11	21	−16.5	−4.2	45.4	0.66	0.54	−13.8	6.0	6.4
2012G	298	30	7	−14	24	−15.6	−2.9	39.8	0.79	0.64	−12.8	6.8	7.2

Note. A subset of the results is shown for the period 1997–2018, for modeling driven by ERA5[E] and GC‐Net[G] AWS data, as well as the extreme 2012 melt year. Melt refers to the total annual surface melting and “net melt” refers to the net surface melting minus refreezing, accounting for meltwater freeze‐thaw cycles. This is the actual surface drawdown associated with summer melting. All energy fluxes have units W m−2, temperatures are in °C, PDD has units °C d, melt is in m w.e., and the maximum summer thaw and meltwater infiltration depths, zt and zw, are in m.

Figure 8 plots the modeled subsurface temperature evolution over this period of GC‐Net data. Mean annual air temperature also had no trend at DYE‐2 from 1997 to 2018, but 10‐m and 20‐m firn temperatures increased over this period, with statistically significant trends of +0.7 and +0.8°C per decade, respectively. The mean air temperature at the GC‐Net station is −17.6 ± 1.2°C, while the firn model predicts mean annual temperatures of −16.3 ± 1.0°C and −16.3 ± 0.8°C at 10‐ and 20‐m depths, with latent heat release from meltwater refreezing causing an average annual warming of 1.3°C over the period 1997–2018. Changes at 20 m lag those at 10 m by ∼one year. The exceptional 2012 melt season caused a warming of ∼4°C at 10‐m depth over the period 2012–2015 (Figure 8b), but the firn recovered from this by 2018, with 10‐m temperatures returning to −16.1°C. Figure 8c plots the envelope of firn temperatures as a function of depth over the 22‐year simulation, calculated from the mean daily temperatures. The maximum envelope illustrates the temperate firn development to 6.8 m depth in 2012, while the red dashed line, from 2016, represents a more typical summer.

Figure 8

(a) Subsurface temperature evolution in the 20‐m firn layer, 1997–2018. (b) Close‐up of the upper 12 m of firn from 2010 to 2018, highlighting the anomalous 2012 melt season. (c) Mean, minimum, and maximum mean daily temperatures in the upper 12 m of snow and firn. Solid lines are for the entire period, 1997–2018, and dashed values are for calendar year 2016, which is a typical year for this period. The T max line in (c) corresponds to 2012.

Measurements from the firn cores provide some validation of the modeled firn density and ice content. Examining the modeled 20‐m firn profile in 2016, the year of firn‐core acquisition, the total ice content and mean firn density and in the model are 6.1 m and 693 kg m−3. These are within 5% of the measured values in the 20.3‐m firn core, a total ice content of 6.4 m and a mean density of 665 kg m−3. These values agree within the uncertainties of the measurements and the modeling, but background firn density (i.e., firn densification from Equation 4) may be overestimated in the model, since bulk density is too high even with a 5% underestimate in ice content.

ERA5 climate reanalyses provide a longer‐term perspective on firn evolution at DYE‐2.

Modeled air temperatures, firn temperatures, and melt conditions from 1950 to 2020 are plotted in Figures 9 and 10. The average summer net melt over this 71‐year period is 0.16 ± 0.12 m w.e., with average meltwater infiltration and thaw depths of 1.9 ± 1.3 m and 1.6 ± 1.3 m, respectively (Table 5). Average GC‐Net and ERA5 results are similar over the common period in the simulations, 1997–2018 (Table 5), but there are some interannual differences. The 2012 extreme melt event was also extreme in the ERA5‐forced model results, but was slightly muted relative to the GC‐Net‐driven simulations. Mean summer temperature in summer 2012 is lower in the ERA5 forcing relative to the GC‐Net data, −4.2°C versus −2.9°C (Table 5). This contributes to lower modeled melt totals, but the ERA‐driven meltwater infiltration depth in 2012 still extends to 6 m depth, warming 10‐m firn temperatures by 2.6°C.

Figure 9

Modeled firn evolution at DYE‐2, 1950–2020, forced by the ERA5 climate reanalyses. (a) Summer air (T JJA) and snow‐surface temperatures (T sJJA). (b) Net annual surface melt. (c) Mean annual air (T ann) and firn temperatures (T 2, T 10, T 20) at 2‐, 10‐ and 20‐m depths. (d) Maximum depth of the annual wetting and thaw fronts.

Figure 10

(a) Subsurface temperature evolution driven by ERA5 climate reanalyses in the 20‐m firn layer, 1950–2020. (b) Close‐up of the upper 12 m of firn from 1965 to 1985 and 2000–2020, highlighting the firn warming effects of meltwater infiltration to 6–7 m depth in 1968, 1978, 2012, and 2019.

There is a statistically significant trend in mean annual air temperature over the ERA5 period, +0.14°C per decade, but mean summer air temperature, firn temperatures, annual melt, and meltwater infiltration depth have no trend over the full period, 1950–2020 (Figure 9). In contrast, there are significant positive trends over the more recent period 1990–2020. Summer (JJA) and mean annual air temperatures have trends of +0.19 and +0.50°C per decade over this period, net melt increased by +0.03 m w.e. per decade, and meltwater infiltration depth increased by +0.33 m per decade. These combined influences produce a statistically significant trend of +1.0°C per decade in 10‐m firn temperatures, that is, a 3°C increase over the last three decades. The firn warming is partially driven by the extreme melt years of 2012 and 2019, with modeled surface melting of more than 0.5 m w.e. and meltwater infiltration to greater than 6‐m depth in each of these summers (Figures 9b and 9d).

The deep meltwater infiltration events of 2012 and 2019 were exceptional but are not isolated or unprecedented in the ERA5‐driven simulations. The summers of 1968 and 1978 featured similar summer temperatures, melt totals, and firn warming impacts (Figures 9 and 10). The 1968 melt event resulted in a 10‐m firn warming of ∼4°C, exceeding the firn warming from 2012 (Figure 9c). The summer of 1968 was the warmest of the ERA5 period, with a mean JJA temperature of −3.1°C and modeled meltwater infiltration to a depth of 6.8 m. The large impact on firn temperatures may have been preconditioned by a strong melt season with meltwater infiltration to 4.8 m depth in 1966. The 2019 melt season was similar to that of 1968, with a predicted 10‐m firn temperature anomaly of ∼4°C that can be expected to persist for several years. Within the model, the amount of seasonal meltwater is the primary driver of interannual variability in the depths of summer thaw and wetting fronts (Figures 9 and 10). The linear correlation coefficient between net summer melt and meltwater infiltration depth is r = 0.97. Winter snow accumulation and temperature (i.e., snow and firn cold content) have some influence, but at DYE‐2 there is still ample cold content and pore space to refreeze all of the summer meltwater, so seasonal meltwater infiltration depth is primarily a function of the amount of meltwater that is available as a latent heat source.

Discussion

Assessment of TDR Methods in Polar Firn

TDR Precision and Accuracy

Prior to the onset of the melt season and over a time period of weeks, changes in density should be minimal and dielectric permittivity should be constant for dry snow and firn. This was true for measurements during the month of May at the study site, with no trends but considerable “white noise” in the dielectric permittivity recorded by the TDR sensors. The eight probes at site B had a mean value and standard deviation of ε  = 2.37 ± 0.05 for the month of May. From Equation 3, δε  = ±0.05 equates to a liquid water content of δθ  = ±0.2%. These dry‐snow data provide an estimate of the TDR precision in snow and firn; taking a 2σ envelope for the instrumental noise, we infer that the changes in liquid water content of less than 0.4% fall within the instrumental noise and cannot reliably be detected. We infer typical liquid water contents of 1%–2% under temperate conditions at DYE‐2, with maximum values reaching ∼5% (Figures 3 and 4), so the signal to noise ratio is adequate to detect wet snow and track meltwater pulses with depth. Low liquid‐water contents or trace melt events are more challenging to capture with the TDR probes.

As an explicit example of the expected liquid water content at DYE‐2, typical July melt rates at DYE‐2 are 10 mm w.e. d−1 (Figure 2d), with peak hourly values of ∼1 mm w.e. hr−1. Over an hour, this equates to 1% by volume in the 0.1‐m surface layer, so is detectable by the TDR sensors but not with high accuracy with respect to the volume of water. Once water content exceeds the capillary water retention capacity, water percolates vertically with high efficiency in temperate snow; the TDR records indicate that surface melt can infiltrate to depths of 0.5 m or more within the same day (Figures 3c and 3f). A total of 10 mm w.e. of daily melt distributed over the upper 0.5 m of snow is equivalent to 2% by volume. This is consistent with the inferred liquid water contents of at DYE‐2, and supports the promise for TDR methods in tracking meltwater infiltration in polar snow and firn.

Comparison With GPR Measurements

Continuous upward‐looking GPR measurements from DYE‐2 provide an independent estimate of meltwater infiltration during the 2016 melt season (Heilig et al., 2018). The timing and pattern of meltwater evolution in the radar data reflects a similar stepwise evolution of the thawing and wetting fronts through the melt season, with a maximum meltwater infiltration depth on August 10 in both the TDR and radar records. The radar data indicate deeper infiltration than the TDR‐based inferences, to 2.2 versus 1.8 m below the surface. Heilig et al. (2018) also infer maximum liquid water contents of ∼2%, whereas the upper TDR sensors recorded values of up to 5% for short periods. The differences in meltwater infiltration depth and liquid water content are significant, that is, exceeding the uncertainty bounds within each method.

Differences between the two methods may be related to lateral variability in meltwater infiltration. Vertical sensor arrays cannot distinguish between the arrival of a uniform wetting front versus preferential flow pathways, that is, meltwater piping. The two firn pits instrumented with TDR sensors are consistent in the inference of wetting and thaw fronts to ∼1.8 m depth, but the radar system may have detected a meltwater pipe that infiltrated to a greater depth. Given its burial depth of 4.3 m, the upward‐looking GPR system of Heilig et al. (2018) monitored a cone with a horizontal footprint of ∼70 m2 at the surface. In contrast, the intrinsic horizontal resolution of the TDR measurements is the area directly sampled by the TDR probes, 0.03 m2. In sweeping over a larger volume of snow and firn, the radar data integrates meter‐scale variability in the snow water content. If there is significant lateral variability, averaging of water content over a given depth range in the radar data would give lower values for maximum water content, relative to the TDR measurements.

Overall, the two TDR arrays and the radar data give similar results from three discrete points within a 400‐m region at DYE‐2. At each location, meltwater infiltration followed a similar evolution and reached a maximum depth of ∼2 m in summer 2016. The question of spatial variability is an important one for meltwater infiltration and refreezing processes in firn (van As et al., 2016), and it is important to understand the extent and significance of lateral variability in meltwater infiltration for the evolution of firn temperature, stratigraphy, and density. This cannot be evaluated from our experimental design, but it would be possible to examine this though a more extensive sensor network that deploys horizontal as well as vertical sensor arrays. Such an approach is recommended for experiments to quantify local‐ to regional‐scale lateral variability.

Dielectric Mixing Models

Mixing models to relate bulk dielectric permittivity to liquid water content in snow and firn are difficult to rigorously validate. Unlike in soils, where it is straightforward to evaporate off the liquid water and compare wet and dry masses, liquid water content in snow or firn is difficult to measure independently. The thermistor and TDR records during melt‐freeze cycles provide some control on the TDR meltwater volume estimates. As an example, the abrupt subsurface warming to 0.5 m depth on June 22 and 23 (see Section 3.2) was driven by an estimated 13.6 mm w.e. of melt. Refreezing of 13.6 mm w.e. of meltwater would release 4.6 × 106 J of latent heat. The upper 0.5 m of snow had a mean density of ∼340 kg m−3 and an average initial temperature of −4.2°C on June 22. Using a specific heat capacity of 2,090 J kg−1°C−1, 1.5 × 106 J are needed to bring the upper 0.5 m of snow to the melting point, so the available meltwater/latent heat release readily accounts for this abrupt warming event. The cold content in the upper 0.5 of the snowpack is sufficient to refreeze ∼33% of the available melt (4.5 mm w.e.), leaving 9.1 mm w.e. to be stored as liquid water, which equates to an average liquid water content of 1.8% in the upper 0.5 m. The measured values of liquid water content in the upper 0.5 m averaged 1.7% for the period 08:00–20:00 on June 23 (Figures 4a and 4b), which is equivalent to the thermodynamic estimate, within the uncertainities of the observations. We note that this is an open system, however, with conductive, sensible heat, and outgoing longwave radiation fluxes acting to transfer energy from the near‐surface snowpack to the atmosphere through much of this period. This energy sink may have enabled more than 4.5 mm w.e. of meltwater to refreeze.

Different empirical equations have been used to estimate snow water content using dielectric and TDR methods (e.g., Denoth, 1994; Denoth et al., 1984; Schneebeli et al., 1997; Techel & Pielmeier, 2011). In lieu of empirical relations, mixing models in the general form of Equation 2 have also been applied with different exponents in firn radar studies. Heilig et al. (2015, 2018) use an exponent β = 0.5 in a three‐component mixing model, as described in detail in Schmid et al. (2014). Alternatively, numerous radar studies have adopted a mixing model after Looyenga (1965) with an exponent of 1/3 rather than 1/2 (e.g., Christianson et al., 2015; Macheret et al., 1993; Murray et al., 2000; Van Pelt et al., 2014). Samimi and Marshall (2017) found Equation 3 with β = 0.5 to work well with TDR sensors in snow, based on validation against a Denoth (1994) capacitance plate and modeling of snow melt (hence, expected liquid water content). We therefore limit our experiments to this mixing model, which also permits a more direct comparison with the results of Heilig et al. (2018).

The optimal mixing model may depend on the measurement technique, as instruments sample different length scales and orientations of the media. TDR probes sample a relatively small radius of influence, ca. 5 cm, in a two‐dimensional sense, parallel to the surface and therefore within any horizontal layering. Effects of layering and anisotropic crystal orientation will therefore be minimal relative to the three‐dimensional sounding that is intrinsic to radar wave propagation. Bayesian mixing models, which include estimates of uncertainty, could be an interesting avenue to explore in future work to examine the optimal mixing model.

The empirical expression of Denoth (1994) relates snow dielectric permittivity to bulk density, ρ , and liquid water content following ε  = 1 + a 1 ρ  + a 2 ρ 2 + a 3 θ  + a 4 θ 2, where a 1 = 0.00192, a 2 = 4.4 × 10−7, a 3 = 18.7, a 4 = 45, and ρ has units kg m−3. Estimates of liquid water content from Equation 3 are less sensitive to changes in bulk dielectric permittivity than the equation of Denoth (1994). For reference values of θ  = 0.02 (2%) and ρ  = 500 kg m−3, the Denoth (1994) equation gives ∂θ /∂ε  = 0.049 (4.9%) and the mixing model in our study gives ∂θ /∂ε  = 0.040 (4.0%). For a given measured increase in ε , estimates of θ from Equation 3 are therefore ∼20% less than would be inferred from the Denoth (1994) equations. Given this conservative estimate, the higher peak values of liquid water content in our data, compared with Heilig et al. (2018), are not likely to be an artefact of the mixing model.

Effects of Changes in Snow and Firn Density

The dielectric permittivity of snow and firn increases with density as well as water content (Denoth et al., 1984; Lundberg, 1997; Schneebeli et al., 1997). An increase in snow density of 100 kg m−3 is equivalent to a liquid water content of ∼1% (Lundberg, 1997), so densification over the melt season may contribute to some of the bulk dielectric permittivity signal that we interpret as liquid water. We can assess the magnitude of this effect for dry firn densification for the two deepest TDR probes at site A. Firn at these depths remained frozen and unaffected by meltwater through the full 2016 melt season. These sensors recorded a gradual increase in dielectric permittivity over the melt season. Averaged over the period September 1–16, the end‐of‐summer values of ε increased above their May baseline values by Δε  = +0.07 and Δε  = +0.03 at 2.8 and 3.7 m depth, respectively. In the mixing model for liquid water content (Equation 3), this is equivalent to apparent liquid water contents of 0.29% and 0.12%, which is less than the instrumental noise (0.4%).

Greater densification is expected near the surface in association with seasonal snow settling and refrozen meltwater. We partially account for these densification effects by resetting the reference dielectric permittivity for each TDR sensor after major refreezing cycles in the snow and firn, that is, when θ drops to ∼0 and temperatures fall below 0°C. Refrozen ice layers will increase the bulk density, but the baseline values of ε are renormalized for these effects. This restoration also eliminates the effects of potential changes in probe coupling with the snow and firn (e.g., air gaps) which may result from meltwater effects on the probes (e.g., Lundberg, 1997). Additional studies are needed to better understand the potential impacts of ice layer formation and freeze‐thaw cycles on coupling of TDR probes and how this may affect the measured dielectric permittivity.

Uncertainties Associated With the Experimental Design

There is some risk that water flow into the firn pit may have affected the observations, as this is a disturbed environment and ice layers were broken during the excavation. We cannot preclude this, but we see no evidence of unnatural influences from the firn pits. The thermal and hydrological evolution was similar at our two sites and in the radar data, so disturbances from the excavated pits would have had to proceed in the same way at each observation site. We also dug to 4.5‐m depth, removing all ice layers that could have acted as barriers to meltwater infiltration, but infiltration was limited to the upper ∼2 m. Surface gradients at DYE‐2 are less than 0.5° and melt rates are low, with no observations of ponded or flowing liquid water at the surface at this site; meltwater infiltration is primarily via vertical percolation, with no expectation of water flowing into the disturbed pits. The TDR probes are also inserted 30 cm into the pit walls, where they are positioned to register natural vertical infiltration of meltwater from above, rather than liquid water conditions within the firn pits.

Surface Melt Estimates

We estimate a total summer melt of 0.18 ± 0.04 m w.e. based on the snow‐surface height and density data at our two firn pits and estimates of snow densification. Energy balance modeling for the summer 2016 melt season gives 0.16 m w.e. with the reference model and 0.21 m w.e. of melt with the CLp7 parameter settings, which give a better fit to the observed meltwater infiltration and firn temperatures. For the CLp7 model, 74% of meltwater refreezes in the seasonal snow (0.16 m w.e.) and 26% (0.05 m w.e.) refreezes in the underlying firn, from 0.8 to 2‐m depth below the surface. Heilig et al. (2018) estimated a total summer melt of 0.15 m w.e. in summer 2016, with 60% of the meltwater (0.09 m w.e.) refreezing in the seasonal snowpack and 40% (0.06 m w.e.) percolating into the underlying firn. Our melt values are higher than those of Heilig et al. (2018), but the estimates of infiltration and refreezing in the firn agree well.

Changes in near‐surface firn ice content from 2016 to 2017 provide an independent estimate of the total melt in summer 2016. Firn cores were acquired each April and were logged at 5‐cm intervals for density and ice layers. There is lateral variability in meltwater infiltration and refreezing on scales of 1 m (e.g., Dunse et al., 2008), which makes it difficult to compare the detailed stratigraphy of firn cores from different locations or years. Changes in total ice content in the upper firn column over one year nevertheless provide an independent estimate of summer meltwater. The TDR and thermistor records indicate that the summer 2016 thaw front reached a depth of 1.8 m at our sites. In April, 2016 (before the onset of melt), the upper 1.8 m of the firn core at site A had a total ice content of 0.15 ± 0.02 m. Based on the snow‐surface height data (Figure 2b), ablation in summer 2016 reduced this 1.8 m of snow and firn to a layer ∼1.55 m thick. When we acquired a new firn core at this site in April 2017, 0.84 m of fresh snow had accumulated on the 2016 summer melt surface. The upper 2.4 m of this core (0.84 m of fresh snow and 1.56 m of firn) contained 0.35 ± 0.04 m of ice, indicating an increase in ice content of 0.2 ± 0.04 m over the year, or 0.18 ± 0.04 m w.e. Within uncertainties, this agrees with our melt estimates from the UDG records and both the reference and CLp7 models. This assessment assumes that all melt from summer 2016 refroze in visible ice layers, which may not be the case for capillary water that refreezes within the pore space. There may be some refrozen meltwater that we don't detect due to this, so estimates of total summer melt based on the ice stratigraphy represent a lower bound.

Firn cores at DYE‐2 also provide longer‐term context for the 2016 observations and the firn modeling from 1997 to 2018. The 20.3‐m core acquired in 2016 had a total ice content of 6.4 m and an average density of 665 kg m−3, giving an estimated total of 13.5 m w.e. Assuming an average accumulation rate of 0.36 m w.e. yr−1 (Mosley‐Thompson et al., 2002), this represents ∼38 years of accumulation, representing the period 1978–2016. If all of the annual meltwater refreezes somewhere in the firn column (i.e., assuming there is no runoff from the site), the 6.4 m of total ice content equates to an average melt/refreeze rate of 0.15 m w.e. yr−1. Similar to the argument above, this estimate represents a lower bound as it only accounts for the visible ice layers and lenses identified in the core stratigraphy. Modeled melt rates averaged 0.14 m w.e. yr−1 from 1978 to 2016 and 0.16 m w.e. yr−1 from 1950 to 2020 (Table 5). These values are consistent with the expectations from the observed firn‐core ice extent.

Firn Hydrology at DYE‐2

Meltwater Infiltration

The observed meltwater infiltration depth in summer 2016, 1.8 m, is similar to the mean maximum meltwater infiltration depth in the ERA5‐forced model run, 1.9 ± 1.3 m. The summer 2016 melt season therefore appears to have been relatively normal for the period 1950–2020 (Table 5). Given the average accumulation rate of 0.36 m w.e. yr−1, this means that summer meltwater at DYE‐2 typically infiltrates the seasonal snow and firn from the previous 2–3 years, with 100% of the surface melt retained through refreezing. In heavy melt seasons such as 2012 and 2019, we model meltwater infiltration to depths of 6–7 m in the firn, through 10–12 years of net accumulation.

The TDR data and model experiments allow a closer analysis of meltwater infiltration and firn hydrological processes at DYE‐2. There were four marked abrupt warming/latent heat release events in the 2016 melt season, on June 10–11, June 22–23, July 19è20, and August 9–11 (Figure 3). Average modeled melt rates on these days were 10.0 mm w.e. d−1, compared with a JJA average of 2.2 mm w.e. d−1. Meltwater was confined to the upper ∼0.3 m of the snowpack in the first of these events, with each subsequent event propagating deeper into the firn. Subsurface warming in each case was associated with high production of meltwater and limited near‐surface refreezing, which drove meltwater infiltration. On average for the summer, 49% of total meltwater refroze in the near‐surface (the upper 0.2 m), but near‐surface refreezing fell to 14% during the four melt events. The July 19–20 event was particularly strong, with average melt rates of 17.0 mm w.e. d−1 but none of this meltwater refreezing in the nea ‐surface (i.e., 100% percolating to below 0.2 m depth).

Meltwater infiltration can be limited by ice layers, if these are able to act as low‐permeability barriers to gravitational drainage (e.g., Bezeau et al., 2013; MacFerrin et al., 2019). Samimi et al. (2020) reported effective meltwater penetration through continuous ice layers up to 0.12 m thick at DYE‐2 in summer 2016. No ice layers were thicker than this in the upper 2 m of the snow and firn in summer 2016, so we cannot assess the effective permeability of thick ice slabs, which are common at lower elevations in the Greenland percolation zone (MacFerrin et al., 2019). The 0.12‐m ice layer was located from 0.96 to 1.08 m depth at site B, and we instrumented this firn pit immediately above and below this layer. Temperatures below 1.1 m remained sub‐zero until the last major meltwater infiltration event on August 9, at which point the entire snow/firn column at site B developed wet, temperate conditions (Figure 3e). Average water content from 1.1 to 1.6 m depth increased from 0.5% on August 7%–% to 1.6% from August 9–12. The ice layer may have limited meltwater percolation prior to August 9, helping to maintain frozen conditions in the underlying firn. Alternatively, this may have been due to limited meltwater supply. Average modeled melt rates from July 21 to August 8 were 3.8 mm w.e. d−1, then tripled to 11.0 mm w.e. d−1 from August 9–11. This final meltwater pulse provided sufficient latent heat to thaw the firn below 1 m, supporting the hydrological breakthrough to a depth of 1.8 m.

Infiltration Velocity

The wetting‐front or infiltration velocity, v , can be characterized from the speed of infiltration during the four main meltwater pulses. As an example, the subsurface warming event on June 22–23, described in Section 4.1.3, involved meltwater infiltration to 0.5 m depth over a period of ∼12 h, indicating a meltwater infiltration velocity of order 10−5 m s−1. Infiltration velocities reach ∼10−4 m s−1 in some instances, for example, meltwater percolating 0.3 m within one hour. An analysis of meltwater arrival times at the TDR probes during all of the infiltration/abrupt warming events gives a range of 1–12 h for the meltwater to travel 0.2–0.4 m distance between probes in the upper 1.8 m of each pit. This equates to infiltration velocities v from 5 × 10−6 to 8 × 10−5 m s−1. The lowest value, 5 × 10−6 m s−1, corresponds to meltwater transmission through the 0.12‐m thick ice layer at site B. The TDR sensor at 1.1‐m depth first registered meltwater at 03:00 on July 20, 12 h after it was first detected by the TDR sensor at 0.9‐m depth. Liquid water content below 1.1 depth remained low until the August 9 meltwater pulse, so this ice layer may have suppressed infiltration.

Differences in the infiltration velocites may be due to a combination of factors, including preferential flow pathways and delays associated with wetting and thawing of the snow and firn. The initial pulse of meltwater can be held through capillary water retention before additional meltwater can drain. Meltwater can also hit thermal barriers to penetration before it can infiltrate to greater depths, that is, limits to infiltration associated with refreezing in cold snow and firn. The thermal and wetting fronts codeveloped in summer 2016 at DYE‐2, so the infiltration velocity implicitly includes delays due to the effects of both capillary water retention and meltwater refreezing. Once the upper ∼1.2 m of snow/firn was temperate through the main melt season from mid‐July to mid‐August, there was no evidence of diurnal cycles in liquid water content. This indicates that once snow was wetted and temperate, meltwater was able to drain as quickly as it is produced. Maximum modeled daily melt rates in summer 2016 are 2 × 10−4 m s−1.

The propagation speed of the wetting front, v , is related to the Darcy velocity, q , via Equation 16, as a function of the saturation and porosity: q  = v S θ = v θ . For typical liquid water contents of θ  = 0.02 above the wetting front, q  ∼ 10−7–2 × 10−6 m s−1. If capillary forces are neglected, q  ∼ k , so these observations provide a rough observational estimate of the hydraulic conductivity in Equation 14. This is a minimum estimate, as capillary tension is a negative force, partially offsetting the gravitational force in the hydraulic potential; with a reduced potential gradient, higher values of k would be needed to explain the observed infiltration velocities. Refined, direct estimates of k should be possible from the TDR methods where capillary forces (i.e., grain size, pore structure, etc.) are directly measured or modeled, and this would be a productive avenue to pursue in further research. Within our modeling, hydraulic conductivity is a function of porosity and liquid water content (saturation), per Equation 15, so it varies with depth and time. As an example for active melt days, the average modeled value in the wetted snow and firn (the upper 1.4 m) on July 22–23, 2016 was k  = 0.9 × 10−5 m s−1, with a maximum value of ∼2 × 10−5 m s−1 for our parameter settings. This is characteristic of hydraulic conductivity values for firn reported by Fountain and Walder (1998).

Model Sensitivity to Hydrological Parameters

Model results are insensitive to changes in our reference value of the effective hydraulic conductivity between 10−3 and 10−7 m s−1, exceeding the range of variations inferred from the infiltration velocities. This indicates that k is adequately constrained in the model for the conditions at DYE‐2, as long as delays associated with capillary water retention and thermal barriers to infiltration are accounted for. In contrast, modeled meltwater infiltration depths are highly sensitive to the treatment of irreducible water content, θ (Figure 6), as identified in previous studies (Reijmer et al., 2012; Verjans et al., 2019). Higher values of capillary water retention reduce meltwater infiltration into the firn and reduce the net melt, since meltwater that is retained near the surface experiences multiple freeze‐thaw cycles during the summer melt season. Our reference model values of θ are typical of the irreducible water content reported in snow hydrology studies, ∼3% by volume (Colbeck, 1974; Coléou & Lesaffre, 1998), although our optimized models have values closer to 2%. Accurate quantification of capillary water retention and its interaction with preferential infiltration pathways require further constraint in firn models. Targeted TDR experiments in polar snow and firn could help to examine this.

Meteorological Conditions During the Melt Events

Almost half of the total melt in summer 2016 was generated during the 9 days represented by the four main melt events. These days were characterized by warm temperatures (average of −0.4°C, vs. a JJA average of −4.6°C) with minimal overnight cooling (Figure 2a). The net radiation, sensible heat flux, and latent heat flux were all above average during these events, giving an average net energy, Q , of 40 Wm−2, compared with a JJA mean of 5 Wm−2. This was driven by high values of incoming longwave radiation, Q ↓: an average of 292 Wm−2, compared with the mean JJA value of 240 Wm−2. Incoming shortwave radiation was anomalously low during these events, indicative of overcast conditions that maintained high values of Q ↓ for a period of two to three days during each melt event, suppressing overnight refreezing. The continuous supply of meltwater and the suppression of near‐surface cooling/refreezing facilitated the meltwater infiltration and firn warming events. The regional‐scale meteorological control of these events is consistent with the similar, stepwise temporal evolution of meltwater infiltration and subsurface warming seen at our two pits and at the radar site.

Long‐Term Firn Evolution at DYE‐2

Decadal‐scale firn warming is expected in Greenland in response to increases in air temperature and summer melting since the 1990s, which have driven strong reductions in surface mass balance (Hanna et al., 2021; Mouginot et al., 2019; Shepherd et al., 2020). This is supported by available firn temperature measurements and modeling at DYE‐2, with previous analyses indicating a 10‐m warming trend of ∼1°C per decade for the period 1998–2017 (Vandecrux, Fausto, et al., 2020). Modeled 10‐m and 20‐m firn temperature trends in our study equal +0.7 per decade over this period (Figure 9c), increasing to +1.0°C per decade for the period 1997–2020, due to the influence of the 2019 melt season. These trends are statistically significant, in contrast with the lack of a significant trend in surface air temperature or meltwater production over this period.

The time series is short for evaluation of temporal trends at DYE‐2 and results are sensitive to initial conditions (the model spin‐up) and the precise time window. The latter is important due to the large interannual variability; the positive trend in the model is partly due to the strong melt seasons of 2012 and 2019, which caused warm anomalies of 3°C–4°C at 10‐m depth. Firn temperatures at 10‐ and 20‐m depth fluctuate with time as a consequence of latent heat release, in contrast to a purely conductive environment, where 10‐m temperatures are representative of the mean annual surface temperature. In our modeling, the 10‐m firn temperature reaches a maximum value of −12.8°C in early 2013 (Figure 8b), in response to the latent heat release from refreezing of the summer 2012 meltwater. A series of moderate melt seasons from 2013 to 2018 allowed 10‐m temperatures to regress towards their 1997–2018 mean. Vandecrux, Fausto, et al. (2020) also documented the transient excursion associated with the 2012 melt event at DYE‐2, with a temporary loss of firn cold content but a restoration to more normal conditions by 2017. Firn temperature evolution at DYE‐2 over the last two decades is dominated by the extreme 2012 and 2019 melt seasons, which strongly influence the 10‐m temperature trend.

Air temperature records from coastal and ice‐sheet stations indicate a warming trend in Greenland from 1991 to 2018, but with no significant temperature trend over the period 2001–2019 (Hanna et al., 2021). Warming in the early 2000s was offset by several cooler years in the period 2013–2018. The model results forced by ERA5 data are consistent with this, with significant increases in air and firn temperatures since 1990 but no trends in air temperature or modeled surface melt in the GC‐Net period, 1997–2018. The extreme melt seasons of 2019 and 2012 are interesting anomalies within this period of relative climate stability at DYE‐2. These years featured two of the three deepest meltwater infiltration events of at DYE‐2 since 1950. Strong melt seasons with similar impacts on 10‐m firn temperatures also occurred in 1966, 1968, and 1978, with modeled meltwater infiltration and firn warming in 1968 exceeding that of 2012. Intermittent heavy melt seasons that cause transient firn warming are therefore not historically unprecedented. It requires several years for firn temperatures to recover from such events, however, so there are strong cumulative impacts from two such events in the last decade.

Broader Implications for the Greenland Ice Sheet

The TDR methods and the coupled model of firn thermodynamics and hydrology in this study are applicable to the entire Greenland accumulation area (i.e., about 80% of the ice sheet). While the inner accumulation area of Greenland seldom experiences melt at present, this interior region is expected to experience increased melt and evolve into a percolation area in future decades (e.g., van den Broeke et al., 2016). DYE‐2 has low surface slopes and modest summer melt totals, and is expected to be representative of much of the accumulation area in future decades. Most of Greenland's accumulation area is low sloping, so meltwater infiltration is well approximated by the vertical (gravitational) flow assumption in the model. Ongoing model improvements are needed, such as physically based treatments of snow/firn metamorphism, grain size development, capillary water retention, and preferential infiltration pathways. Field studies designed to improve understanding of these processes would support the calibration and validation of such models, and TDR methods appear promising to this end.

Meltwater percolation and refreezing in the supraglacial snow are also important processes in the ablation zone, reducing and delaying ablation and runoff. The model is applicable in this setting, configured as a multi‐layer snowpack overlying glacier ice (e.g., Samimi & Marshall, 2017). However, more complete hydrological models are needed in the ablation area and the lower percolation area, where horizontal routing of meltwater must be explicitly modeled. More consideration is also needed for meltwater infiltration processes in regions with thick ice slabs (MacFerrin et al., 2019). The permeability of these ice layers is uncertain, and firn hydrological models need to consider horizontal flow pathways where meltwater spreads out on contact with low‐permeability ice layers. Some of this will refreeze and some is likely to contribute to ice sheet runoff. Finally, the model is relevant to simulation of meltwater infiltration in deep firn, that is, for modeling the recharge of firn aquifers (Forster et al., 2014), but three‐dimensional flow also needs to be considered to simulate meltwater drainage in firn aquifer environments.

Conclusions

Vertical arrays of TDR probes show promise as a method to directly track meltwater infiltration and liquid water content in firn. When combined with thermistors, these provide a coherent picture of the coupled firn hydrological and thermodynamic processes at DYE‐2. Meltwater infiltration in summer 2016 was associated with four main pulses of meltwater and latent heat release and was similar at our two measurement sites. The maximum depth of meltwater infiltration, 1.8 m, coincided with the thaw front and is close to the mean annual meltwater infiltration depth modeled at DYE‐2 over the period 1950–2020, 1.9 ± 1.3 m.

A coupled model of the firn thermodynamic and hydrological evolution using the standard parametrization of irreducible water content by Coléou and Lesaffre (1998) underestimates meltwater infiltration and firn temperatures. Through optimization, we find that a 30% decrease of the irreducible water content gives the best match between simulations and observations through the 2016 melt season. Simulations with the optimized irreducible water parameterization provide insight into the longer‐term firn evolution at DYE‐2 over the multi‐decadal timescale of ERA5 climate reanalyses, 1950–2020. There is a significant long‐term warming signal at DYE‐2 for mean annual temperatures, +0.14°C per decade (dec−1), but summer (JJA) and firn temperatures show no significant long‐term trends in the ERA5‐forced reconstructions from 1950 to 2020. As seen in coastal and ice sheet temperature records, however (Hanna et al., 2012, 2021), there are significant warming trends in Greenland since the early 1990s. Over the period 1990–2020, increases in summer air temperature (+0.19°C dec−1) and meltwater production (+33 mm w.e. dec−1) drove an increase in summer meltwater infiltration depth of +0.36 m dec−1 in our simulations. The resulting trend in 10‐m firn temperature is +1.0°C dec−1, i.e., a 3°C increase in 10‐m firn temperature at DYE‐2 since 1990.

Firn temperature, ice content, and density at DYE‐2 are strongly affected by extreme melt seasons such as 2012 and 2019. High amounts of melt in these two summers, ∼4σ above normal, drove meltwater infiltration to depths of 6–7 m in the model, causing a 10‐m firn warming of 3–4°C that can persist for several years. Firn temperatures from 10 to 20 m depth were almost restored to their pre‐2012 conditions when the 2019 melt event produced similar impacts. These events are not unprecedented in the 71‐year historical reconstruction at DYE‐2, but there are strong cumulative impacts due to two such extreme melt events in a single decade. These extreme melt seasons increase firn temperature, density, and ice content, with implications for meltwater retention in the upper Greenland percolation zone. The simulations suggest that changes in firn structure and meltwater retention capacity may be more sensitive to the frequency of extreme melt seasons than the background climate warming trend.

Further understanding of firn processes, refinement of firn models, and representation of these processes in large‐scale ice sheet models are all needed to support improved understanding of Greenland Ice Sheet response to climate change. TDR is a promising technique to examine firn hydrological processes and calibrate firn models. Future application of TDR at locations with greater meltwater production and ice layer thickness could further increase our understanding of meltwater infiltration in polar firn. These data may help to refine our model treatments of capillary water retention and preferential flow pathways. Two‐dimensional networks of TDR probes and thermistors could provide insight into the heterogeneity of meltwater infiltration and refreezing processes.

The atmospheric warming, increased meltwater production, and evolving firn structure at DYE‐2 are characteristic of most of the percolation area on the Greenland Ice Sheet. As melt advances to higher elevations on the ice sheet in future decades, areas previously untouched by melt may regularly see meltwater infiltration and ice‐layer development similar to the conditions at DYE‐2. Realistic treatments of these processes will be increasingly important for projections of Greenland Ice Sheet surface mass balance and Greenland's contributions to global sea level rise.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.