A customized automatic weather station (AWS) was developed for use in this project. The station housed a suite of meteorological sensors (Table 1), including a four component radiometer to measure each of the radiation fluxes, and a tipping bucket rain gauge to observe precipitation (and calculate Q_R). The temperature profile of the surface layer of ice was measured for the calculation of Q_G using two thermistor arrays drilled into an initial depth of 4 m, with thermistor spacing of 0.25 m in the upper 1 m, and 1 m in the lower 3 m. A sonic ranger was used to observe changes in surface height due to ablation or snowfall. To directly measure the turbulent heat fluxes, the EC method was employed. A combined sensor, comprising a sonic anemometer and an open path infrared gas analyser was selected for this campaign. With this design, the anemometer and analyser components instantaneously measure the same sample volume, minimizing temporal and spatial errors due to the physical separation of the sensors. The sampling frequency of the EC measurements was 20 Hz, with the values averaged over 30 min blocks when calculating the magnitude and direction of the turbulent fluxes.

In the atmospheric surface layer, significant gradients with height of temperature, humidity, and wind speed are typical. Therefore, it is important that meteorological sensors are maintained at a constant height above the surface if comparable measurements over time are required. The EC sensors are particularly sensitive to height changes, as increasing the height of the EC sensor may change the footprint from which the turbulence observed by the sensor is received (Burba, 2013). The sensor will also be at an increased risk of being positioned above the surface layer during strongly stable conditions, meaning that the sensor would no longer be observing eddies with characteristics of the glacier surface fluxes (Aubinet, 2008). To avoid these errors, the main meteorological sensors (including the EC system) were installed on a custom-built four-legged “quadpod,” which sits on the surface of the glacier (rather than being drilled in place), and lowers with the surface during the melt season (Figure 1). This ensured that the sensors maintained a constant height above the surface (wind monitor height z_u: 2.7 m, EC and temperature/humidity sensors height z: 2.0 m; radiometer height: 1.6 m). In order to monitor and correct for any tilting of the station that may occur over the season, an inclinometer was fitted to the quadpod to record its tilt angle in the x and y planes. The sonic ranger, along with the rain gauge, was mounted on a separate mast that was drilled into the ice, so that relative changes in glacier surface height could be measured. Observations were saved in 1 min averages (1 min totals for rainfall; 30 min samples for ice thermistors), while the eddy covariance data was stored in raw (20 Hz) format. A time lapse camera was also installed onsite, recording eight images per day of the surface and atmospheric conditions, and the general position and order of the AWS.

The described AWS was deployed during the summer of 2014 on Nordic Glacier (51°26′N, 117°42′W) in the Purcell Mountains in eastern British Columbia (BC), Canada. Nordic is a small (~5 km²), north-facing glacier, with an elevation range from 2,000 to 2,900 m a.s.l. (meters above sea level), approximately (Figure 1). The glacier is within the basin of the Columbia River; an important aquatic environment, and water source for irrigation and hydroelectric power production in BC, and the United States (Columbia Basin Trust, 2017). Previous studies within this basin indicate that glacier runoff may contribute up to a third of total streamflow during summer months (Jost et al., 2012). The AWS was installed in the ablation zone of the glacier (2,200 m a.s.l.) from July 11th to August 28th 2014, with the first and last days of measurements excluded from the study, to avoid contamination of the data due to installation activity. The study period therefore runs for 47 days from 00:00 on July 12th, day of year (DOY) 193, to 23:30 on August 27th (DOY 239). The station site was selected in the area with the lowest surface slope angle (8°), and with the most uniform upwind footprint, to help minimize the corrections required when processing the radiation and EC data. The station was orientated so that the radiometer was located on the south side of the structure, in order to avoid shadowing of the sensor. The main axis of the EC instrument was positioned to point upslope into the direction of the prevailing wind (downslope flow), to minimize flow distortion due to air flow through the station structure. The AWS performed reliably over the study period, with the quadpod remaining stable throughout (tilting < 2°), the power system maintaining an adequate supply, and the data from all sensors (including the camera) being consistently logged without any evident systematic errors. One of the main concerns prior to deployment, was the reliability of an unattended EC system in a harsh environment for an extended period of time. Calibration of the sensor before and after deployment showed negligible drift in its measurements over the season, while the signal strength of the gas analyser (excluding rain events) decreased from 99 to 97%, staying well above the manufacturer recommended quality cut off point of 70% (Campbell Scientific Inc., 2013).

In the first stage of the model, the observed energy fluxes are converted into 30 min averages and summed to obtain averages of Q_M in Wm⁻² (Equation 1). The general convention in glacier SEB studies is that fluxes resulting in a gain in energy of the surface have a positive sign. The maximum temperature a snow/ice surface can be warmed to is 0°C. When the glacier surface is at this temperature, any positive energy available in the Q_M term is used by the model to generate melting (Hock, 2005):

where M is the average melt per second (for a given 30 min average of Q_M) in m of water equivalent (m w.e. s⁻¹), and L_f is the latent heat of fusion (0.334 MJ kg⁻¹). The unit of m w.e. expresses the thickness of the layer melted (per m² of area) when the snow/ice is converted to water. The surface temperature of the glacier is estimated in the model using the observed (corrected) longwave radiation data (Equation 9). Sublimation (in m w.e. s⁻¹) is also accounted for in the model. A negative/positive value for Q_L (i.e., a loss/gain of latent heat energy) indicates mass loss/gain through sublimation/resublimation and can be accounted for as follows:

where |Q_L| denotes the absolute value. The values for melt and sublimation, converted from averages to 30 min totals, are summed together to calculate surface ablation (a). The performance of the SEB model is evaluated by comparing the values it generates with the data observed by the sonic ranger. Modeled ablation is converted into surface height changes (z_gM) by multiplying the half hourly ablation totals by the ratio of the snow or ice density (ρ_s/i) to ρ_w:

A thin layer of packed snow covered the site surface for the first 3 days of observations, and a constant value of 400 kg m⁻³ was used for ρ_s. For the remaining period, the ice surface was assumed to have a constant density of 900 kg m⁻³ (Arnold et al., 2006; Gillet and Cullen, 2011; Van As, 2011). These surface height values can then be compared directly with the smoothed values from the sonic ranger. Accumulation due to snowfall and deposition is also added to the modeled surface height, so as to more accurately follow the observed values. During the single, brief period of snowfall (0.001 m w.e. approx.), a density of 80 kg m⁻³ was assumed for the fresh snow; a low density value to account for the thin and incomplete coverage of snowfall over the study site.

The performances of bulk method parameterizations of the turbulent heat fluxes (Equations 4 and 5) were evaluated through comparisons with the observed fluxes from the EC system. Three of the bulk transfer coefficient forms most commonly used in glacier studies have been tested here. The first form assumes neutral atmospheric stability (excludes a stability function), and a logarithmic wind speed profile (e.g., Conway and Cullen, 2013):

where k is the von Kármán constant (0.4). The second form includes a function for stability based on the bulk Richardson number, R_ib (Webb, 1970; Brock et al., 2010):

where

ϕ=(1-5Rib)2 for stable conditions (positive R_ib), and ϕ=(1-16Rib)0.75 for unstable conditions (negative R_ib), and g is acceleration due to gravity (9.81 ms⁻²). The third form uses stability functions based on the similarity theory of Monin and Obukhov (1954), and is arguably the most frequently used stability scheme in glacier SEB models (Arnold et al., 2006; Giesen et al., 2014):

where ψm(zuL) and ψh(zL) are the vertically integrated stability functions for momentum and heat, with the stability function for water vapor, ψe(zL), assumed equal to ψh(zL). L is the Obukhov length:

where u_* is the friction velocity (a function of shear or the vertical flux of horizontal momentum). Conceptually, L can be thought of as the height above the surface where buoyant forces equal shear forces in the production of turbulence (Stull, 1988). As L contains Q_H, an iterative scheme is used when solving Equations (4) and (22), following Munro (1989). The parameter zuL is used to determine stability conditions, with zuL < 0 implying unstable, zuL > 0 implying stable, and zuL = 0 implying neutral. In this study, the functions proposed by Beljaars and Holtslag (1991) have been used in Equation (21) for stable conditions:

with a = 1, b = 0.667, c = 5, and d = 0.35. For unstable conditions, the methods of Dyer (1974) have been implemented:

where y=(1-16zuL)14 for Equation (25), and y=(1-16zL)14 for Equation (26).

The meteorological conditions remained mild throughout the study period (Figure 3), with a mean T_z of 8.9°C and consistently positive air temperatures (max: 16.1°C; min: 0.6°C). Cloud coverage was generally low, as verified by camera imagery, resulting in high net shortwave radiation fluxes over the season (mean: 130.6 Wm⁻²) and low net longwave values (mean: −9.2 Wm⁻²). The surface temperature T_s, was estimated to be at melting point for >99% of the study period, with the temperature below 0°C for a total of 5 h (min: −0.3°C). Precipitation was confined to three frontal systems and sporadic convective showers (total: 293.7 mm), with a trace amount of snowfall on July 24th (DOY 205). The mean Q_R value was 1.7 Wm⁻² over the entire period, with a maximum 30 min average of 339.2 Wm⁻² during an individual event with a 34 mm/h rainfall rate (DOY 215). Winds were generally light and steady, with a mean of 3.5 m s⁻¹ and a maximum 30 min average of 7.2 m s⁻¹. The mean wind direction was 204°, which accords with the downslope air flow on the glacier. The wind remained within ±45° of this direction for 91% of the observation period, and within ±90°For 97% of the period, with no evident diurnal cycle. The positioning of the EC sensor within the prevaling wind window (main axis: 225°) ensured that, for the vast majority of observations, air flow was sampled by the EC system prior to being distorted by the station. Q_H remained positive for the entire observation period, with a mean value of 59.4 Wm⁻². Q_L fluctuated about zero, with a mean value of 9.0 Wm⁻². The ice thermistors recorded a small temperature gradient in the near-surface ice layer (~0.9°C over 4 m). As the thermistor cables gradually melted out over the season, the depth range of ice temperature measurements reduced. From the available data, little change is observed in the temperature profile over time, with a slight warming trend commencing a few days after the removal of the snow layer (DOY 196). With such a small temperature gradient, Q_G remained small, with a mean value of -0.7 Wm⁻².

Error analysis performed on each of the observed fluxes returned daily average error values substantially smaller than the fluxes themselves (Table 2). The obtained mean uncertainties for shortwave and longwave fluxes of 1.9 and 3.3% are in line with those estimated in previous studies using similar sensors (Van den Broeke et al., 2004). It should be noted that there remains unquantified uncertainty in S_i and Q_R due to the previously described empirical corrections applied to the data. The combined effect of the errors in each flux were propagated to obtain the error in Q_M (mean relative error: ±9.9%). The uncertainity in the surface height data due to the required smoothing was calculated to be ±0.001 m (±71.4%) and ±0.008 m (±13.8%) for mean 30 min and daily values, respectively.

The ablation values obtained from the SEB model show good agreement with those recorded by the sonic ranger (Figure 4). Both modeled and observed data returned a net surface ablation of 2.74 m over the 47 days. Examining model performance in estimating 30 min changes in surface height, the model showed a correlation (r) of 0.75 with the ranger data, and a root-mean-square error (RMSE) of 0.001 m/30 min. On a daily time scale, model performance was considerably improved (r: 0.86; RMSE: 0.02 m/day). Net radiation (Q_N) was the largest contributor to the melt energy, providing 65.2% over the study period (Figure 5A), with the longwave fluxes effecting a net loss of energy from the surface (−7.6% of Q_N), as shown in Figure 5B. Q_H was the second largest contributor with 29.7%, followed by Q_L with 4.5%. Due to the low magnitude and frequency of negative Q_L fluxes, little sublimation was observed during the season (0.1% of total ablation), while total resublimation (during positive Q_L periods) was estimated at 0.014 m w.e. Providing just 0.9% of the total melt energy, Q_R was observed to be the smallest positive contributor over the season, while Q_G was found to be a small surface energy sink (−0.4%). The day with the greatest observed ablation (in terms of m w.e) was July 29th (DOY 210), with 0.1 m w.e. of loss. This corresponds with the day of maximum net energy received at the surface, with a mean flux value of 288.8 Wm⁻² observed (55% above the average daily value for the study).

As previously described (Equations 27–29), the EC data obtained on Nordic Glacier were used to calculate a dataset of site-specific roughness length values (ec_z). Atmospheric conditions were predominately stable over the glacier (89% of study, as determined by zLec) with few periods of neutral stability. As a result, application of the stability filter (along with the other filters described in Section Roughness Lengths) returned 69, 56, and 11 30-min periods for z_0v, z_0t and z_0q, respectively, when the values were deemed representative of the study surface. The mean log values of the roughness lengths (± 1 standard deviation) for these periods were 10^{−2.3 ± 0.9} m, 10^{−4.6 ± 1.1} m, and 10^{−6.0 ± 0.8} m, respectively. In Figure 6, the ec_z values are compared with those estimated from the most commonly used roughness length schemes from existing glacier SEB studies: “equal” (z_0v = z_0t = z_0q); “1/100” (z_0v/100 = z_0t = z_0q); the Andreas (1987) surface renewal method for z_0t (mean: 10^−3.4 m) and z_0q (mean: 10^−3.2 m). For each of these three methods, z_0v is assumed to be 10⁻³ m for ice and 10⁻⁴ m for snow cover (as discussed in Section Roughness Lengths). Based on camera and albedo observations, the surface was assumed to transition from snow-covered to ice-covered after the first 4 days (DOY 197). Separating the EC-derived z_0v values into periods with a snow or ice surface produced mean values of 10^{−3.8 ± 1.7} m and 10^{−2.2 ± 0.7} m, respectively.

Three forms of the bulk transfer coefficient (Equations 18–26), and the previously discussed range of observed and assumed roughness length values, were implemented to obtain 12 sets of bulk-estimated turbulent heat fluxes (see Table 3 and Figure 7). Where the ec_z data were used, the means of the available roughness length values were implemented (separating z_0v for snow and ice). Evaluating the bulk estimates of the fluxes against the EC observations was only carried out for periods with ideal EC measurement conditions, namely those that passed the filters used when calculating the roughness lengths (Section Roughness Lengths), excluding the stability filter. Therefore, 691 and 380 periods (31 and 17% of study period) were useable for Q_H and Q_L comparisons, respectively. The flux values returned by each bulk method varied greatly depending on the paired roughness scheme. For the ec_z values, C_log and C_MO show significant correlation with the EC-observed fluxes, in particular, for Q_L (r = 0.89, p < 0.01). Weaker correlation exists between C_bR and the observed fluxes, particularly for Q_H, in addition to substantial underestimation of Q_H and Q_L (MBE: −29.1% and −25.5%). Conversely, C_log and C_MO, substantially overestimate both fluxes. For the range of assumed roughness schemes, C_log, when coupled with the 1/100 roughness scheme, produces flux values closest to those observed.

The strong performance of the model in replicating the recorded ablation rates gives a high level of confidence in the successful capture and closure of the SEB by the employed observation and correction methods. The trend of the modeled surface height over the 47 days (Figure 4A) shows no systematic bias to under or over-estimate when compared with the observed values, resulting in matching net values at the end of the study (−2.74 m). While good performance is also noted in the daily and 30 min ablation estimates, there are a number of pronounced periods where modeled values substantially differ from observations. Examining the daily timescale (Figure 4C), the days with the largest model differences are DOY 193 (+0.036 m), DOY 210 (−0.03 m), DOY 226 (+0.053 m), and DOY 234 (−0.034 m). In most cases, however, these differences can be attributed to observation error. Accidental compression of the snow underneath the sonic ranger during installation prior to DOY 193 resulted in a higher snow density and smaller observed height changes, while a review of camera imagery detected a clear tilting of the sonic ranger mast on both DOY 226 and 234, with the mast tilting forward and upslope on DOY 226 (sensor closer to the surface) and backward and downslope on DOY 234 (sensor further from the surface). DOY 210 was the day with the largest observed ablation and mean surface energy flux (288.8 Wm⁻²), but recorded 30% more ablation than DOY 197, which received just 1% less energy (285.5 Wm⁻²). Camera footage from DOY 210 may indicate the enhancement of a small meltwater stream within the footprint of the sonic ranger, which may have increased mass removal during this period, but the camera angle prevents definitive confirmation of this. Evaluating the SEB model performance with confirmed ranger error days removed (DOY 193, 226, 234) results in a marked improvement with r: 0.95 and RMSE: 0.01 m/day.

The SEB model offers good performance at 30 min resolution (r: 0.78; RMSE: 0.001 m/30 min, with aforementioned error days removed). The most prominent anomaly in the comparison between observed and modeled data at this timescale is the increase in surface height frequently observed by the ranger during night hours (negative values in Figure 4B). With only one period of trace snowfall during the study, accumulation does not account for these increases (mean increase of 0.003 m/night). A possible explanation for this behavior is the settling of the ranger mast into the warm surface ice; a signal that only becomes visible in the height data in the evening when surface ablation is reduced. With reduced heat and meltwater input into the borehole as the night progresses, this process slows, and at dawn, its signal is once again masked by ablation. This may also explain why the anomaly appears to be strongest in the later stages of the season (maximum increase of 0.03 m/night), when the base of the mast is closer to the surface, and more exposed to surface heat and under greater tilt leverage. An array of sonic rangers around the observation site would be a useful addition to future studies by providing additional data sources when sensor tilting or noise are obscuring the surface height signal from one. It would also help to constrain the likely range of ablation values at the glacier site where substantial variation in ablation can exist over small spatial scales.

In general, the partitioning of the observed SEB is in line with the values determined in existing studies on other mid-latitude glaciers (e.g., Hock and Holmgren, 1996; Greuell and Smeets, 2001). Q_N dominated melt energy, with its magnitude and variability primarily influenced by the net shortwave radiation (S_N), as shown in Figure 5B. Low cloud amount (reduced L_i fluxes) and a steadily melting surface (large magnitude L_o fluxes) resulted in a negative net longwave radiation flux (L_N), on average. The relatively large energy contribution provided by Q_H was maintained by the strong and consistently positive temperature gradients between the air and surface, and persistent downslope winds. Q_R provided a negligible contribution (<1%) to the total melt energy over the recorded period. However, over daily and sub-daily timescales, Q_R was observed to have a considerable influence on SEB and ablation during heavy rainfall. On July 19th (DOY 200), persistent heavy rain (daily total: 109.8 mm) resulted in a relatively large Q_R contribution of 20.1% to the daily energy flux (Figure 5A). There are previously mentioned uncertainties in the Q_R values arising from the observed rainfall rates, and the assumption made in Equation (13) that drop temperature is in equilibrium with air temperature. However, the general findings from this study show that, on daily and sub-daily timescales, Q_R has the potential to be a significant flux during extreme rain events; a phenomenon projected to become more frequent in a warmer climate (IPCC, 2013). Despite extensive cloud cover on DOY 200, daily ablation (0.058 m w.e.) was above the seasonal average (0.049 m w.e./day), with the non-radiative fluxes dominating Q_M (57.6%), highlighting the importance of accurate observation or estimation of non-radiative fluxes on glaciers (Fausto et al., 2016).

The initial step in evaluating the bulk method and its uncertainties was to compare roughness lengths calculated from the observed EC data with those commonly assumed for glacier surfaces. Despite extensive filtering, the scatter of the EC-derived roughness values was large over time (Figure 6) due to a dependence on a range of mean and turbulent variables in their calculation, as noted in previous studies (Andreas et al., 2010). The mean z_0v value from ec_z over the first 4 day period (10^−3.8 m) was found to be 1.6 orders of magnitude smaller than for the remaining 41 days (10^−2.2 m). The area upwind from the station (i.e., turbulent flux footprint area) was predominately snow covered during this initial period, and may have presented a smoother surface to the flow than the subsequent ice surface. This observed value for snow is similar to the assumed z_0v of 10⁻⁴ m, but there is close to an order of magnitude difference between the observed and assumed (10⁻³ m) values for the ice surface. A z_0v value of 10^−2.2 m is in the upper range of existing observations on glacier ice (Brock et al., 2006), and in line with measurements over rough, heavily gullied ice surfaces (Smeets et al., 1999). This relatively large roughness value may indicate the influence of the meltwater-enhanced channel network observed to develop on the surface of Nordic glacier during the study.

Rather than being equal to z_0v, the observed scalar roughness lengths were found to be several orders of magnitude smaller, as has been noted in previous studies (Beljaars and Holtslag, 1991; Smeets et al., 1998). The assumption that the scalar values are equal to each other was not observed at this site, with mean z_0q (10^{−6.0 ± 0.8} m) 1.4 orders of magnitude smaller than mean z_0t (10^{−4.6 ± 1.1} m). However, only 11 data points were available for z_0q after filtering, and there was substantial variability and overlapping of both sets of scalar values (Figure 8A). Substituting the ec_z values for z_0v into the tested roughness length schemes, the equal method overestimates both of the scalar roughness lengths, in some cases, by up to 4 orders of magnitude. The 1/100 method gives a reasonable mean value for z_0t while overestimating z_0q. The observed data for z_0t/z_0v show general agreement with the curve of the Andreas (1987) model (Figure 8B), but with a tendency for z_0t to be underestimated by the model (mean: 10^−5.3 m). The Andreas model was developed over snow and sea ice, which normally present a smoother surface to the flow than glacier ice. Smeets and van den Broeke (2008) proposed an adaption to this model for rougher ice surfaces (z_0v > 10⁻³ m). This scheme was also applied to the observed z_0v data, resulting in a calculated mean scalar roughness length value of 10^−3.2 m. Therefore, the observed z_0t values for this study lie mostly between both schemes. z_0q is overestimated by both methods, but shows best agreement with the Andreas model, with a modeled mean value of 10^−4.9 m (Figure 8C). Relations between mean meteorological variables (e.g., air temperature, relative humidity) and the scalar roughness lengths have been proposed in previous studies (Calanca, 2001; Park et al., 2010). However, such relationships were not evident in this data set. For the roughness length results in general, the limited number of data points available for analysis due to the persistent stable atmospheric conditions makes it difficult to draw definitive conclusions. An expanded study over more than one location, and with a greater number of measurements during neutral conditions, would be beneficial.

Assuming the ec_z values to be representative of the glacier surface, applying these roughness lengths to the three tested bulk transfer coefficients allows for an assessment of the performance of each stability scheme in returning flux values similar to those observed (Table 3). With no stability functions, C_log does not account for the suppression of turbulence that would theoretically be expected in the predominately stable surface layer at the site. The iterative scheme used to determine L for the C_MO method (Equation 22) was found to overestimate when compared with EC-derived values of L_ec (Figure 9B). As a result, stability is underestimated, and the suppression applied by the stability functions is insufficient. Similar overestimation of L has been observed in a recent comparison on a glacier surface by Radić et al. (2017). For daily mean fluxes (with only ideal EC measurement periods included), the performance of C_MO is significantly improved, with correlation values for Q_H and Q_L of 0.88 and 0.95, respectively (Figure 10). As mentioned, there was substantial observed variability around the mean ec_z values used in the bulk calculations. To test the sensitivity of the modeled fluxes to this variability, Monte Carlo analysis was applied to the utilized roughness lengths. For each 30 min C_MO flux calculation, 10,000 runs were performed where each of the three roughness lengths were selected at random from normally distributed populations with the means and standard deviations of the corresponding ec_z values. The standard deviation of the daily means from these runs are presented as error bars on the modeled flux values in Figure 10. Allowing for this sensitivity to roughness length input, periods of flux overestimation remain, being most pronounced during days with combined conditions of strong stability, large near-surface air temperature gradients, and relatively strong downslope winds, indicative of katabatic flow. Throughout the study period, the strengthening of the near-surface temperature gradient was observed to correspond with a strengthening of the downslope airflow that dominated the observed wind direction (Figure 9A). The possible development during very stable conditions of a low-level wind maximum near measurement height may have decoupled turbulence generation from the surface, and invalidated the assumptions of Monin-Obukhov similarity theory, such as a constant flux layer related to near surface gradients (Van der Avoird and Duynkerke, 1999; Denby and Greuell, 2000; Grisogono and Oerlemans, 2001; Parmhed et al., 2004). Simultaneous wind and momentum flux profile measurements at the study site would be required to confirm the presence and height of a low-level wind maximum, and to identify it as the cause of the error postulated above.

As the majority of studies on glaciers do not have access to direct observations of turbulence or roughness lengths, the performances of the bulk methods were also evaluated using existing assumed roughness values and schemes. As noted, the combination of C_log and the 1/100 roughness scheme provides the best performing assumed bulk method. C_log may have been expected to overestimate the fluxes in the predominately stable atmospheric conditions during the study. However, implementing a roughness scheme with a mean z_0v substantially smaller than the actual value determined for the surface reduces the magnitude of the fluxes calculated by this method. The small underestimation of Q_H can be attributed to implementing a slightly smaller mean z_0t than was observed, and conversely for Q_L, with a slightly larger mean z_0q. C_bR displays the poorest performance of all bulk methods; its stability function excessively suppressing both fluxes, and yielding low correlation values, particularly for Q_H. Q_L is well modeled by C_bR when coupled with the equal roughness scheme; the excessive suppression being balanced by a large z_0q value. While showing general agreement that conditions were predominately stable (Figure 9C), the bulk Richardson number (mean = 0.06 ± 0.02) does not correlate with the observed stability parameter zLec from the EC data (r = 0.09). C_MO returns better performance than C_bR. Q_H is relatively well modeled when coupled with the Andreas roughness model; the small assumed z_0v compensated by an overestimated z_0t, and insufficient turbulence suppression by C_MO. A similar situation using the 1/100 method results in a well modeled Q_L by C_MO.

The implementation of direct observations for all eight fluxes in the melt model returned results that effectively closed the SEB of the study site. However, as complete flux observations on glacier surfaces are rarely available, the effects of simplifying and parameterizing aspects of the SEB are also considered. Some of the most basic assumptions and values used in SEB studies, namely, negligible Q_R and Q_G, and constant T_s (0°C) and L_o (315.6 Wm⁻², with unity emissivity), were applied to the model, and found to have little effect on performance when applied individually or in combination (r = 0.94; RMSE = 0.01 m/day, for daily surface height change). This result is specific to this field study, and may differ significantly if applied to sites other than continental, temperate glaciers. It should also be noted that due to the discussed adjustments made to the L_o data, the assumption of constant T_s has already been partially applied to the model. The impact of parameterizing the turbulent heat fluxes was also assessed, and the choice of bulk method used was found to have a major influence on closure performance. Figure 11A displays the daily deviation from Q_M (obtained through the fully observed SEB) for simplified SEB models with a range of previously discussed bulk methods implemented (in addition to the basic assumptions above). The returned daily model values can vary by ±40%, depending on the scheme utilized, resulting in differences of up to ±0.5 m from the sonic ranger observed (or full SEB modeled) cumulative surface melt over the study period (Figure 11B).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.