Surface energy fluxes on Chilean glaciers: measurements and models

The surface energy fluxes of glaciers determine surface melt, and their adequate parametrization is one of the keys to a successful prediction of future glacier mass balance and freshwater discharge. Chile hosts glaciers in a large range of latitudes under contrasting climatic settings: from 18^∘ S in the Atacama Desert to 55^∘ S on Tierra del Fuego. Using three different methods, we computed surface energy fluxes for five glaciers which represent the main glaciological zones of Chile. We found the main energy sources for surface melt change from the Central Andes, where the net shortwave radiation is driving the melt, to Patagonia, where the turbulent fluxes are an important source of energy. We inferred higher surface melt rates for Patagonian glaciers as compared to the glaciers of the Central Andes due to a higher contribution of the turbulent sensible heat flux, less negative net longwave radiation and a positive contribution of the turbulent latent heat flux. The variability in the atmospheric emissivity was high and not able to be explained exclusively by the variability in the inferred cloud cover. The influence of the stability correction and the roughness length on the magnitude of the turbulent fluxes in the different climate settings was examined. We conclude that, when working towards physical melt models, it is not sufficient to use the observed melt as a measure of model performance; the model parametrizations of individual components of the energy balance have to be validated individually against measurements.

1 Introduction

Glaciers are retreating and thinning in nearly all parts of the planet, and it is expected that these processes are going to continue under the projections of global warming (IPCC, 2019). For mountain glaciers melt is mostly determined by the energy exchange with the atmosphere at their surfaces. The processes leading to this exchange of energy are complex and depend on the detailed (micro)climate on the glacier. Classical empirical melt models like for example degree day models (Braithwaite, 1995a) are being replaced more and more by more complex models which try to quantify the detailed physical processes that govern the energy exchange at the glacier surface. These kind of models are sometimes called “physical melt models” or “physically based models” (Pellicciotti et al., 2008).

In Chile, the only glacier with a climatologically relevant long-term record of surface mass balance is Echaurren Norte Glacier near to Santiago de Chile, which has been monitored since 1975 (WGMS, 2017; Masiokas et al., 2016; Farías-Barahona et al., 2019). Echaurren Norte Glacier (33.5^∘ S) has a general negative trend in its cumulative surface mass balance (−0.48 m w.e. yr⁻¹ in 1976–2017) but also showed stable phases in the 1980s and the first decade of the 21st century (Masiokas et al., 2016; WGMS, 2017; Farías-Barahona et al., 2019). The variations in the surface mass balance of this glacier can be mostly explained by variations in precipitation in the region (Masiokas et al., 2016; Farías-Barahona et al., 2019). In the semiarid Pascua Lama region (29^∘ S) several small glaciers have been monitored since 2003 (Rabatel et al., 2011). These glaciers also show mostly negative surface mass balance and are losing area (Rabatel et al., 2011). During the monitoring period, the limited accumulation of snow was not able to balance out the ablation which was dominated by sublimation (MacDonell et al., 2013). In the Chilean Lake District, Mocho Glacier has been monitored since 2003 (Rivera et al., 2005). Here, a very high interannual variability in the surface mass balance was observed (Schaefer et al., 2017). But, on average, the annual surface mass balance was negative which coincides with the observed areal losses (Rivera et al., 2005).

Figure 1 Middle: glaciological zones in Chile according to Lliboutry (1998) and locations of the studied glaciers: (a) Bello Glacier, (b) Pirámide Glacier, (c) San Francisco Glacier, (d) Mocho Glacier, (e) Exploradores Glacier (background: SRTM), (f) Tyndall Glacier. Background images courtesy of DGA and USGS.

Energy balance studies have been realized in Chile on different glaciers; in the semiarid Andes, MacDonell et al. (2013) quantified in detail the drivers of ablation processes on Guanaco Glacier (29^∘ S). They found that the net shortwave radiation is the main source and the net longwave radiation and turbulent flux of latent heat are the main sinks of energy at the surface of Guanaco Glacier (MacDonell et al., 2013). Due to the low temperatures at this high-elevation site (5324 m a.s.l.), they found that sublimation dominated the surface ablation and that surface melt contributed only during summer. Pellicciotti et al. (2008) and Ayala et al. (2017b) studied the surface energy balance during summer at Juncal Norte Glacier in the Central Andes (33^∘ S, near Santiago de Chile). Similarly to MacDonell et al. (2013) they found that the net shortwave radiation is the main source and that the net longwave radiation and turbulent flux of latent heat are the main sinks of energy. Similar results concerning the influence of the different components of the surface energy balance were obtained by Ayala et al. (2017a), who analyzed meteorological data collected on six glaciers in the semiarid Andes of north-central Chile at elevations spanning 3127 to 5324 m a.s.l.

Brock et al. (2007) studied the surface energy balance of bare snow and tephra-covered ice on Pichillancahue-Turbio Glacier (39.5^∘ S) on Villarrica volcano in the Chilean Lake District during two summers. They found a strong reduction in surface melt on the tephra-covered part of the glacier and a change in sign of the turbulent flux of latent energy to a source due to the higher vapor pressure caused by a more humid atmosphere as compared to the northern and central part of Chile. In southernmost Chile, Schneider et al. (2007) studied the energy balance in the ablation area of Lengua Glacier, which is an outlet glacier of Gran Campo Nevado ice cap (53^∘ S). They found that during February to April 2000, due to the high air temperatures and the high wind speeds, turbulent flux of sensible heat was the main source of melt energy for the glacier surface.

In a comparative study of the surface energy balance of glaciers at different latitudes, Sicart et al. (2008) found that the net shortwave radiation is driving the glacier melt at the tropical Zongo Glacier but that at Storglaciären in northern Sweden the turbulent fluxes of sensible heat and latent heat dominated the melt patterns.

In this study we analyze data from automatic weather stations (AWSs) installed on six glaciers distributed in the different glaciological zones of Chile (Fig. 1), five of them being equipped and maintained by the Unit of Glaciology and Snow of the Chilean Water Directorate (UGN-DGA; UChile, 2012; Geoestudios, 2013; CEAZA, 2015). Using the meteorological observations as input, we compare different ways to compute the glacier surface energy balance. We use direct measurements of the radiative fluxes at the glacier surface and two models that are freely available: the spreadsheet-based point surface energy balance model (EB model) developed by Brock and Arnold (2000) and the COupled Snowpack and Ice surface energy and MAss balance model (COSIMA) (Huintjes et al., 2015b, a).

Instead of validating the ability of the energy balance calculations to adequately predict melt rates, in this study we want to test their ability to reproduce the individual energy fluxes. We want to emphasize the differences between the model parametrizations and their ability to reproduce the directly measured radiative fluxes at the glacier surfaces. We also compare different parametrizations for the turbulent fluxes of sensible and latent heat and discuss the influence of stability corrections and roughness lengths.

2 Sites

The projections of future changes in climate depend on the different climatological and glaciological zones. This is why a detailed analysis of the processes that determine the energy exchange at the surface of the glaciers in the different climatological zones is necessary to be able to make reliable predictions of future surface mass balance and meltwater discharge of Chilean glaciers.

Chile's climate is strongly determined by the Pacific anticyclone and the Andes range which acts as a natural barrier (Fuenzalida-Ponce, 1971; Garreaud, 2009). Large climate differences are observed due to the large north–south extent of the territory (4000 km, 17^∘30–55^∘ S). Despite the different classifications of subglaciological zones (Lliboutry, 1998; Masiokas et al., 2009; Barcaza et al., 2017; Braun et al., 2019; Dussaillant et al., 2019), most authors agree that there is a transition from Dry Andes to Wet Andes at around 35^∘ S (Fig. 1). The Central Andes of Chile (31–35^∘ S) are characterized by a Mediterranean climate, with dry conditions during summer. For the period 1979–2006 Falvey and Garreaud (2009) observed a cooling at the coast and a considerable temperature increase of +0.25 ^∘C per decade inland in the Maipo River catchment in the Central Andes. Precipitation in this area is highly variable and predominantly occurs during winter (Falvey and Garreaud, 2007) controlled by El Niño–Southern Oscillation (ENSO) and the southeast Pacific anticyclone (Montecinos and Aceituno, 2003). Between 2010 and 2015 a megadrought was observed in the Central Andes (Boisier et al., 2016; Garreaud et al., 2017).

In the northern part of the Wet Andes (35–45^∘ S), known in Chile as the Lake District, the elevation range steadily decreases and wetter climatic conditions are predominant. A general decrease in precipitation in the region was observed during the 20th century (Bown et al., 2007; González-Reyes and Muñoz, 2013). The southern part of the Wet Andes, Patagonia, is characterized by a hyperhumid climate (Garreaud, 2018), where the largest glacierized areas in the Southern Hemisphere outside Antarctica can be found. These hyperhumid conditions were recently interrupted by a severe drought during 2016 with a precipitation decrease of more than 50 % (Garreaud, 2018). Under these different climatic settings, Chile hosts the majority of glaciers in South America (more than 80 % of the area), which have been mostly thinning and retreating in the last few decades (e.g.,  Braun et al., 2019; Dussaillant et al., 2019).

In the Central Andes, San Francisco (1.5 km²) and Bello (4.2 km²) glaciers are mountain glaciers which are partially debris-covered at their termini and Pirámide Glacier (4.4 km²) is an almost completely debris-covered glacier (Fig. 1). On San Francisco and Bello glaciers the AWSs were installed over bare ice, and at Pirámide Glacier they were installed over debris cover. Mocho Glacier is part of the ice cap (14 km²) which covers the Mocho-Choshuenco volcanic complex, located in the Lake District (Schaefer et al., 2017). Exploradores Glacier (83.8 km²) is located on the northern margin of the Northern Patagonian Ice Field with a prominent portion of debris cover at its tongue. Recently, at the glacier's front, several lateral lakes have developed and some calving activity has been observed. Finally, Tyndall Glacier (309.8 km²) is one of the large glaciers in the southeastern part of the Southern Patagonian Ice Field. Tyndall Glacier terminates in Geike Lake, where it experiences additional mass losses through calving. All glacier areas are from Barcaza et al. (2017).

Due to the installation of AWSs on several glaciers in the country by the UGN-DGA, detailed meteorological observations from glaciers in the different glaciological zones are now available (UChile, 2012; Geoestudios, 2013; CEAZA, 2015). In Table 1 we present the detailed locations of the AWSs used for this study and some relevant glacier parameters.

(Ayala et al., 2016)(Ayala et al., 2016)(Carrasco et al., 2008)(Schaefer et al., 2017)(Schaefer et al., 2013)(Schaefer et al., 2015)

Table 1 Study period and geographical information of the glaciers and automatic weather stations. ELA is equilibrium line altitude.

Because of their higher relevance for melt modeling, we focused our analysis on summer periods. AWSs of the UGN-DGA have a data record of several years, but during several summers some of the sensors did not work well. In Table 1 we show the selected summer period for every station. For Tyndall Glacier two summers were analyzed. On Mocho Glacier an AWS was installed during only a 50 d period in summer 2006. Figure 2 shows photos from the AWSs installed on the glaciers Bello, Exploradores and Tyndall.

Figure 2 AWS on Bello Glacier (a), Exploradores Glacier (b) and Tyndall Glacier (c).

3 Methods

In this contribution we want to focus on the six most important energy fluxes which normally determine the melt energy available at the glacier surface: the incoming solar radiation (SW_in), the reflected solar radiation (SW_out), the incoming atmospheric longwave radiation (LW_in), the longwave radiation emitted by the glacier surface (LW_out), and turbulent fluxes of sensible heat (SH) and latent heat (LH). We will compare three methods to compute these surface energy fluxes for the selected sites. Before presenting the three methods in detail, we describe the input data for our energy balance calculation in the following section.

3.1 Input data

In Table 2 we present the sensors used at the Mocho AWS in 2006 and the instruments used at the DGA stations on the other glaciers. The main difference in the installation is the CNR4 sensor, which is installed on the DGA stations and provides detailed measurements of all radiative fluxes, while on Mocho Glacier, SW_in, SW_out and the net all-wave radiation is measured. At Mocho AWS mean values of the data were recorded every 15 min, which were resampled to hourly data for the energy balance calculations. At the DGA stations hourly means are recorded and transmitted by a satellite connection. Data at missing hours were interpolated by taking the mean value of the hour before and after the missing one. For Bello Glacier the time resolution of the acquired data changed from hourly to 3 hourly on 20 March 2015. Hourly data were generated using the linear interp1 MATLAB interpolation scheme.

Table 2 Sensors employed at the different AWSs.

3.2 Reference Database

We call this first method the Reference Database, since in this approach direct measurements of the first four fluxes (the radiative fluxes) are used (with the exception of the Mocho AWS where the net longwave radiative flux is inferred from the incoming solar radiation, the reflected solar radiation and the overall net radiative flux). However, on several glaciers the measured mean outgoing longwave radiative fluxes were higher than the ones expected for a blackbody at 0 ^∘C. Therefore we decided to bias correct the measured longwave radiative fluxes in a way that the measured outgoing longwave radiative fluxes in the afternoon correspond to a melting surface at 0 ^∘C. From this calibration of the signal we obtained different correction factors for the longwave radiative fluxes which were applied to both incoming and outgoing longwave radiative fluxes (Table 4).

The turbulent fluxes of sensible and latent heat were calculated according to formulas derived in Cuffey and Paterson (2010). The bulk aerodynamic approach is employed and the following important assumptions are made.

1. The eddy diffusivity for heat has the same value as the eddy diffusivity for water vapor and the eddy viscosity.

2. The shear stress in the first few meters of the atmosphere above the glacier surface is constant.

3. The wind velocity, temperature and water vapor pressure have logarithmic profiles with the same scaling length z₀.

Using these assumptions, the following expression for the turbulent flux of sensible heat can be derived:

$$\begin{array}{}\text{(1)}& \mathrm{SH}=c_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{a}}{C}^{\ast }\left(z\right)U\left(z\right)\left[T\left(z\right)-T\left(s\right)\right],\end{array}$$

where c_a is the specific heat of air at constant pressure which was assumed to be constant at 1.01 kJ (kg K)⁻¹, ρ_a is the air density, U(z) is the wind speed measured at the height z above the surface, T(z) is the air temperature at the height z of the sensor and T(s) is the temperature of the glacier–atmosphere interface.

The dimensionless number C^∗(z) is a proportionality constant called the transfer coefficient. If the above assumptions are fulfilled, it should depend on the measurement height of the sensors of wind velocity and temperature z (2 m in our case) and the roughness length z₀ according to the following expression:

$$\begin{array}{}\text{(2)}& C^{\ast }\left(z\right)={\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{{\ln }^{\mathrm{2}}(z/z_{\mathrm{0}})}},\end{array}$$

where κ is the von Kármán constant, which has an approximate value of 0.4. In practice however the roughness and scaling length z₀ is variable in space and time (Brock et al., 2006). There exist several recommendations in the literature of values for C^∗(z) that have produced satisfying results (Cuffey and Paterson, 2010), which give C^∗(z) the quality of a tuning parameter rather than a physical constant. In the results section we present the results obtained by using an intermediate roughness length of z₀=0.5 mm for all the glaciers. According to Table 5.4 in Cuffey and Paterson (2010), this corresponds to a value between smooth ice and ice in the ablation zone and is also inside the range recommended for new and polar snow. When looking at the glacier surfaces in Fig. 2, we can note that the roughness length probably varies from glacier to glacier, and in Table 3 we present how C^∗(2 m) should vary for typical values of z₀.

Table 3 Variation in the transfer coefficient C^∗(2 m) for typical values of the roughness length z₀.

Assumption 1 is normally fulfilled for a neutral atmosphere, an atmosphere which exhibits a temperature lapse rate equal to the dry adiabatic lapse rate of 9.8 ^∘C km⁻¹. Over a glacier surface, the temperature gradient is often inverted (temperature increases with elevation), especially during summer when the air temperature is positive. This stable layering of air masses reduces the vertical exchange especially for low wind speeds, and we apply a stability correction to Eq. (1) which depends on the bulk Richardson number $Ri=\frac{gT\left(\mathrm{2}\mathrm{m}\right)\mathrm{2}\mathrm{m}}{\left(T\right(\mathrm{2}\mathrm{m})+\mathrm{273.15})U^{\mathrm{2}}}$, where g is gravitational acceleration and T(2 m) has to be taken in degrees Celsius. The correction factor is smaller than 1 for Ri>0.01 (small wind speeds) and approaches zero for Ri=0.2 in the same way as it is implemented in the COSIMA model (Huintjes et al., 2015a).

Using the same arguments from above and assuming the turbulent flux of latent heat to be proportional to the difference in the concentration of water vapor at the glacier surface and the air layer above it, Cuffey and Paterson (2010) derive the following expression:

$$\begin{array}{}\text{(3)}& \mathrm{LH}=\mathrm{0.622}{\mathit{\rho }}_{\mathrm{a}}{L}_{\mathrm{v}}{C}^{\ast }U\left(z\right)\left[P_{\text{vap}}\left(z\right)-P_{\text{vap}}\left(s\right)\right]/P_{\mathrm{a}}.\end{array}$$

Here, P_vap(z) and P_vap(s) are the water vapor pressure at the elevation z=2 m above the glacier and at its surface, respectively; P_a is the air pressure; and L_v is the latent heat of vaporization. The water vapor pressure at z=2 m depends on the (measured) relative humidity and the saturation water vapor pressure P_vap,sat which depends on the air temperature. In Fig. 3a we show measurements of the saturation water vapor pressure at different temperatures (Lide, 2004) and the graphs of several parametrizations as a function of the air temperature, found in the literature (Bolton, 1980; Cuffey and Paterson, 2010; Huintjes et al., 2015a).

Figure 3 (a) Different parametrizations of the saturation vapor pressure as a function of the air temperature found in the literature; we chose the parametrization of Bolton (1980). (b) Emissivity of the atmosphere as a function of the cloud cover as implemented in the EB model and COSIMA.

We decided to use the parametrization proposed in Bolton (1980), since it agrees best with the measurements:

$$\begin{array}{}\text{(4)}& P_{\text{vap,sat}}\left(T\right)=\mathrm{6.112}\exp \left({\displaystyle \frac{\mathrm{17.67}T}{T+\mathrm{243.5}}}\right),\end{array}$$

where P_vap,sat is in hectopascals and the air temperature T is in degrees Celsius. It is assumed that at the glacier surface the water vapor pressure is equal to the saturation vapor pressure at the surface temperature T(s). The turbulent flux of latent heat is corrected in the same way for stability conditions found over glacier surfaces as the turbulent flux of sensible heat (see above).

3.3 EB model

In the spreadsheet-based energy balance model developed by Brock and Arnold (2000) the incoming solar radiation, 2 m air temperature, wind speed and water vapor pressure are the meteorological input variables. The fixed input parameters are the latitude, longitude and elevation of the station; the aspect and slope; the albedo α; and the roughness of the surface z₀. The net shortwave radiation is calculated by multiplying the sum of the direct and diffuse incoming solar radiation by 1−α. The incoming direct and diffuse incoming solar radiation depend on the measured incoming solar radiation SW_in and the glacier's surface slope and aspect at the AWS. However, in our study, these parametrizations produced erroneous values for the net shortwave radiation SW_net in the late afternoon for several glaciers. This is why we computed SW_net from the measured incoming solar radiation SW_in by multiplying it with 1−α:

$$\begin{array}{}\text{(5)}& {\mathrm{SW}}_{\text{net}}=(\mathrm{1}-\mathit{\alpha }){\mathrm{SW}}_{\text{in}},\end{array}$$

where the albedo α is assumed to be a constant, which depends on the characteristics of the glacier surface. This parametrization should exactly agree with the parametrizations proposed in Brock and Arnold (2000) for flat surfaces (zero slope), which should be a very good approximation, since the AWSs are normally placed on flat terrain.

The net longwave radiation is computed by assuming that the snow–ice surface irradiates thermal radiation of a blackbody at 273.15 K (0 ^∘C) which is 315.6 W m⁻² according to the Stefan–Boltzmann law of thermal radiation. This value is subtracted from the incoming longwave radiation from the atmosphere which is computed with the Stefan–Boltzmann law as well, using T(2 m) and an atmospheric emissivity which is a function of the cloud cover. Brock and Arnold (2000) use a parametrization of the atmospheric emissivity ε which increases linearly as a function of the cloudiness n:

$$\begin{array}{}\text{(6)}& {\mathit{\epsilon }}_{\text{(EB)}}(n,T)=(\mathrm{1}+\mathrm{0.26}n){\mathit{\epsilon }}_{\text{cs}}\left(T\right),\end{array}$$

where ε_cs(T) is the clear-sky emissivity which depends on the air temperature T: ε_cs(T)=0.00877T^0.788 (T in Kelvin). The green line in Fig. 3b shows the graph of ε_(EB)(n,T) at 5 ^∘C. The cloudiness is inferred by comparing the theoretically site-specific clear-sky incoming solar radiation with the measured incoming solar radiation.

Similar to the Reference Database, in the EB model the turbulent fluxes of latent and sensible heat are calculated by expressions derived from the bulk aerodynamic method (Brock and Arnold, 2000) according to Eqs. (1) and (3). However, in this model the transport coefficients for the sensible heat flux $C_{\mathrm{EB}\mathrm{1}}^{\ast }$ and latent heat flux $C_{\mathrm{EB}\mathrm{2}}^{\ast }$ have a more complex form and depend not only on the roughness length z₀ but also on the Monin–Obukhov length scale L and the scaling lengths of temperature z_T and humidity z_H (Brock and Arnold, 2000):

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{EB}\mathrm{1}}^{\ast }={\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{\left(\ln \right(z/z_{\mathrm{0}})+\mathrm{5}z/L)\left(\ln \right(z/z_T)+\mathrm{5}z/L)}},\text{(8)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{EB}\mathrm{2}}^{\ast }={\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{\left(\ln \right(z/z_{\mathrm{0}})+\mathrm{5}z/L)\left(\ln \right(z/z_{\mathrm{H}})+\mathrm{5}z/L)}}.\end{array}$$

z_H and z_T are calculated as a function of z₀ and the roughness Reynolds number (Brock and Arnold, 2000; Andreas, 1987). Since normally $z_{\mathrm{H}}<z_T<z_{\mathrm{0}}$ (Ben W. Brock, personal communication, 2018) and L>0, $C_{\mathrm{EB}\mathrm{1}}^{\ast }$ and $C_{\mathrm{EB}\mathrm{2}}^{\ast }$ are smaller than C^∗ and $C_{\mathrm{EB}\mathrm{2}}^{\ast }<C_{\mathrm{EB}\mathrm{1}}^{\ast }$.

3.4 COSIMA

The COupled Snowpack and Ice surface energy and MAss balance model (COSIMA) was developed at RWTH Aachen University (Huintjes et al., 2015a, b) and combines a surface energy balance model with a multilayer subsurface snow and ice model to compute glacier mass balance (Huintjes et al., 2015a). In this work we want to focus on how COSIMA models the six dominant energy fluxes at the glacier surface. The input parameters for COSIMA are the incoming solar radiation (SW_in), the 2 m air temperature (T), the relative humidity (RH), the wind speed (U), the solid precipitation(P_s), the initial snow height, the air pressure (P) and the cloud cover (n; Huintjes et al., 2015a). The daily mean cloud cover over the glacier was estimated by comparing the measured SW_in with the theoretical, site-specific clear-sky radiation computed by a code developed by Corripio (2003). The cloud cover was determined from this cloud transmissivity τ_cl by solving the equation proposed in Greuell et al. (1997):

$$\begin{array}{}\text{(9)}& {\mathit{\tau }}_{\text{cl}}=\mathrm{1}-\mathrm{0.233}n-\mathrm{0.415}{n}^{\mathrm{2}}.\end{array}$$

The net solar radiation is calculated using Eq. (5). In contrast to the EB model, the albedo is variable and depends on the time since the last snowfall t_snow and the thickness h of the snow or firn layer on top of the glacier ice:

$$\begin{array}{}\text{(10)}& \mathit{\alpha }={\mathit{\alpha }}_{\text{snow}}+\left({\mathit{\alpha }}_{\text{ice}}-{\mathit{\alpha }}_{\text{snow}}\right)\exp (-h/d^{*}),\end{array}$$

where α_ice and d^* are constants and

$$\begin{array}{}\text{(11)}& {\mathit{\alpha }}_{\text{snow}}={\mathit{\alpha }}_{\text{firn}}+\left({\mathit{\alpha }}_{\text{frsnow}}-{\mathit{\alpha }}_{\text{firn}}\right)\exp (t_{\text{snow}}/t^{*}),\end{array}$$

with α_firn, α_frsnow and t^* being constants as well. This parametrization of snow albedo in COSIMA was tested at Mocho Glacier, where the glacier surface was covered by snow or firn during the observation period and precipitation data from a nearby automatic weather station were available. On the other glaciers an ice surface and a constant albedo of 0.3 was assumed.

The longwave radiative fluxes are computed using the Stefan–Boltzmann law of thermal radiation as well. The snow–ice surface is considered a blackbody. However, in contrast to the EB model in COSIMA the snow–ice surface temperature is variable depending on the heat fluxes at the glacier–atmosphere interface. Similar to the EB model the emissivity of the atmosphere is modeled as a function of the cloudiness (n) using the following expression:

$$\begin{array}{}\text{(12)}& {\mathit{\epsilon }}_{\text{(COS)}}(n,T,P_{\mathrm{vap}})={\mathit{\epsilon }}_{\mathrm{cs}}(T,P_{\mathrm{vap}})(\mathrm{1}-n^{\mathrm{2}})+\mathrm{0.984}{n}^{\mathrm{2}},\end{array}$$

where the emissivity of the clear sky depends on the air temperature and the water vapor pressure according to the following expression ${\mathit{\epsilon }}_{\text{cs}}(T,P_{\mathrm{vap}})=\mathrm{0.23}+\mathrm{0.433}(P_{\mathrm{vap}}/T{)}^{\mathrm{1}/\mathrm{8}}$, where T is in kelvins and P_vap in pascals. The red lines in Fig. 3b shows the variation in ε_(COS) as a function of the cloud cover at 5 ^∘C and assuming different relative humidities.

The turbulent flux of sensible heat in COSIMA SH_(COS) is calculated using Eq. (1) using the modeled surface temperature on the glacier surface. The same stability correction based on the bulk Richardson number Ri described in Sect. 3.2 is applied here to account for the reduced vertical exchange of air masses in stable conditions (Braithwaite, 1995b).

The turbulent flux of latent heat is calculated by the following expression:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathrm{LH}}_{\text{(COS)}}=& \mathrm{0.622}{\mathit{\rho }}_{\mathrm{a}}{L}_{\mathrm{v}}{C}^{\ast }U\left(z\right)\\ & \left[{\displaystyle \frac{P_{\text{vap}}\left(z\right)}{P_{\mathrm{a}}-P_{\text{vap,sat}}\left(z\right)}}-{\displaystyle \frac{P_{\text{vap}}\left(s\right)}{P_{\mathrm{a}}-P_{\text{vap,sat}}\left(s\right)}}\right].\end{array}\end{array}$$

Since the air pressure P_a is normally much higher than P_vap,sat, this formula should give very similar results as Eq. (3). This expression is multiplied by the same correction factor as SH (see Sect. 3.2). Concerning the parametrization of P_vap,sat as a function of temperature, we decided to replace the original parametrization in COSIMA (red line in Fig. 3a) by the parametrization proposed by Bolton (1980), i.e., Eq. (4), since it agrees best with the measurements. The original parametrization of COSIMA underestimates the water vapor pressure at positive temperatures and therefore underestimates LH for positive air temperatures, which are measured during summer at the AWSs (see below). COSIMA models additional energy fluxes like heat fluxes inside the snow or ice, but for a better comparison of the computed melt rates by the different methods, these fluxes were not considered in this contribution. The modeled heat flux inside the snow–ice with COSIMA depended on the initial temperature distribution inside the snow–ice and was maximum at San Francisco Glacier where it was 3 % of the sum of the modeled fluxes which are considered in this study.

4 Results

4.1 Microclimatic conditions on the glacier surface

In Table 4 we show averages of relevant climatic and glacier surface properties during summer for the six studied glaciers which are ordered according to their latitude from north to south.

Table 4 Mean values of relevant meteorological and glacier surface data during the study periods. For the longwave radiative fluxes, the bias-corrected value is indicated as well as the originally measured one in parentheses.

Variability in conditions is observed not only between the different glaciological regions but also inside each region. In the Central Andes the AWSs installed on Bello and Pirámide glaciers receive considerably more incoming solar radiation SW_in than the AWS on San Francisco Glacier which receives shade from the Mirador del Morado peak in the morning hours (see Figs. 1c and A1 in the Appendix) The mean albedo of the surface was calculated by two methods: firstly by calculating for every day the quotient of the daily sum of outgoing solar radiation divided by the daily sum of incoming solar radiation and taking the average of these values (${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\text{daily}}$) and secondly by simply dividing the mean outgoing solar radiation by the mean incoming solar radiation over the study period (${\mathit{\alpha }}_{\stackrel{\mathrm{‾}}{{\mathrm{SW}}_{\text{out}}}/\stackrel{\mathrm{‾}}{{\mathrm{SW}}_{\text{in}}}}$ ). Both methods give similar results. The heavily debris-covered Pirámide Glacier shows very low albedo. Bello Glacier shows an albedo expected for debris-rich ice, and San Francisco Glacier shows an albedo which can be associated with clean ice (Cuffey and Paterson, 2010). The incoming longwave radiation LW_in is lower on Bello Glacier which can be explained by the lower atmospheric temperature due to its higher elevation (see Table 1). The outgoing longwave radiation LW_out is highest for Pirámide glacier, whose surface is heated considerably in the afternoons (see Fig. 4f and below). No bias correction could be applied to these data since we do not know the real surface temperature of the debris which covers the glacier surface. Bello and San Francisco glaciers show bias-corrected mean values of LW_out of about 300 W m⁻². In the Central Andes wind speed is highest for Pirámide Glacier and the relative humidity is very similar for the three glaciers (around 40 %). At Mocho Glacier in the Lake District, SW_in is slightly lower than for Bello and Pirámide glaciers and the albedo is much higher, which is explained by the fact that on Mocho Glacier the AWS was installed near the equilibrium line altitude (ELA) and snow or firn were covering the glacier surface during the observation period (see Fig. 9 and Sect. 5.2). Both wind speed and relative humidity were clearly higher on Mocho Glacier in comparison with the glaciers of the Central Andes.

For the glaciers of the Patagonian Andes SW_in is clearly lower than for the glaciers of the other regions. This can be explained by their latitudinal dependency, due to the higher absorption of the solar radiation in the more humid and cloudy atmosphere in the Wet Andes and due to the fact that the glaciers in Patagonia are located at lower elevations. The albedo is higher for the clean Tyndall Glacier as compared to the partly debris-covered Exploradores Glacier (Fig. 2). LW_in is highest for Exploradores Glacier where also the highest relative humidity RH is observed. At Tyndall Glacier (the bias-corrected) LW_out is very near to the expected value for a melting ice surface (315.6 W m⁻²) in both years. For Exploradores Glacier it is slightly lower. Mean air temperature on Exploradores Glacier was similar to the one observed at Pirámide and San Francisco glaciers in the Central Andes and a bit lower at Tyndall Glacier. Measured wind speed was higher on Tyndall Glacier.

In order to study the daily cycle of the climatic variables on the glaciers, we calculated the average values which were measured at every hour of the day during the measurement period, presented in Fig. 4.

Figure 4 Averages of meteorological variables at the different hours of the day during the study periods on the six studied glaciers: (a) temperature T, (b) relative humidity RH, (c) water vapor pressure P_vap, (d) wind velocity U, (e) incoming longwave radiation LW_in, (f) outgoing longwave radiation LW_out. Time is local time.

As expected, the air temperature shows a daily cycle for most of the glaciers, with maximum temperatures in the afternoon and minima in the early morning hours. However this daily cycle is much less pronounced for Tyndall and Mocho glaciers. The relative humidity decreases during the daytime for most of the glaciers. The water vapor pressure shows a maximum during the late afternoon, when the air temperature is still elevated and the humidity is increasing. The wind speed shows a very pronounced increase during the daytime at Pirámide Glacier, when its debris-covered surface is heated (Fig. 4f). At Exploradores Glaciers an increase in wind speed during the afternoon is observed as well. The incoming longwave radiation shows a maximum during the daytime, when the atmospheric temperature is highest. The outgoing longwave radiation, which is emitted by the glacier surface has a very distinct maximum during the afternoon for Pirámide Glacier (increase of more than 100 W m⁻²), which means that the rock-covered glacier surface warms during the daytime. Bello, San Francisco and Exploradores glaciers also experience a maximum emission of longwave radiation in the afternoon, whereas Tyndall Glacier shows a constant rate of emission of longwave radiation, which indicates a constant surface temperature during summer.

4.2 Average energy balance and melt

Since both the EB model and COSIMA are models designed to compute the surface energy balance over snow and ice surfaces, we exclude the heavily debris-covered Pirámide Glacier from the analysis of the surface energy fluxes. In Fig. 5 and Table 5 we present the mean energy fluxes and inferred melt rates for the other five glaciers using the three different methods.

Figure 5 Mean modeled and measured energy fluxes and melt during summer for five glaciers: left bar, Reference Database; middle bar, EB model; right bar, COSIMA.

Table 5 Mean values of the computed surface energy fluxes and melt rates during the study periods on the five studied glaciers using the three different methods. For the Reference Database, turbulent fluxes and melt rates calculated without using stability corrections are indicated in parentheses.

Figure 6 Daily modeled and measured energy fluxes and inferred daily melt rates during summer 2015 for Bello Glacier using the three methods.

In Fig. 5 the three columns per glacier correspond to the Reference Database, EB model and COSIMA from left to right. For Tyndall Glacier only the results for the summer season 2016 are shown. The net energy flux towards the glacier was converted into daily melt rates in ice equivalent using an ice density of 917 kg m⁻³. The mean pattern of the energy fluxes changes from the Central Andes to Patagonia. The net shortwave radiation decreases from north to south. The net longwave radiation is negative in the Central Andes and near to zero in Patagonia. The sensible heat flux is a more important source of energy in Patagonia. The latent energy flux changes sign from a sink of energy in the Central Andes to a source of energy in Patagonia. This means that in Patagonia water vapor condensates at the surface of the glaciers, which generates heat for additional melt. There are differences in the prediction of the energy fluxes on the glacier surfaces between the different methods which will be discussed in detail in the next section. The predicted melt rates for the specific study points (locations of the AWSs) in the ablation area of the glaciers are higher for the Patagonian glaciers as compared to the glaciers of the Central Andes.

Figure 7 Daily modeled and measured energy fluxes and melt during summer 2006 for Mocho Glacier using the three methods.

4.3 Daily energy balance and melt

In Figs. 6, 7 and 8, we present the computed daily energy fluxes and melt rates for Bello Glacier, Mocho Glacier and Tyndall Glacier (2016), respectively.

Figure 8 Daily modeled and measured energy fluxes and melt during summer 2016 for Tyndall Glacier using the three methods.

On Bello Glacier the melt rates are clearly modulated by the net shortwave radiation: on days with reduced net shortwave radiation (due to the presence of clouds), melt rates show minima (Fig. 6). On Mocho Glacier this picture changes: high melt rates are instead associated with low net solar radiation, high contributions of the sensible heat flux and positive values of the latent heat flux (see for example 17 February or 8 March in Fig. 7). Similar results can be observed for Tyndall Glacier: peaks in melt rates are associated with a low contribution of the net shortwave radiation and high contributions of the turbulent heat fluxes (15 February or 20 March in Fig. 8).

Statistics of the comparison of the modeled and measured energy fluxes and the modeled melt rates by the three methods are presented in Table A1. A comparison of the hourly modeled and measured energy fluxes during the first three days of February on Bello Glacier and Tyndall Glacier are presented in Fig. A2.

5 Discussion

5.1 Microclimatic conditions on the glacier surface

The systematic variation in several meteorological variables between the different parts of the Chilean Andes crucially determines the importance of the different energy exchange processes at the glacier surfaces (Table 4). SW_in differs around 100 W m⁻² from the Central Andes (considering Bello and Pirámide glaciers) to the Patagonian Andes. The difference in SW_in between the two glaciers in the Patagonian Andes is only 9 W m⁻², although they have opposite exposition. The difference between two summers on Tyndall Glacier is only 4 W m⁻². Another clear trend form north to south is found for the relative humidity. The values measured in the Patagonian Andes are double the values obtained for the Central Andes. This influences the latent fluxes: in the Central Andes moisture is transported away from the glacier surfaces, while in Patagonia moisture is transported towards the glacier surfaces. Although the mean air temperature measured over Pirámide, San Francisco and Exploradores glaciers were very similar, the incoming longwave radiation was much higher (>75 W m⁻²) at Exploradores Glacier. This can be explained by a higher emissivity of the atmosphere due to a higher relative humidity of the air and due to more presence of clouds in the humid conditions of Patagonia.

Considering the variability in the data from the two summers measured on Tyndall Glacier, we can state that the glacier climate was similar in both summers. Especially mean LW_in, LW_out, RH and U were nearly identical. SW_in and T were slightly lower in 2015 as compared to 2016, and the surface albedo was higher in 2015. The mean values of T, RH and U measured on Tyndall Glacier during the summers 2015 and 2016 were also very similar to the mean values measured by Takeuchi et al. during December 1993 (Takeuchi et al., 1999).

Figure 9 Measured and modeled daily albedo on Mocho Glacier and precipitation registered in Puerto Fuy during summer 2006. The albedo value used in the EB model corresponds to the mean value of the measured albedo. In COSIMA the following parameters have been chosen: α_frsnow=0.8, α_firn=0.5, time constant $t^{*}=\mathrm{2}$ d, snow depth constant $d^{*}=\mathrm{8}$ cm (see Eqs. 10 and 11).

5.2 Parametrizations of the surface energy fluxes

The net shortwave radiation is an important source of energy for all the glaciers. According to Eq. (5) it is determined by the incoming shortwave radiation and the albedo of the surface. Albedo of snow and ice surfaces are very variable and depend on grain size and form, liquid water content, impurities, and other factors (Wiscombe and Warren, 1980; Warren and Wiscombe, 1980; Cuffey and Paterson, 2010). Generally, fresh snow has the highest albedo which decreases in time when snow grains grow and the snow eventually becomes dirty. COSIMA tries to reproduce this albedo aging effect by introducing a snow albedo which exponentially decreases in time (Eq. 11). In Fig. 9 we show the comparison between the measured daily albedo on Mocho Glacier during February and March 2006 and the predictions of the COSIMA model. In comparison to other studies that used similar albedo parametrizations (Mölg et al., 2012; Huintjes et al., 2015b), we only slightly reduced the albedo of fresh snow and firn (by 0.05) but significantly reduced the time constant t^* by a factor of 3. Most of the measured increases in the albedo can be associated with precipitation events registered at the nearby automatic weather station in Puerto Fuy (Fig. 9). COSIMA is able to capture these increases; however, the measured increases in the surface albedo are much more variable than the ones obtained from the model. A drawback of this comparison is certainly that we do not know the exact amount of snow falling on the glacier but deduce it from the liquid precipitation measured at an automatic weather station in the valley. The faster reduction in the albedo after snowfall and the associated lower value of the time constant t^* indicate a faster snow metamorphism on Mocho Glacier during the observation period, probably due to higher ambient temperatures in comparison to the high-altitude sites studied in Mölg et al. (2012) and Huintjes et al. (2015b). Considering the statistics of the modeled daily net shortwave radiation (Table A1), it is important to note that the slightly lower root-mean-square deviation for the EB model was obtained by the mean measured albedo, while in COSIMA a standard albedo of ice of 0.3 was applied.

The longwave radiative fluxes make important contributions to the energy exchange at the glacier surface. Since snow and ice emit approximately like blackbodies (ε=1) in this part of the electromagnetic spectrum and the atmosphere mostly shows emissivity smaller than 1, the longwave radiation balance is often negative (even at positive ambient temperatures). However in humid Patagonia the measured net longwave radiation is often small (Fig. 5 – upper panel; left bars; glaciers Mocho, Exploradores and Tyndall). In Fig. 10 we show the “measured” daily emissivity calculated by inverting the Stefan–Boltzmann law:

$$\begin{array}{}\text{(14)}& {\mathit{\epsilon }}_{\text{measured}}={\displaystyle \frac{{\mathrm{LW}}_{\text{in}}}{\mathit{\sigma }{T}^{\mathrm{4}}}},\end{array}$$

where σ denotes the Stefan–Boltzmann constant, as a function of the daily cloudiness for the glaciers Bello, San Francisco, Exploradores and Tyndall. For a direct comparison, we show the parametrizations of the models in the plots of every glacier.

Figure 10 Measured emissivity of the atmosphere as a function of the cloudiness. The data points correspond to the daily emissivity obtained from Eq. (14) by using daily means of LW_in and T plotted against the inferred daily cloudiness values. The solid colored straight lines correspond to linear fits to the data of every glacier and the discontinuous black lines correspond to the model parametrizations (same as Fig. 3b).

Generally we can note that the data of the measured emissivity show considerable scatter around their trend line. This may indicate that the variability in the measured emissivity cannot only be explained by the variability in the cloudiness. The relative humidity and other factors will probably influence the emissivity of the atmosphere. However the big scatter of the measured data might also be associated with the uncertainties in the determination of the cloud cover and the emissivity of the atmosphere. The cloud cover data were obtained from specific parametrizations of the transmissivity of the atmosphere as a function of the cloudiness. However these parametrizations are not unique (see for example Oerlemans, 2001). Also the measured emissivity is associated with some uncertainty: it is not clear if the temperature measured at 2 m over the glacier surface is representative of the temperature of the atmosphere which emits longwave radiation towards the glacier surface. We can recognize that the model parametrizations underestimate the emissivity of the atmosphere at Exploradores Glacier for low-cloudiness conditions. This underestimation of the emissivity of the atmosphere leads to an underestimation of the net longwave radiative balance on clear days. At San Francisco Glacier the model parametrizations overestimate the emissivity of the atmosphere for cloudy conditions. Considering the statistics of the modeled daily net longwave radiation (Table A1), we can recognize that the correlation is much lower than in the case of the modeled net shortwave radiation. Especially the negative correlations in the case of the EB model for the glaciers of the Central Andes indicate that here this model is not able to reproduce the measured daily variations in this flux. Reasons for this could be an unsuccessful determination of the cloudiness conditions by this method and the fact that during nighttime a constant cloudiness was used.

The variability in the modeled turbulent fluxes is very similar in all three methods (see Table A1). This is expected since the formula to compute these fluxes have very similar aspects: in all approaches the sensible heat flux is mainly driven by the temperature difference in the glacier surface and the atmosphere at 2 m elevations and wind speed and the latent heat flux is driven by the difference in the water vapor pressure at the glacier surface and the atmosphere at 2 m elevations and wind speed. The mean values of the turbulent fluxes computed in the EB model are lower than the one obtained by the other two methods (Table 5 and Fig. 5). This is because the EB model assumes a glacier surface at 0 ^∘C which reduces both the temperature difference between glacier surface and the overlying air layer and the difference in water vapor content between both. This causes the melt rates modeled by the EB model to be generally lower (Fig. 5). The modeled turbulent fluxes in the Reference Database and by COSIMA are very similar, except for Mocho Glacier where a glacier surface at 0 ^∘C had to be assumed in the Reference Database, due to the lack of data on the outgoing longwave radiation.

In Table 5, for the Reference Database, we also present in parentheses values obtained for the turbulent fluxes without applying a stability correction and the resulting inferred melt rates. Important differences can be noted especially for the glaciers where the mean wind speed is moderate (Exploradores, San Francisco and Bello). The strongest influence of the stability correction on the melt rate is for Exploradores Glacier, because both turbulent fluxes have the same sign. The stability correction for the sensible heat flux and the latent heat flux perfectly cancel out at Bello Glacier. The influence of the surface roughness on the turbulent fluxes has a similar character: for Bello Glacier, doubling z₀ will change the net energy flux by only 1 W m⁻² and using a roughness length of 10 mm (20 times higher) will change it by 5 W m⁻². At the Patagonian glaciers, however, a higher surface roughness would have a much higher impact since both turbulent fluxes have the same sign.

5.3 Melt rates

Mean melt rates ranging from 2.9 to 4.4 cm d⁻¹ of ice equivalent (cm i.e. d⁻¹) for the Dry Andes and from 3.3 to 6.9 cm i.e. d⁻¹ for the Wet Andes were predicted by the different approaches to quantifying energy fluxes on the surface of five glaciers (Table 5). These values lie within the range of observed melt rates during summer on these glaciers or glaciers with similar climatic conditions.

In the Wet Andes Schaefer et al. (2017) measured ablation rates of 2.6 and 3.2 cm i.e. d⁻¹ in the summers of 2009/10 and 2010/11 and measured and inferred rates of 3.6 and 3.7 cm i.e. d⁻¹ in the summers of 2011/12 and 2012/13 on Mocho Glacier at the same location where the AWS was installed in 2006. This indicates that the Reference Database and COSIMA may overestimate the melt at this location. Here, we have to take into account that for Mocho Glacier net longwave radiation was inferred by subtracting the net shortwave fluxes from the net all-wave radiation, measured with the NR Lite, which is a lower-precision instrument as compared to the newer sensors installed in the CNR4. At Tyndall Glacier in the period of November 2012 to May 2013 an average ablation rate of 3.8 cm i.e. d⁻¹ was observed at two stakes near to the location of the AWS. Considering that this period also includes spring and autumn months, when melt rates should be lower, this is in good agreement with the values of 3.3 to 5.1 cm i.e. d⁻¹ predicted by the different approaches presented in this work. The modeled melt rates for the Wet Andes are also in agreement with the melt rates of 4–5.5 cm i.e. d⁻¹ observed during summer on Perito Moreno Glacier on the Southern Patagonian Ice Field (Stuefer et al., 2007) and the 4–8 cm i.e. d⁻¹ observed at Nef Glacier (Schaefer et al., 2013), both at elevations of about 500 m a.s.l. From January to March 2015 an average ablation rate of 9.1 cm i.e. d⁻¹ was measured at a stake network installed on Grey Glacier at an elevation range of 260 to 380 m a.s.l. This indicates that the high melt rates modeled for Exploradores Glacier seem to have a realistic magnitude for a low-elevation site on a glacier in the Patagonian Andes, although the different climate conditions at Grey Glacier make a direct comparison difficult.

In the Dry Andes ablation was measured at a stake network on Bello Glacier during summers 2013/14 and 2014/15 (CEAZA, 2015). A high variability in ablation rates in space and time was obtained. Several of the ablation rates inferred for stakes near the AWS were of similar size to the ones predicted by the methods presented in this paper. Analyzing the signal of an ultrasonic sensor installed on an ablation gate next to the AWS of San Francisco Glacier, we found a surface lowering of 4.9 cm d⁻¹ during December 2015. Assuming snowmelt in this period of the year with a density of 500 kg m⁻³, this yields a rate of 2.7 cm i.e. d⁻¹, which is in good agreement with the melt rates inferred in this contribution.

When comparing modeled melt rates between the three methods presented in this study, we can state that there exist considerable differences between the three approaches with a root-mean-square deviation of the daily modeled melt rates of around 1 cm i.e. d⁻¹ (see Table A1).

5.4 Implications for glaciological zones in Chile

Our results suggest a transition of energy sources for surface melt from the Central Andes to Patagonia. The energy sources obtained for Mocho Glacier at 40^∘ S are more similar to the ones observed at the Patagonian glaciers than to the ones observed for the glaciers of the Central Andes, where SW_in is the dominating source of energy for melt (Figs. 5 and 6). The greater importance of the turbulent flux of sensible heat as an energy source available to produce surface melt at Mocho Glacier probably contributes to the observed strong dependency of its annual mass balance on the annual mean temperatures measured on a nunatak of the glacier (Scheiter, 2016; Schaefer et al., 2017). This picture changes strongly in the Central Andes, where a very low dependency of the annual mass balance of Echaurren Norte Glacier on the annual mean temperature observed at the El Yeso reservoir was reported (Masiokas et al., 2016; Carrasco, 2018; Farías-Barahona et al., 2019). When comparing the importance of the energy sources between the two summers modeled for Tyndall Glacier (Table 5), we can state that in 2015 SW_net was slightly lower than in 2016 due to the lower SW_in and the higher surface albedo observed in that year (Table 4). However the overall pattern of energy sources (and sinks) is very similar for both years.

Figure 11 Melt rate in the different zones as a function of temperature (a) and as a function of the relative elevation difference from the ELA (b).

In Fig. 11a we plot the predicted melt rates by the Reference Database against the mean air temperature for the six modeling periods (two periods for Tyndall Glacier) and in Fig. 11b against the relative elevation difference in the AWS with respect to the equilibrium line altitude (ELA), which we define as (ELA−ElevationAWS)∕(elevation span glacier). For the glaciers of the Wet Andes we can see a clear increase in the melt rates as a function of the mean temperature. For the modeled glaciers of the Central Andes, however, this trend does not exist. In both regions the melt rates increase with the relative elevation difference from the ELA; however at a similar relative elevation difference, the glaciers of the Wet Andes show clearly higher melt rates.

All these results confirm the general division of the Chilean Andes into the Dry Andes and Wet Andes, with the division being located at approximately 35^∘ S (Lliboutry, 1998).

5.5 Implications for physical melt modeling and transferability of parametrizations

The capacity of the models to reproduce the measured radiative fluxes can still be improved. The albedo aging effect implemented in COSIMA is a big improvement compared to constant albedo parametrizations for snow, firn and ice surfaces. However the parameters of this aging formula seem to vary strongly from one site to another, and a good calibration of this formula seems to be necessary. Regarding the predictions of the net longwave radiative fluxes on the glacier surface, the parametrization of the emissivity of the atmosphere is crucial. The tested models cannot reproduce the variability in the emissivity as a function of the cloudiness for all the glaciers. Especially at Exploradores Glacier the clear-sky emissivity is underestimated by the models. This is because the parametrizations used are fits to data that were obtained in different climatic conditions. Therefore, these parametrizations are not physical and cannot be simply transferred to other sites where the conditions are different.

Different parametrizations for the turbulent fluxes of sensible and latent heat were compared in this study. The transfer coefficient depends directly on the roughness length z₀ (and z_T and z_H in the case of EB model), which therefore crucially determines the magnitude of the turbulent fluxes. However these roughness lengths are constant neither in time nor in space, which makes them very difficult to determine. A common practice is to chose the roughness length in a way that the modeled melt rates agree with the measured ones. However, this exercise does not make these formulas very adequate predictors of melt rates in other situations and on other glacier surfaces. More direct measurements of turbulent fluxes over glacier surfaces (for example using the eddy-covariance technique; Cullen et al., 2007; Litt et al., 2015) are necessary to find physical parametrizations of these fluxes.

Generally, there is still a strong need for measurements of the energy exchange processes over glacier surfaces, and we think that coordinated efforts of governmental agencies, such as the Glaciology and Snow Unit of the Chilean Water Directorate in our case, can make important contributions. If we want to work towards physical melt models, then we have to test the capacity of the models to reproduce the different physical processes that take place at the glacier surface.

6 Conclusions

Performing an extended study of surface energy fluxes during summer on five Chilean glaciers on a north–south transect and under strongly varying climate settings we reached the following conclusions:

• The contribution of the different surface energy fluxes over glacier surfaces change from the Central Andes towards the Patagonian Andes. The net shortwave radiation as a main source of energy in the Central Andes loses importance further south, where the turbulent fluxes of sensible and latent heat provide more energy for melt. The net longwave radiation changes from a strong sink of energy for the glaciers of the Central Andes to a net zero contribution in the Patagonian Andes.

• For the glaciers of the Wet Andes a clear dependency of the modeled melt rates on the mean air temperature was observed. This dependency did not exist for the studied glaciers in the Central Andes.

• The inferred melt rates increased with the relative elevation difference from the ELA in both studied glaciological zones. The modeled melt rates were higher for the Patagonian Andes than for the Central Andes.

• Mocho Glacier in the Chilean Lake District shows similar patterns of surface energy fluxes to the glaciers in the Patagonian Andes.

• The models underestimated the emissivity of the clear-sky atmosphere at Exploradores Glacier, an extremely humid place in the Wet Andes.

• From our study it is difficult to infer which parametrization of the turbulent fluxes is the most appropriate one. More detailed studies on this topic are necessary, which include direct measurements of these fluxes.

• To develop or improve physical melt models, we have to validate every single model parametrization against data and cannot judge the model's performance only by the final output. In these highly parametrized models the effect of physically wrong parametrizations might be canceled out and the final result might be satisfying without reproducing the individual physical processes well.

• Openly shared codes are the best way to improve physical models, since everyone can test the individual parametrizations against their own data and adjust or improve them accordingly. This is the preferred way to improve physical parametrizations as opposed to large chains of models which try to model physical processes without validating intermediate model results.

Appendix A

Figure A1 Incoming shortwave radiation at the AWSs installed on Bello, Pirámide and San Francisco glaciers during March 2016.

Table A1 Statistical comparison of the modeled daily energy fluxes and melt by the EB model and COSIMA with the Reference Database. Pearson correlation coefficient (r) and root-mean-square deviation (RMSD) are indicated.

Figure A2 Comparison of modeled and measured hourly energy fluxes on glaciers Bello (a) and Tyndall 2016 (b) for the first 3 d of February. Time is local time.