Disappearance of the last tropical glaciers in the Western Pacific Warm Pool (Papua, Indonesia) appears imminent

Ice Thickness Reduction.

Changes in ice thickness since June 2010 were determined by measuring the accumulation stake during subsequent visits to the drill site, which revealed progressive exposure of PVC sections due to surface ice loss (SI Appendix, Fig. S1 B–D). During the first site visit on November 3, 2015, 2 pipe segments and 1.49 m of rope were observed lying on the surface, indicating that the ice had thinned a total of ∼5 m at a rate of ∼1.05 m/y since 2010 (Fig. 2A and SI Appendix, Fig. S1B). During a monitoring visit on May 31, 2016, 4 pipe segments were exposed at the surface (Fig. 2A and SI Appendix, Fig. S1C), indicating an additional ice thinning of ∼4.26 m within just over 6 mo. This was likely due to the effects of the very strong 2015–2016 El Niño. The most recent measurement on November 23, 2016, revealed an additional 1.43 m of ice thinning since May 2016, confirming that, although the annual rate of thinning had slowed, it remained 2.7 times greater than the average rate from 2010 to 2015 (Fig. 2A and SI Appendix, Fig. S1D). This thinning accompanied the progressive contraction of the glacier surface area resulting in the separation of the ENF into 2 ice masses (Fig. 2 B–E). These measurements illustrate the overall rapid rate of ice thinning as well as the particularly strong impact of the temperature increase and precipitation decrease associated with the 2015–2016 El Niño.

Ice Coverage Loss.

The 2 remaining large glaciers near Puncak Jaya (Fig. 3A) experienced a dramatic loss of ice coverage from 2002 to 2018 that intensified during the 2015–2016 El Niño (Fig. 3 B and C). The new data from 2015, 2016, and 2018 are integrated with those from previous observations (17) to show that the total glacier area near Puncak Jaya has decreased linearly from ∼1850 to 2018 (Fig. 3D). From 2002 to 2015, the ice cover diminished by 1.45 km², equivalent to a loss rate of 0.11 km² per y. During the 2015–2016 El Niño, an additional 0.11 km² of ice cover was lost from 2015 to 2016. In March 2018, the ice cover was only 0.458 ± 0.036 km² or a loss of 0.09 km² since 2016. During the period from 2002 to 2018, the glacier retreated at a rapid but steady pace (Fig. 3D, Inset).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Changes of total ice area and impact of the 2015–2016 El Niño on Papua glaciers. (A) Landsat-5 satellite image of the glaciated area on October 9, 2009. The ENF and the Carstensz Glacier are outlined in black boxes. (B) Planet’s satellite image of ENF taken on May 16, 2016, showing outlines of ice extent on June 11, 2002 (as in ref. 6); March 16, 2015; May 16, 2016; and March 28, 2018. (C) Planet’s satellite image of Carstensz Glacier taken on May 16, 2016, showing outlines of ice extent on same dates as in B. (D) The changes of total glacier area near Puncak Jaya from ∼1850 to 2018. The Inset highlights the ice area changes from 2002 to 2018. Color-coded triangles depict the coverage outlined in B and C; the gray circles and uncertainty bars are derived from other studies (17). Additional information on the satellite images in B and C are provided in SI Appendix, Table S1.

The closest available instrumental records of temperature and precipitation are located 50 to 450 m below the minimum glacierized elevation (Grasberg [GRS], Dispatch Tower [DISP], and Alpine [ALP] stations, Fig. 1B) and cover the period 1997–2016 (SI Appendix, Table S2). Although these records have a number of data gaps, they indicate a near-zero probability that average daily temperatures have fallen below the freezing point since 1997 (SI Appendix, Fig. S2). This assumption is supported by calculating freezing levels from the average temperature at the nearest station (ALP) using 2 realistic lapse rates, i.e., the mean lapse rate between GRS-DISP and GRS-ALP (4.6 °C/km) (SI Appendix, Fig. S3), and the environmental lapse rate (6.5 °C/km), to identify the remaining ice surface that lies above the freezing line. Coupling these freezing levels with 3 different hypsometric curves of the glacierized area indicate that between 0 and 10% of the remaining ice surface lies above the freezing line (SI Appendix, Fig. S4). Thus, even with steep lapse rates, this hypsometric distribution implies that prevailing temperatures are generally above freezing over the glaciers, which promotes glacier melting. Furthermore, accumulation during these warm conditions is restricted to a few hours during the night and only over the highest elevations where frozen precipitation covers a small fraction of the glacierized area (SI Appendix, Fig. S5).

Glacier Mass Balance Model and Shrinkage Projections.

The rapid rate of area loss and dramatic ∼5.4-fold increase in the rate of thinning in 2015–2016 present a scenario of unstoppable glacier demise. Depending on the interannual variability of temperature and precipitation over the next few years, the remaining ice fields of Papua, Indonesia, are very likely to disappear during the next decade. To provide a quantitative basis for this scenario, a glacier mass balance model was implemented by combining the available instrumental climatic records and dynamical/statistical downscaling of future climate model projections for the region under different climate scenarios (SI Appendix, Figs. S6 and S7 and Tables S3 and S4). We performed simulations of glacier volumetric change for the period 2017–2030 using a spatially distributed mass balance model. The model computes the difference between ablation and accumulation at daily time steps and transforms that change into volume relative to initial volume. Surface area reported in the Randolph Glacier Inventory (RGI) (21) is utilized to compute an initial volume of 0.13 km³ according to an area–volume scaling (22, 23). Ablation is calculated from daily positive temperatures while accumulation utilizes temperature thresholds to determine the amount of solid precipitation falling over the glacierized area. The model was evaluated and calibrated to reproduce the surface thinning based on the available stake measurements by using optimal combinations of lapse rates derived from the closest instrumental records of temperature and precipitation. Simulation of volumetric changes was driven by climate model projections over the Southeast Asia region (Southeast Asia Regional Climate Downscaling [SEACLID]/Coordinated Regional Climate Downscaling Experiment–Southeast Asia [CORDEX-SEA] project; http://www.ukm.edu.my/seaclid-cordex). The model predicts that under practically all scenarios glacier shrinkage will continue unabated, leading to total ice loss no later than 2026 (SI Appendix, Fig. S8).

Ice Core Records.

The ice cores recovered in 2010 from the glaciers near Puncak Jaya, the first to bedrock, provide both an absolute constraint on ice thickness and proxy evidence to examine the role of atmospheric and sea surface warming and other potential mechanisms driving glacier loss in the region. To further understand the mechanisms driving the recently observed retreat of these glaciers, analyses were conducted on the Papua ice cores to reconstruct the past tropical Pacific climate variability. The 2 longest cores drilled on the ENF (D1, 32.13 m, and D1B, 31.25 m long; located ∼5 m apart) were analyzed for δ¹⁸O, δD, and concentrations of insoluble dust and major ions (17). The D2 (26.19 m) core was not analyzed as the drilling was terminated when vertical ice layers, suggesting contorted ice flow, were encountered. Basal melting was confirmed by the 0 °C temperature at the bottom of the D1 and D1B boreholes and by meltwater encountered 80 cm above the bottom of the D2 borehole. Reproducibility of the stable isotopes in both cores is illustrated (Fig. 4 A and B and SI Appendix, Table S5). The presence of the 1964 tritium (³H) peak is detectable in a record of Papua precipitation at the Global Network of Isotopes in Precipitation [GNIP (24)] Jayapura station from 1957 to 1991 (SI Appendix, Fig. S9 A and B). This tritium peak is also recorded at 23.4-m depth in the D1 core, although its value (2.98 ± 0.42 TU) (Fig. 5A) is reduced due to radioactive decay over the 47 intervening years. This serves as the sole chronological marker for the timescale reconstruction. Modest precipitation seasonality in the Papua highlands (13) and post-depositional processes minimize seasonal aerosol variations in the cores (Fig. 4 C–E). Dating of the ice core records was challenging due to 1) this lack of seasonal variability, 2) the presence of only one chronological marker, and 3) melting of the glacier at both the surface and the base. However, construction of post-1964 chronologies for both cores was attempted by δ¹⁸O reference matching with the NINO3 extended reconstructed sea surface temperatures [ERSST (25)]. The rationale for the matching is based on a positive, albeit weak, correlation between rainfall δ¹⁸O at GNIP stations in Papua and NINO3 ERSST (SI Appendix, Figs. S10 and S11 and Table S6). The resulting timescales demonstrate that the ENF cores encompass at most 47 y (1964–2010 CE) in the top 23.4 m (Fig. 5B). The correlations between the 2 ice cores’ δ¹⁸O values and the NINO3 sea surface temperatures (SSTs) are statistically robust with r = 0.58 (P < 0.001) for the D1 core and r = 0.52 (P < 0.001) for the D1B core. The correlations are just slightly stronger for the 1964–1999 period with r = 0.61 (P < 0.001) for the D1 core and r = 0.53 (P < 0.001) for the D1B core. However, ablation at the surface of the glacier, which was observed during the 2010 drilling operation and subsequently documented by stake measurements (Fig. 2), may have truncated the top of the time series. The very high concentrations of soluble and insoluble aerosols in the upper ∼4 m of the ice cores (Fig. 4 C–E) illustrate the extent of this ablation, which may have accelerated in the early 21st century (∼2005) (Fig. 5 C–E). Therefore, this climate record since 1964 should not be regarded as precise at an annual resolution. Since melting is also occurring at the base of the glacier, the bottom ages of the D1 and D1B ice cores cannot be fixed, and thus timescale reconstruction for the cores was not attempted below the depth of the ³H signal.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Ice core records. The (A) δ¹⁸O and (B) δD in the D1 (black) and D1B (red) cores show high reproducibility and positive trends toward the surface. Comparison of (C) mineral dust, (D) chloride, and (E) potassium concentrations among the ENF cores shows reproducibility and increased concentrations toward the top as the stable isotope variability decreases and the trend of isotopic enrichment increases (gray shading).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Annual ice core records since 1964. (A) The δ¹⁸O and ³H data from 68 samples in D1 core (red closed circles with error bars). The 1964 horizon is marked by the ³H peak (2.98 ± 0.42 TU) at 23.4-m depth. (B) The D1 and D1B δ¹⁸O time series are compared, and 12 points in D1 (green closed circles) are paired with corresponding points in 13-mo running means of NINO3 ERSST (24) after assigning the 1964 ³H peak at 23.4 m and assuming the top layer is May 2010. The correlations between the 2 ice cores’ δ¹⁸O values and the NINO3 SSTs are shown. (C) The logged values of annual averages of dust concentration (size of >0.63 μm) in the D1 and D1B cores post-1964. (D) As in C, but for chloride concentrations. (E) As in C, but for potassium concentrations.

An independent check on the timescale is possible by comparing the core D1 tritium data (Fig. 6A) with the record of annual tritium in precipitation from the GNIP stations in Jayapura, Indonesia, and Madang, Papua New Guinea (SI Appendix, Fig. S9A), from 1957 to 1991 (Fig. 6B). Unfortunately, tritium data from these stations do not exist after 1991. However, the comparison of these records over 35 y is comparable with the ice core δ¹⁸O matching with the NINO3 SST time series.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Comparison of core D1 and GNIP station ³H records. (A) The ³H record from core D1 is shown in depth next to (B) the annual ³H in precipitation collected from GNIP stations in Jayapura, Indonesia (blue curve), and Madang, Papua New Guinea (red curve), from 1957 to 1991 (Bottom). The green curve depicts the average ³H values from the 2 stations where the records overlap (1974–1982). The thin lines connect similar features in the ice core and precipitation curves. Note that the ³H concentration in precipitation collected at the stations has been calculated using the decay half-life of 12.3 y in order to show that the value of the peak in the sample collected in 1964 is compatible with that in the ice core record.

In the Tropics, precipitation δ¹⁸O is negatively correlated with precipitation amount at lower altitudes and in coastal regions through the amount effect (26, 27). However, this relationship is much more complex, since stable isotopes in tropical precipitation are affected by source conditions such as SSTs, along with large-scale atmospheric circulation, transport pathways, cloud-top temperature and cloud-top pressure in convective clouds, and prevailing low-latitude upper tropospheric temperatures (10, 13, 28⇓⇓⇓⇓–33). The precipitation type also appears to influence δ¹⁸O, i.e., a higher proportion of convective versus stratiform precipitation results in isotopic enrichment (34). Moreover, convective precipitation is more responsive to increasing temperature than stratiform precipitation (35). The influences of coupled ocean–atmosphere processes on stable isotopes are particularly pronounced during ENSO events, which is evident in the correlation between SSTs in the equatorial Pacific and δ¹⁸O in rainfall at Jayapura (SI Appendix, Fig. S10), which facilitated the development of the ENF ice core timescales.

Since the middle of the 20th century, the increasing tropical Pacific sea surface and upper level temperatures have been augmented by the effects of El Niño (36). The ENF ice core δ¹⁸O records show increasing trends (SI Appendix, Fig. S12A), with an average slope of 0.018 ± 0.008‰ per year (P < 0.03), which reflects the trends in global and regional surface and upper level (550 mb) temperature anomalies (SI Appendix, Fig. S12 B and C). However, from the mid-20th century to 2010, seasonal trends in precipitation anomalies were constant (SI Appendix, Fig. S12D). This also implies that the relationship between δ¹⁸O and the amount effect, which is obvious on a seasonal basis in the Tropics, is indirect on intraannual timescales. Precipitation amount and temperature are components of the complex atmospheric and oceanic processes that control the stable isotopes in the high-altitude precipitation and hence in the glacier ice.

Although the ENF ice core records contain less than a century of climate history, they serve as the only ice core-derived, high-altitude climate reconstruction from the WPWP. These ice core records not only demonstrate the deteriorating condition of these few surviving glaciers, but also provide information on the changing conditions of the atmospheric and oceanic conditions that will drive the loss of the only tropical ice fields between the Himalayas and the Andes over the next few years, especially if the loss is augmented by strong El Niños.

The glaciers near Puncak Jaya are remnants of glaciers that have existed for ∼5,000 y (37, 38), and at the current rate of ice loss these glaciers will disappear in a very short time. The precise ice coverage over several years provides a reasonable determination of the variations in the rate of ice loss between El Niño and non-El Niño years.

The climate records from the ENF ice cores establish that the current glaciers have undergone significant mass loss, are melting from the surface and at the base, and have been doing so for many years, and that the rate of melting and hence the destruction of information regarding interannual climate variability is increasing. There is a strong probability that the ENF now contains fewer years than it did in 2010 as time (ice) is being removed by ablation from the bottom and the top. Along with the Furtwängler Glacier on the top of Kilimanjaro (5), these cores illustrate what happens to a climate record contained in a glacier as it melts away.