1. Introduction
More than 99% of the world’s tropical glaciers are concentrated in the Andes Mountains. Approximately 70% of their total area is located in Peru, 20% in Bolivia and the rest in Ecuador, Colombia and Venezuela (Kaser, Reference Kaser1999). The Peruvian and Bolivian glaciers are situated in the outer tropics, with Bolivian glaciers marking the southernmost boundary (16°S). In mid- and high-latitude regions, seasonal temperature variation is significant. Glacial mass accumulation predominantly occurs during the winter, while ablation is concentrated in summer. A small temperature rise during the winter does not substantially affect precipitation and has a limited impact on accumulation. In the outer tropics, there are two distinct seasons: the wet season (austral summer) and dry season (austral winter). In contrast to the significant seasonal variations in humidity and precipitation, temperature remains relatively constant throughout the year. Hence, glacial accumulation and ablation occur simultaneously during the wet season, while both processes are relatively low during the dry season (Sicart and others, Reference Sicart, Ribstein, Francou, Pouyaud and Condom2007; Sicart and others, Reference Sicart, Hock, Ribstein, Litt and Ramirez2011). During the wet season, precipitation frequently occurs near 0°C. Even a slight temperature increase during this period can shift precipitation from snowfall to rainfall in the ablation area, leading to a significant reduction in accumulation. Tropical glaciers are therefore considered highly sensitive to global warming (Schauwecker and others, Reference Schauwecker2017). Gorin and others (Reference Gorin2024) implied that recent tropical glacier retreat is unprecedented in the Holocene, identifying the Andes as a critical region for understanding the changing state of the cryosphere.
A global study on water supply from mountains to downstream regions (Immerzeel and others, Reference Immerzeel2020) suggested that South American mountains near 16°S latitude are prominent water towers, and downstream societies are vulnerable due to population and economic growth. Meltwater from glaciers in the Bolivian Andes is a primary water resource for the Bolivian capital, La Paz, and its neighboring city, El Alto, which together have a combined population of over two million (Soruco and others, Reference Soruco2015). This region is semi-arid, with annual precipitation averaging only ∼500 mm. Hence, glacier meltwater is crucial to meet water demand, especially during the dry season. Moreover, the population of El Alto doubled from ∼400 000 to ∼800 000 between 1992 and 2012, leading to a continuous increase in water consumption. Bolivian tropical glaciers are also highly sensitive to global warming, with most glaciers rapidly shrinking, fragmenting and disappearing (Sorg and others, Reference Sorg, Bolch, Stoffel, Solomina and Beniston2012; Paul and others, Reference Paul2013). Morizawa and others (Reference Morizawa, Asaoka, Kazama and Gunawardhana2013) reported that the Condoriri Glacier and the fragmented glaciers surrounding it, located upstream of La Paz and El Alto, lost 41% of their area between 1988 and 2010. If this retreat rate persists, their role as a water supplier could be lost by 2035.
Estimating glacier meltwater is essential for sustainable water resource management, particularly for domestic water and irrigation water downstream of glaciers (Immerzeel and others, Reference Immerzeel2020; Qin and others, Reference Qin2020). Degree-day models (DDM) rely on an empirical relationship between air temperature and melt rate (Hock, Reference Hock2003). Air temperature is one of the key factors affecting sensible heat transfer from the atmosphere in the physical process driving snow and ice melt. Additionally, air temperature influences the surface energy balance through atmospheric longwave radiation. DDM are thus considered to provide sufficiently accurate melt rate estimates for most practical applications, despite their simplicity (Ohmura, Reference Ohmura2001). Since net shortwave solar radiation is also a dominant source of melt energy on most glaciers, enhanced temperature-index (ETI) models extend DDM by incorporating solar radiation as an additional explanatory variable in combination with air temperature (Pellicciotti and others, Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005). These empirical models have been widely applied to estimate glacier melt and runoff at both global scales (Hirabayashi and others, Reference Hirabayashi, Döll and Kanae2010; Hirabayashi and others, Reference Hirabayashi, Zang, Watanabe, Koirala and Kanae2013; Immerzeel and others, Reference Immerzeel, Pellicciotti and Bierkens2013; Kraaijenbrink and others, Reference Kraaijenbrink, Stigter, Yao and Immerzeel2021) and catchment scales (Fuchs and others, Reference Fuchs, Asaoka and Kazama2016; Kinouchi and others, Reference Kinouchi, Nakajima, Mendoza, Fuchs and Asaoka2019; Khanal and others, Reference Khanal, Lutz, Kraaijenbrink, Hurk, Yao and Immerzeel2021).
Melting is driven by energy exchange between the atmosphere and the ice or glacier surface. Energy balance models estimate melt energy based on the surface energy balance and are suitable for investigating the specific drivers of melt depending on climatic zones. However, these models require a variety of meteorological input variables and are not always suitable for estimating glacier melt and runoff at the catchment scale due to limited data availability (Sicart and others, Reference Sicart, Hock, Ribstein, Litt and Ramirez2011; Fyffe and others, Reference Fyffe2021). On the other hand, a better understanding of the energy balance can improve empirical melt models, which are widely used for estimating melt rates over large areas. Empirical melt models have previously been validated based on energy balance analyses for glaciers located in mid- and high-latitude regions (Gabbi and others, Reference Gabbi, Carenzo, Pellicciotti, Bauder and Funk2014; Yang and others, Reference Yang2011; Carenzo and others, Reference Carenzo, Pellicciotti, Mabillard, Reid and Brock2016; Litt and others, Reference Litt2019). In the Andes, there are only a few energy balance studies of tropical glaciers in Peru (Winkler and others, Reference Winkler, Juen, Mölg, Wagnon, Gómez and Kaser2009; Aubry-Wake and others, Reference Aubry-Wake, Zéphir, Baraer, Mckenzie and Mark2018) and Bolivia (Wangon and others, Reference Wagnon, Ribstein, Kaser and Berton1999; Sicart and others, Reference Sicart, Wagnon and Ribstein2005; Sicart and others, Reference Sicart, Hock, Ribstein, Litt and Ramirez2011). Fyffe and others (Reference Fyffe2021) applied a physically based energy balance melt model to five on-glacier sites located between 8°S and 15°S in Peru. Their results showed that net shortwave radiation was dominant in the energy balance during the dry season, but melt energy was reduced due to losses associated with net longwave radiation and latent heat fluxes. At three of the sites, the wet season snowpack was discontinuous, and melt rates were highly variable, closely related to the precipitation phase. The southern boundary of tropical glaciers is situated at ∼16°S in Bolivia, where the glacier surface energy balance has been intensively measured only at the Zongo Glacier. Sicart and others (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005) evaluated the energy fluxes on the Zongo Glacier and showed that melt energy was high during the progressive development (from September to January) toward the core wet season. Incoming shortwave radiation was identified as the main source of energy, but seasonal changes in melt energy were driven by variations in longwave radiation. Longwave radiation was thus identified as a key factor driving the melt of tropical glaciers.
A better understanding and quantification of the seasonal energy balance of tropical glaciers can improve empirical glacier melt modeling at larger scales and advance water resource assessments. This study investigates the energy balance characteristics of a tropical glacier using meteorological data measured at the Huayna-Potosí Glacier in Bolivia (16°S) and a multilayer snowpack heat balance model (Yamazaki, Reference Yamazaki2001). The dominant meteorological variables suitable for empirical melt models are identified based on results from the energy balance analysis. Finally, existing empirical melt models for tropical glaciers are validated, and a more suitable empirical melt model, based on the energy balance processes, is recommended for future applications.
2. Methods
2.1. Study site
The Huayna-Potosí Glacier (Figure 1) is located on the western slope of the Huayna-Potosí peak in the Cordillera Real. The peak, with a summit elevation of 6088 m, is situated ∼25 km north of the Bolivian capital, La Paz. The glacier spans an elevation range between 4900 and 6600 m, with over 57% of its area lying above 5100 m (Funaki and Asaoka, Reference Funaki and Asaoka2016). Streamflow from the Huayna-Potosí Glacier supports several wetlands (Navas and others, Reference Navas2024) and glacial lakes along its course (Huang and Kinouchi, Reference Huang and Kinouchi2024). The Tuni Reservoir withdraws water from streams originating from the Huayna-Potosí Glacier via an open channel during the wet season. A constant water volume is delivered from the Tuni Reservoir to a water treatment plant in El Alto at a discharge rate of 1.0 m3 s−1. In 2011, an automatic weather station (AWS) was installed on the Huayna-Potosí Glacier (16.26428°S, 68.31931°W, 5150 m a.s.l). The AWS records downward and upward shortwave radiation, downward and upward longwave radiation, air temperature, humidity, wind speed and direction, snow depth and precipitation at 10 min intervals.
Figure 1. Huayna-Potosí Glacier. Location of the Huayna-Potosí glacier in the tropical Andes (a) and detailed map of the study site showing the AWS location and surface elevation contours (b).
2.2. Application of a multilayer snow energy balance model to a glacier
The Tohoku Snow Multilayer Model (TSMM) is a one-dimensional land surface model based on a multilayer energy balance approach for snow and ice cover (Yamazaki, Reference Yamazaki2001). It divides the snow and ice into vertical layers, each with a thickness of 0.02 m, and calculates parameters such as temperature, freezing and melting times, density, water content and the depth of hoar formation. These calculations are based on processes including the discrimination between rain and snow, changes in albedo, penetration of solar radiation, surface energy balance, vertical heat conduction and vertical moisture transfer. The discrimination between rain and snow is calculated based on the wet bulb temperature derived from air temperature and air pressure. Changes in albedo are applied specifically to snow surfaces. Penetration of solar radiation into the snow layers is also considered. These parameterizations were originally developed and described by Yamazaki (Reference Yamazaki2001). These calculations are conducted separately for each snow and ice layer, thereby allowing for vertical profiles of temperature, density and water content to be resolved throughout the snowpack and glacier ice column.
In applying the TSMM to the Huayna-Potosí Glacier, we considered the albedo, density and thermal conductivity of both the glacier ice and snow. TSMM was originally developed to model snowpack dynamics and has been successfully applied to snow-covered environments. This study extended its application to include glacier ice, enabling the calculation of energy exchange at both snow and ice surfaces. To facilitate this extension, the initial conditions for the simulation were assumed to be a vertical 8-m column of ice with no snow, starting in July 2012. The ice column was assigned a surface albedo of 0.4, and the density of each layer was set to 900 kg m−3 (Asaoka and others, Reference Asaoka, Yamazaki, Miyata, Kazama and Ramirez2013). The thermal conductivity of snow and ice was calculated using an empirical formula originally proposed by Devaux (1933), as cited in Sturm and others (Reference Sturm, Holmgren, König and Morris1997) and implemented in Yamazaki (Reference Yamazaki2001).
In the TSMM, the energy balance at the snow/ice surface without thickness was expressed as follows:
\begin{equation}r\left(1-A\right)S\downarrow+\epsilon\left(L\downarrow-\sigma T_s^4\right)-H-\ell E-G_0=M\end{equation} where 
$r$ (dimensionless) is the solar radiation absorptivity at the surface, 
$A$ (dimensionless) is the surface albedo, 
$S \downarrow $ (W m−2) is the downward shortwave radiation, 
$L \downarrow $ (W m−2) is the downward longwave radiation, 
${T_s}$ (K) is the surface temperature, 
$\sigma $ (W m−2 K−4) is the Stefan–Boltzmann constant, 
$ \epsilon $ (dimensionless) is the snow/ice emissivity (set to 0.97, based on Kondo and Yamazawa (Reference Kondo and Yamazawa1986)). 
$H$ (W m−2) is the sensible heat flux (positive when directed from the surface to the atmosphere), 
$\ell E$ (W m−2) is the latent heat flux (positive when directed from the surface to the atmosphere) and 
${G_0}$ (W m−2) is the conductive heat just beneath the surface (positive when directed downward from the surface into the snow or ice). 
$M$ (W m−2) is the melt energy (M > 0 indicates melting, and M < 0 indicates freezing). The sensible and latent heat fluxes were calculated using the bulk aerodynamic method, following the implementation by Yamazaki (Reference Yamazaki2001). The sensible heat flux is based on the temperature gradient near the surface, while the latent heat flux is derived from the specific humidity gradient. The surface solar absorptivity was estimated based on the albedo using the method described in Yamazaki (Reference Yamazaki2001), and is expressed as follows:
\begin{equation}r=\frac{A_{max}-A}{A_{max}-A_{min}}\end{equation}where Amax = 0.88 and Amin = 0.5. When the albedo was less than 0.5, r = 1.0, implying no penetration of solar radiation into the glacier under these conditions. The original TSMM calculated albedo attenuation after snowfall using an empirical equation (Yamazaki and others, Reference Yamazaki, Taguchi and Kondo1994). However, in this study, the daily albedo was calculated from measured downward and upward shortwave radiation to reduce uncertainty in the energy balance calculation.
To drive the energy exchange simulation, temperature, humidity, wind speed, downward shortwave radiation, upward shortwave radiation and downward longwave radiation measured at the AWS on the Huayna-Potosí Glacier were input into the TSMM. Meteorological data were recorded at 10-min intervals. For the TSMM simulation, the data were averaged to 30-min intervals and then used as input, while the model calculations were performed every 150 s. Specific humidity was calculated from the observed relative humidity assuming an atmospheric pressure of 550 hPa. Due to the large uncertainty in the precipitation recorded by the rain gauge, precipitation input into the TSMM was reconstructed based on changes in snow depth. Specifically, the increase in snow depth over each 30-min interval, measured by the ultrasonic snow depth sensor, was multiplied by the snowfall density (Sicart and others, Reference Sicart, Ribstein, Chazarin and Berthier2002) to estimate snowfall. The snowfall density was assumed to be 150 kg m−3 after performing four simulations with densities of 100, 150, 200 and 250 kg m−3. The simulation with a snowfall density of 150 kg m−3 produced the most reasonable results, which were in close agreement with values reported by Sicart and others (Reference Sicart, Ribstein, Chazarin and Berthier2002) for the Zongo Glacier.
The simulation period for the TSMM energy balance analysis was from 1 June 2012, to 1 July 2013. The initial condition assumed the absence of snow cover on the surface. The first month of the simulation period was used as a warm-up period, and simulation results from July 2012 onward were used for analysis. The calculation time step was set to 150 s.
2.3 Empirical melt modeling
We derived simple empirical melt models through multiple linear regression using meteorological variables selected based on their correlation with melt energy estimated by the energy balance model (TSMM). Although these models include variables commonly used in traditional DDM and ETI models, such as air temperature and shortwave radiation, the structure of our regression models differs from those classical forms. In contrast to DDM and ETI models, which typically involve threshold temperatures and empirical melt factors (e.g. Hock, Reference Hock1999; Pellicciotti and others, Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005), our empirical models take the form of standard linear regressions directly fitted to melt energy values produced by the TSMM. Table 4 presents the regression equations used and their performance in comparison to the energy balance model estimates.
To evaluate the empirical estimation of melt, we identified the meteorological variables most strongly correlated with the melt energy obtained from the TSMM simulations. The DDM uses air temperature as an explanatory variable. The ETI model incorporates not only air temperature but also solar radiation and albedo. In the ETI model, albedo 
$\alpha $ and solar radiation 
$S \downarrow $ are treated as 
$\left( {1 - {\alpha }} \right)S \downarrow $, which represents net shortwave radiation. In several previous applications of the ETI model to catchments or larger scales, a sub-model for estimating albedo (Oerlemans and Knap, Reference Oerlemans and Knap1998; Yamazaki, Reference Yamazaki2001) was incorporated due to the limited availability of albedo measurements on glaciers (Fuchs and others, Reference Fuchs, Asaoka and Kazama2016; Kinoichi and others, Reference Hock2019). In the present study, surface albedo was calculated from measured upward and downward shortwave radiation and directly used as model input.
We analyzed the correlation between melt energy and the following meteorological variables: air temperature, relative humidity, specific humidity, shortwave radiation and net shortwave radiation. While longwave radiation is an important component of the surface energy balance, it is not routinely measured in many glacierized regions. In contrast, air temperature and humidity are more commonly observed, making them practical substitutes for capturing the influence of atmospheric longwave radiation in empirical models. The statistical significance of the correlation analysis was evaluated using p-values. Additionally, multiple regression analyses were conducted using variables that were strongly correlated with melt energy. Since most applications of empirical melt models aim to estimate melt rates on a daily or longer timescale, daily data were used in this analysis. This approach allows for flexibility in variable selection based on local energy balance conditions and the availability of meteorological observations.
3. Meteorological conditions
Sicart and others (Reference Sicart, Hock, Ribstein, Litt and Ramirez2011) observed the energy balance and meteorological conditions at the Zongo Glacier, located on the eastern slope of the Huayna-Potosí peak. Their results showed a distinct seasonality in precipitation and cloud cover, with a pronounced dry season (May to August) and a single wet season (September to April). From September to December, a progressive development toward the core wet season was observed.
Our meteorological data, including temperature, humidity, solar radiation, albedo and snow depth, measured at the AWS from June 2012 to July 2013, exhibit a similar seasonality (Figure 2). The annual mean temperature from August 2012 to July 2013 was −0.6°C. The monthly mean temperature reached a maximum of 1.5°C in November and a minimum of −2.3°C in June. Temperature increased from September onward, marking the beginning of the rainy season. Temperatures remained slightly higher between November and March than during the rest of the year and decreased from April onward, marking the onset of the dry season.
Figure 2. Daily meteorological conditions measured at the AWS on the Huayna-Potosí Glacier. Panels (a)–(d) show daily mean values of air temperature, relative humidity, solar radiation and surface albedo, respectively. Panel (e) shows snow depth recorded as a daily instantaneous value at 00:00.
Relative humidity exhibited significant day-to-day variability; however, clear seasonal patterns allowed the identification of distinct wet and dry seasons. The dry season, from May to September, exhibited relative humidity fluctuating between 20% and 80%. From October to November, relative humidity steadily increased, and the wet season, from December to March, was characterized by relative humidity frequently exceeding 80%.
Solar radiation was strongly influenced by cloud cover. August and September were predominantly characterized by clear skies, resulting in the highest monthly solar radiation. In October and November, during the transition toward the core wet season, clear skies alternated with periods of high cloud cover, leading to highly variable solar radiation. From December to March, solar radiation was consistently low due to persistent high cloud cover.
Albedo depended on surface conditions such as snow cover, firn and glacial ice. The albedo of fresh snow was high but decreased over time in the absence of new snowfall. Furthermore, the albedo of firn and glacial ice with prolonged exposure was significantly lower (Sicart and others, Reference Sicart, Wagnon and Ribstein2005). During the dry season from May to August, albedo ranged from 0.4 to 0.8. It rapidly increased to ∼0.8 following snowfall events but quickly decreased to 0.5 in the absence of fresh snow and when snow transformed into firn. From September to November, albedo fluctuated between 0.6 and 0.8, gradually increasing as snowfall became more frequent toward the peak of the wet season. From December to February, it remained between 0.7 and 0.8 due to frequent fresh snowfall. Starting in March, albedo ranged from 0.5 to 0.8, with rapid increases following fresh snowfalls but decreased to lower values during relatively long periods without snowfall. The albedo in July 2013 was ∼0.2 higher than in July 2012, likely due to lingering snow cover from the wet season.
Snow depth remained nearly constant during the dry season but began to increase steadily from September onward. It peaked in November and decreased from December onward. Snow depth in July 2013 was higher than in July 2012, suggesting that some snow cover accumulated during the wet season persisted. This observation was consistent with the measured albedo.
4. Results
4.1. Application of the TSMM on the glacier
To validate the applicability of the TSMM to the Huayna-Potosí Glacier, we compared calculated surface temperatures with measured surface temperatures (Figure 3). Surface temperature is a direct result of the surface energy balance, determined by the net radiative fluxes, turbulent heat fluxes (sensible and latent) and subsurface heat conduction. Therefore, accurate reproduction of diurnal surface temperature variations implies that the model reasonably represents the surface energy balance. Results are shown for August, representing the dry season with the lowest monthly air temperature, and November, representing the wet season with the highest monthly air temperature. Measured surface temperatures were estimated using the Stefan-Boltzmann law and upward longwave radiation measurements (Jung and others, Reference Jung, Lee and Zo2020). The root-mean-square error (RMSE) was 2.47°C for the year (August to July), 3.23°C in August and 1.50°C in November. The timing of the minimum and maximum surface temperatures in the diurnal cycles showed good agreement with observations, and the diurnal variations were successfully reproduced. In August, the maximum surface temperatures and their timing were well represented, but minimum surface temperatures were overestimated by 2–3°C. In November, the diurnal variation closely matched the minimum daily surface temperatures, but the maximum temperatures were underestimated by 1–2°C. While panels (a) and (b) illustrate short-term variations during representative months, panel (c) provides a broader view of seasonal trends by showing daily average surface temperatures throughout the entire analysis period (August 2012 to July 2013). Panel (c) thus allows assessment of the seasonal consistency between observed and modeled surface temperatures and of the model’s potential biases over time. The TSMM generally reproduced the seasonal pattern well, with slight deviations observed during certain months, particularly from May to July, when the modeled surface temperatures were slightly higher than observed values. This seasonal validation demonstrates that the model can reasonably capture the glacier surface temperature behavior throughout the year, supporting its applicability for simulating energy balance in tropical glaciers. The surface temperature was derived from upward longwave radiation, which may be influenced by atmospheric conditions or sensor characteristics. Therefore, values slightly exceeding 0°C do not necessarily indicate physically unrealistic surface temperatures, but rather may include limitations in longwave radiation-based temperature estimation (Sicart and others, Reference Sicart, Hock, Ribstein, Litt and Ramirez2011).
Figure 3. Comparison of estimated (black) and measured (blue) glacier surface temperatures. (a) Hourly surface temperature variation for August 2012 (dry season), (b) hourly variation for November 2012 (transition to wet season) and (c) daily average surface temperatures from August 2012 to July 2013, showing seasonal trends. Measured surface temperatures were derived from upward longwave radiation using the Stefan–Boltzmann law.
4.2. Seasonal variations in energy balance
Figure 4 shows the seasonal variations in monthly mean values of radiation components, energy fluxes, melt energy and meteorological conditions. The radiation components were measured by the AWS, while the sensible heat flux, latent heat flux and melt energy were estimated using the TSMM. The downward shortwave radiation exhibited a small seasonal variation, with a difference of ∼65 W m−2 between the monthly maximum (246 W m−2 in August) and the monthly minimum (181 W m−2 in March).
Figure 4. Seasonal variations in monthly mean radiation, energy fluxes, melt energy and meteorological conditions observed at the Huayna–Potosí Glacier. (a) Radiation components, (b) net radiation and heat fluxes, (c) melt energy, (d) air and surface temperature and (e) specific humidity of the air and surface.
In general, in the Southern Hemisphere, solar radiation at the top of the atmosphere reaches a maximum during the austral summer and a minimum during the austral winter. However, at the surface, the relationship between the wet and dry seasons is reversed compared to that at the top of the atmosphere due to the effect of cloud cover. Downward shortwave radiation was lower during the wet season than during the dry season because of more frequent cloud cover during the former.
The seasonal variation of upward shortwave radiation was opposite to that of downward shortwave radiation, with a maximum during the wet season and a minimum during the dry season. From December to February, 75–80% of the downward shortwave radiation was reflected due to the high albedo of fresh snow. From July to October, ∼55% of the downward shortwave radiation was reflected. This value represents the lower limit of snow albedo, as the albedo of seasonal snow diminishes over time. During this period, the surface was likely covered by snow that had persisted from the wet season. As a result of these seasonal changes in downward and upward shortwave radiation, net shortwave radiation reached a minimum of ∼40 W m−2 from December to February.
Downward longwave radiation reached a maximum of 296 W m−2 in December and a minimum of 208 W m−2 in August, with a seasonal difference of 88 W m−2. Because downward longwave radiation is emitted by atmospheric vapor, it depends on both humidity and air temperature (Sicart and others, Reference Sicart, Wagnon and Ribstein2005). Even at the same temperature, longwave emission is higher under wet conditions than under dry conditions. This increase in downward longwave radiation is primarily attributed to higher atmospheric moisture and cloud cover during the wet season, which enhance atmospheric longwave emission toward the surface.
Upward longwave radiation, which is emitted from the snow and glacier surfaces, showed a seasonal maximum during the wet season and a minimum during the dry season, with a range of ∼40 W m−2. This seasonal variation was smaller than that of net shortwave radiation. Throughout the year, upward longwave radiation consistently exceeded downward longwave radiation, resulting in negative net longwave radiation. During the wet season, however, persistent cloud cover and high atmospheric humidity increased downward longwave radiation, particularly at night. This suppressed nocturnal radiative cooling and maintained surface temperatures near the melting point throughout much of the diurnal cycle.
During the dry season, however, surface temperatures dropped significantly due to nocturnal cooling, as discussed in Section 4.3. Consequently, the temperature difference between the atmosphere and the surface was larger during the dry season and smaller during the wet season. The seasonal variation in sensible heat flux was primarily driven by this temperature difference between the air and the surface.
Net radiation at the surface reached a seasonal minimum of ∼17–27 W m−2 from December to March, reflecting reduced shortwave input and increased longwave emission during the wet season. A maximum of ∼60 W m−2 occurred in October, prior to the peak of the wet season. Air temperature and surface temperature were both higher during the wet season (austral summer) than during the dry season (austral winter). The sensible heat flux from air to the surface was strongest during the dry season, particularly from June to August, when nocturnal radiative cooling caused surface temperatures to drop well below air temperature, enhancing the temperature gradient and heat transfer.
Both the specific humidity of the atmosphere and the specific humidity at the surface were highest during the wet season and lowest during the dry season. Throughout the year, the specific humidity of the air was lower than that of the surface. This difference results in an upward transport of latent heat via sublimation (Kaser, Reference Kaser2001; Winkler and others, Reference Winkler, Juen, Mölg, Wagnon, Gómez and Kaser2009). In general, the specific humidity at the surface corresponds to saturated conditions and depends on surface temperature. Thus, it was higher during the wet season, similar to the surface temperature. In contrast, air specific humidity depends on both atmospheric moisture and temperature. It was higher during the wet season than the dry season, with a larger seasonal variation compared to surface specific humidity.
The difference in specific humidity between the air and the surface was smaller during the wet season than the dry season. Consequently, the latent heat flux was also lower during the wet season. In this study, both the sensible and latent heat fluxes are defined as positive when directed from the surface to the atmosphere. Annually, sensible heat flux and latent heat flux nearly canceled each other out, implying that most of the net radiation was consumed by conductive heat into the snow and ice. Additionally, part of the net radiation was allocated to melting snow and ice. This melting was highest from October to November, during the period of progressive development toward the core wet season. Although the average albedo in November was moderately high (∼0.61), daily data show several sharp drops in albedo, likely due to transient ice exposure. These short-term events may have locally enhanced net shortwave radiation and melt energy during clear-sky periods, contributing to the seasonal melt peak.
4.3. Diurnal cycles in melt energy
Figure 5 illustrates the diurnal cycles of downward shortwave radiation, downward longwave radiation, melt energy and surface temperature. August 2012, November 2012 and February 2013 are presented as representative examples of the dry season, the period of progressive development toward the core wet season and the core wet season, respectively.
Figure 5. Diurnal cycles in energy fluxes and surface temperature in August, November and February. (a) Downward shortwave radiation, downward longwave radiation and melt energy. (b) Surface temperature.
During the dry season, cloud cover was low and stable. Downward shortwave radiation followed regular diurnal cycles, peaking at ∼1000 W m−2 around noon. In contrast, downward longwave radiation showed no significant periodicity, fluctuating ∼200 W m−2 throughout the day. The low variability in downward longwave radiation during the dry season reflects the persistently clear-sky conditions. Upward longwave radiation consistently exceeded downward longwave radiation, resulting in negative net radiation at night due to the absence of downward shortwave radiation. Consequently, the surface temperature dropped to a minimum of ∼−12°C in the early morning. During the daytime, downward shortwave radiation increased, and the net radiative flux shifted from negative to positive. Most of the net radiative flux was balanced by heat conduction into the snow, causing the surface temperature to rise to 0°C. However, the duration of time when the surface temperature reached 0°C during the day was short, with the maximum energy flux associated with melting reaching ∼100 W m−2 around noon.
In November, during the progression toward the core wet season, cloud cover was greater than during the dry season. The downward longwave radiative flux was stable and relatively high, with a monthly average of ∼269 W m−2 in November. The surface temperature drop at night due to clear-sky conditions followed a pattern similar to that of the dry season, but the minimum surface temperature was higher. This was due to the increased downward flux of longwave radiation, which resulted in a higher net radiative flux compared to the dry season. The duration of time when the surface temperature reached 0°C during the day was longer than in the dry season, and a greater portion of the net radiation was allocated to the melting of snow and ice. However, nighttime cooling caused the surface temperature to drop below 0°C, leading to a repeated pattern of daytime melting and nighttime freezing.
During the core wet season, clear-sky conditions were rare, and cloud cover was frequent. As a result, the downward longwave radiation was stable, with a monthly average of approximately 295 W m−2 in February. The frequency of surface temperatures dropping below 0°C at night was low, and the long periods when surface temperatures remained at 0°C led to extended melting periods. However, because downward shortwave radiation was lower than in November, the total heat allocated to melting was also lower than during November.
4.4. Identification of meteorological variables in empirical estimates of melt energy
Table 1 presents the correlation coefficients between melt energy estimated by the energy balance model and five key meteorological variables: air temperature, specific humidity, relative humidity, downward shortwave radiation and net shortwave radiation. To more specifically assess the meteorological controls on melt processes, we also analyzed correlations using a subset of data where melt energy exceeded 0.1 W m−2. This threshold was used to exclude periods of surface freezing (M < 0) and near-zero melt conditions, thereby isolating actual melting events. Among these, downward shortwave radiation showed the strongest correlation with melt energy (R = 0.62, p < 0.001), followed by net shortwave radiation (R = 0.58, p < 0.001) and air temperature (R = 0.43, p < 0.001).
Table 1. Correlation coefficients between melt energy and individual explanatory variables
* Note:p < .05, **p < .01, ***p < .001.
Additionally, for reference, the table includes correlation coefficients for three radiation-related variables: net radiation, downward longwave radiation and net longwave radiation. Among these, net radiation exhibited the highest correlation with melt energy (R = 0.86, p < 0.001), whereas downward and net longwave radiation showed relatively weak correlations.
In contrast, specific humidity and relative humidity did not show significant correlations with melt energy when analyzed over the entire dataset. However, in the analysis of 3 month intervals, both variables demonstrated statistically significant correlations with melt energy during the April–June period, with higher correlation coefficients than in any other interval.
Table 2 presents the results of a multiple correlation analysis using two explanatory variables for melt energy. The highest multiple correlation coefficient was observed for the combination of net shortwave radiation and specific humidity (R = 0.85, p < 0.001 for both net shortwave radiation and specific humidity), followed closely by the combination of net shortwave radiation and relative humidity (R = 0.82, p < 0.001 for both variables).
Table 2. Multiple correlation coefficients for melt energy with two explanatory variables
Notes: The top row in each cell represents the analysis using the total dataset, while the bottom row represents the analysis using data where melt energy > 0.1 W m−2. * p < .05, **p < .01, ***p < .001. No variables in this table reached significance at the *p < .05 level.
Table 3 shows the results of the multiple correlation analysis using three explanatory variables for melt energy. The combination of net shortwave radiation, specific humidity and air temperature resulted in the highest multiple correlation coefficient (R = 0.85, p < 0.001 for both net shortwave radiation and specific humidity, and p > 0.05 for air temperature), followed by the combination of net shortwave radiation, relative humidity and air temperature (R = 0.83, p < 0.001 for both net shortwave radiation and relative humidity, and p > 0.05 for air temperature). Across most time periods, net shortwave radiation showed slightly stronger correlations with melt energy than downward shortwave radiation. This highlights the importance of accounting for surface albedo, which influences the net shortwave component and thus the actual energy available for melt.
Table 3. Multiple correlation coefficients for melt energy with three explanatory variables
Notes: The top row in each cell represents the analysis using the total dataset, while the bottom row represents the analysis using data where melt energy > 0.1 W m−2. *p < .05, **p < .01, ***p < .001. No variables in this table reached significance at the **p < .01 level.
4.5. Comparison and validation of empirical melt estimates
This paper proposes an empirical melt model for tropical glaciers based on a linear regression approach, in which net shortwave radiation and specific humidity are used as explanatory variables. While the functional structure follows that of commonly used empirical models, such as the DDM and ETI models, this study derives model coefficients by regressing against melt energy estimated from an energy balance model (TSMM). Figure 6 illustrates the seasonal variations in melt energy derived from the energy balance model and three empirical melt models: the DDM, the ETI model and the proposed model. The regression coefficients for these models were determined through regression analysis based on energy balance model outputs (Table 4). The RMSE values were 7.1 W m−2 for the DDM, 4.9 W m−2 for the ETI model and 2.1 W m−2 for the proposed model, indicating that the proposed model had the smallest RMSE and provided the most accurate estimates.
Figure 6. Comparison of empirical melt models (degree-day, enhanced temperature-index and proposed) with the energy balance model. The proposed model uses net shortwave radiation and specific humidity as explanatory variables. Explanatory variables and RMSE values for each model are summarized in Table 4.
Table 4. Coefficients and RMSE values for the empirical melt models
Notes: Coefficients were derived by regressing melt energy estimated by the energy balance model on the selected explanatory variables. The formulations adopt a consistent linear structure for comparison, rather than their original formulations with temperature thresholds or nonlinear radiation terms. The proposed model incorporates net shortwave radiation and specific humidity as explanatory variables.
Regression coefficients 
$\textstyle a_1$, 
${a_2}$, and 
$b$ are derived from regression analysis based on energy balance model outputs. Coefficients marked as dimensionless do not have physical units. RMSE represents the root-mean-square error (W m−2) between the empirical and energy balance model estimates.
5. Discussion
5.1. Application and validation of the TSMM
As demonstrated in Figure 3, the TSMM successfully captured the diurnal cycles of surface temperature throughout the year. Given that surface temperature fluctuations reflect the surface energy balance, this supports the reliability of the TSMM in simulating glacier energy exchanges. Although originally developed for seasonal snowpacks, the TSMM was extended in this study to represent glacier conditions by incorporating physical properties specific to ice, such as albedo, density and thermal conductivity, as described in Section 2.2. One of the key advantages of the TSMM is that it not only resolves the surface energy balance but also accounts for vertical heat conduction and solar radiation penetration, which affect surface temperature and associated fluxes such as sensible heat and longwave radiation. The extended TSMM thus offers a physically consistent framework for assessing glacier energy balance in the outer tropics, where surface temperature undergoes significant diurnal variation. This methodological refinement enhances our ability to interpret melt processes in tropical glacier environments.
Some discrepancies were found when comparing TSMM simulations with observed surface temperatures, particularly during nighttime in the dry season and daytime in the wet season. The surface temperature was derived from upward longwave radiation, which may be influenced by atmospheric conditions or sensor characteristics. Therefore, values slightly exceeding 0°C do not necessarily indicate physically unrealistic surface temperatures, but rather may include limitations in longwave radiation measurement-based temperature estimation (Sicart and others, Reference Sicart, Hock, Ribstein, Litt and Ramirez2011). Considering the limitations of longwave-based surface temperature estimates—particularly under nighttime conditions with strong radiative cooling—the modeled temperature may better reflect actual surface conditions during these periods.
5.2. Characteristics of surface energy balance
The energy balance characteristics observed at Huayna-Potosí Glacier align closely with those reported for the nearby Zongo Glacier in Bolivia. Both glaciers, located in the outer tropics, exhibit melt regimes strongly governed by net radiation. Sicart and others (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and Corripio2005) first highlighted this radiative dominance at the Zongo Glacier, and more recently, Autin and others, Reference Autin, Sicart, Rabatel, Soruco and Hock2022) showed that the transition season—particularly October and November—plays a central role in controlling the annual surface mass balance. During this period, clear-sky conditions and moderately low albedo allow shortwave radiation to dominate, while enhanced downward longwave radiation from moist and cloudy skies also contributes significantly to melt energy.
At Huayna-Potosí Glacier, the seasonal pattern reflects these same dynamics. Most of the net radiation was consumed by conductive heat flux into the snow and ice, with the remaining energy used for melting. Net radiation increased during October and November, contributing to a seasonal peak in melt energy prior to the onset of the core wet season. Although downward shortwave radiation decreased after December, upward shortwave radiation increased due to the influence of fresh snow and high albedo. These seasonal behaviors align with observations at other tropical glaciers such as the Artesonraju Glacier in Peru (Gurgiser and others, Reference Gurgiser, Marzeion, Nicholson, Ortner and Ortner2013), where net radiation similarly drives melt while turbulent heat fluxes play a limited role due to small seasonal temperature variations.
Diurnal-scale analysis using the TSMM provided further insight into how surface energy balance processes differ between seasons. During the dry season, clear skies and low atmospheric humidity resulted in reduced downward longwave radiation at night, while surface longwave emission remained high. This radiative imbalance caused significant nighttime cooling of the glacier surface. Although solar radiation was strong during the day, much of the net energy was consumed by conductive flux to raise surface temperatures to 0°C, limiting the duration of melt. In contrast, during the wet season, persistent cloud cover and increased atmospheric moisture enhanced downward longwave radiation, reducing nocturnal cooling and allowing the surface to remain at 0°C for longer periods. These effects were successfully captured by the TSMM, which integrates surface energy balance, vertical heat conduction and solar radiation penetration.
Broader comparisons across the tropical Andes further support the relevance—and limitations—of radiation-driven melt models. Fyffe and others (Reference Fyffe2021) investigated glaciers in both northern and southern Peru (8°–15°S), and found that southern Peruvian glaciers—closer in latitude and climate to the Huayna-Potosí Glacier—exhibited similar seasonal melt behavior, with enhanced melt energy during the dry-to-wet season transition. Ayala and others (Reference Ayala, Pellicciotti, MacDonell, McPhee and Burlando2017) examined glacier energy balance across north-central Chile (29°–34°S). In the northern part of their study region (∼29°S), which is climatologically close to the Bolivian Andes, low-elevation glaciers were dominated by net shortwave radiation–driven melt. However, at higher elevations, sublimation became the dominant process due to cold, dry and windy conditions. These findings demonstrate both the spatial coherence and the environmental limits of applying net-radiation–based melt models across the outer tropical Andes.
5.3. Key meteorological variables for empirical melt models
As shown in Table 1, downward shortwave radiation (or net shortwave radiation) exhibited the strongest and most significant correlation with melt energy, while air temperature showed a moderate yet statistically significant correlation. This indicates that downward shortwave radiation (or net shortwave radiation) is the most influential variable for empirical melt estimates of tropical glaciers. Neither specific humidity nor relative humidity exhibited significant correlations with melt energy when analyzed across the entire dataset. However, in the analysis of 3 month intervals and for melt energy exceeding 0.1 W m−2, specific humidity and relative humidity demonstrated statistically significant correlations with melt energy during the April–June period. During this time, the correlation coefficients for downward shortwave radiation, net shortwave radiation and air temperature were comparatively lower. This suggests that longwave emissions from atmospheric vapor contribute to melt energy by maintaining surface temperatures near 0°C and creating conditions conducive to melting.
In addition, net radiation, downward longwave radiation and net longwave radiation were evaluated as potential explanatory variables. Among them, net radiation exhibited the highest correlation with melt energy (R = 0.86, p < 0.001), whereas downward and net longwave radiation showed relatively low correlations. However, due to the limited availability of direct longwave radiation measurements in many glacierized regions, these variables were not incorporated into the empirical melt model. Instead, net shortwave radiation and specific humidity were selected as explanatory variables, as they are more commonly observed and can effectively represent the radiative energy that drives melt.
The low seasonal variation in air temperature, the dominant contribution of downward shortwave radiation to melt energy outside of January to March, and the significant contribution of downward longwave radiation during January to March are characteristic of the meteorological conditions and energy balance processes of tropical glaciers. These distinctive characteristics collectively reduce the applicability of empirical melt models for tropical glaciers that rely on a single explanatory variable, such as air temperature, as illustrated in Figure 6 (La Frenierre and Mark, Reference Frenierre and Mark2014; Juen and others, Reference Juen, Kaser and Georges2007).
The strong correlations observed between downward longwave radiation and both specific and relative humidity (see Table 1) suggest that humidity-based variables can effectively serve as proxies for atmospheric longwave radiation in empirical melt models. Specific humidity also showed a strong correlation with relative humidity (R = 0.96), indicating that either variable can be effective, depending on data availability. In particular, relative humidity explained most of the variability in downward longwave radiation, supporting its potential utility when specific humidity data are unavailable.
Statistical analysis showed that downward longwave radiation was strongly correlated with relative humidity (R = 0.96, p < 0.001), whereas air temperature showed only a moderate correlation (R = 0.50, p < 0.001). A multiple regression using both variables yielded a correlation coefficient of R = 0.96; however, only relative humidity was statistically significant. Specific humidity was also strongly correlated with relative humidity (R = 0.96), indicating that either variable can serve as a proxy for downward longwave radiation, depending on data availability.
Furthermore, when used in empirical melt models, relative humidity combined with net shortwave radiation produced a similar multiple correlation coefficient for melt energy as specific humidity. This suggests that relative humidity could serve as a practical alternative to specific humidity in model applications, especially where direct measurements of specific humidity are unavailable.
As shown in Table 3, the multiple correlation analysis with three explanatory variables revealed that the combination of net shortwave radiation, specific humidity and air temperature produced the highest multi correlation coefficient. However, this value was nearly identical to the coefficient obtained using only two explanatory variables: net shortwave radiation and specific humidity. Furthermore, in the regression analysis, net shortwave radiation and specific humidity were statistically significant explanatory variables, while air temperature was not. This finding highlights the relatively low contribution of air temperature to melt energy in tropical glaciers in the outer tropics.
Based on these findings, this paper proposes that the combination of net shortwave radiation and specific humidity represents the most suitable set of explanatory variables for empirical melt estimates of tropical glaciers. The proposed empirical melt model is best suited for tropical glacier environments where seasonal variations in air temperature are small, but fluctuations in atmospheric moisture and cloud cover are substantial. Under such conditions, radiative fluxes—particularly longwave radiation—play a dominant role in melt energy, and relative humidity can serve as a meaningful proxy. In contrast, in mid-latitude or temperate regions with greater seasonal variability in air temperature, the contribution of sensible heat flux tends to increase. For example, Ayala and others (Reference Ayala, Pellicciotti, MacDonell, McPhee and Burlando2017) observed that the role of sensible heat becomes more pronounced in the southern Andes of Chile as one moves away from the tropics. This suggests that while the model may be transferable to some environments outside the tropics, its applicability may be limited in regions where turbulent fluxes exert a stronger influence on the energy balance.
5.4. Empirical melt estimates of tropical glaciers
The ETI model, an extension of the DDM, is widely used for practical melt estimation with two explanatory variables: net shortwave radiation (or downward shortwave radiation) and air temperature. In Figure 6, the ETI model exhibited seasonal variations similar to those of the proposed model but underestimated melt energy in October, when melt energy was at its peak, and overestimated it after March. These findings highlight the importance of incorporating specific humidity as an explanatory variable in empirical melt models, particularly for regions with minimal seasonal temperature variation where air temperature alone fails to adequately capture melt energy variability.
The ETI model appears to underestimate the contribution of atmospheric longwave radiation in October and shows heightened sensitivity to air temperature after March. Fuchs and others (Reference Fuchs, Asaoka and Kazama2016), who applied the ETI model to the Zongo Glacier watershed, reported that the model’s performance was inadequate under cloudy conditions. The results of this study suggest potential improvements to the ETI model by incorporating vapor-related information, which is closely linked to atmospheric longwave radiation. This demonstrates the critical role of atmospheric vapor in enhancing the accuracy of empirical melt models for tropical glaciers.
The proposed empirical melt model, incorporating net shortwave radiation and specific humidity as explanatory variables, is based on the energy balance processes characteristic of tropical glaciers. While the results of this study are derived from observations on the Huayna-Potosí Glacier near its equilibrium line or the upper boundary of the ablation area, the approach has the potential to be adapted to other tropical glaciers in the outer tropics. Addressing the applicability of the proposed model to other parts of the glacier, including areas near the glacier front, as well as to surrounding glaciers, provides a robust framework for improving melt estimates across a wider range of tropical glaciers. Several global glacier modeling studies, including those conducted within the Glacier Model Intercomparison Project (GlacierMIP; Hock and others, Reference Hock2019), have relied on DDM or ETI models for simulating glacier melt. While these models offer global applicability, they assume that temperature is the dominant control on melt energy. This assumption may be less valid in tropical regions, where radiative and atmospheric moisture-related effects are more influential. As noted in previous regional assessments (e.g. La Juen and others, Reference Juen, Kaser and Georges2007; Frenierre and Mark, Reference Frenierre and Mark2014), temperature-index models may not fully capture the processes governing ablation in tropical glaciers. Our results provide region-specific insights that may help reduce this uncertainty in outer tropical settings.
Additionally, the model’s reliance on specific humidity or relative humidity, both of which can serve as proxies for atmospheric longwave radiation, offers practical advantages for regions where direct measurements of longwave radiation are unavailable. This adaptability improves the model’s utility as a tool for regional water resource management and global-scale glacier melt analyses. Future research should explore the applicability of this approach under extreme climatic conditions, such as those associated with El Niño events. Previous studies have shown that ENSO can significantly alter the energy balance of tropical glaciers by influencing cloud cover, precipitation and air temperature (Wagnon and others, Reference Wagnon, Ribstein, Francou and Sicart2001; Autin and others, Reference Autin, Sicart, Rabatel, Soruco and Hock2022). Related impacts include increased net radiation, enhanced melt and glacier mass loss, as well as water resource shortages in affected regions (Francou and others, Reference Francou, Vuille, Wagnon, Mendoza and Sicart2003; Morizawa and others, Reference Morizawa, Asaoka, Kazama and Gunawardhana2013; Veettil and others, Reference Veettil, Bremer, de Souza, Bayer Maier and Simoes2016; Panozo Rivero and others, Reference Panozo Rivero, Morais and de Sousa Pereira2020; Canedo-Rosso and others, Reference Canedo-Rosso, Hochrainer-Stigler, Pflug, Condori and Berndtsson2021). Understanding how the proposed model performs under such climate anomalies could help evaluate its robustness and improve its applicability for long-term hydrological and glaciological assessments in the outer tropics.
6. Conclusions
This study investigated the energy balance characteristics of the Huayna-Potosí Glacier, Bolivia (16°S) using meteorological data and the TSMM, which was extended to incorporate glacier-specific properties. The dominant meteorological variables suitable for empirical melt models were identified based on results from the energy balance model. Additionally, existing empirical melt models were validated under outer tropical conditions, and a more suitable empirical model, based on the characteristics of energy balance, was proposed for future applications. The key findings are summarized as follows:
Energy balance characteristics: The results demonstrated that net radiation is the primary driver of melt energy. Seasonal variations showed an increased contribution of downward longwave radiation and a decreased contribution of downward shortwave radiation from the dry season to the wet season. Diurnal variations revealed mechanisms through which atmospheric longwave emissions suppress melting during the dry season and enhance it during the wet season, which underscores the dynamic interaction between atmospheric and surface conditions.
Key variables for empirical melt estimates: Net shortwave radiation and specific humidity (or relative humidity, in cases where specific humidity data are unavailable) were identified as suitable meteorological variables in the empirical melt estimates of tropical glaciers. These variables can serve as proxies for atmospheric longwave radiation. The proposed empirical melt model incorporates net shortwave radiation and specific humidity as explanatory variables. It outperformed widely used DDM and ETI models, demonstrating seasonal variations in melt energy consistent with those estimated by the TSMM.
The insights gained from the energy balance analysis directly informed the development of the proposed empirical melt model, ensuring its alignment with physical processes specific to tropical glaciers. The simplicity and adaptability of the proposed model make it a practical tool for both regional water resource planning and global-scale climate modeling efforts.
While this study focused on a single observation site on the Huayna-Potosí Glacier near its equilibrium line, future research should evaluate the model’s performance across other parts of the glacier, such as near the glacier front, and on other tropical glaciers. Investigating the model’s robustness under extreme climatic events, such as El Niño, is another critical area for further study. The insights gained from this research have the potential to not only improve water resource management at regional scales but also reduce uncertainties in global glacier melt analyses, particularly for tropical glaciers.
Acknowledgements
This study was financially supported by JST-JICA SATREPS, JSPS KAKENHI (17K06587 and 25K07944) and JSPS Bilateral Joint Research Program (JPJSBP120199979 and JPJSBP120259945). The meteorological data analyzed in this study were collected by the GRANDE project, a Japan-Bolivia joint project conducted between 2010 and 2015. We would like to express our gratitude to Sachie Hayakawa at Tohoku University and Hiro Konno at Nihon University for their invaluable assistance.
Competing interests
The authors declare no competing interests.
Open access
