Local Controls on Near‐Surface Glacier Cooling Under Warm Atmospheric Conditions

The near‐surface boundary layer can mediate the response of mountain glaciers to external climate, cooling the overlying air and promoting a density‐driven glacier wind. The fundamental processes are conceptually well understood, though the magnitudes of cooling and presence of glacier winds are poorly quantified in space and time, increasing the forcing uncertainty for melt models. We utilize a new data set of on‐glacier meteorological measurements on three neighboring glaciers in the Swiss Alps to explore their distinct response to regional climate under the extreme 2022 summer. We find that synoptic wind origins and local terrain modifications, not only glacier size, play an important role in the ability of a glacier to cool the near‐surface air. Warm air intrusions from valley or synoptically‐driven winds onto the glacier can occur between ∼19% and 64% of the time and contribute between 3% and 81% of the total sensible heat flux to the surface during warm afternoon hours, depending on the fetch of the glacier flowline and its susceptibility to boundary layer erosion. In the context of extreme summer warmth, indicative of future conditions, the boundary layer cooling (up to 6.5°C cooler than its surroundings) and resultant katabatic wind flow are highly heterogeneous between the study glaciers, highlighting the complex and likely non‐linear response of glaciers to an uncertain future.

• Glacier size and alignment with valley/synoptic wind gradients control the magnitude of near-surface cooling • Valley/synoptic winds can occur between 19% and 64% of the time and contribute between 3% and 81% of total sensible heat to the ice surface • Localized cooling and turbulence in the boundary layer increase the complexity and non-linearity of glacier response to climate

Supporting Information:
Supporting Information may be found in the online version of this article.
significantly cool the intermediate air layer above the glacier due to sensible heat exchange (B.J. Oerlemans & Grisogono, 2002;Van Den Broeke, 1997).Nevertheless, the magnitude of this "glacier cooling effect" is difficult to quantify due to both the uncertainty or insufficiency of long-term off-glacier observations at high elevation and also due to the general lack of on-glacier meteorological data.Regional studies have accounted for this effect using a single offset value or ignored this process through calibration against off-glacier thermal conditions (e.g., Kraaijenbrink et al., 2017;Rounce et al., 2023).However, recent observational evidence has suggested that this cooling effect is non-stationary, and might depend on the stage of glacier evolution, including its size and fragmentation (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).Ignoring such processes that govern the response of glaciers to future warming may hinder our understanding of glacier mass loss and its implications for water resource management.It has been well-established from the available data that a spatially and temporally variable cooling effect exists on mountain glaciers (Ayala et al., 2015;Conway et al., 2021;Shaw et al., 2021;Shea & Moore, 2010;Yang et al., 2022) and that it has non-negligible implications in the calculation of glacier melt rates (Petersen & Pellicciotti, 2011;Shaw et al., 2017).Nevertheless, there remains much to be understood about the variability of near-surface cooling and to what degree local processes dominate its signal.
We seek to address some of these knowledge gaps using data collected as part of a wider objective to parameterize air temperature on mountain glaciers.Here we compare the simultaneous response of different glaciers to the same sub-regional atmospheric conditions.We note that past research has undertaken simultaneous observations on multiple glaciers to address similar questions (Carturan et al., 2015;Shaw et al., 2021;Shea & Moore, 2010).However, unlike this previous work, we have deployed specific observational networks to address wind interactions near the glacier terminus and its influence on local temperature fluctuations in the glacier boundary layer.Our observational period covers an anomalously warm summer of 2022 following a dry winter (Cremona et al., 2023) whereby glacier snow cover was almost completely absent (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).The European Alps experienced their most negative mass balance year since records began (SCNAT, 2022) under climatic conditions which put into question the role of the glacier micro-climate for modifying glacier response in the future.
Accordingly, our aims are: (i) To quantify the degree of near-surface cooling for three distinct glaciers exposed to an extreme summer.
(ii) To identify the role of valley and synoptic winds in modulating on-glacier meteorological conditions.(iii) To dissect the variable role of katabatic winds and near-surface cooling on the glacier surface energy balance.
Supraglacial debris cover has expanded significantly over the last 20 years on Arolla, now covering ∼37% of its 2022 area, extending even to the upper portions of the glacier as englacial debris melts out with increased frequency.For Otemma and Corbassière, debris covers around 2.34 and 1.68 km 2 of the surface, respectively, largely on the glacier terminuses and in the form of medial moraine bands, comprising ∼17% and 10% of the respective total area.
There has been a large body of past research exploring meteorological conditions on Arolla (Ayala et al., 2015;Petersen et al., 2013;Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023;Strasser et al., 2004), identifying potential controls on the Ta variability above the glacier.Otemma and Corbassière have been mostly studied in terms of their mass balance (e.g., Huss et al., 2015), with Otemma Glacier also being the subject of numerous geomorphological and sediment transfer experiments (e.g., Müller et al., 2022).To the authors' knowledge, no prior studies of the near-surface meteorology have been undertaken on Otemma and Corbassière glaciers.

Local Weather Station Data
In this study, the meteorological setup consisted of on-and off-glacier air temperature stations (hereafter "T-Loggers") and low-cost automatic weather stations ("AWS").The T-loggers were free-standing, 2 m tripods that support an Onset TidBit v2 temperature logger (accuracy 0.2°C) housed in an Onset R4, naturally ventilated radiation shield.AWS were DAVIS Vantage Pro 2 wireless stations mounted on free-standing tripods with a sensor height at 2 m (the setup as in Carturan et al. (2015)).The AWS consisted of an artificially aspirated temperature-relative humidity sensor, barometer, tipping bucket precipitation gauge and a mechanical anemometer and wind vane.The sensors sampled data every 5 min and wirelessly transmit data at 30 min intervals to a Weather Envoy data logger, powered by a 6-W solar panel.T-Loggers in pro-glacial domains (Figure 1) additionally possessed a standalone logger with a cup anemometer and vane for exploring the occurrence of valley winds.The cup anemometer for the pro-glacial zone of Corbassière Glacier was damaged and did not provide data for our analysis, though wind direction measurements were still logged.One T-Logger was installed at a higher elevation, off-glacier location at each site, deemed to be outside of the thermal influence of the glaciers (Figure 1).All data were averaged to hourly intervals for analysis and manually inspected for errors.The 2022 data for Arolla Glacier has previously been used in the analysis by Shaw, Buri, McCarthy, Miles, Ayala, and Pellicciotti (2023) exploring the changing boundary layer conditions of that glacier through time .Here we also utilize this data set with a different sub-period of time, but instead focus on the inter-glacier differences of boundary layer development related to the local glacier and terrain characteristics under which it occurs.

Regional Stations and Reanalysis Data
Additional meteorological information were provided by local meteorological data on an hourly timescale from permanent stations of the MeteoSwiss network as well as surface level reanalysis data from the ERA5 (Hersbach et al., 2020) andERA5-Land (Muñoz-sabater et al., 2021) products.Principally, we obtained information about temperature, wind speed/direction and incoming shortwave radiation to classify the local and regional conditions.One ERA5 grid was centered on the entire domain (Figure 1) whereas the closest pixel centroid to the glacier center was examined for ERA5Land

Temperature Corrections and Uncertainty
We conducted an intercomparison of Ta for naturally and artificially aspirated sensors on Arolla and Corbassière glaciers for radiation shields mounted on the same tripod mast (additional sensors were not available for Otemma Glacier).We compared hourly deviations of Ta for conditions known to promote heating errors of naturally ventilated sensors (Carturan et al., 2015;Cullen & Conway, 2015;Georges & Kaser, 2002).Mean (standard deviation of) hourly differences in natural-aspirated temperatures were 0.42 (0.37) and 0.39 (0.44)°C for the AWS on Arolla and Corbassière, respectively.We applied a multiple linear regression of Ta differences with measured wind speeds at those AWS sites and locally derived information of shortwave incoming radiation from MeteoSwiss climate stations in order to correct for heating errors of sensors in naturally ventilated radiation shields.The linear regression model was applied to all T-Logger sites assuming no spatial changes in the incoming shortwave radiation across the glacier and considering an increasing wind speed along the flowline distance of each glacier based upon the observed differences in wind speed between the upper and lower AWS on Otemma Glacier (mean increase of 1.5 −4 m s −1 m).The mean Ta differences at the AWS on Arolla and Corbassière after correction were zero with a standard deviation of differences of 0.26 and 0.33°C (Figures S1 and S2 in Supporting Information S1).
Remaining absolute Ta differences during the daytime are typically <0.4°C (∼90% of all values) and considered to express the remaining uncertainty of the air temperature sensors and the naturally ventilated radiation shields that is not accounted for by this linear regression correction.This range of uncertainty is in line with that proposed by past research (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023;Troxler et al., 2020) and considered to be appropriate for exploring Ta patterns.

Temperature Sensitivity
A key aim of this study is to quantify the magnitude of glacier cooling and its drivers.The glacier cooling effect is quantified as the temporally-averaged bias of the observed above-glacier Ta Obs compared to the equivalent ambient Ta Est at the same elevation (Ta Obs − Ta Est ), where the latter is derived from the gradient (or lapse rate) of Ta from an off-glacier source (Ta OFF ) following: where z Obs and z OFF are the respective elevations of the on-glacier and off-glacier observations and Γ is the temperature gradient.For this study, we utilize the uppermost off-glacier T-Logger (Figure 1) as the source of Ta OFF .Following Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti (2023), we also calculate the on-glacier "temperature sensitivity" (TS) which describes the ratio of changes in Ta Est to the changes in Ta Obs over a given period by: where ß is the intercept of the regression between Ta Est and Ta Obs .A TS of 0.5 indicates that, on average for every 1°C change outside the glacier boundary layer, a 0.5°C change occurs at the near surface.Ta Obs − Ta Est is used to explore temporally varying cooling at the near surface in response to specific prevailing conditions (e.g., synoptic storms or warming-induced slope winds), whereas TS aids the interpretation of a glacier's sensitivity to temperature changes as a function of a broader classification of conditions (e.g., given wind regimes-see Section 5.3).The value of TS suggests a linear scaling between atmospheric warming (Δ Ta Est ) and the relative changes of the observed temperature at near surface (Δ Ta Obs ), which is largely supported by the data of this analysis (Figures S3-S5 in Supporting Information S1) as well as data presented in previous studies (Greuell & Böhm, 1998;Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).
We use the most locally available off-glacier T-Loggers from our study to provide the best estimate of hourly Γ outside the thermal influence of the glacier (black squares in Figure 1).For Corbassière, the lower T-Logger is observed to be influenced by the glacier boundary conditions and so the higher-elevation, off-glacier T-Logger and MeteoSwiss stations are used to estimate Γ.To provide a comparative means of estimating Ta at the glacier elevations, we additionally calculate Γ from ERA5 surface level data over the region (Figure 2) and apply it to the reanalysis temperature, and use all available MeteoSwiss observations to generate hourly variable Γ for extrapolation from the most local MeteoSwiss station to the glacier in question.

Classification and Modeling of Wind Regimes
To explore the drivers of glacier cooling and TS we pay particular attention to the role of the local wind regimes for each glacier.We perform this in two ways: first we examine the dominant direction of the wind following Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti (2023), utilizing an up-glacier wind index (UWI): where DIR indicates the hourly direction of the wind on (or off) glacier taken from the observations (°) and is the up-valley/glacier direction (°).A value of 1 indicates a precisely up-valley or up-glacier wind direction and −1 indicates a precisely down-valley/glacier wind direction.Combining the UWI values with additional meteorological information (see following sub-section) allows us to classify the occurrence of katabatic conditions on the glacier during the summer and explore its variable role in the surface energy balance between glaciers.
We further explore the local off-glacier wind regimes at each glacier by utilizing the operational wind modeler, WindNinja (Forthofer et al., 2014).WindNinja is a wind prediction tool that conserves both mass and momentum and in this study is forced with domain-averaged wind speeds and directions from ERA5-Land reanalysis to model terrain induced wind patterns on a 100 m pixel resolution.We note that this tool is not able to account for the density-induced katabatic winds over glaciers and snow, the principal interest of this work, and all cells are considered to be grass-covered in its current setup.The outputs from WindNinja are exploited to better understand the local modifications of terrain at each site and provide context to the setting of the observed glacier wind regime within its wider surroundings.Specifically, wind patterns derived from this model are compared to the on-and off-glacier wind observations for conditions classified as being strongly katabatic or otherwise disturbed by synoptic gradients (see following sub-section).

Data Sub-Sampling
AWS and T-loggers were installed on the glaciers in early July and removed in mid-September 2022.Due to strong differential ablation, some stations fell during the earlier part of the observation period.In order to make direct comparisons between all glaciers, we principally focus on a 39 day period (11 August to 19 September 2022) during which all stations were stable and recorded Ta at a 2 m height (Figure S6 in Supporting Information S1).Unless otherwise stated, all further analysis is presented for this time period.
We classify warm afternoon hours to explore the relative importance of katabatic conditions on each glacier compared to those defined as being driven by up-valley wind advection or those resulting from synoptic disturbance.We focus upon clear-sky, afternoon hours when katabatic conditions are expected to occur on mid-latitude mountain glaciers (Greuell & Böhm, 1998;J. Oerlemans, 2001J. Oerlemans, , 2010;;Van Den Broeke, 1997).We consider hours of each day ranging from 12:00 to 19:00 where cloudiness (derived from reanalysis) was less than 20% and where the estimated (Ta Est ) freezing level was above the highest point of the glacier, such that the conditions are met for the generation of a down-glacier wind (Shea & Moore, 2010).Provided this mutual condition (∼20% of the hours from the above 39 days), we subdivide the remaining cases into (a) katabatic conditions, with observed on-glacier UWI (Equation 3) <−0.2, where wind speeds are greater than the stall speed of the mechanical cup anemometers (0.5 m s −1 ) and ERA5 gust speeds are below the 50th percentile, (b) up-valley winds, with a UWI > 0.2 and the wind speed conditions as in (a), and; (c) synoptic/disturbed conditions where ERA5 gust speeds are above the 80th percentile of the entire period and/or where −0.2 < UWI < 0.2.

Energy Balance Modeling
We use the glacier component of the land-surface model Tethys-Chloris (Botter et al., 2021;Fatichi et al., 2012Fatichi et al., , 2021;;Fugger et al., 2022;Fyffe et al., 2021;Mastrotheodoros et al., 2020;Paschalis et al., 2022) to model the surface energy balance across the glacier.The model simulates the energy fluxes of the ice/snow surface and subsurface, following: For snow, and: for ice, where Rn (W m −2 ) is the net radiation absorbed by the surface, Qv (W m −2 ) is the energy from precipitation, Qfm (W m −2 ) is the energy gained or released by melting or refreezing within the snowpack, H (W m −2 ) and LE (W m −2 ) are the sensible and latent energy fluxes and G (W m −2 ) is the conductive energy flux from the surface to the subsurface.Conduction of energy within the ice is represented to a depth of 2 m after which it is assumed the ice pack is isothermal.dQ (W m −2 ) is the net energy input to the snow or ice pack.The sign convention is such that fluxes are positive when directed toward the surface.Iterative numerical methods were used to solve the non-linear energy budget equation until convergence for the ice and snow surface temperature, while concurrently computing the mass fluxes resulting from snow and ice melt and sublimation.
Over snow and ice surfaces, the sensible heat flux (H) is calculated as described in Fugger et al. (2022) and roughness values for snow and ice are derived from literature values for the Alps (Brock et al., 2000).We employ a snow albedo scheme following Ding et al. (2017) and a static ice albedo of 0.23, taken from measurements on Arolla Glacier.For a full description of the model components for glacier energy balance, we refer the reader to Fugger et al. (2022).
To force the model, incoming shortwave and longwave radiation are taken directly from ERA5-Land and precipitation and humidity are taken from the lower AWS on each glacier and assumed spatially constant (precipitation is negligible during the modeled period and no snow accumulation is evident).Ta is taken directly from each station and distributed across the total elevation range of each glacier based upon an hourly on-glacier temperature gradient.Wind speed is modeled using its relation to glacier flowline found from the observations on Otemma Glacier (1.5 −4 m s −1 m) and distributed from the lower AWS record for all glaciers.Modeled melt and sensible heat fluxes are calculated from the model in 50 m elevation bands and subdivided into the hours classified as above.Initial snow depths were prescribed as 0.5 m for the upper most cells of Corbassière seen to be snow covered from remote sensing observations (>3,400 m a.s.l.).All other cells were snow free during the model period.For those elevation bands where >50% of the glacier area is covered by debris, we assign a debris cover albedo value of 0.16 (Brock et al., 2010) and a debris thickness value of 0.15 m, taken as the average of measured points along the medial moraine of Arolla Glacier.While this is a model simplification, we consider it acceptable for the short period of analysis and for exploring the role of the glacier boundary layer and the local meteorology associated with it.The model has been demonstrated to work well in past glacier energy balance studies (Fugger et al., 2022;Fyffe et al., 2021) and we show in Figure S7 in Supporting Information S1 that, for comparison to available stake records, the model performs well at representing the surface melt on Arolla Glacier.
For comparison, we also run two alternative model simulations per glacier where (a) the most locally available, off-glacier Ta Est value (Equation 1) is used in place of the measured on-glacier Ta ("COOL"), or (b) where both the local Ta Est and off-glacier wind speed data are used in place of the measured/calculated on-glacier values ("WIND").For the latter, we consider wind speeds from the off-glacier T-Logger and from MeteoSwiss stations in the case of Corbassiere.We consider sensible heat flux and melt differences as COOL-REF and WIND-REF, respectively, where "REF" refers to the reference model scenario using observed on-glacier meteorology.This aids interpretation of the role of the near-surface meteorology in determining the energy balance and melt of the glacier during warm atmospheric conditions.

Glacier Cooling During an Extreme Summer
The observed, on-glacier Ta for all three sites is notably cooler than the equivalent off-glacier Ta at comparable elevations (Figure 2).For the elevation range of the glaciers, mean off-glacier Ta is fairly consistent when comparing the local (off-glacier T-Loggers), regional (MeteoSwiss stations) and reanalysis (ERA5-Land) Ta, and the mean biases of the on-glacier observations relative to their off-glacier counterparts range from −2 to −6°C (Figure 2b).The main exceptions to this are the pro-glacial station for Corbassière which appears to be strongly influenced by the glacier wind layer and thus cooler, and those ERA5-Land pixels classified with a glacier concentration ≥50% (see Figure S8 in Supporting Information S1).
The Ta measurements show that glacier cooling is highly variable in space and between the three glaciers (Figure 2b), especially toward lower elevations, though variability is reduced (within ∼1°C) above 2,800 m a.s.l.The magnitude of this cooling is dictated by the atmospheric temperature and the gradient of temperature above a 0°C ice surface (Figure 3).This "cooling effect" is most notable once ambient (off-glacier) temperatures at the mean glacier elevation exceed 4°C (Figure 3a).Beyond ∼6°C, the rate of cooling on Arolla is smaller than is visible for Otemma and Corbassière (Figure 3a), though the variability of cooling is greater for Otemma (shaded area), largely dictated by differences in the observations near the terminus.On the warmest afternoon of the sub-period (13 September 2022), with a mean Ta Est of 14.3, 14.6, and 15.5°C for Arolla, Otemma, and Corbassière, respectively, mean cooling across all stations was −2.5, −2.8, and −4.6°C and −2.4,−3.1, and −5.6°C for cooling at the terminus.
These patterns are consistent between the different forcing data sets, though the magnitudes of cooling are significantly different using regional station or reanalysis data (Figures 3b and 3c) due to both the representativeness of the Ta at the station/pixel location and the appropriateness of the lapse rate to describe local conditions at the glacier.For example, using MeteoSwiss station data, the calculated cooling effect for Otemma can be up to 2.5°C greater (more negative) or ∼1.5°C more negative when deriving Ta Est from ERA5-Land (Figures 3b and 3c).These data are presented here in order to highlight the challenges of representing local conditions under warming and their potential impact for the assessment of glacier cooling and its functioning.Nevertheless, the principal analysis conducted in this study is utilizing the most locally available Ta and wind data from the off-glacier T-Logger stations to minimize this forcing uncertainty.

Wind Regimes and Near-Surface Cooling
Observations on all glaciers indicate that up-valley winds erode the presence of a katabatic boundary layer and limit the relative Ta changes at the near-surface (temperature sensitivity-"TS"-Figure 4).For Otemma, A pronounced along-flowline cooling is observed for down-glacier (katabatic) wind conditions (Figure 4) and TS is <0.5 at the terminus.For Arolla, however, during down-glacier wind conditions TS is similar to that during up-glacier wind conditions (Figures 4a and 4c), possibly due to the entrainment of warm air over short distances and increasing presence of debris in the upper reaches of the glacier.At Corbassière, upper sites are found to be cooler than lower ones even under down-glacier wind conditions, though the gradient of TS with flowline distance is stronger for the occurrence of up-valley winds on that glacier (Figure 4c).The pattern of TS for Arolla is in contrast to earlier observations (see Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023), whereby the magnitudes of TS in 2001 (Figure 4a) are similar to those at comparable flowline distances on Corbassière in 2022.Representing these data on a normalized flowline further highlights the differences of the along-glacier TS between glaciers (Figure 4b), but demonstrates the clear role of up-valley winds in diminishing the down-glacier cooling toward the terminus of the glacier (Figure 4d).Changes in slope can modify the ratio of sensible heat exchange to adiabatic warming along a glacier flowline under katabatic conditions (Greuell & Böhm, 1998), assuming negligible influences due to phase change, divergence or entrainment.In our study, the average up-glacier slope angle appears to be only weakly related to TS on Arolla and Otemma, but inversely so on Corbassière (symbol sizes in Figure 4).Representing TS as a function of the combined normalized flowline and slope variations (Figure S9 in Supporting Information S1) produces only a vague relation, which is not itself conclusive.For example, the upper observations on Corbassière are located at the base of steep (∼30°) icefalls where a dominance of adiabatic heating over sensible heat exchange would be expected.
Corbassière experiences a more consistent down-glacier wind (up-glacier winds constitute only 18.8% of all hours) compared to Otemma (up-glacier winds for 31.2% of hours) and especially Arolla (63.6% of hours are classified with up-glacier winds).For Arolla, cooling is similar for up and down-glacier winds though the magnitudes of the cooling are often less than on the terminus of Otemma (Figure 5-"Katabatic" case study).Cooling is mostly diminished during up-glacier wind flow at Corbassière (Figures 5 and 6).For Arolla and Otemma, cooling of ∼1-2°C is still common for up-glacier wind conditions (some fetch over debris-free ice is still present), though cooling is stronger under down-glacier (negative UWI) conditions (Figure 5).Under disturbed conditions (e.g., neutral UWI on Corbassière-Figure 6c) the cooling effect is diminished, and strong gust speeds (interpreted from ERA5), also notably reduce any boundary layer cooling (Figure 5-"Disturbed" case study).When ambient temperatures are close to or below freezing, the cooling effect ceases (Figure 5), as is well documented in the literature (e.g., Ayala et al., 2015;Shaw et al., 2017).In such cases the mean cooling effect is typically within the uncertainty range of the sensors (Figure 6).
Utilizing the wind modeler WindNinja, we find an agreement with the patterns observed on-and off-glacier at Arolla and Corbassière under synoptically disturbed conditions (Figures 7a and 7e-timing indicated in Figure 5).For these glaciers, the terminus is orientated directly against the incoming synoptic flow, such that the boundary might easily be eroded.However, observations on Otemma (Figure 5, red arrows in Figure 7c) show a maintained down-glacier flow, implying an intrusion of synoptic flow from the exposed upper glacier areas which is not fully replicated by WindNinja.These wind model results indicate that there are only slack, lee-ward winds from the northwest that flow weakly across the lower glacier.Uncertainty of synoptic wind direction from the reanalysis data is a plausible reason for this discrepancy (Figures S10 and S11 in Supporting Information S1).
Varying the wind direction (∼40° toward the north) permits greater topographic channeling of this wind toward the lower glacier on Otemma and would then agree with the on-and off-glacier observations that show a coincident down-glacier flow and a lack of near-surface cooling (Figure S12 in Supporting Information S1).
Under katabatic conditions (the case study of the 6 September 2022-Figure 5), north facing glaciers (Arolla and Corbassière) align with a weaker synoptic flow from the south which is likely entrained within the boundary layer flow of the glacier (Figures 7b and 7f).On Arolla this represents a warm, southerly air mass with weak flow that likely advects more warm air from recently deglaciated, snow-free slopes and debris covered areas (e.g., Mott et al., 2015).Combined with the lack of sufficient fetch over which sensible heat can be exchanged with the bare ice surface, this may explain why mean cooling (Figures 5 and 6) is not stronger (more negative) and why TS is not smaller (Figure 4) under down-glacier, katabatic conditions.Modeled synoptic winds on Otemma are in fact stronger, though again in disagreement with observed wind directions, even in the glacier forefield (Figure 7d).In this case, synoptic winds are aligned with the valley axis but opposed to the down-glacier flow direction, showing that the presence of a strong glacier wind and significant cooling (as observed) can overcome the boundary layer erosion that might otherwise occur on a smaller glacier.

Glacier Cooling and Energy Balance
We see that valley winds, which have become more common on Arolla as its katabatic system has decayed over the last 20 years (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023), dominate the sensible heat fluxes at this site (79% of the warm afternoon subset-Figure 8a).Katabatic conditions occur alongside some of the warmest on-glacier observations on Arolla (Figure 8b), but often coincide with periods of weaker, southerly wind flow that likely limits total turbulent heat exchange (e.g., Figure 7b), providing around only 19% of the total sensible heat contribution to the glacier.Differences in sensible heat flux contributions on Arolla from katabatic and valley/synoptic winds are significantly different for the period analyzed (two tailed t-test p < 0.01).
The net sensible heat exchange contributions on Otemma are almost evenly divided between valley-winds and katabatic conditions, but this glacier also shows the largest role of synoptically disturbed conditions for all of the three glaciers (Figure 8c-differences of mean H are not significant).Katabatic conditions that generate the strongest cooling have the highest instantaneous sensible heat fluxes to the surface, largely driven by strong winds speeds in the boundary layer (Figures 8c and 8d), even though the total flux contributions for the subset are similar to those for valley-winds (47 vs. 45%) due to coincidence of both warm and windy conditions (Figure 8d).
Sensible heat fluxes on Corbassière are dominated by katabatic conditions (97% of total sensible heat flux to the surface), and even rare-up-valley wind intrusions were also associated with strong cooling (Figures 6 and 8e).
A few strong katabatic events produced instantaneous fluxes exceeding ∼70 W m −2 on Corbassière (Figure 8f), though sensible heat transfer was often similar or less than found on Otemma Glacier.
When also replacing observed on-glacier wind speeds with local off-glacier wind data (right hand panels of Figures 9 and 10), the model simulates relative decreases in sensible heat flux for Arolla (−0.1%) and Otemma especially (−22.2%-strongeron-glacier winds are measured than those off-glacier-Figures 9b and 9d).This is the inverse on average for Corbassière (+23.7%), with the exception of a few very strong on-glacier winds which are not represented by off-glacier conditions.Differences in glacier-wide melt rates also encompass other processes (e.g., antecedent albedo and surface temperature changes) which are not explored here, so interpretation of changes in melt due to the inclusion/exclusion of boundary meteorological information is more challenging.Exclusion of both glacier cooling and on-glacier wind speed results in a balance of increasing (up-valley winds) and decreasing (katabatic winds) modeled melt for the small Arolla Glacier compared to the reference meteorology.Overall increases in modeled melt are seen on Otemma and Corbassière, which is largely dominated by temperature increases with the exclusion of on-glacier meteorology (Figure 10, Table 1).Katabatic winds only act to balance the effect of cooling on melt when strong wind speeds (>4-5 m s −1 ) are developed (Figures 10d  and 10f).Otherwise, the presence of the glacier boundary layer meteorology acts to reduce the total melt during the summer period.

Dissecting the Role of the Glacier Boundary Layer
It has long since been established that the presence of a 0°C snow or ice surface impacts the representation of Ta above a glacier (Greuell & Böhm, 1998) and some glacio-hydrological studies have ignored or accounted for this in a crude manner (e.g., Huss & Hock, 2015;Jouberton et al., 2022;Rounce et al., 2023), either from attempts to represent the physical processes or through bulk optimization of an empirical model.Past observations have indicated that the flowline distance over which air parcels travel and interact with the glacier surface is highly influential in determining the magnitude of this so-called "cooling effect" (Ayala et al., 2015;Greuell & Böhm, 1998;Shea & Moore, 2010).Past studies have also attempted to model this behavior in local-scale energy balance simulations (Ayala et al., 2017;Carturan et al., 2015;Shaw et al., 2017).What is apparent from this study is that a glacier size and its flowline length alone cannot account for the apparent discrepancy between on and off-glacier Ta under the same warm conditions of the region (Figure 4), despite that glacier retreat can promote a loss of glacier boundary layer development (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).
Our results indicate a glacier's susceptibility to local up-valley slope winds (Figures 4-6) and synoptically driven gusts due to orientation and surrounding topography (Figures 5 and 6c) can significantly control the magnitude of the energy fluxes to the glacier surface (between 3% and 81%-Figure 8).
It has previously been hypothesized that wind interactions and varying wind origins may be responsible for the along-glacier variability of Ta during summer months (Ayala et al., 2015;Munro, 2006;Pellicciotti et al., 2008;Shaw et al., 2021).Moreover, it has been appropriately demonstrated through observational evidence (Shaw et al., 2021) and through detailed large-eddy simulations (Goger et al., 2022) that strong synoptic winds can erode the boundary layer cooling evident on mountain glaciers.The observational evidence presented here clearly conforms to this understanding (Figure 5), though the degree to which this erosion occurs can also be determined by the alignment of the glacier with the predominant syntopic gradient (Figure 7).Moreover, the relative warming on a glacier terminus may also be affected by warm air advection from valley winds (Shaw et al., 2021).This has been more recently demonstrated by Shaw, Buri, McCarthy, Miles, Ayala, and Pellicciotti (2023) using multi-year observations records on Arolla Glacier, highlighting the need to address potential non-linearities of glacier sensitivity to atmospheric temperature changes as mountain glaciers inevitably decay.On Hintereisferner (Austrian Alps) it has been demonstrated that the synoptic gradient can align with the down-glacier flow and limit the development of a low-level jet (Goger et al., 2022), but that exposure to strong crosswinds (flow perpendicular to glacier flowline) can also impact the katabatic flow and relative turbulent heat exchange (Mott et al., 2020).
Interestingly, for Otemma Glacier, the valley axis aligns with the occurrence of substantial southerly and southwesterly synoptic flow, which, in the absence of a glacier (i.e., the WindNinja simulations), would promote strong up-valley flows (e.g., Figure 7d).Observations on the glacier indicate that a strong glacier wind is sufficient to counteract this synoptic/valley flow into the glacier forefield, but suggest that under conditions that promote weaker glacier winds (cloudier and milder days), the role of valley winds could further inhibit cooling on the glacier terminus and produce greater sensible heat transfer to the ice surface through the combination of relatively warmer and windier weather (e.g., Figure 9).Moreover, on Otemma, warm synoptic gusts are evidently more common than for other sites, contributing more to the total sensible heat flux (Figure 8c) and suggesting that this glacier may be increasingly susceptible to a decaying boundary layer in the future, as compared to Corbassière, for example.Ultimately, these findings reinforce the need to account for local topographic conditions, notably the valley/glacier alignment to synoptic winds, when attempting to model the cooling of Ta above mountain glaciers (Ayala et al., 2015;Shea & Moore, 2010) because consideration of flowline distance alone over-simplifies inter-glacier differences (Shaw et al., 2021).
Our model results demonstrate that the differences in meteorological observations on-and off-glacier, even from those off-glacier stations very local to the glacier area (Figure 1), can produce substantial differences in heat flux and melt calculations on Alpine glaciers (Figures 9 and 10).It is very clear that ignoring Ta biases will lead to an over-estimation of melt if using off-glacier temperature forcing (e.g., Figure 9e), especially for melt models which do not overcome this forcing uncertainty through bulk calibration (i.e., degree-day models).However, the coincidence of this cooling effect with stronger down-glacier flow (Figure 8) may also act to compensate for melt differences in physically based models.While our analysis is certainly not exhaustive on this topic, we demonstrate that under warm conditions, the presence of a well-developed glacier boundary layer may act as either net gain (Otemma) or net reduction (Corbassière) in total sensible heat supply to the glacier, despite hours in which a strong katabatic wind layer clearly augments the heat supply to the ice at both sites (Figures 9d and 9f).This is clearly not the same for the decaying Arolla Glacier, however (Figure 9b) which, due to the diminishing boundary layer (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023), reduces the overall importance of local forcing uncertainties (Table 1).
As no past meteorological records on Otemma and Corbassière are known to exist, temporal changes of their near-surface boundary layer cannot be directly assessed, though given a space-for-time substitution one could hypothesize that those glaciers might suffer a similar fate as Arolla as they reduce in size, especially when they already appear to be susceptible to synoptic erosion (Otemma-Figures 7 and 8).In general, this re-emphasizes the need to address forcing uncertainties and biases on glaciers and attempt to represent them either in a physical or statistical framework that accounts for the potential non-linearities under future conditions (e.g., Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).

Implications of Glacier Cooling in a Future of Extreme Summers
It is estimated that during the 2022 ablation season, Swiss glaciers lost ∼6% of their total mass (SCNAT, 2022) and summer ablation rates were the most negative since records began, reflecting a combination of a dry antecedent winter and a hot and dry summer (Cremona et al., 2023).During the warmest point of the total instrumental period (18 July-outside of the data sub-period presented), temperatures reached 18°C at 2,600 m a.s.l. at the lower elevations of Corbassière Glacier and ∼17°C around Arolla and Otemma.Corresponding temperatures were 11.5°C meanwhile on the lower Corbassière terminus, indicating a large 6.5°C of cooling due the presence of the glacier.Figure 9 demonstrates the relevance of accounting for this effect for the accurate calculation of sensible heat fluxes and melt.Recent studies have also discovered the potential role of glacier winds to reduce maximum temperatures in pro-glacial regions (Lin et al., 2021;Salerno et al., 2023), therefore extending the hydrological relevance of glacier cooling beyond to the borders of the ice itself.
Future projections of climate indicate not only an increase in global mean temperature conditions, but also a potential increase in the frequency and duration of heatwaves (e.g., Brown, 2020;Schielicke & Pfahl, 2022).Sustained warm periods, even those occurring over the next decades may already invalidate the utility of static parameters describing temperature biases on mid-latitude glaciers.In long-term modeling studies, temperature biases are often calibrated (amongst other parameters) to fit mass balance observations of simpler models (e.g., Rounce et al., 2023) in order to maintain computationally feasible ensemble simulations.However, this approach assumes a stationarity of the underlying physical processes and suggests that the future relationship of glaciers and climate is linear.We know that this is not valid, either due to hysteresis of certain thermo-dynamical processes (e.g., Marshall, 2021), topographical feedbacks (Bolibar et al., 2022) or regime shifts of glacier boundary layers (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).With regards to the latter, we typically lack the knowledge of local meteorology, and therefore parameters describing temperature biases become less physically meaningful, and their values are often derived through an equifinality-affected optimization.Past temperature correction factors have typically been in the range of 0.5-2°C (Bash & Marshall, 2014;Jouberton et al., 2022;Rounce et al., 2023), sometimes having been derived as an average encompassing many conditions at just one location.Evidence from Figure 5 demonstrates how this is at the lower end of temperature biases for our study site during an extreme year, but that this expected bias also scales with the ambient temperature itself (Figures 3 and 6-e.g., Ayala et al., 2015).
As a simple experiment of the relevance of these processes, we model the location of each terminus AWS (Figure 1) for the data subset period (Figure S6 in Supporting Information S1), applying a Ta perturbation for each year until 2099, based upon the results of a pessimistic SSP5-8.5 CMIP6 ensemble for the region.We run two separate simulations (details are given in Supporting Information S1): (a) observed off-glacier wind speed and Ta, with the CMIP6 summer temperature, humidity, wind and precipitation delta factors relative to 2022 applied to each year ("OFF") and; (b) observed off-glacier meteorology as in (a), but where cooling is applied following the relation of cooling to Ta Est as in Figure 3 and wind speeds are adjusted based upon a linear scaling relation against the observed cooling (Figure S13 in Supporting Information S1-"ON").Analysis of these model runs show how the presence of glacier boundary layer on Otemma and Corbassière creates a reduction in total melt in the near-future (10-30 years), though increasing katabatic wind strength associated with enhanced atmospheric warming can heighten total melt toward the end of the century, varying substantially between glaciers (Figure 11).These back-of-envelope calculations clearly ignore other processes relevant to the above-glacier meteorology, most notably the retreat and decay of the glaciers and resultant glacier fetch (note that Arolla will likely disappear long before the end of the century).Nevertheless, such analysis provides a general indication of how these systems are highly dynamic and non-linear, warranting further investigations into their role under a future climate.

Future Research Directions
The varying role of glacier cooling and katabatic winds on glacier energy balance remains underexplored, though this study has highlighted how variable its importance can be within one region.Our study presents an indication of the melt rates and thus significance of boundary layer processes under a short observation period, but is not currently designed to address the long-term perspective of forcing uncertainty stemming from these processes.Future research should address the long-term implications of the changes in glacier boundary layer meteorology (e.g., Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023) and contextualize their relevance against current practices to model long-term glacier evolution and hydrological contributions using simpler model frameworks.
No one study is likely to provide both a fully distributed and continuous series of on-glacier observations to exhaustively answer how along-and across-glacier micrometeorology can affect Ta variability and sensible heat flux.However, past studies have been able to dissect how synoptic conditions influence boundary layer development (Conway et al., 2021;Goger et al., 2022;Mott et al., 2020), how regional climate can influence the magnitude of Ta biases on glaciers (Yang et al., 2022) and how these biases related to boundary layer conditions can vary in space (Ayala et al., 2015;Carturan et al., 2015;Shaw et al., 2021;Shea & Moore, 2010) and time (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023).Therefore, we have already gained an improved process understanding of many elements related to on-glacier temperature estimation, though a generalized methodology to account for these uncertainties in glacier melt models remains to be developed.
Our study is able to quantify the degree to which localized processes, principally those related to wind regimes, control Ta variability and energy flux on distinct glaciers in the same region.By standardizing our instrumental setup as much as possible, we show how variable cooling effects are in space and time under the same predominant climatic conditions (Figures 4 and 5), highlighting how glacier orientation and surrounding topography can strongly govern a glacier's susceptibility to external conditions (Figures 7 and 8).Superimposed on this are the reductions of heat exchanges with ice surfaces due to the changing glacier area (Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023) which likely enhances the relative contribution of heat advection under slope or synoptically driven winds (e.g., Mott et al., 2015) and which can disproportionately impact the temperature sensitivity of smaller glaciers (Figure 4) to prevailing synoptic winds of the same direction (Figures 7b and 7f).
It has been demonstrated that high resolution atmospheric or large eddy simulations can reveal detailed information about the boundary layer processes in question (Goger et al., 2022;Lin et al., 2021;Potter et al., 2018;Sauter & Galos, 2016).Moving forward, it may be possible to link the information derived from these detailed, albeit very short simulations to the climatic conditions under which they occur in order to aid the parameterization of glacier cooling (Figure 3) or temperature sensitivity (Figure 4) from widely available and computationally "light" data sets (e.g., reanalysis).Ongoing work seeks to blend this information with the wealth of previously gathered on-glacier data to generalize the expected biases for glacier melt models.Nevertheless, more modeling exercises are required to understand the implications of the different model complexities in representing these processes and whether they can be reconciled into an adaptive framework of climate downscaling and Ta estimation on glaciers.

Conclusions
Our study combines a consistent summer monitoring campaign of on-glacier meteorology with regional climate data sets to dissect the occurrence and relevance of the boundary layer on three glaciers in the same region of Switzerland.Under the extreme summer of 2022, observed cooling of the near-surface air above the glacier ranged between 1 and 7°C on warm afternoons, varying by the amount of heat in the atmosphere and by the susceptibility of the glacier to erosion by synoptic wind gradients, even under warm temperatures.We emphasize that the fetch over which air parcels travel down-glacier, thus the size of a glacier, alone cannot determine the amount of near-surface cooling that occurs, and that local factors related to wind speeds and orientation relative to syntopic winds can strongly modulate a glacier's response to climate.
By focusing on the deficits of air temperature at the near surface and the role of local off-and on-glacier winds, we demonstrate a heterogeneous role of warm, up-valley winds in determining sensible heat flux to the ice surface (ranging between 3% and 81% of the total under warm atmospheric conditions during afternoon hours).Meanwhile, the presence of katabatic winds can, depending on their overall strength and persistence, result in a net increase (Corbassière) or a net decrease (Otemma) in sensible heat transfer to the glacier surface compared to modeled results using local, off-glacier meteorology.Ignoring the temperature biases due to the boundary layer cooling effect has clear consequences for heat flux and melt calculations, but the presence of stronger glacier winds can partially compensate for this in physically-based model simulations.As glaciers shrink and fragment (Arolla Glacier), the uncertainties related to boundary layer meteorology become less relevant.Nevertheless, we suggest that any modeling framework should consider the potential for non-linearities in glacier response to future climate as glaciers shrink and decay, not only because of relative temperature changes in the boundary layer (relevant to simpler models), but also due to the related changes in density-driven katabatic wind speeds.

Figure 1 .
Figure 1.The location of the study glaciers (upper panel) with on-glacier T-Logger (dots), automatic weather stations (AWS) (stars) and off-glacier T-Logger (squares) locations provided.The location of the domain within Switzerland is indicated by the inset map outline.The lower panel indicates the total flowline distance (m) from the respective glacier headwall or summit, marking the relative locations of the stations.The Digital Elevation Model (DEM) is provided by ASTER GDEM.

Figure 2 .
Figure 2. The study domain with the grid of mean ERA5-Land Ta (°C) for the period of interest (panel a).The mean Ta against elevation is shown in panel (b), comparing raw ERA5-Land (black dots), MeteoSwiss stations (green squares) and study-specific on-(circles) and off-glacier (triangles) stations.Small colored dots indicate the equivalent ERA5Land pixels that intersect each glacier and dashed black and green lines indicate the linear fit (Ta gradient) of all ERA5-Land points and MeteoSwiss stations against elevation, respectively.

Figure 3 .
Figure 3.The observed glacier cooling for intervals of Ta Obs for different data sources.The cooling (Ta Obs − Ta Est ) is an average of all available stations on each glacier (colors) for the interval of Ta Est (x-axis).The shaded areas indicate one standard deviation of Ta Obs − Ta Est for each interval.

Figure 4 .
Figure 4. Temperature sensitivity (TS-Equation 2) along the flowline for warm conditions (daytime TaOFF >50th percentile), sub-divided into (a, b) katabatic conditions (up-glacier wind index [UWI] < −0.2) and (c, d) up-valley conditions (UWI > 0.2).Circle/star sizes are scaled by the mean up-glacier slope angle (derived from the ASTER GDEM) and the inset number indicates the total number of hours in the condition.Panels (b) and (d) show the same results as (a) and (c), though normalizing the stations' distance along flowline against the total glacier length.Values for Arolla in 2001 (see Shaw, Buri, McCarthy, Miles, Ayala, & Pellicciotti, 2023) are added for comparison.Data for 2022 on Arolla Glacier are also published in Shaw, Buri, McCarthy, Miles, Ayala and Pellicciotti (2023), though express small differences in the analysis sub-period.

Figure 5 .
Figure 5. Mean afternoon cooling effect (Ta Obs − Ta Est ) at the lower automatic weather stations on each glacier (center) in comparison to the mean Ta over the region from ERA5-Land (top color bar), mean ERA5-Land wind direction (arrows), mean ERA5 gust wind speed (black line) and on-glacier up-glacier wind index (lower panel).Shaded areas for the cooling effect indicate the variability based upon the cooling at the lowermost (3) stations on each glacier.Dashed boxes highlight two case studies that are explored for subsequent analysis, related to synoptically disturbed and strong katabatic conditions.Afternoon hours are considered from 12:00-19:00.

Figure 6 .
Figure 6.Hourly values of glacier cooling (Ta Obs − Ta Est ) at the lower automatic weather stations AWS on each glacier plotted against the raw value of TaOFF.Scatter points are scaled by the measured on-glacier wind speed at the AWS and colored by the equivalent up-glacier wind index value (such that blue is down-glacier wind).Dashed horizontal lines indicate the uncertainty band of 0.4°C (see main text) within which differences between Ta Obs and Ta Est have a high probability of being due to measurement error.Inset radial plots indicate the orientation of the up-and down-glacier winds, whereby north is directed upward (Scales not added for neatness).

Figure 7 .
Figure 7. Modeled wind vectors derived from WindNinja for the synoptically disturbed (a, c, e) and katabatic (b, d, f) case studies (see Figure 5) for Arolla (a, b), Otemma (c, d) and Corbassière (e, f).Glacier outlines and station locations are depicted, with observed wind speeds and directions for the given hours shown by the red arrows (Wind speed for the off-glacier station of Corbassière is assumed equal to the on-glacier wind speed for illustration purposes).The color scale indicates the estimated temperature using the local Ta OFF and the mean cooling effect is applied over all glacier cells.

Figure 8 .
Figure 8.The calculated, glacier-wide sensible heat flux (H) for warm afternoon hours, plotted against the cooling effect (left panels) and the observed mean Ta (right panels).Hours are classified into katabatic (blue circles), up-valley (red squares) or synoptically disturbed (crosses) conditions (see main text) and scaled by the observed wind speeds.The distributions of H for the three conditions are given by the right-hand boxplots and their contributions to the total H for the selected warm afternoon hours (all values are positive-toward the surface) are shown by the left-hand pie charts.

Figure 9 .
Figure 9. Calculated differences in sensible heat flux (H) between the reference ("REF") model scenario (Figure 8) and (i) a scenario using Ta Est and no cooling effect applied ("COOL"-left panels) or (ii) a scenario where no cooling effect or local glacier wind speed data are provided ("WIND"-right panels).The results are presented as COOL-REF (left) and WIND-REF (right), such that a positive H difference indicates less sensible heat transfer to the glacier surface under the reference scenario.

Figure 10 .
Figure 10.Calculated differences in hourly melt (m w.e.) between the reference ("REF") model scenario (Figure 8) and (i) a scenario using Ta Est and no cooling effect applied ("COOL"-left panels) or (ii) a scenario where no cooling effect or local glacier wind speed data are provided ("WIND"-right panels).The results are presented as COOL-REF (left) and WIND-REF (right), such that a positive melt difference indicates less surface melt under the reference scenario.

Figure 11 .
Figure11.Total melt (11 August-19 September period) with applied perturbations to the near-surface meteorology based upon a delta change of Ta and meteorology from CMIP6 ensemble simulations ("OFF") and with observed linear relationships to mean cooling and enhanced katabatic wind speeds under those warming scenarios ("ON").Model simulations are based upon a simplification of many processes and only suggestive of possible non-linearities in the future of these glaciers (see Supporting Information S1 for details).