Meteorological drivers of melt at two nearby glaciers in the McMurdo Dry Valleys of Antarctica

Abstract

We study the meteorological drivers of melt at two glaciers in Taylor Valley, Antarctica, using 22 years of weather station observations and surface energy fluxes. The glaciers are located only 30 km apart, but have different local climates; Taylor Glacier is generally drier and windier than Commonwealth Glacier, which receives more snowfall due to its proximity to the coast. Commonwealth Glacier shows more inter-annual melt variability, explained by variable albedo due to summer snowfall events. A significant increase in surface melt at Commonwealth Glacier is associated with a decrease in summer minimum albedo. Inter-annual variability in melt at both glaciers is linked to degree-days above freezing during föhn events, occurring more frequently at Taylor Glacier. At Taylor Glacier melt occurs most often with positive air temperatures, but föhn conditions also favour sublimation, which cools the surface and prevents melt for the majority of the positive air temperatures. At Commonwealth Glacier, most of the melt instead occurs with sub-zero air temperatures, driven by strong solar radiative heating. Future melt at Taylor Glacier will likely be more sensitive to changes in föhn events, while Commonwealth Glacier will be impacted more by changes in near coastal weather, where moisture inputs can drive cloud cover, snowfall and change albedo.

1. Introduction

In the largest ice-free region of Antarctica, the McMurdo Dry Valleys, polar glaciers terminate with steep cliffs in rocky and sandy valleys. In summer, melt from these glaciers fills streams that feed the closed-basin perennial-ice covered lakes. Solar radiation is the largest component of the glacial surface energy balance (SEB) in summer, but because summer air temperatures remain close to zero, the residual energy for melt is small (Hoffman and others, 2008; Hofsteenge and others, 2022). Ablation of the glaciers is monitored with a large network of ablation stakes (Gooseff and others, 2022). The changes observed with this network are small and insignificant and therefore the glaciers are considered to be in a quasi-equilibrium state (Fountain and others, 2016a).

Despite the absence of significant trends in glacial mass balance, glacial melt shows strong variability and is of great importance for the ecosystem of the McMurdo Dry Valleys. In the desert climate of this region, the small amounts of snowfall quickly sublimate and therefore glacial meltwater is the main source of streamflow. Microbial mats flourish around the streams and their distribution and abundance depend on streamflow. In high streamflow years, the mats can be damaged by erosion (Cullis and others, 2014; Kohler and others, 2015). Furthermore, increased wetness can impact the stability of microbial communities that are adapted to the dry conditions (Monteiro and others, 2022). Streamflow in high melt years also impacts the clarity, chemistry and biological activity of the lakes (Fountain and others, 2016). Understanding what drives the variability and changes in glacial meltwater is therefore of great importance to predict changes in the McMurdo Dry Valleys ecosystem in a warming world.

With summer air temperatures close to 0°C, small changes in climate are expected to have large impacts on melt. With these glaciers being on the threshold of change, there is a need to identify the drivers of melt to be able to predict future changes in meltwater generation. In this study we construct a 22 year SEB record, which is the longest record in the McMurdo Dry Valleys to date, to study the meteorological drivers of glacial melt in Taylor Valley. We compare the variability and the drivers of melt for two glaciers in Taylor Valley, which despite their proximity differ in local climate and melt regimes. By comparing two contrasting sites within the McMurdo Dry Valleys, our aim is to improve understanding of the different drivers of meltwater production in this unique region. Ultimately, this should help advance understanding of how glaciers within the McMurdo Dry Valleys will respond to climatic changes as a result of global warming.

2. Study site

Commonwealth and Taylor glaciers are located in Taylor Valley, ~30 km apart (Fig. 1). Taylor Valley has several mountain glaciers, some of which reach the valley floor and typically have 20–30 m tall terminating cliffs. Commonwealth Glacier flows down from the Asgard Range, expanding in the relatively flat valley floor as a piedmont glacier. Of all glaciers in the valley, Commonwealth Glacier is located the closest to the shoreline (~3 km). Fast ice in McMurdo Sound forms around March–April and breaks out between January and February, but rarely does so at the mouth of Taylor Valley (Explorers Cove). The more inland located Taylor Glacier is an outlet glacier of the Antarctic ice sheet and terminates into the west lobe of Lake Bonney. Taylor and Commonwealth glaciers were selected for this study since they had the most complete long-term records of observed meteorology collected by the McMurdo Dry Valleys Long-Term Ecological Research (LTER) programme (Gooseff and others, 2022), and also are expected to have a relatively different local climate due to their different proximity to the coast. The glaciers are separated by the Nussbaum Riegel, a ridge in the middle of the valley which acts as a barrier and control on the wind-driven climatic pattern in Taylor Valley (Fountain and others, 1999). We present data from two automatic weather stations (AWSs), which are located in the ablation zone of the glaciers and at comparable elevations: COHM at 282 m a.s.l. on Commonwealth Glacier (Doran and Fountain, 2023a) and TARM at 331 m a.s.l. on Taylor Glacier (Doran and Fountain, 2023e).

Figure 1. (a) Map of Taylor Glacier Meteorological station (TARM) and Commonwealth Glacier Meteorological station (COHM) in Taylor Valley. The windroses in (b) show the frequency distribution of wind direction and wind speed at TARM and COHM throughout the study period. This map was made using the Quantarctica QGIS package collated by the Norwegian Polar Institute.

3. Methods

3.1 Constructing a 22 year meteorological forcing dataset

AWS observations from the McMurdo Dry Valleys LTER programme at COHM and TARM are used to create a 22 year (1999-11-01 until 2022-01-31) meteorological forcing dataset for energy-balance modelling. The AWSs measure incoming and outgoing shortwave and longwave radiation (S _in, S _out, L _in, L _out), wind speed (WS), relative humidity (RH) and air temperature (T _a), surface height change and surface thermal infrared temperature (T _s). Details on the specifications and accuracies of the instruments are described in Gooseff and others (2022) and summarised in Table S1. Daily accumulated incoming and outgoing shortwave radiation are used to calculate accumulated albedo.

3.1.1 Gap filling: wind speed, air temperature and humidity

Wind speed, air temperature and relative humidity observations are missing at COHM for 0.6% of the time and are filled with observations from the nearby Explorers Cove AWS (Fig. 1b, Doran and Fountain, 2023b). At TARM, relative humidity and air temperature are missing for 1.4% and wind speed for 8.0% of the period. Where possible, the gaps are filled with data from Lake Bonney AWS (Doran and Fountain, 2023c) and the remaining 5.0% of wind speed are filled with data from Lake Hoare AWS (Fig. 1b, Doran and Fountain, 2023d). Air temperature from Lake Bonney and Explorers Cove are corrected for elevation difference by using a monthly mean lapse rate based on the Explorers Cove–COHM, or Lake Bonney–TARM pair. Atmospheric pressure is not measured at either of the two stations and is therefore taken from Lake Hoare AWS and corrected for the difference in station elevation.

3.1.2 Precipitation estimation

Precipitation was calculated from sonic ranging observations at both glaciers. A threshold of an increase of 5 mm in the daily surface height was used to detect precipitation (Fountain and others, 2010; Myers and others, 2022). The observed surface height is used to restrict the snow accumulation and to account for snow removal through wind drift that would occur in reality. Surface height observations cover 64% of the study period at COHM, and 77% at TARM. The remaining 36 and 23% of the period are filled by precipitation measurements using a weighing bucket at the nearby Explorers Cove and Lake Bonney AWSs. A threshold of a 0.5 mm w.e. increase in the daily bucket level was used to detect precipitation events and the following day was checked for a net increase in order to avoid false-precipitation events (Fountain and others, 2010; Myers and others, 2022). An average snowfall density of 200 kg m⁻³ was estimated from the difference between the two methods and used to convert the precipitation based on surface height to mm w.e. The yearly average snowfall density calculated using both methods, however, varied between 40 and 440 kg m⁻³, reflecting the difference in precipitation between both methods. A test run without snowfall showed little impact on the SEB and melt generation and therefore our results are not strongly impacted by the two precipitation data sources used. This is because snowfall has the largest impact on the melt generation through albedo, which is already captured by forcing the model with observed albedo in this study.

3.1.3 Gap filling radiation and derivation of cloudiness

At COHM, there are gaps in S _in for 0.9% of the record, with the largest gap being 36 d in a row. These are filled by observations at the nearby Explorers Cove AWS. Similarly at TARM, 1.5% of the S _in observations are missing and filled with observations from Lake Bonney. For the same gaps albedo is kept constant at the last value before the gap. Longwave radiation components are not measured at TARM and therefore L _in is taken from the Lake Bonney AWS and corrected for difference in sky view factor (svf) and difference in T _a using Eqn (1). svf was calculated as 0.82 at TARM and 0.78 at Lake Bonney using a 30 m digital elevation model. The 24.4% of the Lake Bonney L _in record that is missing is filled by parameterising L _in. A longwave radiation model is optimised for Lake Bonney and used to calculate L _in with TARM air temperature, humidity and cloud factor as described in the following section.

The effect of svf and air temperature on L _in can be estimated using:

(1)$$L_{\rm in} = svf \epsilon_{\rm eff}\sigma T_{\rm a}^4 + ( 1-svf) \sigma T_{\rm a}^4,\; $$

where $\epsilon _{\rm eff}$ is the effective sky emissivity, T _a the air temperature (K) and σ the Stefan–Boltzmann constant. To adjust L _in for differences in svf and T _a between TARM and Lake Bonney, $\epsilon _{\rm eff}$ is first derived from observations at Lake Bonney using Eqn (1). L _in at TARM is then estimated using $\epsilon _{\rm eff}$ from Lake Bonney and T _a and svf values from TARM.

During gaps in the Lake Bonney L _in record, L _in at TARM is estimated using measurements of S _in, T _a and RH from TARM and a parameterisation that describes the relationship between cloud effects on L _in and S _in. Here we use the relationship proposed by Conway and others (2015) which links the transmission of S _in to enhancement of L _in through two complementary cloud metrics: shortwave effective cloudiness (N _eff) and longwave equivalent cloudiness (N _ep).

N _eff is calculated using:

(2)$$N_{\rm eff} = ( 1-trc) /k,\; $$

where k is a cloud extinction coefficient that varies with vapour pressure (Conway and others, 2015) and trc is a transmission factor (Kuipers Munneke and others, 2011) calculated as trc = S _in/S _pot. For optimisation, 316 clear-sky days at Lake Bonney were identified from observed incoming shortwave radiation at the Lake Bonney AWS, based on a monotonic increase and decrease in the morning and afternoon. S _pot is optimised on S _in for days within this clear-sky dataset, using a single scatter albedo of aerosols (alphss) of 0.940 and an aerosol extinction coefficient (kaes) of 0.940 (Iqbal, 1983).

N _ep (Kuipers Munneke and others, 2011) is calculated using:

(3)$$N_{\rm ep} = {\epsilon_{\rm eff}-\epsilon_{\rm cs}\over \epsilon_{\rm ov}-\epsilon_{\rm cs}},\; $$

where $\epsilon _{\rm cs}$ and $\epsilon _{\rm ov}$ (set to 0.95) are emissivity during clear-sky and overcast conditions, respectively.

The model by Konzelmann and others (1994) is used to estimate $\epsilon _{\rm cs}$ as a function of temperature and humidity following:

(4)$$\epsilon_{\rm cs} = \epsilon_{\rm ad} + b {e_{\rm a}\over T_{\rm a}} ^{( 1/m) },\; $$

where $\epsilon _{\rm ad}$ is the emissivity for a clear dry atmosphere (–), e _a the atmospheric vapour pressure (Pa) and b and m are parameters to optimise the model. The model was optimised on the clear-sky dataset for $\epsilon _{\rm cs}$, which was found for $\epsilon _{\rm ad}$ = 0.23, b = 0.43 and m = 8.

The weighted daily average of the N _eff at TARM is then used to calculate $\epsilon _{\rm eff}$ during gaps in L _in observations by rearranging Eqn (3). L _in at either TARM or Lake Bonney can then be estimated using the clear-sky emissivity model optimised at Lake Bonney and observations of S _in, T _a, RH at TARM or Lake Bonney, respectively.

We evaluated the parameterisation of L _in using S _in, T _a and RH at Lake Bonney against observed L _in and found a RMSE of 25 W m⁻² and average bias of −3.40 W m⁻². During December and January, S _in is affected the least by shading, therefore errors in the transmission factor are smallest and RMSE and bias in the simulated L _in are reduced to 23 and 0.4 W m⁻², respectively. Between March and October, solar radiation is absent or too small to estimate N _eff. During these months L _in is estimated using a cloud fraction based on RH as described in Liston and Elder (2006) and previously used at TARM by Hoffman and others (2008). Especially during the melt season accurately representing L _in is critical for melt modelling. The method of Conway and others (2015) captures in this period the variability in observed L _in better compared to the model by Liston and Elder (2006), as evident from the correlation of 0.70 vs 0.53.

For the time periods in which observed L _in and S _in are available, N _ep and N _eff are used to describe the effect of clouds on L _in and S _in. The same procedure is applied on the COHM record to estimate N _ep and N _eff.

3.2 Föhn detection

Föhn winds play a significant role in shaping the climate of the McMurdo Dry Valleys. To identify the occurrence of föhn winds in the 22 year period we use the method of Speirs and others (2010), which was later adapted by Steinhoff and others (2014). With this method föhn is detected if the following conditions occur within a 6 h period: an increase in T _a (1$^\circ$C over 1 h, or any measurement above 0$^\circ$C), a decrease in RH (5% over 1 h, or any measurement below 30%), 15 min maximum wind speeds higher than 5.0 m s⁻¹ for at least 80% of the period and coming from directions between 180$^\circ$ and 315$^\circ$ (Steinhoff and others, 2014). Different methods of föhn detection using AWS observations and/or weather model forecasts have been tested at Joyce Glacier in the McMurdo Dry Valleys by Hofsteenge and others (2022). In the current study the Speirs method was chosen because of the availability of AWS observations over 22 years and the ability of the method to distinguish föhn winds from other winds coming from the expected föhn wind direction (e.g. a nighttime down-glacier wind at Taylor Glacier, Fig. 1b). As shown by Hofsteenge and others (2022), the method detects the onset of föhn winds well but might underestimate its duration.

3.3 Statistical analysis

Pearson's correlations are estimated to test linear relationships between various meteorological variables and SEB terms. To test for monotonic temporal trends in the observations and SEB terms, we use the Mann–Kendall trend test. Lastly, SEB model performance is tested using RMSE and bias.

3.4 SEB model

We use the SEB model by Reijmer and Oerlemans (2002) that was further developed by Jakobs and others (2019) and previously implemented in the McMurdo Dry Valleys at Joyce Glacier (Hofsteenge and others, 2022). The model calculates surface melt by closing the surface energy balance described as (defined with positive fluxes towards the surface):

(5)$$Q_{\rm M} = S_{\rm net} + Q_{\rm P} + L_{\rm net} + SH + LH + Q_{\rm G},\; $$

where Q _M is the energy available for melt, which is 0 when the surface temperature T _s < 273.15 K. S _net is the net shortwave radiation, calculated from observed S _out and daily accumulated albedo. Q _P is the penetrating shortwave radiation, calculated as in Brandt and Warren (1993) and as used in van den Broeke and others (2008) and Kuipers Munneke and others (2009). L _net is the net longwave radiation from observed L _in and calculated L _out from the modelled surface temperature and assuming that surface emissivity equals 1. To estimate the turbulent sensible heat (SH) and latent heat flux (LH), the bulk method is used in which the vertical gradients in wind speed, temperature and humidity are estimated between the AWS measurement level and the surface. For this method Monin–Obukhov similarity is assumed, for which a value is needed for the roughness length of momentum (z _0m), heat (z _0h) and humidity (z _0q). z _0m is kept constant as will be further discussed in Section 3.5, and z _0h and z _0q are calculated using the method described in Andreas (1987). The SEB model is coupled with the SOMARS snow model (Reijmer and Hock, 2008) to model the vertical profile of snow/ice temperatures from which the subsurface conductive energy flux (Q _G) is calculated. Subsurface melt occurs when layers with internal temperatures reach 0$^\circ$C, which is possible in the presence of shortwave radiation penetration. Runoff occurs when subsurface meltwater percolates down to the snow–ice interface or when subsurface melt occurs in the ice. Refreezing can occur in snow, which increases the temperature and density.

In this study we define the melt season as December and January, which are the most relevant months from a hydrological perspective, although melt can occur in November and February (Obryk and others, 2020).

3.5 Model parameter sensitivity

The modelled ablation of glaciers in the McMurdo Dry Valleys is found to be most sensitive to parameters governing the partitioning of solar radiation penetration and z _0m that is required for the estimation of the turbulent heat fluxes (Hoffman and others, 2014; Hofsteenge and others, 2022). These parameters have been used to optimise the SEB model for surface ablation (MacDonell and others, 2013b; Hoffman and others, 2014) or surface temperature under non-melt conditions (Hofsteenge and others, 2022). In Figure 2 we present a sensitivity study for 2013–18 for runs with different combinations of values of z _0m for ice and the skin layer thickness (z _rad) that partitions solar radiation between the surface and penetrating the subsurface (in which a smaller value results in more solar radiation partitioned to the subsurface).

Figure 2. Mean surface temperature bias and mean monthly surface ablation bias (mod-obs) during the melt-season for different values of z _rad and z _0m for ice. To calculate the ablation bias, months with net accumulation in summer are ignored. ‘no-pen’ stands for runs where no solar radiation was allowed to penetrate into the subsurface. The red circles show the parameter combinations with minimum absolute bias.

We find on the one hand model runs with less solar radiation penetrating into the subsurface represent observed surface temperatures better (Fig. 2). On the other hand, some solar penetration is needed to accurately model the surface ablation, which is in agreement with Hoffman and others (2014). An interesting difference is that TARM is more sensitive to changes in z _0m, visible from the vertical gradients in Figure 2 for TARM, while COHM is more sensitive to changes in z _rad (horizontal gradients). This likely reflects the more sublimation-dominated mass balance at TARM, which will be further discussed in Section 4. At TARM, a larger z _0m improves the surface temperature but increases the bias in ablation. Such a trade-off between optimising surface temperatures or ablation has previously been found with SEB models over snow (Conway and others, 2018).

Instead of optimising z _0m, we use a z _0m of 1e⁻³ m for ice and 1e⁻⁴ m for snow at both sites. This choice of z _0m for ice is on the same order of magnitude as estimated from two-level wind measurements at Taylor Glacier by Bliss and others (2011) and from unpublished eddy covariance data on Commonwealth Glacier. A smaller value is used for z _0m when the surface is snow covered, as snow is expected to fill in some of the surface roughness features. The value chosen agrees with z _0m estimates over Antarctic snow (van den Broeke and others, 2005). For z _rad a value of 0.005 m is used as in Hofsteenge and others (2022). The parameters are kept constant between the two sites to eliminate differences in model output between the sites depending on parameter choice. The sensitivities shown in Figure 2 reflect the uncertainty in modelling the absolute and relative contribution of melt and sublimation to the mass balance. However, we found that the temporal variability in melt, which is the focus of this study, did not depend strongly on the parameter choices. Section S2 shows that the drivers of melt found in this study are not significantly impacted by the parameter choices.

4. Results

4.1 Model performance

The modelled and observed surface heights at COHM and TARM are shown in Figure 3. Over the 22 year period there was 5 m of ablation at TARM compared to 2 m at COHM. The modelled surface height record captures the net surface height change over the period well. At COHM some deviation between observed and modelled is visible, for example around the summer of 2001/02 where strong ablation is observed. However, the observed ablation is likely overestimated as the stake was found in the bottom of a meltwater channel. It is evident that TARM experiences less snowfall and no accumulation, while COHM shows some accumulation events. Surface temperature performance is given in Figure S1.

Figure 3. Modelled vs observed surface height at TARM (a) and COHM (b). The observed record is based on surface height derived from sonic ranging observations and stake measurements.

In the context of understanding the drivers of meltwater generation, we also evaluate the modelled runoff (from surface + internal melt). Although observed streamflow at gauges is the total runoff of an on-glacier catchment and modelled runoff is expressed in mm w.e. for a point location, comparing the two can provide useful information about whether inter-annual variability and timing of runoff are captured by the model. Taylor Glacier terminates in Lake Bonney and stream gauges that are solely impacted by runoff from the glacier are missing. However, Commonwealth Glacier is drained by several streams, of which Lost Seal and Commonwealth streams are monitored with stream gauges by the McMurdo Dry Valleys LTER programme (Gooseff and McKnight, 2021, 2023). More information on the stream gauge instruments and measurements is given by Gooseff and others (2022). We justify comparing point-based runoff estimates from the SEB model at COHM with observed streamflow as follows: (1) glaciers in the McMurdo Dry Valleys are typically frozen to the bed and therefore no meltwater is generated or transported below the glacier, (2) runoff of the glaciers is generated in the top 25 cm of ice and the lack of cracks and crevasses on Commonwealth Glacier mean that little meltwater is stored internally in the glacier but rather runs off close to the snow–ice boundary until it leaves the glacier as waterfalls from the terminating cliffs (Bergstrom and others, 2021), (3) snowfall that falls off-glacier is very small and is assumed to not contribute to streamflow and (4) there are no deep aquifers that connect with the streams (Gooseff and others, 2011).

The model captures the inter-annual variability in runoff well (Fig. 4a), especially in representing years of low and high streamflow. Not only is the inter-annual variability captured well, but also the timing of the streamflow as can be seen for the 2001/02 ‘flood’ year (Fig. 4b). We acknowledge that lower parts of the glacier, especially the terminating cliffs, likely experience more melt compared to the AWS location. Still, assuming that our point-based runoff estimation is representative for the entire Lost Seal watershed, and using the indicated watershed size of ~2 km² (Fig. 4c), we find a streamflow estimate which is of the same order of magnitude as the observations. This observed agreement supports our decision to use the model to further study the meteorological drivers of melt.

Figure 4. (a) Seasonal modelled runoff at COHM in mm w.e. and streamflow measured at Lost Seal and Commonwealth stream in m⁻³. (b) Daily total runoff for the 2001/02 ‘flood’ season. (c) Streams, discharge gauges and glacier contributing areas based on Bergstrom and others (2021). Note that the upper boundary of the watersheds are cut off by the DEM extent that was used in Bergstrom and others (2021) and that higher elevations may still contribute to the runoff.

4.2 Local climates and SEB

Despite being in the same valley and separated by only 30  km, the local climates of TARM and COHM differ through their physical setting and local wind regimes (Table 1). The close proximity of COHM to the coast means that it receives a relatively large amount of precipitation for the McMurdo Dry Valleys (annual snowfall of 126 mm w.e. at COHM vs 22 mm w.e. at TARM). The typically drier, windier and warmer conditions at TARM are linked to the prevailing westerly wind direction (see wind-rose in Fig. 1b) and occurrence of föhn wind events (Speirs and others, 2013; Steinhoff and others, 2014). The mean air temperature during the melt season is below zero for both glaciers, with TARM 1.4$^\circ$C warmer (Table 1). However, föhn winds in the McMurdo Dry Valleys can raise air temperatures above zero and trigger melt (Doran and others, 2008; Hofsteenge and others, 2022). Föhn winds typically occur more frequently in the up-valley western side of the McMurdo Dry Valleys and last longer, whereas the down-valley eastern side often experiences a sea/ice-breeze during daytime in summer (see easterly winds in wind-rose Fig. 1), which only infrequently gets interrupted by strong föhn winds. This difference in föhn occurrence explains the difference in degree-days above freezing (DDAF) per year, which are on average 5.6 at TARM compared to 1.7 at COHM.

Table 1. Mean SEB terms and atmospheric variables at TARM and COHM over the full year and melt-seasons (December–January) between 2000 and 2021

Since the wind regimes introduce the largest differences in the local climates of the two glaciers, we expect the turbulent heat fluxes to be the most dissimilar of the energy-balance components. SH is on average positive at both sites, but the sensible heating is stronger at TARM. LH is on average more than double the size at TARM, reflecting increased sublimation rates due to the drier and windier conditions. During the melt season, S _net and L _net are the two energy-balance components with the largest absolute magnitudes at both sites (Table 1). S _net is the main energy source and is on average higher at TARM because of a lower albedo, however, the COHM albedo is more variable as will be discussed in Sections 4.3 and 4.5. Solar radiation that penetrates into the subsurface (Q _P) warms the subsurface and delivers energy back to the surface through the conductive heat flux (Q _G). The most important energy loss is longwave radiative cooling, while LH becomes stronger in summer as sublimation increases, with the strongest sublimation found at TARM. The higher air temperatures in summer at TARM result in a net positive SH at TARM compared to a negative SH at COHM (Table 1). Both sites show the typical diurnal cycle in surface layer stability (not shown), with a positive SH at nighttime and a negative SH during daytime that was previously found at another glacier in the McMurdo Dry Valleys (Hofsteenge and others, 2022).

Typical for the low-melt environment of the McMurdo Dry Valleys, we find that Q _M, the energy available for melt, is the smallest SEB term at both sites and on average larger at TARM. The extra energy input from a lower albedo and a positive SH at TARM in summer is mostly used for sublimation (LH in Table 1), but also explains the higher Q _M. TARM experiences therefore more surface melt and sublimation, and with less snowfall compared to COHM, the overall surface ablation over the period is more than double that of COHM (Fig. 3).

4.3 Drivers of inter-annual melt variability

The average Q _M is lower at COHM, but it has a larger inter-annual variability compared to TARM (Figs 5a, b). The time-series at COHM shows a few years with peak melt, but also years that experience almost no melt (Figs 5b and 4a). In contrast, TARM experiences some melt almost each season and has less distinct extreme years (Fig. 5a). The fact that melt at COHM seems to switch on and off is because COHM receives more snowfall, which can still be present after the winter, keeping the albedo of the surface high and limiting melt (Fig. 5f).

Figure 5. Mean melt-season (December–January) energy available for melt (Q _M) and DDAF, föhn occurrence, albedo and melt-season minimum albedo at TARM (a, c, e, g) and COHM (b, d, f, h). Significant trends are shown when p-value <0.05.

At both COHM and TARM, melt peaked in the summer of 2001/02, which became known as the ‘flood year’ in the McMurdo Dry Valleys and had long-lasting effects on the clarity, chemistry and biological activity of the lakes (Fountain and others, 2016). This extreme melt event was linked by Doran and others (2008) to an extreme increase in DDAF driven by anomalous föhn wind activity. We find for both TARM and COHM a significant correlation between the yearly DDAF and Q _M, with the strongest and most significant relation at TARM (Table 2). While the high DDAF of 2001/02 were an important driver of extreme melt, it is also evident from Figures 5a and b that other high melt seasons cannot be explained by high DDAF alone.

Table 2. Pearson's correlation coefficients (R) between melt-season (December–January) mean energy for melt (Q _M) and SEB terms and atmospheric variables

Correlations with p-value ≤0.05 are shown in bold and correlations with p-value ≤0.01 are shown underlined.

At COHM we find a significant correlation of 0.76 between energy for melt and S _net (Table 2). S _in exhibits at both sites an inter-annual variability of 14 W m⁻² through variability in cloud cover (R of −0.74 at TARM and −0.87 COHM with seasonal average N _ep, both significant with p-value <0.001). Still, we do not find a significant correlation between seasonal S _in and energy for melt. The melt variability through S _net is driven by albedo variability (R = −0.8, p-value <0.01) instead. The stronger correlation between melt and albedo compared to DDAF indicates that variability in albedo explains more of the inter-annual variability in melt at COHM than positive air temperatures. We look further into albedo variability in Section 4.5. A strong correlation is also found between seasonal Q _M and LH (Table 2), showing that seasons with strong melt and strong sublimation co-occur. This enhanced melt and sublimation occurs both in summers with low albedo (R of 0.92 between LH and albedo) which increases surface temperatures and therefore melt and sublimation rates.

In contrast to COHM, no significant correlation is found between seasonal S _net and melt at TARM. At TARM high melt years are correlated with an increase in air temperature, DDAF and föhn hours. We therefore find contrasting drivers of melt at the two sites, with albedo and S _net driving most of the inter-annual variability in melt at COHM, while föhn winds that increase air temperatures to above zero control melt of the surface at TARM. Table S4 shows the same correlation analysis used in Table 2 but using three other model parameter combinations.

4.4 Trends in energy available for melt

A significant increase in Q _M was found at COHM over the 22 years (Fig. 5b). The trend in surface melt is significant with p-value <0.05 and independent of parameter choice (Table S3). A similar but insignificant increase was also found in internal melt and total melt per season (Table S3). We do not find any significant trend in average albedo and DDAF per season at COHM that can explain the trend in Q _M. While mean summer albedo does not exhibit any trend, we find a significant trend in the minimum summer albedo over the 22 year period (Fig. 5h) of −0.007 a⁻¹, which is equivalent to a decrease of ~0.15 over 22 years. This trend likely explains the significant increase in Q _M over the period.

No significant trends were found in Q _M or melt components at TARM. We find a significant increase of 0.05 in average summer albedo and 0.001 in summer minimum albedo over the 22 years at TARM (Fig. 5a), however, this does not result in a significant trend in Q _M.

4.5 Albedo variability over the 22 year period

Albedo is an important driver of melt variability at COHM and to a lesser extent at TARM. The observed broadband albedo at both sites throughout the full study period (Fig. 6) shows that albedo variability is much stronger at COHM, where it varies between 0.4 and 0.9. COHM more frequently reaches albedo values up to 0.8, showing the glacier surface is more often still snow covered after winter or snow covered after a summer snowfall event. TARM is snow-free for most of the time, apart from short events in which the albedo rises up to 0.8 until it lowers again as snow is removed through sublimation or drifting snow erosion. The average summer albedo at TARM shows a significant increase over the study period (slope 0.003y ⁻¹ and p-value <0.01, Fig. 5e), which is likely linked to an increase in summer precipitation (slope 0.25 mm w.e. y ⁻¹ and p-value <0.001), although we are cautious to over-interpret trends in precipitation as two different precipitation measurements were used (Section 3.1.2).

Figure 6. Daily accumulated broadband albedo at TARM (a) and COHM (b) over all summers in the study period. Black line shows the 7 d moving mean albedo over the full period. The year indicates the starting year of each summer.

TARM has a distinct minimum albedo ~0.5 that shows little variation between individual years (Figs 6a and 5g). The ice albedo of TARM lies between 0.5 and 0.65 and is similar to Antarctic blue ice (Bintanja and van den Broeke, 1995; Reijmer and others, 2001).

COHM shows a much larger variability in ice albedo, with some seasons reaching values down to 0.4. The large variability in ice albedo at COHM is likely explained by the more snowmelt-dominated mass balance, in which snow and ice melt along with refreezing of meltwater can lower the ice albedo. While no significant trend is found in summer albedo at COHM, we find a significant decrease in the summer minimum albedo (Fig. 5h). This decreasing trend in minimum albedo suggests a lowering of the ice albedo at COHM. In addition to the impact of melt on the ice albedo, COHM is also more exposed to bare ground in the upwind direction during a föhn event (westerly sector of climatological wind rose in Fig. 1b). Deposition of aeolian sediments can lower ice albedo (Bergstrom and others, 2020) and is expected to be strongest in winter during strong föhn events (Gillies and others, 2013). The summer seasons between 2018 and 2021 that are characterised by some of the lowest ice albedos (Fig. 6) had preceding winters with relatively high föhn wind activity (Fig. S2). However, no significant correlations were found between föhn wind occurrence and minimum albedo to confirm that wind-blown sediments are responsible for this observed albedo trend. The albedo observations at the AWS site are also biased towards the cleaner portions of the glacier compared to lower elevations where sediments have more impact on the ice albedo (Bergstrom and others, 2020).

Interestingly, both sites show on average (black line in Fig. 6) an increase in the albedo between December and January. We propose two possible causes of this increase in albedo. Firstly, it may result from an increase in snowfall events during this period. For example, during January at TARM more events are observed in which the albedo for short periods approaches the albedo of snow (~0.8) (Fig. 6a). Secondly, the albedo increase could be attributed to increase in scattering effects between December and January that increase the ice-albedo. This increase in ice-albedo was previously found by Bliss (2011) and Hoffman and others (2014) at Taylor Glacier and by Bergstrom and others (2020) on several glaciers in the Taylor Valley. They attributed this increase to the development of a porous layer in the near-surface due to internal melt and meltwater drainage, which increases the optical scattering and hence albedo.

4.6 Melt regimes

Figure 7 shows the mean SEB components during hours with or without surface melt, and under positive or negative air temperatures. For most of December and January, the surface is not melting at both COHM and TARM (96 and 93% of the time respectively). Under these conditions the air temperatures are too low or albedo or cloud cover too high to bring the ice surface to the melting point (Table 3). Although TARM experiences more positive air temperatures, 66% of the time there is no surface melt despite air temperatures being above freezing. This happens with dry and windy conditions, often during föhn winds (Table 3), enhancing sublimation and preventing the surface from melting. COHM similarly does not experience melt for the majority of the time air temperatures are positive.

Figure 7. Mean SEB components during the melt-season at time steps with or without surface melt and below or above freezing air temperatures at TARM (a) and COHM (b). Along the top are the percentages given of the time during the melt-season that each condition occurred. SEB components are introduced in Eqn (5).

Table 3. Mean variables during hours with or without surface melt and below or above freezing air temperatures at TARM and COHM

The surface melt that does occur at COHM happens most often with negative air temperatures, when S _net is large (Fig. 7b) because of a low ice albedo and/or little cloud cover (Table 3). Melt with positive air temperatures occurs less frequently, but the energy for melt is in those cases on average higher. Melt under these conditions is still driven by S _net, but SH adds an additional source of energy. In contrast, melt occurs more often with positive air temperatures at TARM. Some melt occurs with negative air temperatures when S _net is large and wind speeds low.

In conclusion, at both sites melt only occurs occasionally during the summer. Melt at COHM is solar radiation driven, where albedo or cloud cover are the limiting factors. At TARM solar radiation is still the largest input of energy during melt, but at this site the amount of sublimation is the limiting factor for surface melt. This makes TARM comparable to other glaciers under (semi-) arid conditions, where windiness and dryness limit the surface melt (Hoffman and others, 2008, 2014; MacDonell and others, 2013).

5. Discussion

5.1 Timescales of melt variability and drivers of longer term trends

This study highlights the meteorological drivers of intra-seasonal and inter-annual variability between 1999 and 2022 and is the longest glacial SEB record to date for the McMurdo Dry Valleys. On intra-seasonal timescales we find that melt can switch on and off within the summer season. High air temperatures and clear-sky conditions favour melt at both sites. At COHM, summer snowfall events increase albedo and switch off the melt and streamflow. Snowfall is less common at TARM and variability in S _in, temperature and wind speed explain most of the intra-seasonal variability, where strong wind speeds increase sublimation and limit melt.

Within the season, cloudy days can interrupt melt, but when assessing the inter-annual variability we do not find a statistical link between cloud factors and melt variability. This is surprising given the strong statistical link between variability in longwave cloudiness (N _ep) and S _in. Clouds have previously been eliminated as a driver of inter-annual variability in S _in in the McMurdo Dry Valleys by Obryk and others (2018), who instead found changes in atmospheric sulphur dioxide concentrations from volcanic activity and anthropogenic emissions as the main driver. However, Obryk and others (2018) used annual cloud cover over McMurdo Station and it is expected that large differences in cloud cover occur between Taylor Valley and McMurdo Station as the two locations are almost 100 km apart and experience very different local weather conditions. Furthermore, statistical links with annual cloud cover are expected to be weak as only summer cloud cover impacts the solar radiation variability. We find at both sites that inter-annual variability in annual S _in is explained by variability in longwave equivalent cloudiness (independent of aerosol depth and scattering processes described by Obryk and others (2018)) in the months with the highest solar radiation (November–February), supported by the correlation of −0.65 (p-value 0.001) at COHM and −0.56 (p-value 0.01) at TARM. We therefore suggest that variability in cloud cover is the main driver of year-to-year variability in S _in. Longer term trends beyond the year-to-year variability, such as the increase in S _in between 1991 and 2000 described by Obryk and others (2018), might still be attributed to a decrease in global sulphur dioxide emissions, but more studies are needed to confirm this.

The study period between 1999 and 2021 was chosen based on data availability for SEB modelling, but how representative is the inter-annual variability of this period to longer term changes in streamflow? Previous studies on the hydrology in Taylor Valley have shown a cooling period between 1969 and 2000 and a decreasing trend in lake levels in the McMurdo Dry Valleys, suggesting lower levels of glacial melt (Wlostowski and others, 2016). The period following the high melt summer of 2001/02 shows on average higher streamflow and larger inter-annual variability in which lake levels rise (Cross and others, 2022). The extended SEB record to 2021 presented here confirms that this period of enhanced variability is still ongoing.

The high melt season of 2001/02 that formed the end of a cooling period in the McMurdo Dry Valleys has been described as a catalyst for change in previous literature Fountain and others, 2016a; Levy and others, 2018; Bergstrom and others, 2021), with different mechanisms being identified to explain the transition. Bergstrom and others (2021) compared roughness estimates from two lidar elevation surveys and found increased roughness at the lower ablation zones of Commonwealth and Canada glaciers after the 2001–02 summer and hypothesised that this roughening led to a transition to more melt-dominated glacial ablation. Using the same lidar surveys, but a different approach to calculate surface roughness, Levy and others (2018) instead found an overall glacier smoothing and suggested this was caused by extensive surface melting during the 2001–02 summer, removing roughness features. Fountain and others (2016a) attributed the increased melt after the 2001–02 summer to exposure of glacial sediments after the extensive surface melt. The latter agrees with the current study where albedo is the main driver of variability at COHM. The significant trend in summer minimum albedo over the 22 year period shows a similar feedback to Fountain and others (2016a), which is a lowering of surface albedo leads to enhanced melt.

We find that air temperatures, snowfall and albedo are most important in explaining temporal variability in melt. These drivers account for the inter-annual and intra-seasonal variability in observed runoff in this study, while parameters reflecting glacial properties such as the roughness length have been kept constant. The main drivers we find might be valid for the AWS locations which are relatively flat, however, at lower elevations of the ablation zones of the glaciers, feedbacks stemming from localised climatic effects on surface roughness as discussed by Bergstrom and others (2021) might be more important as will be discussed in Section 5.3.

5.2 The future fate of glaciers with contrasting sensitivity to climate forcing

At both TARM and COHM we find that seasonal melt is correlated significantly with (positive) summer air temperatures, suggesting that small changes in air temperatures can have large impacts on melt. At both sites the occurrence of DDAF are linked with föhn wind activity. Föhn winds in Taylor Valley occur typically with low-pressure over the Ross Sea, which creates a southerly synoptic flow towards the McMurdo Dry Valleys. This air flow is then forced through a mountain gap, accelerated, warmed and channelled into Taylor Valley (Steinhoff and others, 2013). TARM experiences more föhn winds and more frequent positive air temperatures compared to COHM, but also most of the melt occurs with positive air temperatures. Many of the föhn winds do not fully penetrate through the Taylor Valley to impact the down-valley glaciers such as Commonwealth Glacier. To understand how glacial melt in Taylor Valley and the wider region might change in the future, more research is needed to predict changes in the synoptic weather circulation patterns in the Ross Sea and the mesoscale föhn wind activity in the McMurdo Dry Valleys.

We hypothesise that changes in large-scale and synoptic-scale atmospheric circulation will impact the glaciers in different ways: Taylor Glacier is more sensitive to changes in föhn wind activity, while Commonwealth Glacier is more sensitive to changes in snowfall. Moisture input with snowfall is more dominant in the eastern McMurdo Dry Valleys, while snowfall rarely reaches Taylor Glacier. The 22 year albedo record at Commonwealth shows strong variability in albedo and occurrence of snowfall throughout the summer. Both sites show on average an increase in albedo between December and January linked to snowfall events, which might be linked to the break-up of fast ice in the McMurdo Sound that typically occurs at this time of the year (Kim and others, 2018). This may lead to an increase in evaporation sources and precipitation amounts. Future research should therefore study the synoptic drivers of snowfall but also it's moisture sources, to assess how snowfall might change with ocean warming and changing sea ice conditions. Our findings provide evidence to suggest that an increase in snowfall in the McMurdo Dry Valleys could lead to a transition into a lower-melt regime. However, the physical processes controlling longer term changes in snowfall need to be investigated further to determine whether a negative feedback could occur on the eastern glaciers of the McMurdo Dry Valleys.

5.3 Limitations of point-based SEB modelling and uncertainty in ablation contributions

To fully understand the temporal and spatial variations in melt of glaciers in the McMurdo Dry Valleys, it is essential to quantify temporal variability and horizontal gradients in mass and energy balance through high-quality observational data and SEB modelling. The two long-term records of AWS observations at COHM and TARM used in this study therefore have value for identifying trends in forcings and dominant mechanisms driving the glacier SEB at two contrasting sites. Working with SEB models at AWS locations does also have limitations. The glacial AWSs in the McMurdo Dry Valleys are generally located in flatter areas of the glaciers, which may not represent the regions of most significant melt. Previous studies have shown increased melt rates along vertical cliffs (Lewis and others, 1999), in surface channels in the lower ablation zone (Johnston and others, 2005) or debris-covered areas (MacDonell and others, 2013b). Hoffman and others (2016) did an experimental modelling study in which micro-climates in topographic basins were represented in a distributed SEB model by increasing air temperature and lowering wind speeds and albedo compared to AWS observations. These experiments showed that parts of the ablation zone with topographic basins can generate melt rates that are ten times higher than flat glacier surfaces. Therefore, the drivers of melt that are found in this study might not be valid for all parts of the glacier surface. Nonetheless, the agreement in inter-annual variability between observed streamflow from large parts of the ablation zone and modelled runoff at the single AWS point at Commonwealth Glacier (Fig. 4) gives us confidence our results are representative for large parts of the ablation zone.

The parameter sensitivity analysis of Section 3.5 reflects the uncertainty in modelling the contribution of sublimation and subsurface and surface melt to the glacial mass balance. This uncertainty is also reflected in the literature. Considering the full extent of Taylor Glacier, only melt is expected to occur in the lowest part of the glacier tongue in Taylor Valley, and melt was therefore neglected in studies focusing on sublimation of the full ablation zone of Taylor Glacier (Bliss and others, 2011). However, other studies have shown the importance of surface melt to explain the channels on top of the lower part of Taylor Glacier (Johnston and others, 2005). Hoffman and others (2014) instead did a modelling study in which almost all melt occurred in the sub-surface and rarely as surface melt. We find contrasting sensitivities of the two studied glaciers, where surface ablation at COHM was more sensitive to changes in the parameter governing solar radiation penetration, while the sublimation-dominated TARM was more sensitive to the roughness length. There is still a need for direct observations of the turbulent heat fluxes over glaciers in the McMurdo Dry Valleys to better quantify the contributions of sublimation vs melt to ablation.

6. Conclusions

In conclusion, this study has advanced our understanding of the meteorological drivers of melt in the McMurdo Dry Valleys, by examining two contrasting glacier sites over a 22 year period. We show how differences in their local climates determine the drivers of melt variability, suggesting their response to regional atmospheric circulation changes will also be different.

On intra-seasonal timescales, melt happens at COHM most often with sub-zero air temperatures when large S _net can bring the surface to melting point under clear-sky and low albedo conditions. At TARM, melt happens more often with positive air temperatures. Still, for the majority of time when air temperatures are positive, melt is limited at TARM because sublimation is more dominant due to dry and windy conditions, often during föhn wind events.

At both sites, the seasonal amount of positive degree-days, which are strongly linked to föhn wind activity, are an important driver of inter-annual melt variability. At COHM however, S _net and albedo are the main drivers of inter-annual melt variability. Albedo variability is larger at COHM, which experiences more frequent snowfall, but also has a stronger varying ice albedo. We find a significant increase in surface melt over the 22 years at COHM, which correlates with a significant decrease in summer minimum albedo. No significant trend in melt was found at TARM.

Given these findings, we hypothesise that changes in atmospheric circulation might impact glaciers within the McMurdo Dry Valleys in different ways: the further inland glaciers such as Taylor Glacier being more influenced by föhn wind activity and the more coastal glaciers such as Commonwealth Glacier more responsive to changes in snowfall.