A method for monitoring glacial loss and temperature variation using satellite observations: Case study of Pico de Orizaba and Iztaccíhuatl (Mexico)

ABSTRACT

Throughout the world, tropical glaciers are rapidly receding and some are at risk of complete loss within the next several decades. It is therefore important to monitor these glacial regions to better understand the factors affecting glacier loss and risks to nearby communities. Here, we provide an update on the summit glacier extents of Pico de Orizaba (19.03°N, 97.27°W) and Iztaccíhuatl (19.18°N, 98.64°W) in central Mexico, reporting areas every 1–11 years between 2001 and 2019 (depending on the availability of high-quality imagery). Glacier extents are derived from multispectral satellite imagery (Landsat 7–8 and the European Space Agency’s Sentinel-2) using a semi-automated mapping method that takes the ratio of the near-infrared (or panchromatic) and shortwave infrared bands and applies appropriate threshold(s) for glacier identification. We also use thermal band imagery from Landsat 7–8 to estimate surface temperatures of both the glaciers and the surrounding terrain to provide a more comprehensive understanding of the summit environment. We find that glacier retreat has continued on both summits, and surface temperatures (even over glacier ice) can be at or near the melting point of water during some parts of the year, particularly on lower-elevation Iztaccíhuatl, suggesting a risk of continued ice loss in the future.

KEYWORDS:

• Glacier retreat
• surface temperature
• Landsat
• Mexico

Introduction

Glaciers can be found even in the Tropics (±20° latitude) at sufficiently high elevations. These represent significant sources of water for nearby communities but can also pose potential hazards. Equatorial glaciers have been a focus of studies characterizing glacial loss (Paul et al. 2004) and hazard analyses (Kääb et al. 2005), where events such as glacier outbursts, enhanced runoff, slope instability, mass movement, and interactions between ice and volcanic activity (e.g., lahars) may threaten downslope populations. Low-latitude glacial extents are restricted to the highest mountain regimes of Asia (Indonesia), Africa (e.g., Mt. Kilimanjaro in Tanzania), South America (including Ecuador, Peru, and Colombia), and North America (Mexico). Worldwide, alpine glaciers are particularly susceptible to temperature changes associated with climate change, because they are above the treeline and are subject to unobstructed influxes of solar radiation and positive feedbacks associated with their high albedos (Haeberli and Beniston 1998). In the case of tropical glaciers, changes in both air humidity and temperature (not just temperature alone) are thought to be drivers of glacier recession (Kaser 1999). In both mid- and low-latitude alpine regions, retreat has been observed in recent decades (e.g., Francou et al. 2005; Cullen et al. 2006; Soruco et al. 2009; Thompson et al. 2009, 2011), and it is possible that some of these glaciers will be lost completely within the next several decades (Zemp et al. 2006).

Here, we focus on two glacierized summits in south-central Mexico’s Cordillera Neovolcánica: Pico de Orizaba and Iztaccíhuatl (Figure 1). The Cordillera Neovolcánica spans the southern edge of the North American plate and contains the highest peaks in Mexico, all dormant or active stratovolcanoes. The glaciers in this region are unique because they are the only permanent ice masses found near the 20°N latitude band. Although not strictly characterized as tropical glaciers because they lie outside of the oscillation range of the Intertropical Convergence Zone, Mexico’s glaciers are located within the latitudinal bounds of the Tropics and meet the thermal criteria for the Tropics (within the range where diurnal temperature variations exceed annual temperature variations; Kaser and Osmaston 2002). Due to the circulation of arid westerly winds during northern winter/spring and the northeasterly trade winds that bring in moist air during northern summer/fall, these glaciers do have seasonal patterns similar to those in some parts of the Tropics, with distinct humid and dry seasons (Kaser and Osmaston 2002). It is therefore plausible that these seasonal variations in atmospheric humidity, in addition to temperature, may play a role in driving the processes affecting these glaciers today. Some work has indicated that the glaciers in Mexico may not last beyond 2050 (Delgado-Granados et al. 2007), making this a critical region for monitoring glacier losses in real time; many other mountain regions are suffering from similar glacier retreat (e.g., Paul et al. 2004). Previous mapping efforts in these regions used aerial photography, fieldwork, orbital remote sensing, or a combination of these. The most recent updates are largely performed with orbital multispectral imagery using band ratio calculations that are shown to provide a rapid, robust means of identifying glacier boundaries (e.g., Schneider et al. 2008; Cortés-Ramos and Delgado-Granados 2012).

Figure 1. Map of the locations of the two glacierized summits in central Mexico: Pico de Orizaba (top right) and Iztaccíhuatl (lower left).

Pico de Orizaba (19.03°N, 97.27°W), at 5,636 m, is the highest point in Mexico and contains the country’s largest glacier. It is also notable for having one of the highest treelines in the world (Cruz-Kuri, McKay, and Navarro-González 2001). Previous mapping of the glaciers on Pico de Orizaba reported a total areal extent of 9.5 km² in 1958 (Lorenzo 1964; White 1978); however, cartographic inconsistencies noted by Cortés-Ramos and Delgado-Granados (2015) indicated that the mid-twentieth-century areal estimates were, in fact, significantly overestimated. By correcting these errors via rescaling of the glacier outlines to align with the INEGI (2002) topographic map, an updated estimate of the 1958 glacier coverage was found to be 2.24 km² (Cortés-Ramos and Delgado-Granados 2015). In addition, Cortés-Ramos and Delgado-Granados (2015) constructed an outline of the glacier extent in 1975 based on aerial photography and found that the summit glaciers had lost approximately 24 percent of its area between 1958 and 1975. Historically, there were seven outlet glaciers extending from the Gran Glaciar Norte ice cap, along with Glaciar Oriental, a detached mountain glacier extending east from the summit (Lorenzo 1964; White 1978). The Chichimeco glacier tongue, which extended north from the ice cap, lost significant area between 1959 and 1975 (Cortés-Ramos and Delgado-Granados 2015) and was completely gone by 1988 (Palacios and Vázquez-Selem 1996). Extending northwest from Glaciar Norte was Glaciar de Jamapa (Lorenzo 1964; White 1978), which similarly experienced fairly significant losses by 1975 (Cortés-Ramos and Delgado-Granados 2015) and was no longer attached to the larger ice cap by the mid-1990s, limited in extent to two patches of ice downslope from the rest of the glacier (Palacios and Vázquez-Selem 1996). Several additional glacial tongues extended west and southwest from Glaciar Norte (Glaciar Occidental and Suroccidental, respectively) and experienced substantial losses between 1959 and 1996 (Palacios and Vázquez-Selem 1996). In 2001, the area of Glaciar Occidental was ~0.0895 km², but by 2007 it has lost 74 percent of its area (Cortés-Ramos and Delgado-Granados 2012). Due to this retreat, most of the ice mass that remained on Pico de Orizaba by 1996 was limited to the ice cap, extending to a minimum elevation of 5,037 m (Palacios and Vázquez-Selem 1996). In 2001, the area of the Glaciar Norte ice cap was 0.93 km², which shrank by approximately one third between 2001 and 2007 (Cortés-Ramos and Delgado-Granados 2013). Glaciar Oriental, which had a minimum altitude of 5,070 m and an area of roughly 0.118 km² in 1958 (Lorenzo 1964; Palacios and Vázquez-Selem 1996; Cortés-Ramos and Delgado-Granados 2015), had remained more constant in size, with little apparent area loss between 1958 and 1988 (Palacios and Vázquez-Selem 1996). In 1996, Glaciar Oriental had a minimum altitude of 5,100 m (Palacios and Vázquez-Selem 1996). Between 2002 and 2007, only approximately 13 percent of its glacial area was lost (Cortés-Ramos and Delgado-Granados 2012).

The summit of Iztaccíhuatl (19.18°N, 98.64°W) once had a large glacial area, covering an area of ~6.4 km² at the peak of the Little Ice Age (around 1850) based on mapping of the associated moraines (Schneider et al. [2008], using the geologic maps of Nixon [1989]; INEGI 1995). The earliest detailed map of glacier extent at Iztaccíhuatl combined aerial photographs from 1945 with fieldwork likely conducted in 1958 to obtain glacier coverage of 1.21 km² across twelve named glaciers and two firn fields (Lorenzo 1964; White 1978). These glaciers spanned several subpeaks of Iztaccíhuatl; from northwest to southeast, relevant subpeaks are named La Cabeza, El Pecho, and La Panza. Several smaller ice masses were present in 1959 (Lorenzo 1964; White 1978): Glaciar de la Cabeza, a thin, 200-m-long ice mass on the northwest side of the La Cabeza summit (1959 area = 0.0144 km²); the Nororiental firn field southeast of the El Pecho summit (1959 area = 0.025 km²); Glaciar de San Augustin (1959 area = 0.0113 km²), southeast of La Panza; and Atzintli (1959 area = 0.0575 km²), which prior to 1918 was connected to Glaciar de Ayoloco, part of the La Panza glacier system (White 1978). However, many of these smaller masses have disappeared completely since the observations of Lorenzo (1964). The San Augustin glacier was lost by 1982, the glaciers on and around the La Cabeza summit were gone by 1994, and the Atzintli glacier had receded completely by 2004 (Schneider et al. 2008). Over time, retreat has occurred such that only two glacierized regions remain, concentrated at the El Pecho and La Panza subpeaks. The total glacier area at the Iztaccíhuatl summit in 2007 was 0.273 km² (Schneider et al. 2008), where El Pecho and La Panza had glacial extents of 0.14 and 0.133 km², respectively. This reflects a loss of about 80 percent of the glacier area between the 1959 maps of Lorenzo (1964) and 2007. Additionally, using the contours of the early glacier mapping and estimating the surface elevations using ground Global Positioning System measurements and orbital digital terrain models (from the Shuttle Radar Topography Mission, Shuttle Radar Topography Mission, and Terra/Advanced Spaceborne Thermal Emission and Reflection Radiometer [ASTER] Channel 3 images), it was observed that thinning of the remaining glaciers had occurred between 1959 and 2000, where El Pecho lost ~70 m in elevation and La Panza lost ~30 m (Schneider et al. 2008). Based on the degree of deglaciation observed, it was further surmised that total deglaciation was likely within several decades, perhaps as early as 2020 (Schneider et al. 2008).

Historically, the northern slope of Popocatépetl (19.02°N, 98.63°W), ~17 km south of Iztaccíhuatl, was also glaciated. In 1958, these glaciers covered an area of 0.729 km² (Lorenzo 1964), comprising a firn field with three northward-extending glaciers (White 1978). However, these glaciers experienced rapid retreat through the latter half of the 1990s into the early 2000s. Much of this loss was driven by an increase in volcanic activity at the summit (Delgado-Granados et al. 2007), which began in 1994 with a series of explosive volcanic events and, after 1996, alternating effusive and explosive events that accelerated ice retreat. The deposition of tephra on the glacial ice likely played a key role in increasing melting (White 2002; Julio-Miranda et al. 2008). Through the use of pairs of aerial photographs and photogrammetric techniques, Andrés et al. (2007) estimated the 3D ice loss at several time points, where almost 4 × 10⁶ m³ of volumetric loss was recorded between November 1997 and December 2002. By 2001, although ice masses were present on part of the summit, the ice no longer flowed and the glaciers were thus declared extinct (Delgado-Granados et al. 2007). Although the increase in volcanic activity (which has largely continued to the present day) is thought to be the primary driver for the glacier extinction at Popocatépetl, it has been hypothesized that climate change and air pollution due to its proximity to Mexico City may have played a small role (Delgado-Granados 2007).

Though neither Pico de Orizaba nor Iztaccíhuatl is currently erupting, they are not presently classified as extinct. Iztaccíhuatl is considered dormant, because there is no historic record of significant volcanic activity and much of the lavas at the summit are thought to be Pleistocene in age (Nixon 1989). Pico de Orizaba has experienced several ash eruptions and some fumarolic activity in recorded history (between 1530 and 1870 C.E.; Carrasco-Nùñez and Rose 1995; De la Cruz-Reyna and Carrasco-Nùñez 2002) and is considered active but presently dormant (Macías 2007). Therefore, future volcanic activity cannot be ruled out as a potential factor affecting the remaining glaciers in the upcoming decades.

Regions of temperate climate surround both Pico de Orizaba and Iztaccíhuatl, although at higher elevations (>4,300 m on Pico de Orizaba and >4,000 m on Iztaccíhuatl) tundra and ice cap climates predominate (White 1978). At lower elevations surrounding Pico de Orizaba, it is typically drier during the winter and rainier between late spring and early autumn (White 1978). Near Iztaccíhuatl, rain is more prevalent year-round at lower elevations, although winter is still typically the driest season (White 1978). Though long-term climate records do not appear to be available at either summit (within the tundra/ice cap climate regimes), weather stations are present at Sierra Negra (7 km southwest of Pico de Orizaba and over a kilometer lower in elevation). This summit area, home to the Large Millimeter Telescope and the High-Altitude Water Cherenkov Gamma Ray Observatory, is classified as a tundra environment (Carramiñana et al. 2007), and the annual average temperatures from 2000–2008 were at or just above 0°C (Carrasco et al. 2009).

Here, we aim to update the glacier monitoring of these summit regions over the last two decades, with a focus on variations since the 2007 update of Pico de Orizaba (Cortés-Ramos and Delgado-Granados 2013) and the 2008 update of Iztaccíhuatl (Schneider et al. 2008). We use the method of band ratios in order to provide a more robust understanding of the current and recent glacier changes in these areas. Some recent mapping of these regions was performed using the Normalized Difference Snow Index and a range of different Landsat and Sentinel-2 images over 15-year intervals (Veettil and Wang 2018); however, studies have shown that using band ratios, along with an additional threshold to help identify ice in shadows for Sentinel-2 images, provides more accurate results (Paul et al. 2016). We further explore higher-resolution temporal variations, comparing outlines derived from images collected 1–6 years apart (depending on the availability of high-quality, low-snow images). In addition, we use thermal-band Landsat images to estimate surface temperatures in order to obtain a more comprehensive understanding of the ground conditions within the summit areas.

Methods

Multispectral satellite imagery over the regions of interest was obtained from the U.S. Geological Survey (USGS) using the GloVis tool (glovis.usgs.gov; U.S. Department of the Interior/U.S. Geological Survey 2005). Several different instruments were used in this analysis:

1. Enhanced Thematic Imager Plus (ETM+) on Landsat 7. Visible, near-infrared (NIR), and short-wave infrared (SWIR) bands have a resolution of 30 m, the thermal infrared band resolution is 60 m (resampled to 30 m), and the panchromatic band resolution is 15 m. In operation since 1999, permitting a longer temporal observation range but limited by the 31 May 2003 failure of the Scan Line Corrector that led to significant gaps in the data coverage of images collected after this date.

2. Operational Land Imager (OLI) and Thermal Infrared Sensor on Landsat 8. The visible, NIR, and SWIR band resolution is 30 m, and the panchromatic band resolution is 15 m. The thermal infrared band resolution is 100 m (resampled to 30 m). In operation since 2013.

3. Multispectral Imager (MSI) on Sentinel-2 (European Space Agency). The visible and NIR band resolution is 10 m, and SWIR resolution is 20 m. Sentinel-2A has been in operation since 2015; Sentinel-2B has been in operation since 2017. Permits slightly higher-resolution observations but only in recent years.

For all instruments, cloud-free or nearly cloud-free images were identified, and only those where the summit region was not obscured were used for further analysis. Due to the frequency of re-imaging inherent to these orbital platforms (e.g., each Landsat orbiter covers the Earth’s surface every 16 days), it was possible to obtain many cloud-free images at each site, although, as expected, most of these images fall within the drier months except as noted.

All mapping was performed using the QGIS open-source software suite (http://qgis.osgeo.org; QGIS Development Team 2018), with some data processing steps performed with GRASS GIS  (http://grass.osgeo.org; Grass Development Team 2018). The two primary tasks performed were semi-automated glacier mapping and surface temperature estimation (using thermal infrared imagery). In addition, the ASTER GDEM2 tiles over the study areas were obtained from the National Aeronautics and Space Administration’s (NASA) Earthdata (earthdata.nasa.gov) and used to produce elevation profiles in each summit area and to identify the lowest-elevation glacier extents at several different time steps. In general, the spatial resolution of the ASTER GDEM2 is ~75 m/pixel, and the vertical root mean square error (RMSE) is 8.68 m (ASTER GDEM Validation Team 2011). However, accuracy errors vary by location and by the number of images used to generate the GDEM at each pixel, where RMSEs are higher when fewer images are used (less than ten) but stabilize at scene counts above approximately fifteen (ASTER GDEM Validation Team 2011). Without obtaining higher resolution topography, we cannot validate the local RMSE but will make approximations of the accuracy based on scene counts.

Semi-automated glacier mapping

Because ice and snow are highly reflective in the red-infrared to NIR parts of the spectrum (Dozier 1989) and highly absorptive in the SWIR, taking the ratio of these bands has proven to be a highly robust method for identifying the boundaries of glaciers (Albert 2002; Paul et al. 2002; Paul and Kääb 2005). Historically, data from ASTER were one of the primary sources for this type of glacier mapping (e.g., Kääb et al. 2002; Svoboda and Paul 2009; Gjermundsen et al. 2011). However, with the failure of the onboard SWIR detectors in 2008, ASTER scenes collected after that point cannot be used for glacier mapping by this method. We therefore focus on Landsat data, supplemented with some more recent imagery from the European Space Agency’s Sentinel-2, in this study. The panchromatic bands of both ETM+ and OLI, which include both the red-infrared and NIR wavelengths, can be used in lieu of the narrower bands to obtain higher-resolution glacier boundaries (Paul et al. 2016). Due to the resolution of the instruments of interest (15 m at best using the panchromatic band of the Landsat imagers, 10 m for Sentinel-2), variations must be on the order of tens of meters in order to be detected using this method. Though month-to-month and/or interannual variations rarely meet this threshold, it is nevertheless a viable method for detecting glacier variations over time spans of several years to decades.

In all cases, the lower-resolution SWIR band was resampled such that both bands used have an equivalent resolution. The Raster Calculator was used to take the ratio of the NIR or panchromatic band to the SWIR band using raw digital numbers, as recommended by Paul et al. (2002). A threshold ratio value for glacier classification was identified manually, using a false-color composite (NIR–red–green) for comparison. Typical glacier ratio thresholds are between two and three (Paul et al. [2016] and references therein). For Sentinel-2 images, an additional blue band threshold was applied, because this has been shown to improve the classification of ice in shadows (Paul et al. 2016). Typical threshold values for this parameter are between 500 and 2,000 (Paul et al. 2016). The threshold(s) were then used to reclassify the raster and convert it to vector polygons. Polygons were manually edited down on an as-needed basis using visual interpretation, and glacier area extents were calculated. Glacier extents from Landsat imaging instruments (ETM+ and OLI) can be directly compared, although comparisons between outlines derived from Landsat and MSI on Sentinel-2 are more challenging due to inaccuracies in co-registration (Storey et al. 2016) between the two platforms. Nevertheless, qualitative interpretations can still be made in some cases (with potential errors noted throughout).

We only report glacial extents derived from images that contain the least visible snow. However, because ice and snow have similar spectral signatures (Dozier 1989), it is possible that there is some snow contamination in the glacier boundary mapping. Changes in the glacier extent over time were noted, along with instances of possible/likely snow contamination, and are described in the Discussion section.

Surface temperature estimation

To estimate the surface temperature throughout the summit region, thermal band data from Landsat 7 and 8 can be used, along with several assumptions regarding surface materials and properties. The basic equation for calculating land surface temperature (LST) is

(1) $LST=\frac{k_2}{ln\left(1+\frac{k_1}{L}\right)},$(1)

where k₁ and k₂ are constants based on the thermal band (found in image metadata) and L is the surface radiance of the thermal band (USGS 2018). Surface radiance is calculated by first converting the raw digital numbers (Q_cal) to at-sensor radiance (L_λ):

(2) $L_{\lambda }=M_L\times Q_{cal}+A_L,$(2)

where M_L and A_L are the multiplicative and additive constants for the thermal band (found in image metadata; USGS 2018). A simple atmospheric correction is performed using the tool developed by Barsi et al. 2003, 2005), available at atmcorr.gsfc.nasa.gov  (Barsi et al. 2003). This tool identifies appropriate atmospheric profiles (of pressure, air temperature, and water vapor) from the National Centers for Environmental Prediction based on the image date, time, and location and applies a radiative transfer model to calculate the applicable values for transmissivity (τ) and upwelling and downwelling radiance (L_u and L_d, respectively). A calculation of the corrected, at-surface thermal radiance (L) is performed as follows (Barsi, Barker, and Schott 2003; Lo Vecchio et al. 2018):

(3) $L=\frac{L_{\lambda }-L_d\times \left(1-\epsilon \right)\times \tau -L_u}{\tau \times \epsilon }.$(3)

This atmospheric correction and surface radiance calculation method is expected to predict surface temperatures within ±2°C where surface emissivity (ε) is known. In this study, in situ measurements of emissivity are not available. We instead use simple landscape classification assignments, where typical emissivity values for each surface type are assumed. For the glacierized and snow-covered regions, extents are derived using the methods described in the previous subsection, and a constant emissivity of 0.98 is assigned (Rees 1990). There is some variability in the emissivity of snow and different types of ice (between 0.98 and 0.995), which introduces uncertainties of around ±0.4°C when calculating the brightness temperature (Rees 1993). For all other surface types, we use the Normalized Difference Vegetation Index (NDVI), which has a known relationship to surface emissivity (Caselles and Sobrino 1989). NDVI is often used in the classification of vegetation from remote sensing data sets (Tucker and Sellers 1986) and is calculated from the red-infrared and NIR reflectance using the following equation:

(4) $NDVI=\frac{NIR-Red}{NIR+Red}.$(4)

In order to more accurately obtain the NDVI values at the surface, it is important to perform the calculation in Equation (4) using atmospherically and topographically corrected surface reflectance rasters (Moreira et al. 2016). To obtain these for both the NIR and red channels, we first convert the raw digital number values to at-sensor reflectance using appropriate constants (USGS 2018) and perform a simple Dark Object Subtraction (DOS) correction using the GRASS GIS function i.landsat.toar. DOS is a purely image-based atmospheric method that infers a uniform correction where the darkest pixel in a scene is assumed to have a reflectance of zero. The GRASS GIS function i.topo.corr was then used to generate an illumination model based on (1) the scene’s topography from the ASTER GDEM2 tile over the region of interest and (2) local illumination conditions at the time of image acquisition (solar zenith and azimuth, found in image metadata). This illumination model and the GRASS GIS function i.topo.corr were then used to perform a topographic correction on the DOS-corrected surface reflectance data.

From the NDVI map, and using threshold values for the NDVI of fully vegetated (NDVI_v) and bare soil (NDVI_s), we can calculate a proportion of vegetation (PV), or the fraction of vegetation per pixel (Carlson and Ripley 1997):

(5) $PV={\left(\frac{NDVI-NDV{I}_s}{NDV{I}_v-NDV{I}_s}\right)}^2.$(5)

All pixels with NDVI values less than NDVI_s, which was set to 0.2, are assumed to be 100 percent bare soil (PV = 0), whereas all pixels with NDVI values greater than NDVI_v, set to 0.5, are assumed to be fully (PV = 1) vegetated (Sobrino, Jiménez-Muñoz, and Paolini 2004). Constant emissivity values are assumed for each of these surfaces using reasonable values for each material type. For mixed pixels, where the vegetation proportion is between 0.2 and 0.5, the emissivity (ε) is calculated as follows (Caselles and Sobrino 1989; Sobrino, Jiménez-Muñoz, and Paolini 2004):

(6) $\epsilon ={\epsilon }_v\times PV+{\epsilon }_s\times \left(1-PV\right),$(6)

where ε_v and ε_s are the emissivity values for pure vegetation and soil, respectively assigned as 0.973 and 0.996 (Avdan and Jovanovska 2016). It is important to note that the atmospheric correction performed on the thermal band data is expected to estimate surface temperatures within ±2°C, where emissivity values are known (Barsi et al. 2005). Because emissivity values for each land-type classification span a range and we do not have field measurements of this parameter, the values used are only approximations and may result in additional errors in surface temperature estimation. In the case of snow/ice and vegetation, with narrower ranges of typical emissivity, this error is on the order of around ±0.5°C. For bare soil, with much more variation in typical emissivity values and a stronger dependence on soil moisture (Mira et al. 2017), errors could be on the order of up to around 5°C. Because the output land surface temperature maps produce realistic results (see next section), we conclude that they provide a reasonable estimate of the actual surface temperature.

Results

Semi-automated glacier mapping

Pico de Orizaba, Glaciar Norte

At Glaciar Norte on Pico de Orizaba, reasonably snow-free and cloud-free images were obtained between 2001 and 2019. Most of these images were collected during the dry season (November through February), although a couple of August images were cloud-free and appeared to have very little snow cover (2013 and 2015) and were included to allow for better temporal coverage of glacier changes. Figure 2 shows the derived glacial boundaries, as well as the glacier areas at most of these time points. Due to gaps in the image data for the Landsat 7, areas could not be derived for any ETM+ images after May 2003. Although there are slight orthorectification differences between Landsat and Sentinel-2 (causing a small shift in the derived glacier outlines), we include an outline from the 23 February 2018 MSI image for comparison. A recently obtained outline from a Landsat 8 image obtained on 11 March 2019 was also included in Figure 2.

Figure 2. Boundaries of Glaciar Norte (Pico de Orizaba) derived from mostly snow-free images since 2001 using Landsat 7 (ETM+), Landsat 8 (OLI), and Sentinel-2 (MSI). Outlines are overlain on a false-color composite of the Landsat 8 image from March 2019. Note that there is a slight shift in the georeferencing between Landsat and Sentinel-2, so the 2018 outline is somewhat offset. White boxes highlight the regions of greatest change since 2001.

The extent of Glaciar Norte has generally decreased over time, where an areal loss of approximately 30 percent was observed between 2001 and 2019. The largest region of retreat was located on the northeastern edge of the glacier (Figure 2a). Between 2001 and 2017, a total loss of 0.18 km² was observed in this region of the summit. In the decade between 2007 and 2017, approximately 100 m of linear retreat was observed in this location, accounting for an area loss of approximately 0.066 km².

We also find that glacier changes have occurred to the west of the summit crater, where there was a larger extent observed in 2001 as compared to 2019 (Figure 2b). There are some small increases in apparent glacier extent in this region (e.g., in February 2017), which are likely due to classification errors where snow adjacent to the glacier is indistinguishable from ice. However, the larger overall change from 2001 to 2019, a difference in area of 0.064 km², may be related to the final retreat of the occidental/suroccidental glacier tongues.

Pico de Orizaba, Glaciar Oriental

Because Glaciar Oriental is rather small, it is not always detectable in Landsat images due to a range of challenges, including obvious snow adjacent to the glacier, small clouds that obscure observations, difficulties identifying boundaries at the resolution of orbital imagery, and, in the case of ETM+ images, losing all or a significant part of the glacier to data gaps. We use the extents in the Randolph Glacier Inventory to help identify the glacier location (RGI Consortium 2017). It should be noted that though the latest inventory released is from 2017, the outlines pertaining to this particular region have not been updated and still reflect the approximate extent in 1999. Moreover, the Randolph Glacier Inventory (RGI) outlines are rather generalized and do not do a good job of characterizing the more subtle surface variations that are apparent when analyzing these features in detail. Due to the size of this glacier, it was possible to identify Landsat 7 ETM+ images where data gaps do not interfere with mapping the extent, so we exclusively use images that yield complete outlines and valid areal extents for this analysis.

Figure 3 shows the outline of Glaciar Oriental on Pico de Orizaba. We find that the area of Glaciar Oriental as of 2001 was 0.116 km², similar in extent to the outline found in the RGI (black line in Figure 3; RGI Consortium 2017) and roughly consistent with the observations of Cortés-Ramos and Delgado-Granados (2015), who found that the area of Glaciar Oriental did not change significantly from the observations of Lorenzo (1964). However, over the last several decades, we have recorded additional losses to the areal extent of Glaciar Oriental: by 2018, the areal extent measured was 50 percent smaller than that which was observed in 2001, where the bulk of the ice loss occurred from 2013 to 2018. Furthermore, the most recent outlines indicate that ice is currently only present in the southern part of the area that was glacierized in 2001. There is a north-facing cliff (around tens of meters high as estimated from the ASTER GDEM2) that bisects this region; ice only remains above this cliff, on the higher-elevation southern surface and above ~5,315 m in altitude.

Figure 3. Changes in Glaciar Oriental on Pico de Orizaba between 2001 and 2018, compared to the outline in the RGI from 2009. Most outlines are derived from images from ETM+ on Landsat 7 with the exception of 23 February 2018, from MSI on Sentinel-2. The background image is a color composite of the image from February 2018.

Iztaccíhuatl (El Pecho and La Panza)

Figure 4 shows the glacier outlines at the summit of Iztaccíhuatl from 2003 to 2018 using the best cloud- and snow-free images from Landsat 7 and 8 and Sentinel-2 and includes the areal extents of both the El Pecho (northern) and La Panza (southern) glacial regions. The Sentinel-2 outline from 2017 is slightly shifted south and east of the Landsat-derived outlines, which accounts for some of the apparent deviations. Much of the glacier retreat observed at Iztaccíhuatl has been on the northern/western part of El Pecho (Figure 4a), constituting ~0.024 km² of areal ice/snow loss between 2003 and 2014, with an additional loss of ~0.017 km² distributed around the entire glacier between 2014 and 2018 (for a total areal loss of ~0.034 km²). At La Panza, there has also been some ice loss in the northwestern (Figure 4b) and northeastern (Figure 4c) parts of the glacier, as well as a gap in the middle of the glacier that account for a total area difference of 0.078 km² between 2003 and 2018. It is also possible that some ice has been misclassified due to debris cover or other darkening at La Panza, particularly in the 2014 image (which appears to have the lowest areal extent in Figure 4), because although this region looks somewhat darker in 2014, it once again has a band ratio consistent with ice in the subsequent images from MSI in 2017 and OLI in 2018. It is also notable that there is a small ice detection south of La Panza derived from the earliest images in Figure 4, consistent with the location of the Atzintli Glacier, which is no longer present (Schneider et al. 2008). In 2018, we observe that 80 percent of the 2001 El Pecho area remains, whereas La Panza lost approximately 52.5 percent of its areal extent since 2003.

Figure 4. Glacier outlines at the summit of Iztaccíhuatl derived from Landsat images collected between 2003 and 2018 and a Sentinel-2 image from 2017. White boxes highlight some of the observed areas of retreat. Outlines are overlain on the OLI image from January 2018. *Two different ETM+ images from 14 November and 30 November 2005 are included to obtain an improved outline across data gaps. However, we do not use this merged outline to obtain the areal extent at this time point.

Surface temperature estimation

Surface temperature maps were derived using the methods described in the previous section for all noncloudy Landsat 7 and 8 images over both summit areas, including images where the summit is snow covered. Note that all surface temperatures are rough estimates, due to both the approximation of surface properties like emissivity and the intrinsic errors in the atmospheric correction of thermal data (Barsi et al. 2005). In addition, all Landsat images are collected at a similar time of day (around 5 p.m. GMT/11 a.m. local time), and many images collected during the wet season include clouds over the summit. Therefore, all temperature estimates described here only offer rough approximations and do not represent all possible times/variations.

Pico de Orizaba

Figures 5a–5f show the land surface temperature maps for a selection of cloud-free January images collected since 2002. Because these images are from the peak of the dry season at Pico de Orizaba, they tend to have the least visible snow and cloud contamination. The temperatures over Glaciar Norte itself are generally well below freezing, as is the surrounding soil at the summit (although it is important to note that, due to the inherent variability in soil emissivity and the lack of ground observations, the temperature of bare soil is subject to higher errors than that of ice or vegetation). We also show some of the LST maps collected during the wet season (Figures 5g–5l), although there were fewer usable (cloud-free) images available during these months. Ice surface temperatures during the wet season are 5°C–10°C warmer than those seen during the dry season but are generally well below the freezing point of water. Snow-covered soil, abundant in some of the wet season images (most notably in July 2014 [Figure 5j] and April 2016 [Figure 5k]) also tend to be below freezing for tens to hundreds of meters downslope from the glacier. However, the image from 27 May 2018 indicates that the ice surface temperatures on the northern part of the glacier are at or near the melting point of water. We also note that the region northeast of the glacier generally tends to have the warmest surface temperatures over the summit area, reaching about 15°C during several months.

Figure 5. Land surface temperature maps (for temperatures less than 4°C) at Pico de Orizaba, derived from a selection of Landsat images (a)–(f) during the dry season in January and (g)–(l) during the wet season (April–October) on Pico de Orizaba. Glacier outlines shown are from Figure 2, where the outline from the temporally closest snow-free image was used for each time point: (a) 28 January 2001; (b) 7 January 2005; (c), (d), (g), (h) 19 February 2009; (e), (i) 17 August 2013; (j) 18 December 2014; (f), (l) 11 January 2018; and (k) 29 November 2016. At all January time points shown, the ice temperatures are colder than −8°C. In many cases, the surface temperatures adjacent to the glacier are also below freezing. Many of the wet season images contain extensive snow cover, with the exception of (g) and (l).

In order to better understand surface temperature variations throughout the year, we collected the average temperature profile for each month using as many cloud-free images as possible across all available years from all Landsat 7 and 8 LST maps. Profiles were taken crosswise over Glaciar Norte to span the region of retreat shown in Figure 2, as well as down the center of the glacier, extending northwest from the summit crater. Each monthly profile is the result of averaging over two to twelve observations depending on image availability and data gaps (ETM+), where most dry season months were composed of more than seven profiles and most wet season months were averaged over less than five profiles (and months with only one image available were excluded). Resulting average monthly temperature profiles are shown in Figures 6 and 7, respectively. For both of these plots, regions with a white background are on the glacier itself (based on the glacier extent mapped in 2018) and the regions that are shaded gray are bare ground adjacent to the glacier. Each figure also includes the elevation derived from the ASTER GDEM2 along the profile. We find that the surface temperature of the glacier itself is generally below freezing along both profiles, with some seasonal variability. In the wet season (approximately April through October), ice temperatures are warmer by about 5°C–10°C compared to the peak of the winter dry season (November through February). In addition, we observe average monthly soil temperatures that are 15°C–20°C warmer on the northeastern side of profile A to A′ compared to the southwestern part of the profile (Figure 6).

Figure 6. (Top) Monthly average surface temperatures across Pico de Orizaba’s Glaciar Norte and the (bottom) associated elevation profile. Regions shaded in gray are adjacent to the glacier (using the boundaries mapped in 2018). Inset map shows the location of the profile (superposed on a 2018 OLI image), and the dashed cyan line demarcates the freezing point of water. The standard error for each profile is within a few degrees (see Supplemental Figure 1 for error bands).

Figure 7. (Top) Monthly average surface temperature profiles of Glaciar Norte on Pico de Orizaba, along a profile that extends northwest from the summit crater, with the (bottom) associated elevation profile. The gray-shaded region is the soil beyond the edge of the glacier, and the dashed horizontal cyan line represents the freezing temperature of water. The standard error for each profile is generally plus or minus a few degrees (see Supplemental Figure 2 for error bands).

At Glaciar Oriental, we observe that surface temperatures over parts of the glacierized region remain below freezing year-round. However, both the northern and southern sides of the profile shown in Figure 8 peak at reasonably high temperatures, particularly during the spring and summer months (March–July). As discussed previously, we found that in 2001, the glacier area was roughly equivalent to the extent in the RGI (black outline in inset of Figure 8; RGI Consortium 2017). However, by 2013, the only region containing ice was the southern side of the small, north-facing cliff that bisects the RGI extent, as indicated by the shaded regions in Figure 8 (where the vertical dashed line represents the cliff location). The cliff height as measured using the ASTER GDEM2 is on the order of 10–30 m, fairly close to the vertical resolution and dwarfed by the overall elevation changes across profile C to C′ (bottom plot in Figure 8). Surface temperatures on the northern side of the cliff, with overall average monthly temperatures that peak at 5°C–10°C during the wet season at the lowest elevations, still tend to drop below freezing at the base of the cliff, likely due to shadowing effects. However, in images collected after 2013, no glacier ice is identified in this region, even when using an additional threshold in the blue band to better characterize ice in shadow (as recommended by Paul et al. 2016). The southern side of Glaciar Oriental also reaches temperatures that appear to be above the freezing point of water during several dry-season months, although this region remained glacierized in the images we analyzed. Note that, because there are fewer good images over smaller Glaciar Oriental, the standard error of these average profiles is fairly high (up to ±5°C; Supplemental Figure 3) in addition to the errors inherent to this method of estimating surface temperatures. Furthermore, it is probable that at the edge of the glacier extent, pixels are a mix of ice and bare soil. These uncertainties likely account for the apparent higher (above melting) temperatures at the southernmost edge of profile C to C′, where the ice itself is probably closer to 0°C. It is possible that this region of the glacier is thinning due to these warmer temperatures and may be at risk for additional areal retreat in the near future.

Figure 8. Average monthly surface temperature profiles across Glaciar Oriental on Pico de Orizaba, spanning the outline from the RGI (black outline in inset). Gray-shaded regions do not contain glaciers (each boundary is labeled with the year at which the glacier ended at that point along the profile). The dashed vertical gray line represents the ridge that bisects the region, where the glacier is presently limited in extent to south of the ridge. The bottom plot shows the elevation along the profile from the ASTER GDEM2 (see Supplemental Figure 3 for error bands).

Iztaccíhuatl

In general, surface temperatures at the subpeaks of Iztaccíhuatl are warmer than those at Pico de Orizaba. We similarly average over all available LST profiles by month to obtain monthly averages using all available cloud-free Landsat images between 1999 and 2018, where there are more total usable profiles between November and February (more than seven) compared to those in March and May (less than four). There are two or less images available for April and June through October, so these months are excluded. We see a trend similar to that of Pico de Orizaba, where the average temperatures are generally higher in the spring months than during the winter. In addition, the temperatures of the bare ground surrounding the glaciers is significantly higher than at Pico de Orizaba (averaging >15°C year-round and as high as 25°C in March and April). The only non-icy region that tends to have surface temperatures below 10°C is the northernmost part of profile D to D′ (Figure 9), which is relatively higher in elevation and experiences some topographic shadowing. The surface temperatures of the glaciers are markedly lower than the surrounding regions, consistently at or just below the freezing point of water throughout the year. Apparent ice temperatures above 0°C at the glacier edges are likely attributable to errors in temperature calculation and the influenced of mixed materials in the pixels at the ice boundaries. Winter (dry season) surface temperatures also appear to be slightly higher (by a few degrees) at the lower-elevation La Panza summit compared to the El Pecho summit. The La Panza summit has experienced more extensive ice loss over the last decade (see section “Iztaccíhuatl (El Pecho and La Panza)”), although the differences in ice surface temperatures over each glacier are similar in magnitude to the error bars (Supplemental Figure 4) and are therefore not considered significant. It is plausible that these glaciers are experiencing thinning in addition to the areal losses observed, and they may also be at risk of continued areal loss due to the fact that their surface temperatures are very close to the melting point of water throughout the year.

Figure 9. Surface temperature profiles, averaged by month since 2002, spanning the two remaining glaciers on Iztaccíhuatl. Gray-shaded regions do not contain glaciers, whereas the regions with white backgrounds cover (from left to right) La Panza and El Pecho. The bottom plot shows the variation in elevation along the track from the ASTER GDEM2 (for error bands, see Supplemental Figure 4).

Discussion

At Pico de Orizaba, we find that most of the areal ice loss is concentrated in the northern/eastern side of the Glaciar Norte ice cap, as well as just west of the summit crater. In 2001, the downslope remnants of Glaciar de Jamapa observed by Palacios and Vázquez-Selem (1996) were no longer apparent, consistent with their estimate that the two ice patches would not survive more than 5–10 years. We also observed changes in Glaciar Oriental, on the eastern part of the Pico de Orizaba summit.

Unsurprisingly, the location with the most consistent retreat over the last two decades (the northeast side of Glaciar Norte) tended to have higher surface temperatures as estimated from the Landsat thermal band imagery. Figure 10 shows four selected glacier outlines (from 2001, 2007, 2013, and 2018), as well as profiles of surface temperatures across the glacier derived from the associated thermal band data (plotted above the profile’s elevation from the ASTER GDEM2). We find that the surface temperatures within this region of concentrated retreat are significantly warmer than the soil on the opposite side of the glacier (by about 15°C–20°C). The southwestern side of the glacier is ~30 m higher in elevation and is topographically shadowed, allowing the surface temperatures adjacent to the glacier to persist at temperature below freezing for several hundred meters in all images observed. To what extent regional air pollution is affecting the retreat of Glaciar Norte is unknown. Atmospheric sampling of nitrogen dioxide (NO₂), carbon monoxide (CO), and sulfur dioxide (SO₂) performed on the northwest of Pico de Orizaba (at 3,300 m.a.s.l.) shows a daily pattern that suggests the influence of pollution generated by anthropogenic processes from the Mexico City basin (Márquez et al. 2005). The transport of air pollutants outside the Mexico City basin exhibit three wind circulation patterns (de Foy et al. 2005). For one in particular, referred to as O3-South (because it corresponds to days when the ozone is highest in the south of Mexico City), the plume leaves the basin mainly to the south and fans out and impacts both the Pacific coast to the west and the Gulf coast to the east. The air pollutants that are transported to the east show some looping around the Pico de Orizaba volcano due east, with northward transport along the flanks of the Sierra Madre Oriental (de Foy et al. 2006). This atmospheric transport of air pollutants from the Mexico City basin is correlated with the retreat of the glacier, and it is possible that pollution plays some role in glacier loss. However, we cannot definitively determine causation with the present data. The altitudinal retreat of the glacier in this northeastern summit area (as measured using the ASTER GDEM2) between 2001 and 2018 was approximately 40 m, from 5,020 to 5,058 m.a.s.l. Over the Pico de Orizaba summit region, approximately 4 percent of ASTER GDEM2 pixels had scene counts of less than ten, and aproximately 85 percent of pixels had scene counts greater than fifteen, so we expect that the RMSEs are reasonably low (≤9 m). Because the measured elevation changes are on the order of around four times larger than the typical RSME, this difference is well above the expected limit of detection using the ASTER GDEM2 and provides a reasonable estimate of altitudinal retreat.

Figure 10. (Top) Retreat of Glaciar Norte at Pico de Orizaba and (bottom) associated surface temperature and elevation profiles from A to A′. Glacier outlines were obtained from ETM+ images collected around the same time of year (dashed parts of the outlines reflect gaps in image coverage) and are overlain on a false-color composite from January 2018. The dashed horizontal cyan line in the temperature profile plot highlights the freezing point of water, and the vertical dashed lines show the retreat of the glacier edge at each time point. All yellow lines correspond to 2018; all orange lines correspond to 2013; all purple lines correspond to 2007; and all black lines correspond to 2001. The bottom part of the profile plot shows the elevation along the profile.

Smaller Glaciar Oriental on Pico de Orizaba was also found to have experienced retreat in recent decades, in contrast to previous observations. Between 1958 and 2007, the area of Glaciar Oriental did not vary significantly compared to the other glaciers on Pico de Orizaba (Palacios and Vázquez-Selem 1996; Cortés-Ramos and Delgado-Granados 2012). The glacier extent that we measured in 2001 was roughly consistent with previous observations (Cortés-Ramos and Delgado-Granados 2015). However, areal ice loss appeared more rapid recently, especially in the last 5 years. By 2018, we observe that Glaciar Oriental lost approximately 50 percent of its 2001 area, where most of the remaining ice mass was limited to the higher-elevation, southern side of a cliff that extends ENE from the summit crater. It is possible that Glaciar Oriental has experienced thinning/surface lowering between 1958 and 2007, a phenomenon observed in many glacier systems throughout the world (e.g., Berthier et al. 2004; Paul and Haeberli 2008; Malecki, 2017; Dehecq et al. 2019). However, with the lack of elevation data with sufficiently high resolution, we cannot confirm this hypothesis. We further find that the surface temperatures are higher both in the region where retreat has occurred and at the southern margin of the remaining glacier area, indicating that it is possible that this region is at risk for future ice losses despite its relatively higher elevation. In 2001, we find that the lowest elevation glacial extent at Glaciar Oriental was 5,228 m, north of the cliff. By 2018, the glacier only extends to 5,293 m, at the southeasternmost edge of the glacier.

At Iztaccíhuatl, it was previously estimated that the glaciers may not persist beyond 2020 (Schneider et al. 2008). As of 2018, these glaciers were somewhat smaller in extent but not necessarily consistent with total loss within the next 2 years. Areal losses of 0.032 and 0.077 km² were observed at El Pecho and La Panza between 2003 and 2018, where lower-elevation La Panza has suffered the most extensive loss. Compared to the rate of areal shrinking found between 1958 and 2007 (Schneider et al. 2008), the rate of loss that we observe over the last 15 years appears to be somewhat lower. However, it is important to note that these observations only account for area changes, not volume loss. Because volume change requires high-resolution topography, it is impossible to characterize such variations with existing orbital data and would likely require repeated field campaigns in the future to fully characterize the ice loss in 3D space. The surface temperatures at Iztaccíhuatl are largely found to be greater than those found at Pico de Orizaba, over both the glaciers themselves and especially over the surrounding bare soil at the summit. This is likely due in large part to lower elevation, where the peak of Iztaccíhuatl is approximately 500 m lower than Pico de Orizaba. Due to the warmer ice temperatures that approach the melting point of water, it is plausible that volume loss is actively taking place, at least during part of the year.

Over the summit region of Iztaccíhuatl, 68 percent of pixels in the ASTER GDEM2 had scene counts less than fifteen, suggesting higher overall RMSEs. We are therefore less confident in the accuracy of elevation changes at Iztaccíhuatl (for pixels with scene counts less than fifteen, typical RMSE ranges between 10 and 16 m). The observed change in elevation of La Panza between 2003 and 2018 was ~70 m, along the western side of the glacier. At El Pecho, over the same time span, approximately 45 m of altitudinal retreat was noted on the north side of the glacier (from 5,192 m to 5,238 m).

Conclusions

It is clear that glacier retreat is taking place at both Pico de Orizaba and Iztaccíhuatl. Semi-automated mapping using orbital multispectral imagery reveals two primary concentrations of retreat on Pico de Orizaba’s Glaciar Norte ice cap: the northeastern glacier face, which has retreated approximately 100 m in horizontal distance (and ~40 m in elevation) between 2007 and 2018, and west of the summit crater, where the occidental and suroccidental glaciers once extended down to 4,930 m and presently extends to only 5,310 m. We also observe retreat of Glaciar Oriental, on the eastern summit of Pico de Orizaba, which had shown little area change between 1958 and 2007 but lost approximately 50 percent of its areal extent between 2003 and 2018. At the lower-elevation summits of Iztaccíhuatl, we find that glaciers at El Pecho and La Panza have retreated by approximately 20 percent and 50 percent, respectively, between 2003 and 2018, although it is plausible that these glaciers will persist beyond 2020 (the time frame for total glacier loss estimated by Schneider et al. 2008). Fairly high temperatures over ice (at/near 0°C) and bare soil (up to ~15°C adjacent to glaciers) are noted, especially near/over parts of Glaciar Oriental, El Pecho, and La Panza, suggesting that these regions may be particularly susceptible to future ice loss. Continued monitoring of these sites is important, because these glaciers are part of a dwindling population of equatorial glaciers and are the only permanent ice masses at this latitude across the Northern Hemisphere. Among other things, the retreat of such glaciers could cause the loss of microbial diversity that is endemic to these environments (Callegan et al. 2008). Future in situ observations would be helpful to better quantify the factors that are impacting modern ice losses in this unique environment and to identify higher-resolution changes in ice volume over time.

Acknowledgments

All Landsat 7, Landsat 8, and Sentinel-2 imagery used in this study are courtesy of the U.S. Geological Survey, obtained through the GloVis tool (glovis.usgs.gov). ASTER GDEM2 is a product of METI and NASA and was acquired from the Land Processes Distributed Active Archive Center (LP-DAAC) using the Earthdata tool (earthdata.nasa.gov).

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

This work was supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Universities Space Research Association under contract with NASA.