Modelling spatial patterns of near-surface air temperature over a decade of melt seasons on McCall Glacier, Alaska

Abstract

We examine the spatial patterns of near-surface air temperature (T_a) over a melting glacier using a multi-annual dataset from McCall Glacier, Alaska. The dataset consists of a 10-year (2005–2014) meteorological record along the glacier centreline up to an upper glacier cirque, spanning an elevation difference of 900 m. We test the validity of on-glacier linear lapse rates, and a model that calculates T_a based on the influence of katabatic winds and other heat sources along the glacier flow line. During the coldest hours of each summer (10% of time), average lapse rates across the entire glacier range from −4.7 to −6.7°C km⁻¹, with a strong relationship between T_a and elevation (R² > 0.7). During warm conditions, T_a shows more complex, non-linear patterns that are better explained by the flow line-dependent model, reducing errors by up to 0.5°C compared with linear lapse rates, although more uncertainty might be associated with these observations due to occasionally poor sensor ventilation. We conclude that T_a spatial distribution can vary significantly from year to year, and from one glacier section to another. Importantly, extrapolations using linear lapse rates from the ablation zone might lead to large underestimations of T_a on the upper glacier areas.

1. Introduction

Glacier mass balance is controlled in large measure by the energy balance at the glacier surface during the melt season. The fluxes of energy exchanged at the glacier–atmosphere interface are controlled by the local climate, which can differ significantly from that at the same elevations outside of the glacier boundary layer. Therefore, the use of climate data, namely air temperature, extrapolated from off-glacier sites has been increasingly questioned for energy balance modelling (Greuell and Böhm, 1998; Shea and Moore, 2010; Petersen and others, 2013; Carturan and others, 2015; Shaw and others, 2017), especially for the calculation of the turbulent fluxes (Braithwaite, 1995; Shea and Moore, 2010; Sauter and Galos, 2016).

The air mass modification induced by a melting glacier, defined as the cooling of a relatively warm air mass flowing over a colder surface, affects the atmospheric layer above the glacier in a manner that compromises the use of off-glacier meteorological variables as the representative of those in the glacier boundary layer (van den Broeke, 1997; Oerlemans, 2001), especially for air temperatures. The presence of a sloped, relatively large surface at 0°C (i.e. a melting glacier) affects both vertical variations of temperature in the column of atmosphere directly above the surface and variations along the glacier flow lines (distance from an upslope summit or ridge to a given point, Shea and Moore, 2010). The cooling from below causes the formation of a stable layer which should preclude further turbulent vertical mixing due to extreme stability associated with buoyancy forces (Treidl, 1970; van den Broeke, 1997; Oerlemans and Grisogono, 2002). Furthermore, as the ambient temperature increases, the heat-transfer coefficient is lowered due to increasing stability. The change in the near-surface air density along the sloping glacier surface is therefore the driving force for downglacier katabatic flow (van den Broeke, 1997; Greuell and others, 1997). Katabatic winds over valley glaciers are normally shallow, and generate a well-defined low-level jet at only few meters above the surface (Oerlemans and Grisogono, 2002). They act as a ‘heat pump’ for the glacier surface, as they generate turbulence that can then disrupt the stable stratification and thus bring heat to the surface (Klok and others, 2005).

In this context, assuming a linear change in air temperature with terrain elevation is questionable, but has nonetheless been widely applied in the studies of glacier energy balance (Arnold and others, 1996, 2006; Nolin and others, 2010). Greuell and Böhm (1998) suggested an analytical thermodynamic glacier-wind model (hereafter indicated as GB model) to calculate temperature distributions along a glacier dominated by katabatic conditions assuming that the temperature of air parcels travelling along the glacier flow line is governed by two main processes: (i) adiabatic heating; and (ii) exchange of sensible heat with the underlying surface. This pattern (exemplified in Fig. 1) usually is in contrast with what is prescribed by linear lapse rates, i.e. colder temperatures at higher elevations and warmer ones on the glacier ‘tongues’ (here in reference to the lower glacier ablation zone). Despite its relevance and insights, the model had not been tested until recently (Petersen and others, 2013; Ayala and others, 2015; Carturan and others, 2015; Shaw and others, 2017). This is partly because air temperature records at multiple locations over glaciers are difficult to obtain and therefore very few datasets existed to support or disprove the theory and model.

Fig. 1. Conceptual scheme of the spatial patterns of near-surface air temperature over a melting glacier under warm conditions prescribed by the three approaches discussed in this paper: a linear lapse rate derived from the data collected on the ablation zone, a typical solution of the GB model equations (Greuell and Böhm, 1998), and the modified GB (ModGB) model (Ayala and others, 2015).

Petersen and others (2013) tested the model with a dataset of 2 m air temperature measurements at seven locations on Haut Glacier d'Arolla and found that the model was only a small improvement over linear lapse rates in its original formulation. The authors thus suggested an approach where the ratio of the boundary layer height to the bulk heat transfer coefficient, a tuning parameter in the model, was allowed to vary along the glacier flow line. Using the Haut Glacier d'Arolla data and additional datasets from Place Glacier (Canada) and Juncal Norte Glacier (Chile), Ayala and others (2015) showed that, during warm off-glacier conditions, a pattern of decreasing air temperature along the flow line derived from the GB equations (one of the possible solutions shown in figure 7 in Greuell and Böhm, 1998) was evident only up to a given distance from the upper section of the glacier, after which, air temperature increased on the lower sections. They attributed the observed pattern to additional heating sources over well-confined glacier tongues (van den Broeke, 1997; Greuell and Böhm, 1998; Shea and Moore, 2010) and modified the model formulation to account for this additional, non-adiabatic warming shown by the data, which resulted in a non-linear relation between elevation and air temperature.

The modified model (indicated in that publication as ModGB model) has a key relevance for melt modelling purposes, as linear lapse rates reconstructed from observations from the glacier tongue only would be very steep and result in cold air temperatures and reduced melting at high elevations (Ayala and others, 2015). Carturan and others (2015) found that the original GB model was the best alternative temperature distribution method to a linear lapse rate for modelling mass balance of very small, fragmenting glaciers in the Ortles-Cevedale range, Italy, though with evidence of additional adiabatic warming towards the glacier tongue. The modified model of Ayala and others (2015) has been tested only once to date, by Shaw and others (2017) on the small Tsanteleina Glacier in the Italian Alps during one summer, where the model provided an improvement over linear lapse rates for warm conditions favouring the onset of katabatic flow. However, the short flow line length and consequent weak effect of the glacier boundary layer on Tsanteleina Glacier emphasised the difficulties in the transferability of ModGB parameters and that work alone could not provide conclusive evidence about the future utility of the model.

Here, we analyse a decade-long, distributed dataset of near-surface air temperature (T _a) from McCall Glacier (69.3°N, 143.8°W), a ~7 km long glacier in the eastern Brooks Range of northeast Alaska. This Arctic environment is significantly different from those in which the ModGB model has been tested so far. Data spanning 10 years, from 2005 to 2014 (Nolan, 2019), are used to analyse the spatial patterns of air temperature over the glacier and test the applicability of the ModGB model in time. Accordingly, we address the following specific aims:

1. (1) Identify the main patterns of T _a variability during multiple summer ablation seasons (June, July, August) on a high latitude (69.3°N), Arctic valley glacier.

2. (2) Test the ability of linear lapse rates and the non-linear ModGB model to describe distributed T _a with a focus on model robustness and transferability in time.

2. Study site

McCall Glacier is a polythermal valley glacier in the eastern Brooks Range (69°3′N, 143°8′W), northeast Alaska (Klok and others, 2005; Delcourt and others, 2013) (Figs 2a, b). It covers an area of about 6.35 km² with a maximum flow line distance of 7200 m over an elevation range from 1375 to 2635 m above sea level (m a.s.l.) (Fig. 2a). The glacier lies about 100 km south of the Arctic Ocean in a north-facing, steep-sided valley (Nolan and others, 2005). The average slope of the glacier along the flow line is 7.6° and the flow line is uniformly inclined along the line instrumented in this study (Fig. 2c). A hanging glacier, disconnected from McCall Glacier, is located on a small sub-catchment to the east of the main glacier.

Fig. 2. (a) Map of McCall Glacier, showing the distance along the flow line (m) and the elevation contour lines (m a.s.l). Meteorological stations correspond to single temperature loggers (T-loggers) and automatic weather stations (AWS). (b) The location of McCall Glacier in Alaska. (c) Elevation of the stations versus their distance along the flow line with its origin in the upper cirque of the glacier, where T6 is located. A line is fitted to the relation between the distance along the flow line and elevation of the meteorological observations. Flow line distance was calculated using Matlab's TOPOtoolbox (Schwanghart and Kuhn, 2010), and it is only shown for the upper area where T6 is located.

The glacier has a mountain climate with higher precipitation rates (around 500 mm per year) compared to coastal and interior Alaska (Wendler and Ishikawa, 1974a). The Arctic Ocean (100 km to the North) and the Bering Sea (700 km to the West) are the main sources of moisture for precipitation under a dominant southwesterly wind (Klok and others, 2005). A small mass turnover including relatively small ablation rates and small spatial differences in the mass balance are characteristic of the glacier (Wendler and others, 1975; Wendler and Ishikawa, 1974a; Rabus and Echelmeyer, 1998). Mean monthly air temperature ranges from around −30°C (winter) to 5°C (summer) and the ablation periods last on average from June to August (3 months). This period is shorter than those generally observed in mid-latitude or maritime glaciers (Klok and others, 2005). Klok and others (2005) found that ablation during summer season 2004 was larger than previous measurements in the early 1970s, which was attributed to a combination of higher net radiation, a lower albedo and larger turbulent heat fluxes. In their study, they also found a strong dampening effect of the glacier surface on the diurnal variation of meteorological variables over the glacier tongue.

The climate of McCall Glacier is different from the climate of the glaciers where temperature extrapolation approaches and models used in this study have been developed or tested so far: (i) Juncal Norte Glacier in the dry Andes of central Chile (semi-arid, seasonality-driven climate with relatively low humidity and high temperature Masiokas and others, 2006; Ayala and others, 2017); (ii) Haut Glacier d'Arolla (temperate glacier in the Val d'Herens, Valais, Switzerland, dominated by an Alpine climate Arnold and others, 1996); (iii) Place Glacier, located in the southern Coast Mountains of British Columbia, Canada (in between maritime and continental climatic influences Moore and Demuth, 2001; Shea and others, 2004, 2009; Shea and Moore, 2010); and (iv) Tsanteleina Glacier, a small glacier in the temperate climate of the northwestern Italian Alps where the effect of katabatic winds was found to be weaker than for Haut Glacier d'Arolla (Shaw and others, 2017).

3. Data

3.1. Topography

We use a Digital Elevation Model (DEM) of McCall Glacier acquired in summer 2015 at a resolution of 1.1 m using airborne photogrammetry (Nolan and others, 2015). The distance along the flow line is calculated from the DEM using Topotoolbox (Schwanghart and Kuhn, 2010) in MATLAB. We only make use of the flow line descending from the east upper glacier cirque, where the meteorological instrumentation is located. The coordinate system of the flow line follows that defined by Greuell and Böhm (1998), in which the surface is projected over the horizontal plane. Figure 2 shows the topography of the study area (in contours), and the calculated flow line.

3.2. Meteorological data

Since May 2003, various automatic weather stations (AWSs) with different setups have been installed on a transect along McCall Glacier, and on its surroundings. For the observation period considered in this work (1 June 2005 to 31 August 2014) on- and off-glacier AWS provided data for T _a (°C), relative humidity (RH, %), wind speed (m s⁻¹) and direction (°), off-glacier net radiation (W m⁻²), and on-glacier 4-component radiation (W m⁻²). Hereafter, the on- and off-glacier AWSs are referred to as T4 and M1, respectively. In addition to T4 and M1, several AWSs equipped with an air temperature logger and sometimes a few other sensors have been deployed over the glacier at varying heights (from 1 to 3 m due to initial station setup). These AWSs are hereafter referred as T-loggers.

The stations have different tripod designs (free-standing or drilled into the ice) and sensors (Table 1). To avoid issues derived from the use of different sensors and ventilation systems, we only use temperature data measured with Onset (Onset Computer Corp) S-TMB-M002 sensors in Onset shields, with the exception of the off-glacier AWS, where only a Vaisala HMP45C sensor in a Young 12-fin shield was available (Table 1). We use data from six T-loggers, indicated as T1–T6, which are located along the glacier flow line (Fig. 2), selected because they contain the longest records of T _a and have a consistent location over the study period. For comparison with the online published archives of stations (Nolan, 2019), the original nomenclature of each station is given in parentheses in Table 1.

Table 1. Main characteristics of meteorological stations used in this study

Note: Station names in parentheses refer to original station names provided by M. Nolan.

* Years in which measurements for June, July and August are recorded.

† Might have occasionally varied due to snow accumulation.

‡ Only variables used in this study are listed.

§ Young 12-fin shield, Onset shield with undersides painted black, Onset shield with low gills flipped and black interiors with fewer fins total, stock Onset shield.

The data used in this study are summarised in Table 1. All data were logged in 15 min intervals and averaged into hourly values prior to analyses. The coordinates, elevation and location along the glacier flow line of the meteorological stations are specified in Table 2.

Table 2. Location of meteorological stations

4. Methods

4.1. Data preparation

As sensor heights vary from station to station (Table 1) we correct them to the reference 2 m above the surface. For this, we use observations at T4 from 2006, when we have complete measurements of air temperature at a height of 1, 2 and 3 m above the surface for the entire melt season. In so doing, we assume a stable relationship of temperature with vertical heights ≤3 m for different years. We use the T4 measurements to create vertical temperature profiles of near-surface temperature between 0 and 3 m for three off-glacier weather conditions (cold, medium and warm off-glacier ambient temperatures, see next subsection for details). As the maximum temperature over the vertical profile is typically formed at heights >10 m (Munro and Davies, 1978; Oerlemans and Grisogono, 2002; Shea and Moore, 2010), well above our highest sensor, we use linear profiles to derive the air temperature at the 2 m height. This is supported by the lack of a logarithmic profile in the observations at T4, whereby deviations from a linear profile were < 0.3°C for the warmest conditions.

4.2. Data sub-setting

To understand how the on-glacier T _a is influenced by the surrounding meteorological conditions, we analyse the spatial patterns derived from our observations of on-glacier T _a for distinct off-glacier temperature conditions (measured at M1). For this, we select a subset of on-glacier T _a data that correspond to three groups limited by the percentiles 0–10th (P10, cold), 45–55th (P45–55, medium) and 90–100th (P90, warm) recorded at the off-glacier AWS (M1). This means that, if an hourly temperature value at M1 lies in the 90–100th percentile of that melt season, contemporaneous on-glacier T _a at all T-loggers are extracted and placed in the corresponding group. If any station was not recording during a particular hour (due to logger failure, e.g.), we exclude that hour for all stations in our study, so that every timestep used is comparable with another; in this way gaps in any record do not bias our means or analyses.

We use this data sub-setting to test the performance of: (i) linear lapse rates (Section 4.3), and (ii) the non-linear ModGB model of Ayala and others (2015) (Section 4.4). Additionally, we evaluate on-glacier wind at the AWS T4 to determine the dominant wind direction and infer the presence of katabatic winds, the occurrence of which defines the conditions for which the ModGB model was developed (van den Broeke, 1997; Greuell and others, 1997; Greuell and Böhm, 1998). Finally, to understand the inter-annual variability of the observed spatial patterns of air temperature and their possible controls, we explore the relationship of ambient air temperature to incoming shortwave radiation at M1 station and on-glacier albedo at T4.

4.3. Linear lapse rates

We calculate linear lapse rates by a regression of T _a observations against their elevation. The strength of the linear relationship between air temperature and elevation is provided by the coefficient of determination (R ²). Lapse rates are considered positive (negative) when T _a increases (decreases) with elevation (e.g. Minder and others, 2010). To avoid ambiguities, we follow the terminology established by Pepin and Losleben (2002), and followed by several other authors (Chutko and Lamoureux, 2009; Petersen and Pellicciotti, 2011; Petersen and others, 2013; Shaw and others, 2017), where a steep lapse rate has a rapid decrease of temperature with elevation, and a shallow lapse rate is a less negative or positive lapse rate. We also note that we use the term lapse rate as a change in T _a with increasing terrain elevation, which is different from the definition frequently used by the meteorological community, in which lapse rate refers to variations along a vertical profile in the atmosphere (e.g. Gardner and others, 2009).

4.4. ModGB model

The ModGB model (Ayala and others, 2015) is a modification of the Greuell and Böhm (1998) model for on-glacier T _a under the presence of katabatic winds to account for additional heat sources over glacier tongues. Compared to the use of a uniform lapse rate, the original GB model is markedly different in terms of its sensitivity to the temperature outside the glacier boundary layer, as the influence of off-glacier temperature asymptotically decreases along the flow line of the glacier (Greuell and Böhm, 1998). Ayala and others (2015) added a term to the GB model to account for the warming over the glacier tongue that was observed on several glaciers. The resulting equations derived by solving the energy balance of a descending air parcel due to katabatic winds, plus the term including the additional heating over the glacier tongue, are (see Greuell and Böhm, 1998; Ayala and others, 2015 for explicit derivation):

(1)$$\underbrace{\underbrace{T_a\lpar x\rpar = \lpar T_0 - T_{{\rm eq}}\rpar \exp\left(-{x-x_0\over L}\right)+ T_{eq}}_{\rm Original\, GB\, model}+K\left({x-x_0\over L}\right)}_{\rm ModGB\, model}$$

(2)$$L = {H \cos \alpha\over C_H}$$

(3)$$T_{{\rm eq}} = b L$$

(4)$$b = \Gamma_{{\rm d}} \tan \alpha\comma\; $$

where x is the distance along the flow line (m), T ₀ (°C) is the air temperature at the top of the flow line, x ₀, T _eq (°C) is the air temperature at x = ∞ (Eqn (3)), L (m) is the characteristic length scale (Eqn (2)), K (°C) is an empirical factor used to parameterise the influence of additional heat sources on the glacier tongue, b (°C m⁻¹) is the modified dry adiabatic vertical lapse rate (Eqn (4)), α (°) is the mean glacier slope, $C_H\ \lpar \!-\!\rpar$ is the bulk transfer coefficient (approximated as 0.002) (Greuell and Böhm, 1998), Γ_d (°C m⁻¹) the dry adiabatic vertical lapse rate (− 0.0098°C m⁻¹), and H (m) is the katabatic layer height. We derive T ₀ using the environmental lapse rate (ELR = −0.0065°C m⁻¹) to distribute T _a from M1. As the difference in elevation between M1 and T ₀ is 89 m, a 10% change in the lapse rate to derive T ₀ results in a difference of only ±0.05°C.

We use H and K as tuning parameters when fitting the model to the air temperature observations based on the off-glacier temperature sub-setting for each year (Ayala and others, 2015; Shaw and others, 2017). In doing so, we derive parameters which can be used to predict on-glacier air temperatures as a function of off-glacier conditions.

4.5. Inter-annual transferability

To determine the temporal stability of the ModGB parameters and linear lapse rates, we evaluate their inter-annual transferability using a leave-one-out cross-validation procedure. Here, instead of a percentile-based classification in which percentiles boundaries vary from year to year depending on weather conditions (e.g. P90), we test the performance of the models using fixed 1°C bins of off-glacier temperature data (recorded at M1 station). The procedure is as follows: for each bin and year, we fit the ModGB model to the average values of contemporaneous on-glacier T _a measurements at each site, and obtain a set of model parameters (H, K and slopes in the case of the linear lapse rates). The inter-annual median of these parameters are calculated from all other years except the test year then used to predict bin-averaged T _a in that test year. For example, if 2005 was the test year, we apply the median parameters derived from the relevant temperature bin in years 2006–2014 to predict the along-flow line temperature of 2005, and so on.

4.6. Uncertainty of measurements

To provide an estimate of uncertainty in the temperature observations, we conduct an inter-comparison of 2 m air temperatures at sites T4 and T6. Details of these experiments are given in the supplementary material and summarised here for brevity. For all years under analysis, we compare four Onset thermistor measurements in separate naturally-ventilated radiation shields at T4 (Table 1), under warm P90 conditions at M1 and low (< 1 m s⁻¹) wind speeds recorded at T4. We select these conditions because they promote the highest heating errors for naturally ventilated shields. Tests reveal mean absolute differences of < 0.2°C (n = 814) and < 0.3°C for the above conditions with the presence of high snow depths (> 1 m). Propagated manufacturer uncertainty for sensors and radiation shields is 0.4°C for well-ventilated conditions, and given small differences in inter-comparison tests, we consider this as our maximum uncertainty for all sites, with the exception of T6.

The fixed station T6 is subject to specific conditions that likely promote larger heating errors associated with the instrumentation, such as a small fetch (Fig. 2), low glacier winds and deeper snow depths (closer proximity of sensor and snow surface). Based on a comparison with another nearby temperature sensor that was present in some of the study years, we estimate an uncertainty of 0.8°C for this site under P90 conditions.

5. Results

5.1. Meteorological conditions

Figure 3 shows the dominant meteorological conditions and inter-annual variability at McCall Glacier during the study period. Ambient air temperature (°C), wind speed (m s⁻¹) and incoming shortwave radiation (W m⁻²) are shown for the three different off-glacier temperature percentile groups (P10, P45–55 and P90) defined in Section 4.2.

Fig. 3. Average off-glacier meteorological variables recorded at M1 in the melt seasons (June to August) of years 2005–2014. (a) Air temperature, (b) wind speed and (c) incoming shortwave radiation. Colours indicate different off-glacier temperature conditions (T > P90 percentile group: red, T = P45–55: green, T < P10: blue).

Ambient summer air temperatures (as recorded at the off-glacier AWS, M1), averaged over P10, P45–55 and P90 percentile groups, vary between −10 and 10°C, and the inter-annual variability of the lowest temperatures presents the highest amplitude (Fig. 3a). This is mostly explained by the relatively high values of the coldest temperatures recorded in 2007, which is the warmest year in this record. Air temperature in 2013 shows also a very high value during warm conditions, but a low value in the P10 percentile group, suggesting that temperature varied over a large range over the course of that summer. We do not find a trend in the off-glacier temperature records given at M1 station during our period of analysis.

Average wind speed in the percentile groups vary between 1 and 5 m s⁻¹ (Fig. 3b). The inter-annual variability of wind speed is highest for cold conditions, which also show the highest mean values, with 2006 and 2013 showing the highest average wind speeds (> 4 m s⁻¹), associated with low-pressure conditions over the northern coast of Alaska, as it can be inferred from reanalysis composites of the National Oceanographic and Atmospheric Association (https://www.esrl.noaa.gov/psd/data/composites/day/) (not shown). Incoming shortwave radiation relates well to the inter-annual variability of off-glacier air temperature (Fig. 3c), suggesting that clear sky conditions are associated with high air temperatures. During warm conditions, incoming shortwave radiation is lowest in 2006 (370 Wm⁻²), and highest in 2007 (540 Wm⁻²), while for the rest of the years it varies between 370 and 470 Wm⁻² (Fig. 3c). The years with the highest incoming shortwave radiation coincides with the dominance of relatively high-pressure north of the study site (not shown). During medium and cold conditions, there are no significant differences in the observed incoming shortwave radiation.

In Figure 4, we analyse the time distribution of wind speed and direction on the glacier tongue to investigate the presence of katabatic winds. For this, we plot wind roses using T4 wind data during warm (P90) and cold conditions (P10) for two given years: 2006 (relatively cold) and 2007 (relatively warm). Additionally, we calculate the directional constancy (DC) (Bintanja and others, 2014) for both percentile groups. The directional constancy is calculated as:

(5)$${\rm DC}={\lpar \overline{u}^2+\overline{v}^2\rpar ^{1/2}\over \lpar \overline{u^2+v^2}\rpar ^{1/2}}$$

Fig. 4. Wind roses for on-glacier AWS T4 during relatively warm (P90) and cold (P10) off-glacier temperature conditions. The two selected years represent cold (2006) and warm (2007) off-glacier temperatures (recorded at M1 AWS). DC is the directional constancy of each wind rose.

A DC value of 1 means that the wind blows in one direction all the time, whereas a DC value of zero means that wind direction is variable.

Prevailing downglacier winds and a high directional constancy (>0.95) suggest the development of katabatic winds during warm conditions. During cold off-glacier conditions (P10), wind speeds are similar in magnitude to those observed during warm off-glacier conditions, but wind is less directionally consistent (as found on other glaciers by Greuell and others, 1997), meaning that the glacier boundary layer is weak and likely disrupted by synoptic or warmer up-glacier winds controlled by the valley climate system (Oerlemans, 2001). During the warm 2007 year, directional constancy in the P10 percentile group reaches a value of 0.94, which suggests that katabatic winds can develop also during summer cold periods. A relatively large hanging glacier to the southeast of T4 (Fig. 2) could be related to the development of winds from the SE direction (Fig. 4b).

5.2. T_a patterns

To examine how the T _a distribution is affected by processes that develop and strengthen along the length of the glacier, such as katabatic winds (see Introduction), we plot T _a along the glacier flow line for the three off-glacier conditions defined by percentiles groups P10, P45–55 and P90 (in Fig. 5). Percentile ranges vary according to the ambient conditions outside the glacier boundary layer for each year, and are presented in Table 3. For warm conditions, when the downglacier katabatic wind is dominant, T _a follows a downglacier cooling (or no variation) from the top of the flow line towards the mid-flow line (from T6 to T5) and successive warming towards the lowest elevations (from T4 to T1) (Fig. 5 – red shaded areas). This pattern during warm conditions is noticeable for the years 2007, 2008 and 2010–2014, though a more linear temperature relation with the flow line is found for the years 2005, 2006 and 2009. For cold and medium conditions, T _a has a more linear relation with elevation, without thermal inversions on any glacier section.

Fig. 5. Near-surface air temperature (T _a) at different meteorological stations along the flow line of McCall Glacier during cold (blue, P10), average (green, P45-55) and warm (red, P90) off-glacier conditions. Dots represent mean air temperature at a specific station. Standard deviations are shown by shaded colour bounds that are linearly interpolated between the stations.

Table 3. Off-glacier temperature (T _off, measured at M1) ranges for the different percentiles in Figure 5

5.3. Linear lapse rates

Figure 6 shows the variation of average on-glacier linear lapse rates – and their confidence intervals – with off-glacier thermal conditions, in every year of the study period. Lapse rates are calculated for a given hour and averaged over each group of percentiles, in each year. The figure includes average lapse rates calculated using all the on-glacier data (from T1 to T6) and lapse rates calculated using only data collected on the ablation zone (from T1 to T5). During cold conditions (Fig. 6a), average linear lapse rates over the entire glacier (black dots) vary between −4.7 and −6.7°C km⁻¹, whereas lapse rates over the ablation zone (orange dots) are slightly shallower. The differences between lapse rates that include or not the upper glacier area (i.e. T6) are relatively small, and excluding a few years, the coefficient of determination is high (Fig. 6b), suggesting a strong relation between T _a and elevation in every section of the glacier. During medium off-glacier temperature conditions (Fig. 6c), the observed spatial pattern of on-glacier air temperature can be described over every glacier section by linear lapse rates (R ² higher than 0.7) which are typically between −4 and − 5°C km⁻¹, which is shallower than for cold conditions. Similarly to the cold conditions, the values of R ² in medium warm conditions are also high (Fig. 6d). For warm conditions (Fig. 6e), the temperature pattern can no longer be described by linear lapse rates over the entire glacier, resulting in shallow lapse rates (Fig. 6e) with a poor goodness of fit (R ² values as low as 0.28 – Fig. 6f). The P90 lapse rates over the entire glacier demonstrate a trend towards steepening in latter years (Fig. 6), though the strength of these linear relationships with elevation remains relatively poor (<0.7). In opposite, except for a few years, the lapse rates over the ablation zone are actually very steep (up to −7°C km⁻¹), with higher values of R ² (>0.7). With the exception of P90 conditions, all R ² values are significant on the 0.95 level.

Fig. 6. Average temperature linear lapse rates with confidence intervals (a,c,e), and their corresponding coefficients of determination R ² between air temperature and elevation (b,d,f) for the period 2005–2014. The lapse rates are calculated using all T-loggers (from T1 to T6, black colours) and only those on the lower area of the glacier (from T1 to T5, orange colours). Dots represent the average and the bars show $95\%$ confidence intervals.

In Figure 7, we investigate the relation between the strength of the fit of the linear lapse rates during warm off-glacier conditions and the snow cover conditions on the glacier tongue. The plot shows summer time series of the daily average values of the R ² between T _a and elevation, daily albedo at the location of T4 (calculated as daily averages of hourly values that were previously checked for outliers), and the number of hours per day in the P90 percentile group of each year. During the melt season of years 2005, 2006, 2008 and 2009, there were several summer snowfall events, and thus a smaller duration of exposure of bare ice. The number of days with a larger bare ice exposure appears to be larger during the melt season of year 2007 (the warmest observation year - Fig. 3a) and in the period 2010–2013. R ² values are generally above 0.75 over the study period, but frequently drop to values down to 0.25. In every year, P90 hours are not uniformly distributed over the melt seasons, and they tend to group in clusters likely associated to periods of fine weather in mid-summer. While some of the drops in R ² are evidently coincident with the occurrence of P90 hours (see, e.g. the last two groups of P90 hours in 2005, or the medium group in 2007), they are not always clearly connected. In relation to the albedo, we find low and high R ² values during periods of bare ice exposure (see, e.g. the long periods of low albedo in the mid-summer of 2011), but it is difficult to find periods with a highly reflective snow surface and high R ² values.

Fig. 7. Comparison of daily albedo measured at T4 (blue line), temporal distribution of the hours included in the P90 interval (orange bars), and daily average determination coefficient of hourly linear lapse rates (green line) between June and August of each year in the period 2005–2014. X and Y -axes are identical for each subplot.

5.4. ModGB model

Figure 8 shows the comparison between fits of linear lapse rates and the non-linear ModGB model (black and blue lines) to the on-glacier temperature data during warm conditions (P90) for each year, as a function of distance along the flow line. By definition linear lapse rates are calculated against terrain elevation, but given the strong linear relation of the flow line with elevation (see Fig. 2c), lapse rates also show a linear behaviour when plotted against the distance along the flow line. Error bars for each T-logger are 0.8°C for T6 and 0.4°C for the rest of the sensors, as estimated in Section 4.6. The Root Mean Squared Error (RMSE) metric is used to compare the performance of both models. During some years (2005, 2006 and 2009), the application of ModGB results in similar performance for warm conditions as when applying a linear lapse rate, despite not relating clearly to any prevailing meteorological conditions measured at the site (Fig. 3). In the rest of the years, particularly in the warmest ones, the ModGB model represents an improvement over linear lapse rates for warmer conditions in almost all years, with reductions in RMSE varying from 0.2°C (in 2008) up to 0.5°C (in 2014). We note that the uncertainty in measurements at T6 (error bars in Fig. 8 – see supplementary information) could decrease the improvement in RMSE of the ModGB model in comparison with linear lapse rates. This is because a decrease in the air temperature measured at T6 would place that point closer to the linear relation extrapolated from the stations on the ablation zone.

Fig. 8. Fit of linear lapse rates and the ModGB model for the warmest conditions of each year (P90). The lower limit of air temperatures in P90 for each year is shown, together with the RMSE of linear lapse rates and the ModGB model. Black dots show the mean measured temperature at each station with error bars showing the measurement uncertainty (supplementary material) for T6 (light blue) and other stations (orange).

The calibrated parameters H and K for each year are plotted against average T ₀ values in Figure 9, and can additionally be found in Table 4. We discard from Figure 9 the parameters derived for the three years of the analysis period (2005, 2006 and 2009) when the temperature data can be adjusted to a linear relation yielding similar errors. The parameters calibrated for those years are actually outside the range defined by the rest of the years (H > 20 m and K > 7°C). During years when the ModGB model outperforms the linear lapse rates, H, K and the ratio K/L (which represents the magnitude of the temperature increase with the distance along the flow line), are stable and lie between 5.9 and 9.7 m, 3.5 and 6.1°C, and 1.0 and 1.3°C km⁻¹, respectively. We observe a certain increasing relation between K/L and T ₀, but this relation cannot be found for the rest of the parameters. The median values of each parameter correspond to H=7.6 m, K = 4.2°C and K/L = 1.2°C km⁻¹.

Fig. 9. Fitted ModGB model parameters (H, K and K L ⁻¹) as a function of T ₀ temperatures during warm conditions (P90) for the 2005–2014 period, excluding the three years when ModGB shows a performance similar to that of the linear lapse rates. The dashed lines represent the median values (H=7.6 m, K = 4.2°C and K/L = 1.2°C km⁻¹).

Table 4. Fitted ModGB model parameters and T ₀ temperatures for McCall Glacier during P90 conditions of several years (Figure 4)

Note: We highlight the years in which ModGB parameters did not converge to typical values in bold.

In Figure 10 we show a summary of the results of the inter-annual transferability test of ModGB and linear lapse rates. This test is based on a leave-one-out cross-validation procedure using temperature bins of off-glacier data (see Section 4.5 for details). As explained above, the years 2005, 2006 and 2009 are not included in the calculation of median parameters because the ModGB parameters converged to values well outside the range defined by the rest of the years. On the left axis, the figure shows the inter-annual average RMSE of the fitted models (in bars) against T ₀ bins. On the right, the axis shows the variation of the median parameters of the ModGB model. On the selected temperature range, the average RMSE of ModGB is up to 0.5°C smaller than that of the linear lapse rates. While H tends to decrease from 9 to about 6 m with increasing ambient temperatures, K/L slightly increases from 1 to 1.5°C km⁻¹. We identify that 20% of the total hours of data experience off-glacier temperature > 5°C at M1 station.

Fig. 10. Comparison of the inter-annual transferability of linear lapse rates and the ModGB model as a function of the off-glacier air temperature. Left axis shows the average RMSE of each model (in bars) and the right panel shows the ModGB model parameters (in lines).

6. Discussion

6.1. Meteorological conditions and spatial patterns of T_a

Our results provide an overview of the meteorological summer conditions on McCall Glacier in the period 2005–2014, with a special focus on air temperature and its spatial distribution over the glacier. Such decadal meteorological records of on-glacier observations are rare, particularly those with observations that cover entire glaciers (Carturan and others, 2015; Shaw and others, 2017). Although we did not find an evident temporal trend in ambient (off-glacier) air temperature, or other meteorological variables such as wind speed, we found two years, 2007 and 2013, in which the average of the 10% highest air temperatures (P90) recorded at the off-glacier AWS were higher than 10°C (see Fig. 3a), with hourly maxima of 15.1 and 14.7°C, respectively. Average temperature in these years was higher than that in 2004, which was recognised by Klok and others (2005) as a year with very negative surface mass balance, suggesting that at least two years in the decade of observations presented here had large ablation rates as well. Klok and others (2005) also found that katabatic winds are constantly present on McCall Glacier, but that they are relatively weak, and that their daily cycle might be affected by valley winds. Our results from the analysis of percentile categories show that katabatic winds dominate over valley winds during warmer conditions (Fig. 4), and that valley and synoptic winds are present only during colder conditions, when they may disrupt the katabatic layer.

The analysis of air temperature spatial distribution reveals average patterns that are also common in mid-latitude glaciers, such as Haut Glacier d'Arolla, Tsanteleina and Juncal Norte glaciers (Ayala and others, 2015; Shaw and others, 2017). During cold conditions, there is a weak thermal contrast between the glacier surface and the surrounding terrain, and air temperature observations over the glacier surface can be well described by linear lapse rates on both upper and lower glacier sections (Figs 6a, b), which may likely be extrapolated with greater confidence from off-glacier stations, depending on their site positioning (Shaw and others, 2017). These conditions are strengthened when the glacier surface is covered by snow (Fig. 7). As off-glacier air temperature increases, the thermal contrast intensifies, katabatic winds appear (Fig. 4), and the glacier cooling effect emerges as dominant. As a consequence, air temperature decreases from the top of the flow line in a downglacier direction due to the dominant role of turbulent heat exchange of downglacier-moving air with the glacier surface, and rapidly increases over the glacier tongue due to adiabatic heating, heat advection from local winds, terrain irradiance and air entrainment from upper atmospheric layers. Some of these processes on the glacier tongue, such as heat advection from surrounding terrain and longwave radiation from mountain slopes, have also been observed over glaciers under fragmentation (Jiskoot and Mueller, 2012; Carturan and others, 2015), and are evident in detailed numerical simulations of the glacier boundary layer (Sauter and Galos, 2016). However, as suggested by Klok and others (2005), the temperature contrast between the ambient atmosphere and the glacier surface seems less intense than on mid-latitude glaciers, likely due to the absence of hot valley winds that penetrate on glacier valleys, and the relatively weak solar heating of mountain slopes (Wendler and Ishikawa, 1974b). An additional problem can be the disconnection of the glacier boundary layers above the ablation zone and the upper glacier cirque where the uppermost sensor T6 is located. However, downglacier winds from all the upper areas are usually funnelled to the ablation zone, and extensive literature supports this effect (e.g. Munro, 2006; Jiskoot and Mueller, 2012).

As the spatial patterns we found are derived from temporal averages, they are not necessarily representative of the air temperature distribution at a particular time step, when heat advection can disrupt stable structures in the glacier boundary layer. Sauter and Galos (2016) found that heat advection associated with local wind systems causes small-scale variations of sensible heat flux of up to 100 W m⁻², strongly modifying the on-glacier air temperature distribution. They concluded, using computationally-intensive idealised large-eddy simulations of 2 h duration, that the derivation of wind and temperature distribution at a few observation sites is insufficient to capture temperature distributions over an entire glacier (also a conclusion of Shaw and others, 2017), and that only mass-consistent models including representative wind fields are able to characterise the cross-glacier variability and the glacier boundary layer at particular time steps. However, the very short simulation period of such complex models limits a complete understanding of spatial variability of air temperature because they cannot be tested under different conditions. We think that the study of average patterns is useful for process understanding and the derivation of simple models that can be applied for large glacier samples or long time periods. Nevertheless, an in-depth analysis of transitions in meteorological conditions during hourly or daily time steps is still warranted with large datasets such as this, and is planned for future work.

6.2. Model performance

Observed linear lapse rates during cold conditions (Fig. 6) are distributed around the ELR (−6.5°C km⁻¹) and show high coefficients of determination. These results are consistent with those previously derived from temperature records on Haut Glacier d'Arolla (Petersen and others, 2013; Ayala and others, 2015) and Tsanteleina Glacier (Shaw and others, 2017). During average off-glacier temperature conditions, the observed linear lapse rates are more shallow and vary around −4 and −5°C km⁻¹, and reproduce the observed temperature patterns accurately (see linear lapse rates and coefficients of determinations of P45–55 in Figs 6c, d). Interestingly, the inter-annual variability of linear lapse rates during cold conditions is high (about from −4.7 to 6.7°C km⁻¹) (Fig. 6a), which might be linked to different types of cold synoptic conditions, or different spatial extent (or albedo) of the seasonal snow cover. For example, Figure 7 shows that temporal variability in albedo measurements at T4 appears to coincide with fluctuations in the elevation-dependency of T _a, though point-based measurements only at this site are not enough to understand how distributed snow cover might affect the linearity of lapse rates. Future work could leverage long-term observations from remotely sensed data (e.g. MODIS or Landsat) to relate snow line retreat to the spatial variability in along-flow line T _a. This, for example, could be a cause of a steepening in the supposed ‘lapse rate’ for P90 conditions in latter years of this analysis (Fig. 6e).

Our results show that during warm off-glacier conditions the ModGB model represents the observed temperature patterns better than linear lapse rates, reducing RMSE values in seven out of ten of the observed years and up to 0.5°C in certain years (Fig. 8). While this overall performance improvement might seem small at the sites of observations, this finding has a higher relevance considering that temperature is normally extrapolated from one or only a few observations from the lowest glacier sections in mass-balance models. This is exemplified in Figure 11 for year 2011. In this year, the differences between the temperatures prescribed by ModGB and linear lapse rates are large. If air temperature is extrapolated from data of stations T1–T4, T _a on the upper glacier areas is underestimated by up to 3°C. This would lead to large errors in melt and mass balance models, and a significant underestimation of melt, especially when calculated by temperature-index models. On the other hand, if only temperature data from upper stations (T5–T6) are considered, the extrapolated near-surface air temperature on the glacier tongue would be too low (about 2.5°C) compared to the actual measured values as well. These findings support the results from previous studies, such as Ragettli and Pellicciotti (2012) or Pellicciotti and others (2014), and point to the advantage of the new model when used to calculate air temperature for input to melt and mass-balance models from only a limited number of observations. To quantify the consequences for melt estimates of using linear lapse rates calculated only using data from the lowest stations (T1–T4) as opposed to the ModGB approach (as shown in Fig. 11), we perform a simple modelling experiment using a simple temperature-index model with two melt factors. We use an hourly melt factor of 0.1 mm h^−1°C⁻¹ for snow and 0.5 mm h^−1°C⁻¹ for ice to estimate total melt in the summer season of year 2011 at two locations the accumulation (T6) and ablation (T5) zones. With this approach, we find an underestimation of 331 mm w.e. of melt at T5 and 152.2 mm w.e. of melt at T6 if using linear lapse rates from the lowest stations instead of the ModGB model. Even considering the minimum T6 temperature associated with the estimated sensor uncertainty (supplementary information), the minimum underestimation of melt using a lapse rate compared to the ModGB would still be 106 mm w.e., equivalent to almost 32% of the summer melt total at that station.

Fig. 11. ModGB model fit and linear lapse rates for the observations in the melt season 2011. Linear lapse rates are fitted to (i) only lower data (T1–T4), (ii) only upper data (T5–T6) and (iii) all data with starting point T4.

Two key uncertainty sources in our study are (i) the occasionally shorter distance from the air temperature sensors to the surface due to snow accumulation, and (ii) occasionally poor ventilation conditions for the temperature sensors during warm conditions, low wind speed and strong solar radiation. The first uncertainty source could impact the shape of the spatial distribution of air temperature, but it should be less relevant during mid-summer, when long periods of exposure of the bare ice occur (see Fig. 7). Second, the specific conditions for poor sensor ventilation could artificially increase the temperature measured by T6 on the upper glacier area. To address this issue, we have analysed the temperature observations at T6, and compared them to those made in a nearby station that measured only during part of the study period. Based on this analysis (presented in the supplementary material), we estimate an uncertainty of 0.8°C for T _a observations at T6, which is about the double of that in the rest of the T-loggers. A value of 0.8°C could partially affect our conclusions (see error bars in Fig. 8), but the main outcomes of the study are maintained, and are supported by both the theory (Greuell and Böhm, 1998), and evidence from other glaciers (Ayala and others, 2015; Shea and Moore, 2010; Shaw and others, 2017).

6.3. Future application of the model

Our findings point to a general transferability of the ModGB model in time: the elevation-dependent temperature patterns and their disruption under warm ambient conditions are consistent with those of former studies (Ayala and others, 2015; Shaw and others, 2017), and the parameters used to model the non-linear behaviour of T _a under warm ambient conditions (H and K) are also comparable to those studies (Table 4). ModGB model parameters are relatively stable over several years with the exception of some (Fig. 9), in which the lapse rates perform similarly to the ModGB model (2005, 2006 and 2009). For example, Ayala and others (2015) and Shaw and others (2017) found H parameters of 4–8 m and K parameters of 6-8°C above the 8°C ambient temperature threshold (see, e.g. Fig. 7 of Shaw and others, 2017). This would suggest that ModGB could be used to distribute the temperature forcing for mass-balance modelling studies on glaciers, assuming the appropriate conditions for its use are met (i.e. warm off-glacier temperatures and development of katabatic winds). The ModGB approach, similarly to the original GB model and the Shea and Moore (2010) approach, is designed to utilise only off-glacier air temperature data. This is an important advantage given that off-glacier temperature are typically the only available data for glacier energy/mass-balance modelling. We present here a direct comparison to on-glacier data and on-glacier lapse rates to demonstrate the departure from the linearity of distributed along-glacier temperatures for warm ambient conditions. Nonetheless, previous work has already highlighted the problems of off-glacier temperature extrapolation and its impact upon over-estimation of on-glacier temperatures (Carturan and others, 2015; Shaw and others, 2017). Localised processes and the distribution of stations used to construct the ModGB curve of each glacier might present a challenge for spatial transferability of the model between different glaciers. For example, defining the top of the boundary layer (x ₀) can affect the strength of the fitted model curve and its ability to explain the same along-glacier behaviour on other glaciers.

In summary, the model's full transferability needs to be investigated further, and related to both an improved understanding of the glacier boundary layer and deeper insights on how it is affected by synoptic conditions, cloudiness and snow conditions (Fig. 7), and how this interplay shapes the temporal and spatial patterns of air temperature. This calls for a clear need for more observations and detailed studies of the glacier boundary layer to understand the physical meaning of the model parameters, and what controls their value and stability. Multi-annual on-glacier datasets of distributed air temperatures are rare, and it seems imperative that distributed air temperature observations are collected over glaciers of varying size, orientation and setting in order to constrain the limitations to its use in replacing linear lapse rates. It is clearly evident from Figures 5 and 10 (as well as suggestions from the literature) that ModGB is not appropriate for estimating T _a for cooler ambient conditions, and as such, well-defined thresholds for its application are still required.

7. Conclusions

This study presents and analyses a 10-year (2005–2014) distributed dataset of near-surface air temperature observations collected during the summer season on McCall Glacier, north Alaska. This unique dataset is used to study the temporal variations of the spatial distribution of air temperature. We also test the applicability of linear lapse rates and the ModGB model, a model that accounts for more complex variations along the glacier flow line linked to heat advection during warm conditions. The key findings of this study are:

1. (1) During cold and average ambient conditions of each year, linear lapse rates reproduce the observed average patterns of air temperature across the glacier accurately (with a strong relationship between temperature and elevation). These lapse rates are steeper than for warmer conditions and vary between − 6.7 and − 4.7°C km⁻¹. This relatively large inter-annual variability might be linked to different synoptic conditions, or different spatial extent (or albedo) of the seasonal snow cover, though the influence of these factors require more detailed investigation.

2. (2) During warm conditions, the average spatial distribution is poorly described by linear lapse rates in agreement with the literature to date. Instead, air temperature follows a pattern that has been previously observed in mid-latitude mountain glaciers, which is driven by the complex spatial variation in the thermal difference between the glacier surface and that of the overlying atmosphere. Although more uncertainty might be associated with these results due to occasionally poor sensor ventilation, we find that ignoring such deviations from linear temperature extrapolations at high elevations could lead to errors and bias in total modelled melt rates.

We suggest that future field campaigns, such as those presented in this study and those from Petersen and Pellicciotti (2011), Shea and Moore (2010), Petersen and others (2013) and Shaw and others (2017), should be combined with more detailed numerical simulations, such as the large-eddy simulations presented by Sauter and Galos (2016). This would help to improve our current understanding of the glacier boundary layer and its dominant energy fluxes, and aid the development of air temperature distribution models, such as presented here, for long-term mass-balance simulations (Greuell and others, 1997; Shea and Moore, 2010; Ayala and others, 2015).