Quantifying meltwater refreezing along a transect of sites on the Greenland ice sheet

On the Greenland ice sheet, a significant quantity of surface meltwater refreezes within the firn, creating uncertainty in surface mass balance estimates. This refreezing has the potential to buffer seasonal runoff to future increases in melting, but direct measurement of the process remains difficult. We present a method for quantifying refreezing at point locations using in situ firn temperature observations. A time series of sub-hourly firn temperature profiles were collected over the course of two melt seasons from 2007 to 2009 along a transect of 11 sites in the accumulation zone of Greenland. Seasonal changes in temperature profiles combined with heat flux estimates based on high-temporal-resolution temperature gradients enable us to isolate the heat released by refreezing using conservation of energy. Our method is verified from winter data when no refreezing takes place, and uncertainty is estimated using a Monte Carlo technique. While we limit our method to a subsection of firn between depths of 1 and 10 m, our refreezing estimates appear to differ significantly from model-based estimates. Furthermore, results indicate that a significant amount of refreezing takes place at depths greater than 1 m and that lateral migration of meltwater significantly complicates the relationship between total surface melt and total refreezing.


Introduction
The mass balance of the Greenland ice sheet has become increasingly negative over the course of the last decade (Shepherd et al., 2012;Rignot et al., 2011;Vaughan et al., 2014). Recent measurements of outlet glacier discharge indicate that surface mass balance is now the dominant source of mass loss (Enderlin et al., 2014), and regional climate-surface process coupled models show that most of the increases in surface mass loss are due to significant increases in surface melting and runoff beginning in the early 1990s (Ettema et al., 2009). Increases in the areal and temporal extent of surface melt, evident from remote-sensing, support-model-based increases in surface melting. However, the increase in meltwater leaving the ice sheet is not as well constrained.
Refreezing of surface meltwater, as it infiltrates into the underlying cold firn, creates a significant source of uncertainty in both remotely sensed and model-based estimates of Greenland mass balance. Modeling studies have estimated that almost half of the meltwater generated annually in Greenland refreezes (Ettema et al., 2009). However, the complexity of the infiltration process remains difficult to incorporate into snowpack models and infiltration hydrology is usually modeled as a purely uniform process (Greuell and Konzelmann, 1994;Mernild et al., 2010;Fettweis, 2007;Ettema et al., 2009;Bougamont et al., 2005). There have been some efforts to quantify refreezing using a variety of parameterizations typically based on simple energy balance ideas where refreezing is controlled by the total cold content of the firn (Pfeffer et al., 1991;Pfeffer and Humphrey, 1996;Reeh, 1991;Janssens and Huybrechts, 2000;Oerlemans, 1991) (see Reijmer et al., 2012, for an overview). However, these parameterizations remain largely unverified with in situ data (Reijmer et al., 2012). A limited number of studies have attempted direct quantification of refreezing (Wright et al., 2007) or indirect measurements using density (Parry et al., 2007). However, these techniques are resource intensive, non-continuous, and do not adequately capture the complex three-dimensional spatial variability of meltwater refreezing.
The most challenging aspect of quantifying refreezing is that both infiltration and refreezing of meltwater in cold firn is highly heterogeneous in space and time. Infiltration is characterized by a complex network of rapidly developing ice lenses and vertical pipe structures (Pfeffer and Humphrey, 1998). Vertical pipe flow can transport water deeper into the firn than uniform infiltration. Recent field observations in the percolation region of the Greenland ice sheet accumulation zone suggest that meltwater is able to penetrate cold firn to depths greater than 10 m and remain mobile throughout the winter Forster et al., 2014). Piping of meltwater into cold firn occurs without warming the entire profile. The water moves down the pipes with only minimal heating of the surrounding firn, making a locally complex temperature field of both cold and near-melting temperatures. This deep penetration and heterogeneous heating make rudimentary parameterizations based on Pfeffer et al. (1991) questionable.
Here we demonstrate that temperature measurements in the firn provide an alternative method for investigating and quantifying refreezing. Latent heat released during refreezing diffuses through the firn, causing a thermal perturbation in the temperature profile that can be quantified using a conservation of energy approach. Use of thermal data has advantages: installation of thermal sensors is relatively easy with handheld snow drills and logging of thermal measurements is simple and robust. Furthermore, using measured temperatures takes advantage of the diffusive nature of heat conduction which helps reduce the effect of extreme spatial discontinuities inherent to heterogeneous infiltration and refreezing processes. The method is not without its challenges as melting, accumulations, and solar radiation in the near surface make it difficult to account for energy transfers from temperature measurements alone. However, by limiting our domain to firn depths between 1 and 10 m, we are able to apply our method to a transect of melt season thermal profiles on the Greenland ice sheet. The result is the first in situ measurements of refreezing on the Greenland ice sheet that completely span the percolation zone.

Field measurements
Firn temperature data were collected in 2007, 2008, and 2009 from a transect of 11 sites in southwestern Greenland, about 100 km northeast of Jakobshavn Isbrae (Figs. 1 and 5). All of the sites are within the accumulation zone, and annual snow accumulation along the transect is on the order of 1 m (density of 0.35 g cm −3 ) (Parry et al., 2007;Benson, 1962;Hanna et al., 2006). In contrast to accumulation, there is a strong gradient in the degree of summer surface melting in this region that is a result of changing elevation and albedo. This important subregion of the accumulation zone is known as the percolation zone, and the gradient in melt is reflected in the physical and thermal characteristics of the underlying firn (Cuffey and Paterson, 2010). The highest elevation site (Crawford Point) is near the upper edge of the percolation zone where summer surface melting is rare. In contrast, the lowest elevation site (H4) is within 20 km of the equilibrium line altitude (ELA), and enough melting takes place that layers within the firn can, at times, become saturated with water. While extensive refreezing, and thus warming, takes place in the percolation zone, parts of the underlying firn remain subfreezing even at the lowest elevation of H4 . Sensor spacing is 0.25 m from 0 to 5.5 m depths and 0.5 m from 5.5 to 10 m depths. The sensors were installed with reference to the surface at time of installation. After temperature string emplacement, the boreholes were backfilled with fine-grained, cold snow, and our temperature measurements show rapid thermal equilibrium with the surrounding undisturbed firn pack (For details see Harper et al. (2011). During the subsequent year or more of data collection, the surface at most sites showed some net snow accumulation of the order of 0.5 m or less and seasonal ablation of the same magnitude. Thus, quoted sensor depths are not relative to the actual surface at subsequent times. Site density profiles were obtained from 10 m snow cores extracted during temperature string emplacement . Field density measurements were obtained at variable spacings in order to accurately sample observed firn stratigraphy. Since the thermal sensors were regularly spaced, the density data were averaged and re-sampled on a matching 0.25 m spacing.
Not all sites were instrumented at the same time, for the same period, or with exactly the same sampling interval (Table 1). The higher elevation sites from CP to T1 recorded data in the summer and fall of 2007, while the lower sites from T1 to H4 recorded data from the summer of 2008 to the spring of 2009. Two temperature strings located 20 m apart were installed at site T1 in 2008 to investigate lateral heterogeneity of the temperature field (see Fig. 5a). Most of the summer data have a 20-30 min sampling interval. However, power requirements limited the sampling interval of the 2008/2009 winter data to every 8 h.

Method theory
Our approach quantifies refreezing using the change in heat content at each site as measured by the change in the vertical temperature profile over the summer season (Fig. 2). When we assume only vertical temperature gradients (discussed below), the change in the heat content of a section of firn results from the net conduction flux across the top and bottom boundaries plus the advection of latent heat associated with refreezing of infiltration water. Both the change in heat content and the net heat conducted across the boundaries can be estimated from the data. The latent heat released  during refreezing is the difference between the total change in heat content and the net heat conducted across the section boundaries.
Heat from refreezing = change in heat content ( H ) where is used here and in the following discussion to indicate change of a variable over the data time period, typically over the summer melt season. A one-dimensional approach requires that horizontal temperature gradients within the firn are negligible when averaged over the summer melt period despite the horizontally heterogeneous nature of melt infiltration. Fortunately, the dif-fusive nature of heat flow ensures that horizontal gradients in temperature decay rapidly as long as the length scale of the lateral heterogeneity is less than the vertical scale. Furthermore, the gross firn structure, which dominates both thermal and hydrologic properties, is predominantly horizontally layered, and we therefore assume lateral variations in both temperature and infiltration should be stochastically distributed. A random distribution of lateral temperature variations coupled with rapid decay in horizontal temperature gradients implies that the transfer of heat laterally should sum towards zero on a seasonal timescale.
Our assumption that lateral gradients can be ignored is given credence both by analysis of the data from the two adjacent temperature strings located at T1 (Fig. 5a) and by a theoretical scaling analysis. The two profiles at T1 show occasional, significant differences during melt events that last a few days (also see Humphrey et al., 2012), demonstrating that firn temperatures show some local, lateral variability. However, the calculated refreezing quantities (see Table 1) are within measurement error of each other, indicating that the temperature variations are insignificant when averaged over the summer season. The significance of the lateral variations can be investigated theoretically by using the analytical solution of exponential thermal decay of a line heat source (representing a vertical pipe) in a three-dimensional homogeneous firn solid (Carslaw and Jaeger, 1986). As an example of this scaling for 1 m pipe spacing, using appropriate parameters for firn, the analytical solution yields a temporal decay scale of 3 days, which also supports our assumption that horizontal gradients decay over the summer season. When the pipe spacing is larger, on the order of 10 m, then our assumption is not valid over the season. We note that the horizontal spacing of observed pipes was typically considerably less that 10 m . From this we conclude that the lateral temperature variations at each site averages out to a uniform system on seasonal timescales.
The change in heat content over the summer melt season can be quantified from the changes in profile temperature (z = depth from the top of the profiles) using This formulation assumes that density (ρ) and heat capacity (C p ) do not change over time. This is reasonable because, at most of the sites, the input of meltwater is minor compared to the water equivalent of the firn column. This assumption may break down for the lowest sites (H3, H4). Densification due to compaction of the firn is also assumed to be minimal within our seasonal timescale.
The net heat conducted through the boundaries is where q net (t) is the net heat flux as a function of time and is defined by Fourier's law operating at the boundaries of the profile (K is thermal conductivity).
The net heat flux is integrated over the time period corresponding to T in Eq. (1) (see Fig. 2). This time period is typically the summer melt season.

Method implementation
Numerical approximations to Eqs. (1)-(3) are used to calculate refreezing quantities from the observed temperature profile data. We strive to calculate refreezing quantities at each site corresponding to the entire melt season. However, leakage of the water into the data loggers at some of the 2007 sites resulted in sections of unusable data, and we are therefore forced to limit our analysis at some sites to shorter time periods (i.e., site T1-07, see Table 1).
The method is applied to firn depths ranging from 1 to 10 m, and we therefore ignore the data from the upper four sensors. We refer to this subsection of the firn pack as our analysis domain. This domain is deep enough to remain unexposed, as melting, sastrugi migration, and accumulation lead to significant variations in the surface elevation. Furthermore, the influence of solar radiation is greatly reduced below about 0.5 m. Refreezing that occurs above or below the domain remains unaccounted for by this analysis, but, as is shown below, we estimate it captures a majority of total refreezing. Heat content in this domain is assumed to change only from conduction across the region boundaries and advection of heat energy in the form of the phase change of refreezing percolating meltwater. We assume meltwater is at 0 • C and has no additional sensible heat. Firn densities are derived from the borehole core data. Heat capacity is considered to be constant at 2097 J kg −1 • C −1 . The boundary temperature gradients ( dT dz ) in Eq.

The
(3) are approximated using the temperature gradient of the two sensors closest to the 1 and 10 m bounds. Equation (2) is approximated by numerically integrating net heat flux using the trapezoid rule. Figure 2a shows a time series of net heat flux at site H2. High-frequency variations on the order of 0.5 • C are a result of random electronic noise in each temperature measurement. Since this noise is random, the integrated flux derived from the gradients is not biased. Thermal conductivity (K) is calculated as a function of the boundary densities following Schwerdtfeger (1963) (see Appendix).

Error analysis and testing
Since our method differentiates discrete data when the heat flux is calculated, we assume that the largest errors stem from amplification of data noise by differentiation. Additionally, uncertainties in our profile density measurements create the potential for error as the values are utilized to calculate thermal conductivities and are direct inputs into Eq. (1). A Monte Carlo approach is used to estimate how these errors contribute to overall method uncertainty. Random perturbations were added to the observed temperatures and densities, and then the refreezing quantity was recalculated for each case. Since the noise in the data is electronic, the distribution of the perturbations is Gaussian with a standard deviation equal to the uncertainties in the observed data. Temperatures showed noise with an uncertainty conservatively estimated to be 0.5 • C and uncertainty in density is estimated to be 20 kg m −3 . After 1000 Monte Carlo trials, the average and standard deviation of the refreezing is then used to estimate mean refreezing values and corresponding uncertainty.
We test our method using profile temperature data from the winter season. Temperature profiles were chosen from the data in which no surface warming approaching 0 • C occurred and when the quantity of refreezing is assumed to be zero. With no latent heat input, our energy balance should show that the temperature change within the firn is exactly balanced by the heat flow across the boundaries, and errors in the method would show up as spurious melt or refreezing. Multiple tests were performed at all sites except T4 where there are no temperature data outside of the melt season (roughly May-September). In order to verify consistent results, multiple tests were performed with time spans ranging from 1 to 4 months and include data with different sampling frequencies. With the exception of three sites, all tests resulted in refreezing quantities within 1 cm w.e. of zero. Since the two standard deviation uncertainty bounds on all refreezing estimates are on the order of 1.5 cm w.e., these tests confirm the method produces a refreezing value not significantly different from zero in the winter.
At sites H165, H2, and H3, tests showed that the method produced refreezing quantities on the order of −2 cm w.e., indicative of a mismatch between the change in the firn in-ternal temperature structure and the flow of heat across the boundaries. In the winter, a negative refreezing value results from the temperature profile losing more heat than would be expected given the calculated heat flux through the domain boundaries. The heat flux needs to be higher (more heat loss) in order to balance the profile temperature change. The most important parameter in this balance is the firn conductivity at the boundary, and we find that a small increase in conductivity at the upper boundary eliminates the energy imbalance. This increase in conductivity implies that these sites require a corresponding increase in the density of the boundary firn. We found that, in all cases, the mismatch can be eliminated when the densities near 1 m depth are increased by around 200 to 600 kg m −3 . This is a reasonable near-surface seasonal density change, and furthermore, the near-surface measured densities (obtained at the beginning of the melt season) at these sites are unusually low compared to the underlying firn. It should be noted that although the above discussion is somewhat speculative, this same density change applied during the melt season has minimal effect on our calculated refreezing quantities as the melt season temperature gradient near 1 m is often near zero.

Results
The calculated quantity of water refreezing between 1 and 10 m depths at each site is plotted in Fig. 3. The error bars on each value are equal to 2 standard deviations of the variability generated in the Monte Carlo trials. The higher elevation sites show melt results calculated from 2007 temperature data, while the lower elevation sites use data collected from 2008 (Table 1). As might be expected due to increases in melting, refreezing quantities generally increase with decreasing elevation. Sites T2 and T1 have temperature profile data from both 2007 and 2008. In addition, there were two temperature profiles available at site T1, placed about 20 m apart. Refreezing values calculated from these two profiles show no significant difference from one another.
The results show a large difference in overall refreezing magnitudes between 2007 and 2008 that reflect the overall melt conditions experienced during each melt season (2007 was a high melt year, see Fettweis et al., 2011). Unfortunately, data quality problems at some of the sites reduced the time period over which our method could be applied (see Table 1 column 4). With the exception of site T1, refreezing quantities in 2007 correspond to July and August, while all sites in 2008 correspond to June-August. Site T1 in 2007 corresponds to only July. If the T1-07 data also included June, the difference in refreezing trends between 2007 and 2008 would likely have been even more substantial than is shown in Fig. 3.
It is useful to give our results some context by comparing them to two separate and independent analyses. First, a simple positive degree-day melt model (PDD), following Hock  Table 1. (2005), is used to calculate a plausible surface melt range at each site (Fig. 3). In Greenland, empirically determined melt factors (DDF) for snow range from 2.5 to 5 mm day −1 • C −1 (Hock, 2003;Janssens and Huybrechts, 2000;Braithwaite, 1995;Cuffey and Paterson, 2010). Air temperatures at 2 m height from the Crawford Point weather station are used as input to the melt model. Prior to calculating the sum of positive degree days, the raw hourly air temperatures are adjusted for each site using a seasonally variable slope lapse rate given by Hanna et al. (2005). Also, the hourly temperatures are averaged to daily values. The upper and lower bounds correspond to our estimates of maximum and minimum degreeday factors. Note that the PDD model only produces melt; it does not deal with refreezing.
We also compare our results to refreezing quantities output by the regional climate model MAR (Fettweis, 2007;Fettweis et al., 2011;Tedesco et al., 2014) (Fig. 3). MAR has a resolution of 25 km and a time step of 120 s and has been utilized in numerous studies related to modeling surface mass balance on the Greenland ice sheet (for a list see http://www.cryocity.org/papers.html). MAR utilizes the physically based, one dimensional snowpack model CRO-CUS to calculate refreezing to a depth of 15 m. MAR-based estimates of refreezing quantities were determined for each of our sites by summing the daily average values of refreezing output by MAR over the time periods shown in column  Table 1). Profiles with triangles correspond to 23 July 2007  In contrast, the MAR refreezing values corresponding to the 2008 sites are substantially higher than our values, including some values over 10 times higher (Fig. 3). There is also a significant difference in the data trends. MAR refreezing values do not show any clear difference between 2007 and 2008 or a decreasing amount of total refreezing near the ELA. Differences between our refreezing values and the PDD become more significant at lower elevations. Below H165 the refreezing begins to exceed the estimated melting, peaking at site H2, and below that the refreezing apparently becomes much less than the melt.

Discussion
The large discrepancy between our values and that of MAR may be partly a result of refreezing taking place in the upper 1 m of firn. Refreezing in the near surface is not included in our analysis as, unlike the deeper thermistors, the data from thermistors in the first meter are influenced by the effects of ablation, accumulation, and solar radiation, which are complexities that are outside the scope of our simple energy balance approach. Nonetheless, the near-surface data can still be used make qualitative interpretations of the thermal conditions in the upper meter of firn, enabling us to investigate the significance of refreezing in the upper meter. At each site, a daily mean of the average of the temperatures at depths of 0, 0.25, 0.5, and 0.75 m was calculated, and a subset of sites are plotted in Fig. 4 (some sites are omitted for clarity). At all of the sites in 2007 and most of the sites in 2008, warming was sufficient to bring the average temperature in the upper 1 m of snow to 0 • C for almost the entire melt season (Fig. 4a). In some cases, the average temperatures are even above the melting point, indicating that either the sensors are close enough to the surface to be warmed by radiative heating or even exposed by ablation or wind scour. These observations imply that the capacity for refreezing in the upper 1 m of firn at these sites is almost zero for most of the melt season. Consequently, most meltwater generated will infiltrate without refreezing until it reaches the deeper, colder firn within our method domain. Sites T2 and T1 in 2008 are the only sites with near-surface firn temperatures below 0 for a significant part of the melt season. However, the large difference in refreezing values between the two methods (our values vs. MAR) is present at T1, T2 (2008), and the other sites. It is therefore unlikely that refreezing in the upper 1 m of firn is the sole source of the difference in values.
Diurnal melting and refreezing of pooled water at the surface is also unaccounted for using our method. In this situation, refreezing takes place at the surface from radiational cooling without the need for subfreezing firn temperatures. Several melt-freeze cycles could take place before the water finally infiltrates and/or runs off. The cumulative effect could drive the increasing difference between our values and MAR sites H3 and H4. However, the significance of this process is highly uncertain since lower elevations, with fully saturated firn, also have smaller diurnal temperature changes. Furthermore, there was no evidence of extensive surface water present at any of the sites .
A final possibility is that the large difference is not due to physical parameters. The MAR output is representative of 25 km grid cells that have not been downscaled to each site. Furthermore, climate models are known to have inherent bias (Tedesco et al., 2013) and the 1-3 month time period used may be unreasonably short for assessing MAR results.
It may be more appropriate to interpret our results in relation to the PDD melt estimates, and we give two reasons. First, the temperatures used in the PDD model have been adjusted to elevation and may be better for point measurements than the coarse atmospheric model resolution. Second, a PDD melt model does not capture short duration, minor melt events that are unlikely to infiltrate to within our model domain. PDD melt factors are often calibrated using daily observations of ablation stakes to calculate melt (Hock, 2003;Braithwaite, 1995). This type of observation is more sensitive to significant melt events and may even ignore short duration or diurnal melting. Therefore, the PDD melt range shown in Fig. 3 can be interpreted as the total amount of melt that was likely to penetrate more deeply within the firn. Since the refreezing capacity of the near-surface firn is minimal at most sites (Fig. 4), we can therefore utilize the PDD melt range to interpret our results in relation to how much meltwater is infiltrating within the method domain.
In 2007, the refreezing quantities lie within or slightly below the PDD melt range, implying that a significant fraction of meltwater (>50 % in most cases) is infiltrating deep within the firn and that there is sufficient refreezing capacity in this region to capture most of it. The highest site, Crawford Point, has more cold content in the upper 1 m of firn (Fig. 4a) than most of the other sites. This could lead to more refreezing in the near surface and may explain why the total refreezing value is significantly below the PDD melt range. Refreezing quantities at sites T1 and T2 overlap with the PDD melt range in both 2007 and 2008 despite the substantial change in total melting between the 2 years. This shows that the firn has some ability to at least temporarily buffer large changes in melt as the refreezing capacity in this region was not completely eliminated during the 2007 season, or it was able to sufficiently recover over the 2007/2008 winter.
For all sites above H165 (Fig. 5b), there is no significant difference between estimated melt and refreezing values, indicating that all melt produced at the surface appears to refreeze in the upper 10 m of the firn column. In contrast, the in situ refreezing in the lower percolation zone, below H165 (Fig. 5c), cannot be simply described as full refreezing of the predicted melt. We interpret our results to indicate two separate processes occurring in the lower percolation zone. At elevations below H2, the refreezing quantities begin to decrease. Additionally, profile temperatures at the lowest site www.the-cryosphere.net/9/691/2015/ The Cryosphere, 9, 691-701, 2015 H4 are actually colder than either H3 or H2 (Fig. 5c). This is inconsistent with the expected increased melt at lower elevations and therefore must result from meltwater running off rather than refreezing . This region is the location of the runoff limit, where some of the meltwater may ultimately leave the ice sheet and contribute to sea-level rise. The zone encompassing H165, H2, and H3 is more difficult to interpret. The refreezing values in this region may be higher than the calculated melt due to lateral migration of meltwater in the firn. It is also plausible that some of the total refreezing quantity is derived from meltwater generated during a previous melt season that remained unfrozen within the firn in a manner described by Forster et al. (2014). Lastly, it is also possible that the PDD melt model significantly underpredicts melt in this region. It is interesting to note that Ambach (1988) calculated DDF for the ELA elevation of this transect region in Greenland and found a particularly high DDF, which may indicate that it is not realistic to use a single DDF in this region but that the DDF should increase with decreasing elevation. Nevertheless, this transitional region is not fully explained by this comparison with the simple PDD model.

Conclusions
Firn temperature profiles provide an effective means of estimating in situ quantities of meltwater refreezing. Our method treats the firn as a one-dimensional system and the heat conducted into and out of the system through time is tracked using temperature measurements. The latent heat released into the system by refreezing is calculated using conservation of energy, and this value is then converted to a water equivalent.
Application of our method is problematic in the upper meter of the snow profile because of radiative heating and cooling and because of unobserved snow accumulation and ablation during the yearly cycle. Any refreezing in this layer represents an unaccounted error in our melt estimates. Nonetheless, testing of the method using winter temperature measurements, when it is assumed that no refreezing is taking place, verified that, in most cases, our method is robust. Four of the sites required slight tuning of the near-surface densities, but this is reasonable given expected melt season densification of the firn.
The calculated refreezing quantities reveal a transition from complete refreezing of meltwater at higher elevations to eventual runoff of meltwater near an elevation of around 1500 m, up to 40 km inland from the ELA. Even where complete refreezing does occur, a significant portion of the overall refreezing takes place at depths greater than 1 m. This may be a result of piping of meltwater to much greater depths than would otherwise occur by uniform infiltration. Since heterogeneous infiltration is not currently accounted for in either snow hydrological models or simple theoretical parameterization, these in situ refreezing values provide an important source of snow/firn model validation. Finally, our results also give some indication that lateral movement of infiltrated meltwater, in some cases from prior melt seasons, may be significant in this region of Greenland, complicating the classic understanding of percolation zone processes.