# CLIMATE IMPACTS ON GLACIERS AND BIOSPHERE IN FUEGO-PATAGONIA

EDITED BY : Christoph Schneider, Matthias Holger Braun, Marius Schaefer, Jussi Grießinger and Gino Casassa

PUBLISHED IN : Frontiers in Earth Science and Frontiers in Ecology and Evolution

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-758-4 DOI 10.3389/978-2-88963-758-4

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# CLIMATE IMPACTS ON GLACIERS AND BIOSPHERE IN FUEGO-PATAGONIA

Topic Editors:

Christoph Schneider, Humboldt-Universität zu Berlin, Germany Matthias Holger Braun, Friedrich-Alexander-University of Erlangen-Nuremberg, Germany Marius Schaefer, Austral University of Chile, Chile Jussi Grießinger, Friedrich-Alexander-University of Erlangen-Nuremberg, Germany Gino Casassa, Universidad de Magallanes, Punta Arenas, and General Directorate

of Water (DGA), Chile

Citation: Schneider, C., Braun, M. H., Schaefer, M., Grießinger, J., Casassa, G., eds. (2020). Climate Impacts on Glaciers and Biosphere in Fuego-Patagonia. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-758-4

# Table of Contents

Inés Dussaillant, Etienne Berthier and Fanny Brun

*20 Cirque Glacier on South Georgia Shows Centennial Variability Over the Last 7000 Years*

Lea T. Oppedal, Jostein Bakke, Øyvind Paasche, Johannes P. Werner and Willem G. M. van der Bilt

Stephanie S. Weidemann, Tobias Sauter, Rolf Kilian, David Steger, Nicolas Butorovic and Christoph Schneider

*73 Climate and* Nothofagus pumilio *Establishment at Upper Treelines in the Patagonian Andes*

Ana M. Srur, Ricardo Villalba, Milagros Rodríguez-Catón, Mariano M. Amoroso and Eugenia Marcotti

Stephanie S. Weidemann, Tobias Sauter, Philipp Malz, Ricardo Jaña, Jorge Arigony-Neto, Gino Casassa and Christoph Schneider

*121 Snow Cover Change as a Climate Indicator in Brunswick Peninsula, Patagonia*

Francisco Aguirre, Jorge Carrasco, Tobias Sauter, Christoph Schneider, Katherine Gaete, Enrique Garín, Rodrigo Adaros, Nicolas Butorovic, Ricardo Jaña and Gino Casassa

Lukas Langhamer, Tobias Sauter and Georg J. Mayr

# Editorial: Climate Impacts on Glaciers and Biosphere in Fuego-Patagonia

Christoph Schneider 1 \*, Matthias Holger Braun2 , Marius Schaefer 3 , Jussi Grießinger 2 and Gino Casassa4

1 Geography Department, Humboldt-Universität zu Berlin, Berlin, Germany, 2 Institute of Geography, Friedrich-Alexander-University Erlangen-Nuremberg, Erlangen, Germany, 3 Institute of Physics and Mathematics, Austral University of Chile, Valdivia, Chile, 4 Centro de Investigación GAIA Antártica, Universidad de Magallanes, Punta Arenas, Chile

Keywords: Patagonia, Tierra del Fuego, cryosphere, biosphere, climate variability, climate change

**Editorial on the Research Topic**

#### **Climate Impacts on Glaciers and Biosphere in Fuego-Patagonia**

Southernmost South America including Tierra del Fuego and its adjacent islands form the only significant land masses within the mid-latitude Westerlies of the Southern Hemisphere. Since climate variability in the Southern Hemisphere and its impacts are still insufficiently understood, it is essential to improve our understanding of climate impacts and resulting ecological and glacier responses in the region. The Conference and Workshop "Climate Impacts on Glaciers and Biosphere in Fuego-Patagonia1 " in Berlin, Germany, 14–19 July 2017, explored these topics and resulted in the current Frontiers Research Topic. The activities were initiated and linked to the joint research project "Responses of Glaciers, Biosphere and Hydrology to Climate Variability and Climate Change across the Southern Andes" (GABY-VASA) by University of Magallanes and the Instituto Antártico Chileno, both in Punta Arenas, Chile, and the Humboldt-Universität zu Berlin, the Friedrich-Alexander-University Erlangen-Nuremberg and RWTH Aachen University in Germany on dendroclimatology, climatology and glaciology in Southern Patagonia and the Cordillera Darwin, jointly funded by the Chilean Comisión Nacional de Investigación Científica y Tecnológica and the German Federal Ministry of Research and Education.

The regional climate in Southernmost South America is heavily influenced by the proximity to the oceans. This generates rather weak seasonal cycles with cool to cold summers and moderate to cold winters, especially on the western Pacific side. Slightly more pronounced, continental seasonal cycles are observed in the East of the Andes. While annual mean air temperatures across the region are decreasing from North to South precipitation patterns show very pronounced east-west gradients. The distinctive gradients in precipitation are caused by the north-south striking mountain ranges of the Patagonian Andes, and the northwest-southeast stretching mountain chains of the Cordillera Darwin. Both mountain ranges enforce heavy precipitation on the west and southwest exposed flanks by uplift and dry foehn-like conditions on the leesides (e.g., Holmlund and Fuenzalida, 1995; Schneider et al., 2003; Rasmussen et al., 2007) which produces extremely high drying ratios (Escobar et al., 1992; Carrasco et al., 2002; Smith and Evans, 2007). At inter-annual to decadal time scales atmospheric teleconnections such as the El Niño Southern Oscillation (ENSO) (Schneider and Gies, 2004), Southern Annular Mode (SAM), and Pacific Decadal Oscillation (PDO) are influencing spatial and temporal patterns of both, precipitation and air temperature. For example, positive SAM modes (Garreaud, 2009; Weidemann, Sauter, Kilian et al.) and the PDO (Villalba et al., 2003) are associated with higher air temperatures. Langhamer et al. show that the source of precipitation in the Southern Andes also depends on these teleconnections.

1https://www.geographie.hu-berlin.de/en/professorships/climate\_geography/patagonia\_workshop

#### Edited and reviewed by:

Regine Hock, University of Alaska Fairbanks, United States

\*Correspondence: Christoph Schneider<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Cryospheric Sciences, a section of the journal Frontiers in Earth Science

Received: 26 February 2020 Accepted: 16 March 2020 Published: 09 April 2020

#### Citation:

Schneider C, Braun MH, Schaefer M, Grießinger J and Casassa G (2020) Editorial: Climate Impacts on Glaciers and Biosphere in Fuego-Patagonia. Front. Earth Sci. 8:91. doi: 10.3389/feart.2020.00091

An important aspect is that Patagonia and Tierra del Fuego, together with the sub-Antarctic islands are the only regions where direct proximity between Antarctica and land masses north of the Southern Ocean is given. Such linkages are for example explored with investigations by Hebel et al. on the biosphere and Oppedal et al. for the regional glacier history.

The Humboldt or Peru Current transports cold waters along the Pacific coast of South America and the strong Westerlies cause reduced air temperatures combined with heavy precipitation in Southern South America. This setting induces both, unique eco-systems comprised by extensive bogs, evergreen and broadleaf Nothofagus (rain) forests, as well as characteristical elements of the cryosphere including extensive seasonal snow cover (Aguirre et al.). The mountain ranges host major ice bodies such as the Northern and Southern Patagonia Icefields and the Icefield of the Cordillera Darwin all of them with outlet glaciers reaching down as far as sea level. In addition, a multitude of smaller ice caps and glaciers are wide-spread.

On the eastern side of the Andes, the small ocean basins along the Strait of Magallan, the glacially formed gulfs of Seno Skyring and Seno Otway and large glacial lakes such as Lago Argentino and Lago Viedma provide evidence of extended glaciation until the Last Glacial Maximum (e.g., Rabassa et al., 2011). In recent centuries glacier re-advances cumulated in Little Ice Age maximum glacier extents and associated trimlines and moraines mostly dating from the second half of the nineteenth century as for example explored in Meier et al.. The underlying changes in climate can be evidenced by analyzing dendroclimatological records as pinpointed for example for N. pumilio in this Research Topic by Srur et al. and Grießinger et al..

Due to the temperate nature of glaciers associated with its high mass turnover, ice bodies in Patagonia and Tierra del Fuego tend to respond quickly and highly sensitive to atmospheric forcing as shown by analyzing the energy and climatic mass

#### REFERENCES

balance in details in Weidemann, Sauter, Malz et al.. There is overall consensus that anthropogenically induced global climate warming has a major impact on the cryosphere in Fuego-Patagonia. Nonetheless, it is still under debate how precisely recent glacier retreat and glacier thinning (e.g., Braun et al., 2019; Dussaillant et al.) can be explained by the combination of recent climate warming and glacier adjustment to late nineteenth century to twentieth century climate variability. Since there are no direct measurements on the icefields, accumulation can only be deduced from firn cores and numerical modeling providing a wide range of 7–34 m per year (Aristarain and Delmas, 1993; Shiraiwa et al., 2002; Schwikowski et al., 2013; Lenaerts et al., 2014; Schaefer et al., 2015). It is disputable if—despite the high drying ratios—such high accumulation rates are feasible by atmospheric moisture transport. An effort toward solving this issue is undertaken by Langhamer et al. by analyzing sources of precipitation.

#### AUTHOR CONTRIBUTIONS

CS drafted the editorial. MB, JG, MS, and GC subsequently contributed to the text of the editorial. This contents of the editorial was jointly discussed by all authors.

#### FUNDING

The workshop on Climate Impacts on Glaciers and Biosphere in Fuego-Patagonia on which this Research Topic builds was funded by the CONICYT-BMBF project GABY-VASA under Grant Nos. 01DN15007 and 01DN15020. Further support for the workshop by Geo.X, the Research Network for Geosciences in Berlin and Potsdam, Germany and the Embassy of Chile in Berlin, Germany are gladly acknowledged.

Extreme precipitation and climate gradients in Patagonia revealed by highresolution regional atmospheric climate modeling. J. Clim. 27, 4607–4621. doi: 10.1175/JCLI-D-13-00579.1

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Schneider, Braun, Schaefer, Grießinger and Casassa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Geodetic Mass Balance of the Northern Patagonian Icefield from 2000 to 2012 Using Two Independent Methods

#### Inés Dussaillant 1,2 \*, Etienne Berthier 1 and Fanny Brun1,2

1 LEGOS/OMP, Université de Toulouse, CNES, Centre National de la Recherche Scientifique, IRD, UPS, Toulouse, France, 2 IGE, Université Grenoble Alpes, Centre National de la Recherche Scientifique, IRD, Grenoble, France

We compare two independent estimates of the rate of elevation change and geodetic mass balance of the Northern Patagonian Icefield (NPI) between 2000 (3,856 km2 ) and 2012 (3,740 km2 ) from space-borne data. The first is obtained by differencing the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) from February 2000 and a Satellite pour l'Observation de la Terre 5 (SPOT5) DEM from March 2012. The second is deduced by fitting pixel-based linear elevation trends over 118 DEMs calculated from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) stereo images acquired between 2000 and 2012. Both methods lead to similar and strongly negative icefield-wide mass balance rates of −1.02 ± 0.21 and −1.06 ± 0.14 m w.e. yr−1 respectively, which is in agreement with earlier studies. Contrasting glacier responses are observed, with individual glacier mass balance rates ranging from −0.15 to −2.30 m w.e. yr−1 (standard deviation = 0.49 m w.e. yr−1 ; N = 38). For individual glaciers, the two methods agree within error bars, except for small glaciers poorly sampled in the SPOT5 DEM due to clouds. Importantly, our study confirms the lack of penetration of the C-band SRTM radar signal into the NPI snow and firn except for a region above 2,900 m a.s.l. covering <1% of the total area. Ignoring penetration would bias the mass balance by only 0.005 m w.e. yr−1 . A strong advantage of the ASTER method is that it relies only on freely available data and can thus be extended to other glacierized areas.

#### Edited by:

Matthias Holger Braun, University of Erlangen-Nuremberg, Germany

#### Reviewed by:

Marco Möller, University of Bremen, Germany Ninglian Wang, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, China

> \*Correspondence: Inés Dussaillant

<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Cryospheric Sciences, a section of the journal Frontiers in Earth Science

Received: 19 September 2017 Accepted: 22 January 2018 Published: 12 February 2018

#### Citation:

Dussaillant I, Berthier E and Brun F (2018) Geodetic Mass Balance of the Northern Patagonian Icefield from 2000 to 2012 Using Two Independent Methods. Front. Earth Sci. 6:8. doi: 10.3389/feart.2018.00008 Keywords: glacier mass balance, geodetic method, North Patagonian Icefield, radar penetration, Patagonia

## INTRODUCTION

Patagonian glaciers and icefields, including the Northern Patagonian Icefield (NPI), are essential indicators of climate change in the southern hemisphere and also stand as one of the major contributors to recent sea level rise (Gardner et al., 2013; Mernild and Wilson, 2016). Direct glaciological observations in this remote region are still scarce, mainly due to the region's inaccessibility and its harsh climate (Buttstädt et al., 2009; Mernild et al., 2015; Schaefer et al., 2015). In recent decades, remote sensing techniques have enabled the spatial coverage of glaciological observations to be extended. They are particularly useful in remote and inaccessible areas like the NPI, making it possible to observe and measure glacier variations and to estimate mass balance rate at regional scales (Rignot et al., 2003; Rivera et al., 2007; Lopez et al., 2010; Davies and Glasser, 2012; Willis et al., 2012a,b; Dixon and Ambinakudige, 2015; Jaber et al., 2016). Among those studies, mass balance rate estimates are often based on the geodetic method, which allows calculation of glacier volume changes by differentiation of two or more multi-temporal digital elevation models (DEMs) (e.g., Paul et al., 2015; Marzeion et al., 2017). In particular for the NPI, Willis et al. (2012a) measured mass balance using the Shuttle Radar Topography Mission (SRTM) DEM and a collection of freely available Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs from 2001 to 2011. Jaber et al. (2016) also provided a mass balance estimate for 2000–2014 by comparing two interferometric DEMs, SRTM, and TanDEM-X.

The goal of our study is to produce two new and independent geodetic estimates of the NPI mass balance rate for the time period 2000–2012 using, on one hand, SRTM and SPOT5 and, on the other hand, ASTER DEMs. Although the NPI mass balance rate during this time period has already been estimated in two previous studies (Willis et al., 2012a; Jaber et al., 2016) our contribution is justified by:

(i) The use of ASTER DEMs calculated using the Ames Stereo Pipeline (ASP) (Shean et al., 2016). As shown in the supplement of our study, these DEMs have less artifacts than the 14DMO DEM available from LP-DAAC and used in many earlier studies including the one of Willis et al. (2012a) over the NPI.

(ii) The need to verify the lack of penetration of the February 2000 SRTM radar signal. Until now the effect of radar wave penetration into snow and ice has been assumed to be negligible in Patagonia due to the likely wet condition at the glacier surface in summer, favoring surface scattering of the signal (Willis et al., 2012a; Jaber et al., 2013, 2016). However, this assumption has not been verified yet.

(iii) The need to account for changes in ice-covered areas when calculating the glacier mass balance. None of the earlier studies considered these area changes, known to be large in Patagonia (Davies and Glasser, 2012).

### STUDY AREA: NORTHERN PATAGONIAN ICEFIELD, CHILE

The NPI is the second largest temperate ice mass in the southern hemisphere, after the Southern Patagonian Icefield. Located in Chile, it extends from 46◦ 30′ S to 47◦ 30′ S, stretching for almost 125 km north–south (**Figure 1**). It covers an area of ∼3,740 km2 (measured in this study for year 2012) extending from sea level to altitudes higher than 4000 m a.s.l. at the summit of Mount San Valentin. It is composed of 38 glaciers larger than 0.5 km2 . San Quintin is the largest glacier, covering 790 km2 . San Rafael (714 km2 ) is the only tidewater calving glacier of the icefield and the lowest-latitude tidewater glacier as well as one of the fastest-flowing glaciers in the world (Rignot et al., 1996). The largest proportion of glaciers present fresh-water calving fronts (66% of the total area), whilst a smaller proportion (16%) correspond to land-terminating fronts (Rivera et al., 2007).

The climate of the Northern Patagonian region is dominated by the southern hemisphere westerlies that bring moisture from the Pacific, and the equatorial Pacific sea surface temperature that regulates the El Niño Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (Aravena and Luckman, 2009). This later two large scale modes of natural climate variability, together with the southern annular mode (SAM), are key modulators of temperature and precipitation in the southern Andes at interannual and inter-decadal timescales (Gillett et al., 2006; Garreaud et al., 2009). The north–south orientation and high altitude of the Andes chain at this latitude acts as a geographical barrier for the moisture, generating a föhn effect that creates a strong west-to-east gradient of precipitation, with considerably dryer conditions and higher temperatures on the eastern side of the icefield (Fujiyoshi et al., 1987). The average annual precipitation between 1975 and 2011 is 8.03 ± 0.37 m (Schaefer et al., 2013) with little seasonal variability. The 0◦C isotherm altitude is at about 2,000 m a.s.l. during the summer, dropping to 900 m a.s.l. during winter (Barcaza et al., 2009).

### DATA

#### DEMs and Satellite Stereo Images SRTM DEM

The SRTM acquired elevation data between 11 and 22 February 2000 using an interferometric aperture radar sensor (SAR) (Farr et al., 2007). We used the SRTM DEM processed by NASA from the raw C-band radar signal at 1 arcsecond resolution (about 30 m). Due to radar layover or loss of coherence over vegetation and steep terrain, the SRTM DEM may contain data voids that are sometimes filled by interpolation or using an external DEM. For the purposes of DEM comparison, we use the non-void-filled version in which data voids cover 13% of the total NPI glaciated area. By specific elevation bands the percentage of data voids is 3% below 1,000 m a.s.l (total area 1,163 km2 ), 22% between 1,000 and 2,000 m a.s.l. (2,323 km2 ) and 13% above 2,000 m a.s.l. (370 km2 ). A second, void-filled SRTM DEM (version 4, 3 arcsec resolution) produced by the CGIAR-CSI (http://srtm.csi.cgiar.org) is also used to obtain a complete hypsometry of the icefield, necessary for the mass balance calculation.

#### SPOT5-HRS DEM

The SPIRIT project produced DEMs over glaciated zones from stereo-pairs acquired by the SPOT5 high-resolution stereoscopic (HRS) sensor (Korona et al., 2009). A 40 m DEM, acquired on 18 March 2012, provides a complete coverage of the NPI (**Supplementary Figure S1**). The DEM is delivered with a reliability mask and one ortho-image generated from the rear HRS image. Reliability mask values range from 0 to 255, with 0 indicating the least reliable elevations that are thus masked out in the SPOT5 DEM. Data voids cover 16% of the total NPI glaciated area. By specific elevation bands the percentage of data voids is 7% below 1,000 m a.s.l, 22% between 1,000 and 2,000 m a.s.l. and 13% above 2,000 m a.s.l.

#### TERRA–ASTER Stereo Images

Raw stereo images obtained by the ASTER sensor on board the TERRA satellite were ordered free of charge from the

ECHO/REVERB website (REVERB)1 from April 2000 to April 2012. After excluding all images with cloud cover fraction larger than 80%, a total of 118 ASTER Level 1A stereo images intersecting the NPI were selected. For the sake of concision, the list of ASTER images used in this study were not provided but is available upon request to the corresponding author.

#### LANDSAT Images

Four late-summer nearly cloud-free LANDSAT L7 ETM+ orthoimages at 30 m spatial resolution acquired from the USGS Earth Explorer were used to improve glacier and water outlines for years 2000 and 2012 (8 March 2000 and 21 January, 22 February, 26 April 2012).

### Glacier Outlines, Water Bodies, and Ice Divides

Initial outlines for glacier limits and water bodies over the NPI area were obtained from online datasets. The Global Land Ice Measurements from Space (GLIMS) database and the Randolph Glacier Inventory version 5.0 (RGI 5.0) offer access to glacier outlines from previous studies and inventories over the Southern Andes region (Rivera et al., 2007; Davies and Glasser, 2012; Pfeffer et al., 2014; Mouginot and Rignot, 2015). Following visual assessment, RGI5.0 was selected as the best existing outline, corresponding to the ice extent of year 2001 (RGI6.0 was not available at the time of our study). RGI5.0 was not used "as is" but was improved as explained in section Improvement of the Glacier Outlines. Water body outlines were acquired from the USGS Earth Explorer SRTM water body data (earthexplorer.usgs.gov), and the ice divides derived from velocity fields (Mouginot and Rignot, 2015) were used to split the icefield into individual glaciers.

### METHODOLOGY

#### Improvement of the Glacier Outlines

The RGI 5.0 glacier outline was corrected manually using latesummer and cloud-free LANDSAT, ASTER, and SPOT5-HRS images as references to adjust the position of the ice for the years 2000 and 2012. RGI 5.0 ice divides were replaced by the ice divides from Mouginot and Rignot (2015). A visual comparison between our final outlines and best existing previous ones can be observed in **Supplementary Figure S2**. Significant errors and gaps were found in the SRTM water body data, so manual corrections were also performed. Furthermore, a cloud mask was drawn manually to exclude the cloudy pixels over the SPOT5- HRS DEM. No cloud mask was created for ASTER DEMs because the lack of correlation over the clouds leads directly to data voids in the DEMs.

### Glacier Volume Change and Mass Balance

As explained above, two different methods were applied to derive the map of the rate of elevation changes over the NPI: SRTM and SPOT5 DEM differencing (SPOT-SRTM hereafter) and

1REVERB. Available online at: https://reverb.echo.nasa.gov/#utf8=%E2%9C%93& spatial\_map=satellite&spatial\_type=rectangle

linear trend analysis of multiple ASTER DEMs (ASTER\_trend hereafter).

#### Generation of ASTER DEMs

ASP is a NASA open source pipeline that allows for the mass production of DEMs out of high-resolution stereo images (Shean et al., 2016). We generated 118 ASTER-ASP DEMs with a grid spacing of 30 m from the raw ASTER L1A stereo images. ASP processing parameters are available upon request to the authors. Each ASTER DEM has a different footprint and, given the 60 km swath of ASTER, does not cover the entire NPI. The result is a spatially-varying number of ASTER DEMs available over the NPI (**Supplementary Figure S1**).

#### Adjustments of the DEMs

#### **Horizontal adjustment of the DEMs**

The DEMs used in this study are produced independently using different methods and without ground control points such that a horizontal shift is expected to exist between them. 2D coregistration (i.e., horizontal adjustment) should be performed on stable terrain and in regions where DEMs are most reliable (Paul et al., 2015). Stable terrain corresponds to all pixels with valid elevation values after excluding ice, snow, water bodies, and clouds.

The SRTM DEM was used as reference, and all other DEMs (SPOT5 and ASTER) were adjusted to it. For SPOT-SRTM, the horizontal shift is determined by an iterative minimization of the standard deviation of the elevation difference (dh) on stable terrain (Berthier et al., 2007). For ASTER\_trend, an approach based on terrain slope and aspect was used to calculate the horizontal shift (Kääb, 2005; Nuth and Kääb, 2011). Both coregistration methods lead to similar results. The Nuth and Kääb method is preferred to coregister ASTER DEMs as it is numerically faster and thus more suitable for managing large amounts of data (Paul et al., 2015).

#### **Vertical adjustment of the DEMs**

Three different types of vertical biases are quantified on stable terrain and corrected in SPOT5 and ASTER DEMs. We first correct for a simple vertical shift, then biases along and across the satellite tracks and finally we account for elevation difference due to different DEM resolution.

For SPOT5 and ASTER DEMs, the median of the altitude difference with the SRTM DEM over stable terrain is first substracted. To avoid anomalous values that are usually found over steep terrain, only pixels where the slope is <30◦ are considered to estimate this median vertical shift (Berthier et al., 2016). We recognize that this limit is somewhat arbitrary but further justified by the fact that over 90% of the NPI surface has a slope lower than 30◦ .

For SPOT5 and ASTER DEMs, all derived from stereo imagery, the satellite acquisition geometry and its particular trajectory may induce biases in the DEMs along and across the satellite track direction (Berthier et al., 2007; Scherler et al., 2008; Nuth and Kääb, 2011). We rotate the original coordinate system of the SPOT5 and ASTER DEMs tracks to match their specific trajectory (Nuth and Kääb, 2011). The elevation difference (dh) is then computed over stable terrain along and across the satellite track, and corrected using a fifth-order polynomial fit.

We applied a curvature-dependent correction to the SPOT5 and ASTER DEMs. This correction is calculated on the stable terrain by fitting a fifth-order polynomial to the elevation difference with the SRTM DEM plotted as a function of terrain maximum curvature (Gardelle et al., 2012).

After all adjustments, pixels from the SPOT5 and ASTER DEMs which have an absolute elevation difference with SRTM greater than 150 m are considered as outliers and excluded from all subsequent analysis.

### Rate of Elevation Change and Mass Balance Rate Calculation

The maps of elevation change rates (dh/dt) are obtained differently for the two methods. For SPOT-SRTM, a simple DEM difference is computed and divided by their time separation (12.1 yr) to calculate dh/dtSPOT−SRTM. In the case of ASTER\_trend, for each pixel, a first linear trend is fit through all the available elevation values (from all ASTER DEMs and the SRTM DEM) during 2000–2012. The final dh/dtASTER\_trend is calculated from a second fit considering only the ASTER DEMs within the 95% confidence interval of the first regression line and excluding the SRTM DEM to avoid possible bias due to signal penetration. For the majority of the NPI area, the first valid elevation measurements used in this second linear fit come from a DEM acquired in April 2000, and the last measurements from two DEMs acquired in March 2012 and April 2012 (**Supplementary Figure S3**). We did not observe any obvious temporal concentration of ASTER DEMs. Thus, the time-stamp of the dh/dt maps from the two methods is similar, i.e., 2000– 2012.

In the conversion from dh/dt to mass balance, only reliable pixels from glaciated terrain are considered. For ASTER\_trend, we further exclude all pixels for which the uncertainty of the second linear fit is larger than 3 m yr−1 (at the 95% confidence level).

NPI individual glaciers have shown different volume change rates during recent years, with differential responses within glaciers and a strong contrast between the eastern and western margins (Rivera et al., 2007; Willis et al., 2012a). To take into account these local effects, separate calculations were made for each glacier, and then the overall volume change rate and mass balance rate of the NPI were calculated as the area-weighted sum of every glacier. Each individual glacier is divided into 50 m altitude bands using the void-filled version of the SRTM DEM. For each altitude band, we compute the mean dh/dt after excluding all values lying outside of ±3 ∗NMAD (normalized median absolute deviation).This is an efficient way to exclude outliers (Berthier et al., 2004; Gardelle et al., 2013). The band average dh/dt is then weighted by its area to calculate the glacierwide average dh/dt. In the rare cases where no measurement is available for an altitude band, the rate of elevation change is assumed to be 0.

Following Fischer et al. (2015), we considered the varying areas of each individual glacier between 2000 and 2012 to convert volume change rate to mass balance rate. We use a value of 850 ± 60 kg m−3 as the volume-to-mass conversion factor, which is reasonable regarding the duration of our study (>5 years) (Huss, 2013).

#### Penetration of the SRTM C-Band Radar Signal Over Snow

The ASTER\_trend method allows rough estimation of the SRTM penetration depth (Wang and Kääb, 2015; Berthier et al., 2016). By extrapolating dh/dt to the date of acquisition of the SRTM DEM, a reconstructed SRTM DEM (called hereafter SRTMrec) can be obtained. The differencing of SRTMrec-SRTM maps the penetration depth of the C-band radar signal.

#### Error Assessment

Formal uncertainties for mass balance estimates were calculated taking into account six main sources of error: errors on dh/dt, errors in glacier area, error on the density conversion factor, errors due to data voids, errors due to seasonal cycle sampling irregularities and departure from a linear trend.

#### Uncertainty in Rate of Elevation Change σ dh dt

The uncertainty on the elevation change of a pixel is mainly given by the uncertainty inherent in the original DEMs, including the random and systematic errors produced during their coregistration. To assess this error we split the stable terrain dh/dt maps into 3 × 3 tiles as in Berthier et al. (2016). For each tile the median of dh/dt on stable ground is computed to visualize the spatial pattern of biases in dh/dtASTER and dh/dtSPOT−SRTMmaps (**Figure 2**). We then compute the error in dh/dt as the mean of the absolute difference of the median dh/dt for these 9 tiles (**Figure 2**). We obtained an error σdh dt of 0.07 m yr−1 for SPOT-SRTM and 0.06 m yr−1 for ASTER\_trend.

#### Uncertainty in Area (σA)

Glacier outlines were improved manually and contain an error that is inherent to the individual that performed the work. The error in glacier area was assessed by drawing buffers around the glacier outline (Granshaw and Fountain, 2006). A one-pixel buffer is a reasonable error level for glacier outlines produced manually, except for debris-covered glaciers (Paul et al., 2015). As the NPI debris coverage is limited, a one-pixel buffer was applied to the NPI outline, leading to an error of about 5% of the measured area. This value is consistent with Paul et al. (2015) and Pfeffer et al. (2014).

#### Uncertainty in Density Factor σfρ

For density the error of ±60 kg m−3 is obtained from Huss (2013).

#### Uncertainty in Data Gaps (g)

A conservative multiplicative factor of 5 is used for the uncertainties in the unsurveyed areas (Berthier et al., 2014).

#### Uncertainties due to the Uneven Sampling of the Seasonal Cycle

Systematic errors in the dh/dt from ASTER\_trend could come from an irregular sampling of the seasonal elevation change cycle as the ASTER DEMs were acquired at different times of the year. For example, a large systematic error (i.e., too negative mass balance) could occur if ASTER DEMs were acquired mostly in winter at the beginning of the study period and mostly in summer toward the end of the study period. To quantify the magnitude of this potential error, we sample a hypothetical seasonal cycle at the time of acquisition of the ASTER DEMs. The amplitude of the seasonal cycle is poorly known for Patagonia and must be spatially variable, so we prefer to use a normalized value of 1 m, with a minimum at the end of April and a maximum at the end of October. As shown in **Figure 3**, the ASTER temporal sampling is sufficiently random that the regression slope is only −0.0013 m yr−1 . Even if the seasonal elevation change cycle was as large as 10 m (a high value but not unrealistic for the maritime NPI), the resulting systematic error would be 0.013 m yr−1 , well below other sources of uncertainties. This source of error is thus neglected.

#### Uncertainties Due to Departure from a Linear Trend

The assumption of linearity is one of the weaknesses of our ASTER-based method, as in principle glacier elevation change do not follow a linear behavior. To detect potential changes of the ASTER trends (and mass balances) over time, we divided 2000– 2012 into two sub-periods: 2000 to March 2007 (2007.3) and 2007 to May 2012 (2012.5). We did not split the 12-year study period in two sub-periods of equal duration because there was no DEM for year 2006. Conversely, several useful ASTER DEMs were acquired in 2007 (late summer, no clouds, and good spatial coverage of the NPI), thus we decided to cut the period at 2007 and include summer 2007 DEMs in both trends. The NPI-wide mass balance from 2000 to 2007.3 was −1.06 ± 0.20 m w.e. yr−1 , and −0.99 ± 0.20 m w.e. yr−1 from 2007 to 2012.5. Given error bars, mass balance rates for both sub-periods are not different and agree with the estimate for the whole period 2000–2012 (−1.06 ± 0.15 m w.e. yr−1 ). This result justifies the assumption of linearity and demonstrates that, within the error bars of our method, we are not able to detect a temporal change in the NPI mass balance.

We then consider all the significant sources of error to calculate our formal uncertainties on mass balance. First, the uncertainty in glacier rate of volume change σdV dt , assuming that the uncertainty on area is independent of the uncertainty on rate of elevation change, was calculated as:

$$
\sigma\_{\frac{dV}{dt}} = \sqrt{\left(\sigma\_A \cdot \frac{\overline{dh}}{dt}\right)^2 + \left(\sigma\_{\frac{dh}{dt}} \cdot \left(5\mathfrak{g} + \left(1-\mathfrak{g}\right)\right)A\right)^2} \tag{1}
$$

where dh dt is the mean rate of elevation change, σA = 5 100A, and g is the proportion of data gaps from the total area.

Assuming that the uncertainty in rate of volume change is independent of the uncertainty of the density conversion

factor, the uncertainty on geodetic mass balance rate σM˙ , was calculated as:

$$
\sigma\_{\dot{M}} = \sqrt{\left(\sigma\_{\frac{dV}{dt}} \cdot f\_{\rho}\right)^2 + \left(\sigma\_{f\_{\rho}} \cdot \frac{dV}{dt}\right)^2} \tag{2}
$$

where fρ = 850 kg m−3 is the density conversion factor and σfρ = 60 kg m−3 is the uncertainty on the density conversion factor (Huss, 2013).

#### RESULTS

#### Area Changes

The NPI had a total glacierized area of 3,856 ± 211 km2 in 2000, and 3,740 ± 200 km2 in 2012. The total area loss for the period 2000–2012 is 116 km2 (3%), at a rate of −9.7 km2 yr−1 .

Areal changes are heterogeneous and do not follow any obvious spatial pattern. Largest areal changes are observed on glaciers U4 (19%) located in the south-east, and glaciers Circo (12%) and Grosse (10%) located in the north of the NPI. The largest frontal retreats were observed in the tongues of two glaciers in the south-west, Steffens (3.3 km frontal retreat) and HPN1 (2.5 km retreat over its main tongue and even 3.1 km retreat over one of its tributaries), and San Quintin in the northwest (2.5 km frontal retreat). Glaciers Gualas, Grosse and Reicher in the north and glacier Acodado in the south also showed significant frontal retreats of the order of 1 and 2 km between 2000 and 2012.

## Elevation Change Rates and Glacier Mass Balance Rates

#### Entire NPI

Both methods (SPOT-SRTM and ASTER\_trend) resulted in a similar spatial distribution of dh/dt over the NPI (**Figure 4**), with maximum thinning rates localized at low elevations over the glacier tongues. Icefield-wide mass balances rates during 2000– 2012 are also in excellent agreement: −1.02 ± 0.21 m w.e. yr−1 for SPOT-SRTM and −1.06 ± 0.15 m w.e. yr−1 for ASTER\_trend (**Table 1**).

Strongest thinning rates are observed over the outlet glaciers of the western margins of the icefield, with dh/dt rates as negative as −10 m yr−1 for the terminating tongues of Frankel, Benito, HPN1 and Acodado glaciers. The largest differences between the two dh/dt maps are observed over San Rafael glacier, where thinning rates are more negative for SPOT-SRTM, and over the north-east margins. Looking closely at the dh/dt vs. altitude curves for the two methods, a strong thickening is observed above 2,900 m a.s.l. for SPOT-SRTM, whereas a slight thinning

TABLE 1 | NPI-wide volume change rate, thinning rate, and mass balance rate estimates for the 2000–2012 time period for the two methods.

is measured by ASTER\_trend (**Figure 4C**). The anomalously positive values in SPOT-SRTM dh/dt curves at high altitude will be discussed in section Penetration of SRTM Signal Over NPI.

Spatial coverage of valid dh/dt rates is different for the two methods, with 68% of the total NPI area covered in SPOT-SRTM, and a better coverage of 82% in ASTER\_trend. Almost complete coverage was achieved below 1,000 m a.s.l. in both cases. Higher concentrations of data gaps exist in the accumulation zone (>1,000 m a.s.l.) with 35% gaps for SPOT-SRTM vs. 22% for ASTER\_trend. Both the SPOT5 reliability mask and the ASTER error grid produced during the linear fit suggest larger errors in this region. The lack of image texture over the flat and homogeneous snow surface makes correlation of the stereo images either impossible or erroneous. Higher altitudes are better represented in ASTER\_trend, with only 10% of data gaps above 2,000 m a.s.l. vs. 49% for SPOT-SRTM; the latter is mainly due to voids in the SRTM DEM.

#### Individual Glaciers

The coverage of individual glaciers with valid dh/dt measurements is much better for ASTER\_trend, where all glaciers show two-thirds (66%) of their total area with valid data. In the case of SPOT-SRTM, coverage is not as good. Fourteen out of thirty-eight glaciers, corresponding to 77% of the total NPI area, show between one- and two-thirds of their area with valid measurements, while 11 glaciers (6% of the NPI area) show less than one-third (**Figure 5A**).

Individual glacier mass balance rates (M˙ ) for the period 2000–2012 are shown in **Figure 5** and **Supplementary Table S1** for both methods. All glaciers show negative M˙ ASTER\_trend. Contrarily, M˙ SPOT−SRTM for glaciers Grosse, Exploradores, Bayo, and Leones in the north-east part of the NPI and glacier Pissis in the south are slightly positive. However, these five glaciers were poorly sampled due to gaps in the SPOT5 DEM. Both methods agree on highly negative mass balances rates for southwest glaciers. The most negative M coincide in glacier HPN1 ˙ with values of −2.38 ± 0.32 and −2.30 ± 0.24 m w.e. yr−1 for ASTER\_trend and SPOT-SRTM respectively. The M for the ˙ three largest glaciers of the icefield, San Quintin, San Rafael and Steffens, accounting together for 51% of the total area, are strongly negative and in good agreement between methods.

For all glaciers larger than 100 km2 , M from both methods ˙ agree within error bars, and their absolute difference is smaller than 0.3 m w.e. yr−1 (**Figure 6**). Conversely, smaller glaciers tend to show larger disagreement between methods, with M˙ differences up to 1.6 m w.e. yr−1 . The largest differences are observed for glaciers poorly sampled in SPOT-SRTM (green and orange dots).

#### DISCUSSION

#### Area Changes

Our area estimate for 2000 is in good agreement with previous studies, especially with Jaber et al. (2016), suggesting that our error bars are conservative (**Table 2**). Davies and Glasser (2012) obtained a slightly larger value for the NPI area in 2001 because they considered some of the rock nunataks as part of the icefield

altitude; red dots represent SPOT-SRTM and orange dots represent ASTER\_trend. The histograms show the hypsometry of the NPI by 100 m elevation bands (gray bars), and the corresponding measured area from SPOT-SRTM (blue bars) and ASTER\_trend (purple bars). (Note that hypsometry is shown here every 100 m bands. This is to help visualization although the entire processing was made using 50 m bands).

area. Similarly, Rivera et al. (2007) included adjacent small icy areas as part of the icefield. To our knowledge, we provide the first estimate of NPI area for 2012.

Glacier areal changes are not homogeneous through the NPI. In general the largest frontal retreats were observed for glacier tongues situated in the south-west of the icefield (glaciers Steffens and HPN1), also in agreement with previous studies (Rivera et al., 2007; Lopez et al., 2010; Davies and Glasser, 2012; Willis et al., 2012a). No convincing relationship was observed between glacier-specific area loss and mass balance. This conclusion is in agreement with earlier findings for Alaskan glaciers (Arendt et al., 2002).

We accounted for these area changes in our mass balance rate estimates. At the scale of the entire icefield, the effect is minor as

the NPI mass balance rate is changed by 0.01 m w.e. yr−1 if only the 2000 inventory is available. Yet, the mass balance rate can vary by up to 10% for the fastest retreating glaciers.

#### Penetration of SRTM Signal Over NPI

Some positive dh values were found above 2,900 m a.s.l. with the SPOT-SRTM method. These positive values are anomalous as, to our knowledge, no study documented thickening in the upper part of the NPI since 2000. **Figure 7** shows the distribution with altitude of our estimate of the SRTM penetration depth (SRTMrec-SRTM), a side product of the ASTER\_trend method. Positive values are observed above 2,900 m a.s.l. This is consistent with the apparent thickening observed above 2,900 m a.s.l. in the SPOT-SRTM dh/dt and thus suggest possible penetration of the SRTM C-band signal on the NPI snow located above 2,900 m a.s.l.

The SRTMrec-SRTM and the SPOT-SRTM curves in **Figure 7** do not match perfectly over this high-altitude zone. Likewise, the SRTMrec-SRTM curve shows positive values between 1,000 and 2,900 m a.s.l. However, these differences remain within the uncertainties of the SRTM penetration depth, assumed to be roughly ±3.1 m in Berthier et al. (2016), and cannot be further interpreted. We further confirmed that this ±3.1 m is conservative by extrapolating dh/dtASTER\_trend to reconstruct two DEMs at the dates of acquisition of the SPOT5 2012 DEM and of an additional SPOT5 DEM acquired 18 May 2005 (not used in the study for measuring elevation changes but with similar specifications as the SPOT5 DEM 2012). Differencing between SPOT5rec\_2012-SPOT52012 and SPOT5rec\_2005-SPOT52005 yielded mean dh on the NPI of 1.03 m and 0.14 m, respectively, both values well below the 3.1 m assumed error.

If we consider the region above 2,900 m a.s.l., the mean penetration depth is 10 m, a value similar to the 10 m C-band penetration depth found on the dry and cold firn of Greenland and Alaska (Rignot et al., 2001). Averaged over the whole icefield, the mean penetration depth is negligible (0.02 m), in agreement with a recent study in New Zealand (also in the southern hemisphere), where no penetration was found (Wang and Kääb, 2015).

Up to now, SRTM penetration depth has been considered negligible in Patagonia due to almost constant wet conditions at the glacier surface in February 2000 (Jaber et al., 2013). Here we show that penetration is actually occurring at the highest elevations, where temperatures are lower and firn is probably dry even in summer. Nevertheless, only 0.75% of the NPI area (29 km2 ) is found above 2,900 m a.s.l., a fraction small enough not to affect the overall mass balance rate of the icefield (the corresponding bias is <0.01 m w.e. yr−1 ). Therefore, SRTM penetration can be assumed negligible over the NPI. Still, future geodetic mass balance studies in Patagonia using SRTM Cband DEM, should carefully consider the penetration depth, especially when a significant fraction of the glacier lies at high elevations.

The lack of SRTM penetration, except at the highest elevations, is in line with a backscattering analysis using the SRTM swath image data (SRTM-IMGR) in Jaber et al. (2016). They found low backscattering values over the plateau,

indicating the presence of wet snow and consequent insignificant penetration depth. Less negative backscattering values are also observed over the highest zones, suggesting dryer snow where the C-band signal can penetrate.

#### Comparison of the Two Geodetic Methods

Both methods led to similar estimates of the NPI mass balance rate for the 2000–2012 period (**Table 3**).

The SPOT-SRTM method has already been used in the southern Andes (Falaschi et al., 2016) but this method is all the more relevant now that we have verified the negligible penetration of the SRTM C-band signal at this latitude.

Even though individual ASTER DEMs have a higher noise level than the SPOT5 DEM when subtracted from the SRTM DEM (mean standard deviation of the residual for ASTER DEMs: 14 vs. 8 m for SPOT5 DEMs), their larger number led to final errors on stable ground very similar to those obtained with the SPOT-STRM method.

Icefield-wide mass balance error differs slightly between methods (±0.21 m w.e. yr−1 for SPOT-SRTM and ±0.15 m w.e yr−1 for ASTER\_trend). This difference is mainly related to the amount of data gaps, higher in SPOT-SRTM. The distribution of data gaps coincides over the NPI plateau (∼1,200–1,700 m a.s.l.) for both methods. Yet, over some glaciers in the north-east margins of the NPI, specifically Grosse, Bayo, and Exploradores glaciers, less than one-third of their area is covered with the SPOT-SRTM method due to clouds in the SPOT5 images. Coverage is also restricted on small steep glaciers like glacier TABLE 2 | Comparison of the NPI total area in this and previous studies.

Pissis in the south, where a high percentage of data gaps exists over the SRTM DEM. These latter glaciers are those for which the mass balance rate differs most between the two methods. The redundancy of the acquisitions used in the ASTER\_trend method overcomes this issue.

The 2000–2012 M are in agreement within error bars between ˙ the two methods for all glaciers with coverage larger than 66%. Specific cases where M do not agree within error bars correspond ˙ to glaciers poorly sampled in the SPOT-SRTM method due to clouds or SRTM data gaps. The three largest glaciers of the icefield (51% of the total area) present M differences smaller than 0.2 m ˙ w.e. yr−1 .

Availability of the data is also important to consider when comparing the two methods. A strong advantage of the ASTER\_trend method is that it is based on freely and nearly globally available ASTER L1A stereo images. SPOT5 stereo images were only acquired over specific areas and are not yet freely available. Thus, the ASTER\_trend has the potential to be extended to wider areas.

#### Comparison with Previous Estimates

Many studies already estimated the geodetic mass balance rates of the NPI for diverse time periods (Rignot et al., 2003; Aniya, 2007; Rivera et al., 2007; Glasser et al., 2011; Willis et al., 2012a; Jaber et al., 2016). However, only two of these studies examined the same time period as us and can therefore be directly compared with our results (**Table 3**). NPI-wide surface mass balances have also been estimated from atmospheric models (Schaefer et al., 2013; Mernild and Wilson, 2016), but these numbers cannot be compared with ours as they do not account for the mass loss by calving which is an important component of the mass loss for the NPI.

Generally, good agreement is found between our two estimates and previous geodetic estimates (**Table 3**). Although they used ASTER stereo images like us, Willis et al. (2012a) found a smaller mass loss. The difference may result either from the use of different ASTER DEMs, different coregistration methods, the amount of DEMs considered to extract dh/dt (55 in Willis et al., 2012a, 118 in our study) or the slight difference between time periods. We generated our own ASTER DEM using ASP, and they exhibit less artifacts than the DEMs from the 14DMO product used in Willis et al. (2012a) (**Supplementary Figure S4**).

Our mass loss estimates agree with Jaber et al. (2016) who compared the SRTM DEM with a TanDEM-X DEM mosaic. These two DEMs, derived from radar interferometry,

TABLE 3 | Global volume change rates estimates for the NPI from this and previous studies during similar time periods.

are not affected by clouds and provide good coverage of the accumulation zone. Nevertheless, their coverage of the icefield was not complete, as northernmost NPI glaciers were not covered by TanDEM-X DEMs.

Error estimates in Jaber et al. (2016) and Willis et al. (2012a) are three to four times smaller than ours. One reason for our larger errors could be the factor of five that we conservatively allocated to errors over data gaps. However, mass balance differences of typically 0.2 m w.e. yr−1 between our two methods for well-sampled glaciers (**Figure 6**, **Supplementary Table S1**) suggest that our error bars are reasonable and not too conservative.

### CONCLUSIONS

We compared two geodetic methods to obtain the region-wide and glacier-wide mass balance rates of the NPI. The first is based on the standard differencing of two DEMs (SRTM and SPOT5) whereas the second method extracts the rate of surface elevation change by fitting a linear trend to over 118 DEMs derived from ASTER stereo images.

For the entire NPI (about 3,800 km2 ), excellent agreement is found between the two methods. Following previous studies, our analysis confirms the strong mass loss, with an icefield-wide mass balance rate of −1.06 ± 0.15 m w.e. yr−1 (ASTER\_trend). The thinning rate reaches up to 10 m yr−1 for some low-elevation glacier tongues. Importantly, we also verify the lack of significant penetration of the SRTM C-band radar signal in the wet and temperate firn of the icefield except for the highest altitudes where dry snow was present in February 2000 (>2,900 m a.s.l.; <1% of the icefield).

For individual glaciers, larger differences exist between the two methods, especially when the data coverage is more limited. By definition the SPOT5-SRTM estimate does not cover the area where clouds are present in the SPOT5 stereo images. The ASTER\_trend method is less prone to this effect because it relies on numerous (>100) ASTER DEMs. Both methods exhibit data gaps in the upper accumulation areas due to the lack of texture in the optical stereo images. This issue could be solved in the future by using images with higher spatial resolution and better radiometric depth.

Here we have provided insight into the effectiveness, advantages and disadvantages of two methods for estimating the geodetic mass balance rates of Patagonian glaciers and icefields. Both methods lead to coherent and satisfying results, though the free availability of the extensive archive of ASTER stereo images is a strong advantage.

### AUTHOR CONTRIBUTIONS

ID: performed all the DEM analysis with inputs from EB and FB. ID: Led the writing of the paper and all other co-authors contributed to it.

#### ACKNOWLEDGMENTS

ID acknowledges a Ph.D. fellowship from the France Space Agency (CNES) and the Région Occitanie. Etienne Berthier acknowledges support from the CNES and the Programme National de Télédétection Spatiale (PNTS, http://www.insu.cnrs. fr/pnts), grant PNTS-2016-01. ASTER images are courtesy of NASA/METI/AIST/Japan Space systems and U.S./Japan ASTER Science Team. SPOT5 images and DEMs were provided by the SPIRIT project (Korona et al., 2009), funded by CNES. This work could not have been performed without the GLIMS project that allowed populating a vast archive of freely available ASTER stereo images over glaciers. We also acknowledge the Ames Stereo Pipeline developer teams for their help and Vincent Favier for his comments on an earlier version of this manuscript. We are grateful to Ken Moxham for reviewing the English. We thank the editor, Matthias Braun, and two reviewers for their very useful comments on our manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00008/full#supplementary-material

Supplementary Figure S1 | Spatial footprint of the SPOT5-HRS DEM for year 2012 (red line). The color code corresponds to the number of ASTER elevation measurements available for every pixel from year 2000 to 2012 (after excluding elevation values outside the following accepted elevation range: within 150 m from the median of all DEMs on ice, and within 100 m off ice).

### REFERENCES

Supplementary Figure S2 | Comparison of three NPI glacier inventories. (a) The Blue line corresponds to the inventory created and used in this study for year 2000 (see section Glacier Outlines, Water Bodies, and Ice Divides), ice divides are obtained from Mouginot and Rignot (2015). The red and green lines correspond to the RGI 5.0 and Rivera et al. (2007) inventories for year 2001 (recently incorporated in RGI 6.0), respectively. Named glaciers are those showing the largest area differences between inventories. (b) Zoom over the Cachet glacier area where inventories largely disagree on ice divides. Landsat Image for the 8 of mars 2000 as background. (c) Table shows the differences in glacier area between the inventories from Rivera et al. (2007) and the RGI 5.0 compared to the one used in this study, for the 7 glaciers that showed the largest differences. Note that the three largest glaciers of the NPI, San Quintin, San Rafael and Steffens, covering together 51% of the NPI area are included in this table.

Supplementary Figure S3 | Spatial distribution of the minimum (a) and maximum (b) dates of ASTER DEMs used in our study. White gaps correspond to water bodies and locations where the uncertainty of dh/dtASTER\_trend is larger than 3 m yr−1 (at the 95% confidence level).

Supplementary Figure S4 | Comparison of SRTM, ASTER-14DMO and ASTER-ASP DEMs along a transverse profile of San Quintin glacier. (a) Location of the profile including glacier and stable terrain; transitions are presented as blue lines. (b) Elevation profiles extracted from the SRTM (February 2000) ASTER-14DMO and ASTER-ASP (18 March 2012) DEMs. Both ASTER DEMs are derived from the same Level 1A stereo pair. The ASTER-14DMO DEM (used in Willis et al., 2012a) was downloaded from LPDAAC. This profile illustrates the occurrence of numerous artifacts over the glacier surface in the ASTER-14DMO DEM. Conversely, the ASTER-ASP DEM is smoother and follows nicely the undulations of the SRTM DEM with an offset due to glacier thinning between 2000 and 2012. Hence, ASTER-ASP DEMs are more precise and should be preferred to ASTER-DMO DEMs to study glacier elevation changes.

Supplementary Table S1 | Mass balance of NPI individual glaciers from 2000 to 2012 and percentage area measured for SPOT-SRTM and ASTER\_trend methods.

elevation change and velocity products. Remote Sens. Environ. 162, 408–426. doi: 10.1016/j.rse.2013.07.043

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Dussaillant, Berthier and Brun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Cirque Glacier on South Georgia Shows Centennial Variability over the Last 7000 Years

Lea T. Oppedal 1,2 \*, Jostein Bakke1,2, Øyvind Paasche2 , Johannes P. Werner 1,2 and Willem G. M. van der Bilt 1,2

*1 Department of Earth Science, University of Bergen, Bergen, Norway, 2 Bjerknes Centre for Climate Research, Bergen, Norway*

A 7000 year-long cirque glacier reconstruction from South Georgia, based on detailed analysis of fine-grained sediments deposited downstream in a bog and a lake, suggests continued presence during most of the Holocene. Glacier activity is inferred from various sedimentary properties including magnetic susceptibility (MS), dry bulk density (DBD), loss-on-ignition (LOI) and geochemical elements (XRF), and tallied to a set of terminal moraines. The two independently dated sediment records document concurring events of enhanced glacigenic sediment influx to the bog and lake, whereas the upstream moraines afford the opportunity to calculate past Equilibrium Line Altitudes (ELA) which has varied in the order of 70 m altitude. Combined, the records provide new evidence of cirque glacier fluctuations on South Georgia. Based on the onset of peat formation, the study site was deglaciated prior to 9900 ± 250 years ago when Neumayer tidewater glacier retreated up-fjord. Changes in the lake and bog sediment properties indicate that the cirque glacier was close to its maximum Holocene extent between 7200 ± 400 and 4800 ± 200 cal BP, 2700 ± 150 and 2000 ± 200 cal BP, 500 ± 150–300 ± 100 cal BP, and in the Twentieth century (likely 1930s). The glacier fluctuations are largely in-phase with reconstructed Patagonian glaciers, implying that they respond to centennial climate variability possibly connected to corresponding modulations of the Southern Westerly Winds.

Keywords: glacier, distal glacier-fed lake, distal glacier-fed bog, glacial sediments, holocene, neoglacial, Southern Ocean, sub-Antarctica

### INTRODUCTION

The accelerated melting of glaciers in the Southern Ocean (Cook et al., 2010; Scambos et al., 2014; Koppes et al., 2015) manifests rapid climate change over the last decades (Gille, 2002, 2008; Vaughan et al., 2003; Böning et al., 2008). Although anthropogenic forcing probably explains parts of this regional-scaled glacier retreat, strong regional feedback mechanisms still predominate (Thompson et al., 2011; Wang and Cai, 2013). Many of these mechanisms can be traced back to the Southern Westerly Winds, which are closely linked to the Southern Ocean's uptake of CO2 through their effect on upwelling of carbon-rich deep water, and therefore have the potential to affect climate globally (Lenton and Matear, 2007; Le Quéré et al., 2007; Hodgson and Sime, 2010; Sigman et al., 2010). At present, the strength and position of the westerlies are rapidly shifting, but the processes behind these changes remain poorly understood, due to scarcity and of

#### Edited by:

*Matthias Holger Braun, University of Erlangen-Nuremberg, Germany*

#### Reviewed by:

*Bethan Joan Davies, Royal Holloway, University of London, United Kingdom Olga Nikolaevna Solomina, Institute of Geography (RAS), Russia*

> \*Correspondence: *Lea T. Oppedal<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Quaternary Science, Geomorphology and Paleoenvironment, a section of the journal Frontiers in Earth Science*

> Received: *29 September 2017* Accepted: *11 January 2018* Published: *13 February 2018*

#### Citation:

*Oppedal LT, Bakke J, Paasche Ø, Werner JP and van der Bilt WGM (2018) Cirque Glacier on South Georgia Shows Centennial Variability over the Last 7000 Years. Front. Earth Sci. 6:2. doi: 10.3389/feart.2018.00002*

**20**

instrumental time series and paleoclimatic reconstructions (Hodgson and Sime, 2010), including glacier records which can help map out temporal and spatial changes during the Holocene.

Glaciers are established proxies of climate change (e.g., Oerlemans, 2005) and ubiquitous on Southern Ocean islands, but few reconstructions currently exists. Moreover, most work on Holocene glacier variability in the Southern Ocean is restricted to moraine chronologies (Solomina et al., 2015 and references therein). Mapping and dating of moraines provide valuable information of a glacier's extent, size, and position at points in time. However, moraine chronologies do not capture retreat phases or length of events, and cannot exclude the possibility that moraines have been erased by subsequent glacier advances of equal or larger magnitude (Balco, 2009).

The study of sediments deposited in distal glacier-fed lakes provide the opportunity to continuously reconstruct up-valley glacier variability, although not without caveats. This approach (Karlén, 1976; Bakke et al., 2010) is based on quantification of the sediment delivery by glacial erosion beneath temperate glaciers and its relationship to glacier size (Herman et al., 2015), and assumes that para- and extra-glacial material can either be ignored or identified and removed from the final glacier reconstruction, which is not always the case. Either way, larger glaciers produce more erosional products available for downstream transport by meltwater. These fine-grained glacigenic sediments settle out of suspension in downstream lakes and thus builds unremitting archives of variations in glacier size. Here we take advantage of this lake approach in an effort to produce a new and reliable glacier reconstruction from South Georgia.

The Late Pleistocene-Holocene glacier history of South Georgia has previously been reconstructed in some detail (e.g., Clapperton et al., 1989a; Gordon and Timmis, 1992; Bentley et al., 2007; Roberts et al., 2010) suggesting large Neoglacial advances. At present, only two continuous reconstructions of past glacier change exist: Rosqvist and Schuber (2003) provides a sediment record from Block lake (Husvik, **Figure 1**) based on gray-scale density, loss-on-ignition and grain-size spanning from 7400 to 1000 cal BP (**Figure 9C**), while Van Der Bilt et al. (2017) reconstruct the last 1250 years with a record from the Middle Hamberg lake (south of Grytviken, **Figure 1**) based on titanium counts validated against sediment density and grain size (**Figure 9B**). Both studies suggest centennial-scaled variability, although we do not know whether these records are representative for South Georgia. In order to enhance our understanding of Holocene glacier fluctuations on the island, we need to test if these records can be replicated.

In this study, we present a continuous reconstruction of Diamond glacier (unofficially named) on north-central South Georgia (**Figure 1**) covering the last 7200 years. Our glacier reconstruction is based upon measurements of physical, geochemical, and magnetic properties (cf. Bakke et al., 2010) of sediments deposited in the distal glacier-fed Diamond lake as well as the adjacent Diamond bog (both unofficially named). Mapped marginal moraines in front of Diamond glacier provide an independent framework for the past glacier fluctuations, although the timing of these moraines remains unresolved. Together, the lake-, bog- and moraine records demonstrate a consistent Holocene glacial history of Diamond glacier, which compares well with reconstructed glaciers on South Georgia and in Patagonia.

## SETTING

The sub-Antarctic island of South Georgia is about 170 km long and up to 40 km wide, with extensive glacier cover and an axial, alpine mountain range where mountains rise to nearly 3,000 m (**Figure 1**). Regional climate is influenced by the position and intensity of the westerlies, the position of the Polar Front Zone, and the Antarctic winter sea ice limit, which together mediate the temperature and moisture supply to the island (Pendlebury and Barnes-Keoghan, 2007; **Figure 1**). An instrumental temperature record exists for King Edward Point (KEP, c. 2 m altitude), Grytviken (**Figure 2A**). For the period 1905–1982, mean annual temperature was 2.0◦C and mean summer (December–Feburary) temperature was 4.9◦C. Mean annual precipitation was c. 1,500 mm and mean winter precipitation was 144 mm over the same time interval.

Due to a significant leeside effect on the island, glaciation on the northeast coast is restricted to small cirque and plateau glaciers, while tidewater and calving valley and outlet glaciers occupy the valleys and fjords (**Figure 1**; Gordon et al., 2008). According to observations, temperature is the dominant driver of glacier mass balance (Gordon and Timmis, 1992). Cirque and small valley glaciers respond most sensitively to shifts in climate, and have been receding in response to regional warming since the 1950s (Gordon and Timmis, 1992). With a delayed response, the larger valley and tidewater glaciers followed suit after the early 1970s (Gordon and Timmis, 1992), and now more than 97 % of the glaciers on South Georgia are retreating (Cook et al., 2010).

Olsen Valley in Cumberland West Bay, (**Figure 2**; Supplementary Figure 1) offers a favorable setting for reconstructing past cirque glacier variability using sediments deposited in distal glacier-fed Diamond lake and Diamond bog in combination with a sequence of mapped marginal moraines in front of Diamond glacier. Diamond glacier is located in a cirque southwest of the summit Diamond Peak (∼510 m altitude) between c. 360 and 500 m altitude (2011) (54.20◦ S 36.64◦W).  In the lower valley, marginal moraines mark the former glacier front positions of Neumayer tidewater glacier when it advanced into the valley from the west (cf. Bentley et al., 2007).

In the upper reaches of the valley, sediment-laden meltwater, originating from Diamond glacier, is funneled through a welldefined channel before feeding into the gently sloping valley below. The meltwater follows a course of about six kilometers down Olsen Valley before entering Cumberland West Bay (**Figure 2B**). The stream is intersected by an alluvial fan in the central part of the valley, containing both fresh sediments and sediments covered by vegetation, indicating deposition of glacial sediments both today and in the past. The meltwater has cut deep meanders (up to 4 m) into the bog covering the valley floor in the

FIGURE 1 | (A) Overview map showing South Georgia in the Southern Ocean with major circulation patterns such as the Antarctic Polar Front Zone (APFZ), the Southern Westerly Winds (SWW), and the Antarctic Circumpolar Current (ACC), as well as the Antarctic Circle and winter sea ice limit after Scar (2015) and Orsi et al. (1995). (B) Inset of South Georgia showing the current glacier extent after SGGIS (2017) and our study area.

lower valley. Sediment sections along the river show alternating layers of peat and fine-grained sediments in the upper meter, indicating deposition of glacially derived sediments when the river floods its plains (Matthews et al., 2005).

Further downstream, the stream enters another alluvial fan separated from Diamond Lake by a marsh. The meltwater drains into the alluvial fan, which during times of moderate to high runoff has standing water. If the water level increases above c. 0.5 m, water overflows the marsh and drains into Diamond Lake. We find it likely that this occurs every year during the ablation season. The lake (54.24◦ S 36.65◦W, 16 m altitude) covers around 10,000 m2 and is ∼3 m deep. Its catchment covers c.11 km2 , with 3% of that surface currently occupied by the Diamond glacier. Bedrock, consisting of andesitic volcaniclastic greywackes (Tanner, 1981 in Quilty, 2007), is exposed in the steeper part of the valley, which is free of vegetation, and unprotected paraglacial sediments are present in an area of c. 2 km2 downstream of the glacier. The gentle slopes of the lower valley, below c. 150 m altitude, are covered by vegetation which protects the lake from mass-wasting processes. The grass Festuca contracta covers the slightly raised areas in the valley, while vegetation on the bog is dominated by the moss Tortula robusta, with grass and flowering plants present such as Acaena magellanica and Rostkovia magellanica. Following the classification of peat on South Georgia by Smith (1981), Diamond Bog is a soligenous eutrophic mire. It is developed on the floor of a broad valley and traversed by drainage channels which increases runoff and lowers the nutrient levels of the peat. At the coring site, Diamond Bog is 220 cm.

### MATERIALS AND METHODS

Our reconstruction of glacier fluctuations at the Diamond cirque glacier is based on a multi-proxy approach (Bakke et al., 2010), integrating geomorphological field mapping, and physical, magnetic and geochemical analyses on sediments deposited in the distal glacier-fed lake and bog.

#### Geomorphological Mapping

We mapped the geomorphology of the Diamond catchment (**Figure 2B**) based on ESA Sentinal satellite imagery (Google, n.d.)1 combined with field surveys during January 2008 and December 2011. Our main goal was to identify marginal moraines marking former glacier coverage, and to determine past and present pathways of glacial meltwater. The topographic maps of South Georgia were, at the time of our field work, not of a resolution suitable for detailed mapping and positioning; we therefore used a hand-held GPS device (Garmin CSX60) for sites and sampling locations. The average off-set on the positions are in the order of 2–10 m.

#### Lake Survey and Sediment Coring

We mapped the bathymetry of Diamond lake (**Figure 3**) using a Garmin Fishfinder 160x echo sounder combined with a GPS and retrieved the sediment cores from the deepest and most central part of the lake (**Figure 3**, **Table 1**). For the longer piston cores (LDP1, 75 cm and LDP2, 246 cm) we used a modified piston corer able to capture up to six meters of sediments. The surface sediments at each of the two coring sites (LDS1, 34 cm and LDS2, 46 cm), we extracted using a gravity corer (Renberg and Hansson, 2008). In addition, we retrieved three sediment cores from a section of Diamond bog, upstream of the lake, by hammering PVC tubes into the bog (**Figure 2B**, **Table 1**).

### Sample Treatment and Laboratory Analyses

Prior to sub-sampling, we visually logged the lithological characteristics of the sediment cores, and used an ITRAX Xray fluorescence (XRF) core scanner to map out the sediment geochemistry (c.f. Kylander et al., 2011). Together with the geochemical mapping at 200-µm resolution, we obtained optical and radiographic images at the same time (**Figure 6**). The XRF analyses were made using a molybdenum tube with power settings of 30 kV and 35 mA, employing a counting time of 10 s

LDP2 in the deepest part of the lake. Gravity cores LDS1 and LDS2 were taken at the same two locations.

TABLE 1 | Overview of sediment cores.

(NLDP1 = 3179, NLDP2 = 12280, NDB1 = 6771, NDB2 = 10201, NDB3 = 3551).

We scanned the cores for surface Magnetic Susceptibility (MS) at 0.5 cm intervals on a Bartington MS2 point sensor (range = 0.1) (NLDP1 = 132, NLDP2 = 492, NDB1 = 277, NDB2 = 297, NDB3 = 135). Since the gravity cores (LDS1 and LDS2) had already been sub-sampled in the field (every 0.25 cm and 0.5 cm, respectively), we measured the bulk magnetic susceptibility of these samples on a MFK1-FA (sensitivity: 2 × 10−8 SI). We use the MS data to tentatively correlate the four cores from Diamond lake (**Figure 5**). Being the two longest cores from, respectively, the lake and the bog, we selected core LDP2 and core DB2 for further and more detailed analyses.

We followed the protocol for Loss-On-Ignition (LOI) after Dean (1974) and dry-bulk density (DBD) after Bakke et al. (2005) of 1 cm3 dry sub-samples every 0.5 cm (NLDP2 = 492, NDB2 = 408). Using an MFK1-FA (sensitivity: 2 × 10−8 SI), we analyzed the mass specific magnetic susceptibility (χbulk)

1Google (n.d.). Google Map of South Georgia [Online]. Available online at: https:// goo.gl/Dr1pYP (Accessed 2008).

(Thompson and Oldfield, 1986) of discrete samples (6.8 cm3 ) at 1 cm intervals (nLDP2 = 238, nDB2 = 204). Prior to measurements, we freeze-dried the samples to remove the water and thereby reduce chances of post-sampling diagenesis. With the aim of distinguishing between para- and ferromagnetic sediments, we conducted the measurements at both room temperature (293 K) and 77 K (liquid nitrogen), and calculated the paramagnetic ratio (χbulk77K/χbulk293K) for each sample (Lanci and Lowrie, 1997).

The grain size distribution of 20 samples from LDP2 and 5 samples from DB2 were analyzed using a Mastersizer 3000 with an LV dispersion unit. All samples were pre-treated for c. 24 h with 35% hydrogen peroxide to remove the organic material (Vasskog et al., 2016) and sonicated with ultrasound for 60 s. prior to measurement. Bulk sample sizes were adjusted to have an obscuration rate of 10–20%, and measured for 10 s (blue laser) and 20 s (red laser) using refraction index of 2.5, absorption index of 0.01, and stirring rate of 2,500 rpm.

Age-depth relationships were established for the lakeand bog cores with the help of 24 individual samples of terrestrial macrofossils collected from LDP2 and DB2. Samples were submitted to Poznan Radiocarbon Laboratory for AMS radiocarbon dating (**Table 2**). Based on these dates, we constructed age-depth models for each core with the software package Clam 2.2 (Blaauw, 2010). We applied statistical crossmatching techniques to further constrain the age models, and to demonstrate correspondence between the three sediment cores (section Age-Depth Modeling).

#### Equilibrium Line Altitude Calculations

We constrained former glacier extents combining satellite imagery (present) with fairly intact moraines, which were mapped during the two fieldwork seasons (2008 and 2011). We used ArcMap v. 10.4 as the platform to apply the cartographic method (e.g., Ballantyne, 1989) for reconstruction of past glacier extents and surface areas. Based on this, we estimated the present and past Equilibrium Line Altitudes (ELA) of Diamond glacier (**Figure 2B**) using the Accumulation Area Ratio method (AAR) (Porter, 1975) with an accumulation area of 0.6 ± 0.05, which is the standard AAR value used for cirque glaciers. In addition, we used the Area x Altitude method and the Area x Altitude Balance-Ratio method (AABR) (Osmaston, 2005) with a range of balance ratios (BR) (1.5, 1.9, and 2.5) since little is known about the balance ratio in this region. The ELA reported for each moraine is the mean value of the different estimates.

### RESULTS

### Moraines

We mapped nine terminal moraine ridges (M1-M9) in front of Diamond glacier (**Figure 2B**), which we use to reconstruct former glacier extent and calculate the corresponding ELAs (**Figure 4**; Supplementary Table 1). The moraines are free of vegetation, and none of them is ice-cored. We group them into four clusters described in the following.

The outermost moraine cluster (M9-M8) consists of two distinct free-standing continuous ridges (c. 2–3 m high and 300 m) located, c. 40 m apart, at a shelf (c. 170 m altitude) right

*Radiocarbon dates obtained for the Lake Diamond Piston core (LDP2) and Diamond Bog core (DB2).*

before the valley drops steeply into the gently sloping lower Olsen Valley. The ridges are located at both sides of the present-day meltwater stream emerging from Diamond glacier. They consist of diamictic sediments, ranging from clay to gravel. Scattered boulders (10–30 cm in diameter) are partly embedded within or resting on the ridge surfaces. Cross sections of the moraines show an asymmetric profile with a steeper distal slope (c. 30◦ ) and a low-angle proximal slope (c. 10◦ ). The area distal to the ridges has exposed bedrock on both sides of the meltwater stream, which cuts through the moraines in a south westerly direction. Both ridges have undulating long profiles, which indicates partial erosion by meltwater, implying that the drainage pattern was more directly oriented toward the south when the ice front was closer to the terminal moraines. We cannot ascertain the processes of moraine formation, as no exposure in the ridge was present. At the time of moraine formation, we estimate that the ELAs were 64 ± 13 and 70 ± 12 m lower than in 2011, respectively for M8 and M9 (**Figure 4**; Supplementary Table 1).

The intermediate moraine cluster consists of two freestanding ridges (M7-M6), varying from 1 to 3 m in height, that continue for more than 600 m and are only cut by the river from Diamond glacier confined to a canyon (with a depth between 3 and 10 m) distal from the moraine, which funnels the meltwater toward the lower valley. Proximal to the two ridges is another partly preserved ridge (M5) that continue for approximately 550 m. Like the outermost moraines, ridges M7-M5 are composed of

diamictic sediments with grain sizes ranging from clay to gravel, but with fewer boulders embedded and on the ridge surfaces. The cross profile of the ridges shows a steep distal side (c. 36◦ ) and a more gently sloping proximal side (c. 15◦ ). The similarities with the outermost moraine system is striking. At the time M7 and M6 were formed, ELA was 61 ± 13 m lower than in 2011, while the corresponding ELA for M5 is 47 ± 14 m (**Figure 4**; Supplementary Table 1).

The inner moraine cluster is located less than 50 meters proximal to ridge M5. However, the geomorphology of this moraine system is different with numerous small ridges located very close to each other. One of the ridges (M4) is freestanding and slightly higher than the others (4 m high). It can be followed over a distance of 700 m, only crosscut by the present-day meltwater stream. On the north side of M4 more than 15 small (<1 m high) saw-tooth shaped ridges can be mapped over a distance of less than 50 meters (M3). The diamicton in the moraines consists of fine sand, gravel and boulders up to 1 m in diameter. The lack of clay and fine silt indicate presence of water during deposition. Our interpretation is that the moraine cluster nested to M4, including the many smaller moraine ridges, represent an oscillatory glacier front with moraines deposited rapidly by a shrinking glacier. At the time of moraine formation, ELA was 46–47 ± 14 m lower than in 2011 (**Figure 4**; Supplementary Table 1).

The innermost moraine cluster consists of more than 30 individual moraine ridges deposited in a 200-meter-wide area around 320 m altitude. The two most pronounced ridges are marked M2 and M1 on the map. The ridges are between 0.5 and 1.5 meters high and are strongly asymmetrical, with a much steeper (35◦ ) distal face and gentler (20◦ ) proximal face. The geomorphology of the moraines is very similar to what Lukas (2012) describes as annual moraines formed when the glacier advances and acts as a bulldozer on the sediments in front. Such moraines are observed in front of several mountain glaciers on the peninsulas of north eastern South Georgia, documenting glacier retreat during the past century (Gordon et al., 2008). During the time of moraine formation, ELA was 42 ± 15 and 38 ± 15 m lower than in 2011 (**Figure 4**; Supplementary Table 1).

Based on GPS measurements of the glacier front during our two field campaigns, the glacier front retreated approximately 150 m between January 2008 and December 2011. This suggests that Diamond glacier follows the same pattern of rapid glacier retreat observed on South Georgia, in the sub-Antarctic, in Patagonia and indeed, globally (e.g., Cook et al., 2010; Koppes et al., 2015; Zemp et al., 2015).

In the lower valley, we mapped three moraines, covered by grasses and herbs, which outline past glacier front positions of Neumayer tidewater glacier, which advanced into the valley from the west. Parallel to Cumberland West Bay, a prominent lateral moraine (OM3), 1–4 m in height, sloping from 60 to 45 m altitude is deposited superimposed on a 100 m-high and 1 km long bedrock ridge. About 1 km further north, several cones consisting of diamict material rising above the bog likely represent remnants of a marginal moraine system (OM2). A pronounced ridge along the mountainside between Olsen Valley and the next tributary valley to the west marks a younger glacier advance (OM1). This ridge is between 2 and 6 m high and the material in the ridge is a diamicton of fine sand and boulders up to 40 cm in diameter.

#### Lake Sediments

Based on the correlation of the MS stratigraphy of the sediment cores from Diamond lake, (**Figure 5**), we argue that core LDP2 is representative for the lake. Piston core equipment sometimes fails to capture the sediment-water interface, and correlation between LDP2 and the gravity core LDS2 (taken from the same location) suggests that this is the case for LDP2. Possible overlaps between LDP2 and LDS2 were tested using a Bayesian cross matching algorithm by Werner and Tingley (2015) (Supplementary Material; Section Age-Depth Modeling), which suggests that 10 cm of the topmost sediments of LDP2 were likely lost during coring. We sub-divide LDP2 into six

units (A–F) based on visual logging of the core (**Figure 6A**), and further divides unit F into five sub units (F1–F5) based on visual logging, grain size measurements and χbulk. Below, we give a description of each unit along with the sediment parameters LOI, DBD, χbulk, and Ti XRF count rate which we use to detect changes in the minerogenic vs. organic content in the core. The ratio of χbulk measured at 77 K and 293 K (χbulk77K/χbulk293K) is close to 3.8 throughout the core (**Figure 6A**), indicating that paramagnetic minerals dominate the samples (Lanci and Lowrie, 1997). Most values are in fact higher than 3.8, which is theoretically impossible, but can be explained. We ascribe these high values to the fact that actual temperatures of the samples when measured deviated from either 77 or 293 K, as the room might have been too cold. Samples of weak magnetic susceptibility are more prone to this error.

#### Unit A

Unit A (243.5–246 cm) is a matrix-supported diamict, consisting of sub-angular clasts with no apparent orientation embedded in a well consolidated very dark gray matrix of clay to very fine sand (**Figure 9**). The clasts range in size from c. 2–20 mm and include one cobble at the bottom of the unit, measuring c. 80 mm across. The unit is characterized by high DBD, χbulk and Ti values [medians of 0.5 g/cm3 , 9.0 × 10−8 m3 kg−1 and 1,030 counts per second (cps)], reflecting its high minerogenic content. Based on the poor sorting and the sub-angular clasts, we interpret unit A as glacial till deposited during the last glacier retreat from the coring site.

#### Unit B

Unit A grades into unit B (223–243.5 cm), a dense, structureless and highly decomposed very dark brown peat with one laminae of sandy silt at 230 cm. The high organic content of the unit is reflected by high LOI values (median of 80 %) corresponding with low values of DBD, χbulk, and Ti (medians of 0.1 g/cm3 , 2.8 × 10−8 m3 kg−1 and 80 cps). We interpret unit A to represent a peat bog environment with little (glacio-) fluvial influence.

#### Unit C

A sharp erosive bedding plane separates unit B from unit C (221.5–223 cm), a very thin bed of dark gray moderately sorted sand, with sub-rounded grains ranging from c. 1–5 mm, and containing one wood fragment (10 mm long, 5 mm in diameter). The unit is best resolved in the high-resolution XRF counts, with high Ti counts around 1,000 cps, reflecting the minerogenic nature of the sediments. Based on the sediment texture and wood fragment, we interpret the unit to represent a flooding event, causing erosion into unit B. According to our age-depth model (**Figure 7**), the resulting unconformity represents a hiatus of ∼2000 years.

#### Unit D

Unit C has a sharp bedding contact with Unit D (194–221.5 cm), which consists of a dark grayish brown organic-rich silt with grain sizes ranging from clay to medium sand (**Figure 9**). The sediments are well consolidated (median DBD of 0.6 g/cm3 ), and small plant fragments are dispersed throughout the unit (median LOI of 11 %). The high minerogenic content is reflected by high median values of DBD, χbulk (7.9 × 10−8 m3 kg−1 ) and Ti (930 cps) (**Figure 5**). We interpret unit C to represent glacial flour (clay and silt) and reworked sediments (sand) deposited in a lake.

#### Unit E

A sharp undulating non-erosive bedding contact separates unit D from unit E (185.5–194 cm), which is composed of a poorly sorted very coarse dark gray loosely consolidated sand containing sub-angular pebbles up to 30 mm across and no silt or clay particles. The high minerogenic content is reflected by low LOI values (median of 2%) and high values of DBD (median of 1.1 g/cm3 ), Ti (median of 960 cps) and χbulk (one bulk

measurement of 4.0 × 10−7 m3 kg−1 ) (**Figure 5**). We interpret unit E as a slump, formed by mass wasting of sediments from a moraine.

#### Unit F

Unit F (0–185.5 cm) consists of dark brown to dark grayish brown organic-rich sediments with dominant grain sizes ranging from clay to fine sand. We interpret unit F to represent a lake similar as today, into which the meltwater stream carried glacially eroded silt and clay directly from the glacier, but also remobilized coarser grain fractions as well as terrestrial plants in the catchment during times of high runoff. We argue that the input of minerogenic sediments to the lake is positively correlated with glacier size (Section Reconstructing the Glacier). Based on variations in color, minerogenic content, χbulk, grain size distribution, and organic content, we divide unit F into five sub units (F1–F5).

Unit F1 (185.5–182 cm) is a dark grayish brown organicrich silty sand displaying similar characteristics as unit D, but is dominated by sand, rather than silt. It was likely deposited under similar environmental conditions as unit D, but also affected by the preceding mass wasting event (Unit E). One grain size

constrained using Bayesian statistics. Bands in a lighter red and pink shade marks the 50 and 95% confidence intervals, respectively.

measurement at 182 cm shows that very fine sand dominates, and median χbulk is 8.1 × 10−8 m3 kg−1 .

Unit F2 (182–128 cm) consists of consolidated organic-rich sandy silt with a dark grayish brown color, interspersed by lamination of well sorted very fine sand at 181, 153, 148, 146, 137, and 134 cm. Grain size measurements reveal a bimodal distribution with peaks around very fine sand and medium silt (**Figure 9**). The finer fractions indicate input of glacial sediments while the sand is likely remobilized paraglacial sediments from the fan next to the lake. Throughout the unit, χbulk generally decreases. Values are low (median of 6.4 × 10−8 m3 kg−1 ), but increases over the sand layers (**Figure 6A**).

Unit F3 (128–58 cm) is a dark brown organic-rich silty sand with high content of moderately decomposed plant fragments. The grain size distribution of four samples (118, 88, 70, and 65 cm) peaks around very fine sand, but silt is also present. χbulk generally increases through the unit (median of 8.0 x 10−8 m3 kg−1 ). Between 71 and 75 cm, large terrestrial plant fragments (likely the herb Aceana magellanica/tenera) are found at the center of the core, probably representing a flooding event, with the high organic content of the layer reflected in a drop in the DBD and χbulk values (**Figure 6A**). Part of this organic material fell out when splitting the core into two half tubes. The resulting hole caused very low XRF counts (190 cps) as the detector was not close enough to the sediment surface (**Figure 6A**). Terrestrial plant fragments are found at several depths in the unit (e.g., 86, 97, 110, and 120 cm). The terrestrial plant remains and coarse grain sizes indicate a high energy system causing remobilization of sediments and terrestrial organic material in the catchment, and subsequent deposition into the lake. While sand settles in the lake, the finer fractions are transported all the way into Cumberland West Bay.

Unit F4 (58–10 cm) is dark grayish brown organic-rich silt. Grain size measurements at 55, 33, 24 and 14 cm show dominating clay and silt fractions. At 44 and 48 cm, fine sand dominates but the finer fractions are also present. In the lower 24 cm χbulk values are comparable to unit D and F1, but starts to decrease from 5.6 × 10−8 40 cm, reaching around 13 cm with χbulk (5.6 × 10−8 m3 kg−1 ). We interpret this unit to represent lower energy conditions than in unit F3. The presence of clay and silt indicate input of glacial sediments. According to χbulk glacial input is high, comparable to that of unit D, between c. 58 and 38 cm, but decreases through the rest of the unit. Deposition of remobilized sand occurs at some intervals, indicating events or short periods of higher run-off.

Unit F5 (10–0 cm), consists of dark brown organic-rich silty sand (**Figure 10**) with a high content of woody plant fragments, increasing toward the top. The grain size distribution of one sample at 3 cm peaks around very fine sand, similar to unit F3, but the silt (and clay) content is a little higher. χbulk increases through the unit from 5.6 × 10−8 to 7.3 × 10−8 m3 kg−1 . We interpret unit F5 to represent a system of increasingly high runoff, where both clay and silt from the glacier as well as remobilized sediments and terrestrial organic material from the catchment are deposited in the lake.

### Bog Sediments

During the coring of the DB2, we measured compaction due to the friction in the corer to be ∼16 cm, but since we know little about how the compaction was distributed through the core, we do not attempt to correct for it. We sub-divide DB2 into two units (A and B) based on visual inspection of the core (**Figure 6B**).

#### Unit A

Unit A (95–204 cm) is composed of dense, highly decomposed black peat with laminations of clayey silt at 100 and 155 cm. The high organic content is reflected in the high median value of LOI (87 %) and low median values of DBD, χbulk and Ti (0.2 g/cm3 , 1.0 × 10−8 m3 kg−1 and 90 cps) (**Figure 6B**). We interpret unit A to represent a peat bog environment with little influence of glacial meltwater.

#### Unit B

Unit A grades into unit B (0–95 cm), which consists of very dark brown clayey peat interbedded with gray to grayish brown layers of fine minerogenic material and fragments of moss and rootlets. Bedding planes are indistinct, and bed thicknesses vary from <1 mm to 5 cm, with three 3–5 cm thick layers in the upper 25 cm of the unit. Using the XRF count ratio Fe/Ti, we identified 33 minerogenic layers in unit B, which we confirmed by inspection of the sediments. The grain size distribution of four minerogenic layers at 3–5, 6.5–10, 18.5–23, and 73–74 cm, reveal a domination by silt and clay. At the bottom of the unit (94 cm) grain sizes are dominated by very fine to fine sand in addition to a high content of silt and clay. The organic content of the unit is moderately to poorly decomposed, with the degree of decomposition decreasing toward the top of the core which preserves the topsoil. Compared to unit A, the measured sediment parameters show higher variability in Unit B, reflecting the variable content of minerogenic vs. organic material. LOI varies between 4 and 88%, DBD varies between 0.2 g/cm3 and 1.1 g/cm3 , χbulk varies between 3.0 × 10−9 m3 kg−1 and 5.8 × 10−8 m3 kg−1 , and Ti varies between 210 cps and 1460 cps. The paramagnetic ratio is close to 3.8 throughout the unit (**Figures 6B**, **10**), indicating one main sediment source (Section Lake Sediments). This, together with the grain size distribution of the minerogenic layers and the transport of glacial meltwater through the bog in a deep channel, suggests that the minerogenic layers in the bog are composed of glacial sediments. We interpret unit B as a peat bog environment such as today, where glacially eroded sediments are deposited onto the bog when the river overflows its banks. The high content of sand in the bottom of the unit concurs in time with the increased sand input to the lake (Section Reconstructing the Glacier) and infers higher energy in the system.

### Age-Depth Modeling

We established age-depth relationships for LDP2, LDS2, and DB2 (**Figure 7**) in the following steps (for more details on the approach see Supplementary Material):

4) From cross-matching LDP2 and LDS2 by eye (**Figure 5**), we argue that the top sediments of LDP2 are missing. In order to extend LDP2 to the present, we constructed an age model for the undated core LDS2 based solely on the sedimentation rates inferred for LDP2. Using the aforementioned algorithm, we tried to reconcile the possible age models for LDP2/LDS2 (Supplementary Figures 5, 6).

## DISCUSSION

### Reconstructing the Glacier

There is a long tradition for utilizing downstream basins as sedimentary archives for reconstructing glaciers dating back to pioneering work by Karlén (1981). Since these early studies, progress have been made regarding the methodological toolkit invoked in order to track past glacier activity (e.g., Bakke et al., 2010) as well as tools for precise modeling of the age-depth relationship based on 14C dates, 210Pb dates and tephra (Blaauw, 2010) ensuring robust chronologies. One of the challenges with lake-based glacier reconstructions is the necessity of recognizing various catchment processes that potentially can impact downstream sedimentary records (Rubensdotter and Rosqvist, 2005; Bakke and Paasche, 2011). To disentangle the different processes acting on the archives, it is increasingly common to employ multiple proxies (e.g., Bakke et al., 2010; Vasskog et al., 2012). Sediment variables, commonly reflecting the glacial signal include elements such as Si, K, Ca, Ti, Fe, Rb, and Zr (Davies et al., 2015), variations in bulk organic matter as given by loss-on-ignition (e.g., Karlén, 1976), DBD (e.g., Bakke et al., 2005) and grain size (e.g., Matthews and Karlén, 1992), as well as a range of magnetic properties (e.g., Paasche et al., 2007). In this study, we use mass specific magnetic susceptibility (χbulk) as an indicator of change in the suspended glacigenic input to the bog and lake. Magnetic susceptibility is a measure of the magnetic minerals in a sample which is closely related to the concentration of minerogenic material (Thompson et al., 1975; Thompson and Oldfield, 1986). The approach of using χbulk as a proxy for glacier variations is validated by the positive correlation of χbulk with DBD, surface magnetic susceptibility and common detrital elements such as K, Ca, Ti, V, Mn, Fe, and Sr (Supplementary Table 2).

Olsen Valley represents a relative simple system where a single river drains a small upstream cirque glacier, feeding meltwater through the lower valley from a point source. Given the high ablation season temperatures (mean TOct−Apr = 3.6◦C) we assume that the glacier is temperate and also remained temperate during the period covered by our sediment record (last 7200 years). We neither find geomorphological evidence, such as ice-cored moraines, that the glacier was cold based or polythermal during this period, nor are such previously reported from the island. For the sake of comparison, polytemperate glaciers in Scandinavia and Greenland typically have mean ablation temperatures around 1◦C or colder (Østrem et al., 1988; Weidick, 1992). Such a decrease in temperature is highly unlikely for the Holocene on South Georgia, as indicated by the GDGT temperature record by Foster et al. (2016) and the δD temperature record from the Antarctic Peninsula by Mulvaney et al. (2012). The size variations of temperate glaciers are typically positively correlated to the production and evacuation of glacial sediments (Herman et al., 2015).

We regard the minor minerogenic content in the bog during the late Holocene climate optimum (4800 ± 200 to 2700 ± 150 cal BP, section Glacier retreat: 4800 ± 200 to 2700 ± 150 cal BP), when the glacier is inferred to be very small, as evidence that the glacier is the main driver of the sediment variability in the bog (Matthews et al., 2005). The sediment influx to the bog and lake is predominantly associated with the ablation season (summer) when increased runoff causes the river to transport and subsequently deposit glacial flour in the bog and lake. Based on the low mean summer precipitation (mean PJJA = 105 mm), we argue that the additional extra-glacial sedimentary delivery from the other parts of the catchment is small compared to that which stems from glacial meltwater. There exsists no record of extreme precipitaiton events from South Georgia, however such events can cause inwash of inorganic sedimentation to our sites and is a possible error source for our interpretations. Sitting in a gently sloping valley covered by vegetation, we consider mass wasting processes not likely to affect the bog and lake.

While the gently sloping part of the valley is covered by peatland, thereby limiting the possibility of reworking of older sediments (paraglacial activity), glacigenic deposits closer to the glacier, as well as in three fans deposited in the lower valley, are potentially available for erosion and thus re-mobilization (Ballantyne, 2002). This could potentially be a source of error when using the sediments to reconstruct the up-valley glacier. Still, since sedimentation into the bog and lake happens mainly during the ablation season, we argue that re-deposition of paraglacial sediments would only reinforce the glacier signal rather than distort it. The presence of fine sand in the lake (**Figure 9**) indicates that there is indeed a paraglacial component present. These sediments are likely remobilized from the alluvial fan to the west of the lake during periods of high runoff. This is supported by the domination of clay and silt in the bog, which is located up valley from the fan, suggesting low paraglacial input. Sand in the bog is limited to an interval between 93 and 95 cm (c. 2700–2400 cal BP). This interval concurs in time with increased sand input to unit F3 (2700 ± 150–2000 ± 200), suggesting that both lake and bog received paraglacial sediments from up valley sources (likely the alluvial fans) between 2700 ± 150 and 2000 ± 200 cal BP when runoff was at its highest in the last 7200 years.

Seeing as both the bog and lake receive input of glacigenic material during periods of high river discharge, we claim that the two records should co-vary in phase over time. This assumption is confirmed in **Figure 8**, which shows increased χbulk in LDP2 corresponding in time with the 33 silt layers identified in DB2. This relationship serves to validate the lake record and demonstrates the potential in using distal-glacier fed bogs as archives of glacier variability. Other lines of evidence also suggest that the lake and bog capture the same glacial history. The paramagnetic ratios (χbulk77K/χbulk293K) of LDP2 and DB2 indicate that both records are dominated by paramagnetic minerals (Lanci and Lowrie, 1997) (**Figure 11**; section Lake Sediments) and that this doesn't change during the period in question. This shows that the minerogenic material deposited in the lake and bog derives predominantly from one source only, unless one assumes that the signal is the same for the entire catchment, which is unlikely considering the geological processes involved. Further supporting evidence is the overlapping grain size distributions of the bog and lake. Characteristic of glacial flour, we argue that the source of the minerogenic and finegrained sediments is erosion by the Diamond glacier, which, drawing from moraine evidence, was present in the catchment during the Holocene.

The mapped moraines in front of Diamond glacier provide spatial constraints for the glacier fluctuations reconstructed from the sediments, as we arguably detect the documented glacier advances of the LIA and in the 1930s both in the sediment records and in the moraine sequence. Based on previous studies from South Georgia, we argue that M1 was deposited during the 1930s while M2 was deposited during the Little Ice Age (LIA). An extensive glacier advance during the LIA (seventeenth

to ninetieth centuries) is well-documented on the peninsulas of northeastern South Georgia (Gordon and Timmis, 1992). Clapperton et al. (1989b) estimate an average ELA depression of c. 50 m during this period, corresponding with a cooling of 0.5◦C. Within the uncertainty range, this is comparable to the ELA of M2. Proximal to the LIA moraine it is common to find a moraine deposited during the twentieth century (Gordon and Timmis, 1992). This advance was documented by historical evidence (Hayward, 1983), showing that many glaciers advanced in the twentieth century. Cirque and valley glaciers were at its most advanced position in the 1930s, while larger valley and tidewater glaciers reached their maximum glacier extent in the 1970s. Such a glacier advance is also documented for the Hamberg glacier by Van Der Bilt et al. (2017). Furthermore, during the recession phase after the twentieth century advance, many cirque glaciers deposited annual moraines (Gordon and Timmis, 1992), such as the ones observed in the innermost moraine cluster. Thus, Diamond glacier followed a similar pattern to that observed for small glaciers (0.1–4.0 km2 ) on South Georgia during the late Holocene, with a Little Ice Age advance, a period of recession, a twentieth century advance and a recent recession (Gordon and Timmis, 1992). Based on the lake sediments, we argue that M8 and M9 were deposited prior to lake formation 7200 ± 400 years ago. We assume that the glacier would produce a distinctly higher χbulk in the lake at this time (being 30% larger than the glacier which deposited M7). The lake record, however, indicates that the glacier repeatedly reached close to its maximum Holocene position (**Figure 11**).

#### Holocene Glacier Variability in Olsen Valley

In the following we present the reconstructed Holocene glacier variability in Olsen Valley as interpreted from the lake and bog sediment records and linked to the moraine sequence based on previous studies of glacier front positions on South Georgia (section Reconstructing the Glacier). We use χbulk as a proxy for glacier size, except in the topmost part of the lake record where we use κbulk (see section Sample Treatment and Laboratory Analyses). Since we combine two different measures of magnetic susceptibility in the lake, we standardize the χbulk and κbulk results (**Figure 11A**). To account for glacier response time, we apply 30-year bins in the presentation of the glacier reconstruction following Van Der Bilt et al. (2017).

#### Deglaciation of Lower Olsen Valley: <9900 ± 250 cal BP

Based on the onset of peat accumulation in Diamond lake, lower Olsen valley was ice free prior to 9900 ± 250 cal BP. We argue that the diamicton found in Unit A (LDP2) was deposited by the Neumayer Glacier, as moraines OM2 and OM3 outline a glacier advance from the west (**Figure 2B**) (Bentley et al., 2007). This is supported by the slightly lower values of χbulk77K/χbulk293K in Unit A compared to the rest of the record (**Figure 6A**), indicating a different sediment source. A retreat of Neumayer Glacier from Olsen Valley before 9900 ± 250 cal BP agrees with radiocarbon dates of peat and lake bases in the Cumberland Bay and Stromness Bay area (cf. Van Der Putten and Verbruggen, 2005) suggesting that the lower zone (<50 m a.s.l.) was deglaciated in in the first centuries of the Holocene.

Marine evidence connects OM1 and OM2 to the glacier retreat following the Antarctic Cold Reversal on South Georgia (15 170–13 340 cal BP). Hodgson et al. (2014) document a partially preserved moraine in the outer basin of Cumberland West Bay (**Figure 2A**). This moraine is associated to the Antarctic Cold Reversal advance on South Georgia by Graham et al. (2017) who dated the transition from subglacial to ice-proximal sediments in a core taken through the equivalent moraine in Cumberland East Bay to 13,340 cal BP. Another date from the same sediment core suggests nearby glacier sources at the outer basin moraine until c. 10,600 cal BP (Graham et al., 2017). Based on the above dates, we suggest that OM2 and OM3 were deposited between c. 10,600 and 9900 ± 250 cal BP.

#### Glacier Advance: 7200 ± 400 to 4800 ± 200 cal BP

We have evidence of continuous peat formation and dry land conditions until 9350 cal BP, but the unconformity between unit C and D makes it difficult to establish when the lake formed, as we have no account of what happened between 9350 and 7200 ± 400 cal BP. The period between 7200 ± 400 and 4800 ± 200 cal BP (unit D) is characterized by high input of clay and silt (**Figure 9**). Characteristic of glacial flour, these fine-grained sediments were transported in suspension by glacial meltwater directly from the glacier. We interpret the high χbulk values to represent a glacier close to its maximum Holocene position (**Figure 11A**). In line with our findings, evidence from bryophytes and seeds on nearby Tønsberg Peninsula (**Figure 2A**; Van Der Putten et al., 2004) suggests cold and wet conditions between 7000 and 4000 years BP, which would favor glacier a during this period.

#### Glacier Retreat: 4800 ± 200 to 2700 ± 150 cal BP

From 4800 ± 200 cal BP χbulk decreases and remains low until 2700 ± 150 cal BP (Unit F2, **Figure 6A**, **Figure 11A**). In Diamond Bog median LOI is close to 90% (**Figure 6B**), suggesting that the bog was isolated from meltwater input (**Figure 10**). However, the presence of clay and silt in Diamond Lake (**Figure 9**) suggests that the glacier was present in the catchment, and had not completely melted away. Evidence for a contemporaneous Late Holocene climate optimum on South Georgia has been found in several studies, reporting reduced glacier activity (Rosqvist and Schuber, 2003) higher temperatures (Clapperton et al., 1989a; Strother et al., 2015; Foster et al., 2016) and a drier climate (Van Der Putten et al., 2004, 2009). Paleoclimatic records from other sub-Antarctic locations such as the Antarctic Peninsula (Sterken et al., 2012), Livingstone Island (Björck et al., 1991), South Shetland Islands (Björck et al., 1993), and James Ross Island (Mulvaney et al., 2012) demonstrate consistent climate change over a wide region.

#### Glacier Advance: 2700 ± 150 to 2000 ± 200 cal BP

From 2700 ± 150 cal BP we see a sudden shift in the sedimentary regime of both the bog and the lake. Large amounts of sand and silt suddenly starts to wash into the bog (**Figure 6B**). Sand, silt and terrestrial plants are deposited in the lake (Unit F3, **Figure 6A**), and modeled sedimentation rates increase in both records (**Figure 7**). We interpret this shift to mark a sudden transition from warm and dry climate conditions with restricted glaciation to a wetter and colder period favoring glacier advance. The energy in the hydrological system of the valley is high, as indicated by the coarser grain fractions and the input of terrestrial organic material, and the contribution of precipitation vs. glacier advance to the increased χbulk values are difficult to disentangle. However, based on the sensitive response of present cirque glaciers on South Georgia (Van Der Bilt et al., 2017), we find it likely that Diamond glacier advanced within short time after 2700 ± 150 cal BP. The change in the sediment records coincides with a shift in climatic conditions, which is recorded in various proxy records from South Georgia. Clapperton et al. (1989a) dated moraines on north-central South Georgia to around 2000 cal BP. Rosqvist and Schuber (2003) document cirque glacier advance at the same time in Husvik (**Figure 2A**), while evidence from changes in pollen and macrofossils from Tønsberg Peninsula (**Figure 2A**; Van Der Putten et al., 2004, 2009) and Annenkov Island (**Figure 1**; Strother et al., 2015) mark a sudden change in vegetation in response to cooler and wetter conditions from around 2700 cal BP, and Foster et al. (2016) infer cooling on Annenkov Island based on glycerol dialkyl glycerol tetraethers (GDGTs) (**Figure 11E**). The shift to a cool and wet climate around 2700 cal BP also falls in line with evidence of a major glacier advance prior to 2600 cal BP in the south-central Andes (Grosjean et al., 1998). On James Ross Island, temperatures starts to decrease after 3000 cal BP (Mulvaney et al., 2012; **Figure 11F**), while evidence of peat growth initiation since 2800 cal BP suggests warmer summers in the western Antarctic Peninsula (Yu et al., 2016), demonstrating the Antarctic dipole (Yuan, 2004).

#### Glacier Retreat: 2000 ± 200 to 1000 ± 200 cal BP

Clay and silt dominate the sediments deposited in Diamond bog and Diamond lake between 2000 ± 200 cal BP and 1000 ± 200 cal BP clay (Unit F4, **Figures 6A**, **Figure 11A**). The decrease in relative sand content indicates lower runoff, and a gradual decrease in χbulk suggests glacier retreat back to

mid-Holocene conditions (**Figure 11A**). This coincides with cirque glacier retreat in neighboring Husvik (**Figure 2A**), as demonstrated by Rosqvist and Schuber (2003) (**Figure 11C**) and higher temperatures at the western Antarctic Peninsula inferred from changes in peat growth (Yu et al., 2016).

#### Glacier Advance: 1000 ± 200 to 300 ± 100 cal BP

From 1000 ± 200 to 600 ± 100 cal BP there is a gradual increase in χbulk (**Figure 11A**), which we interpret as glacier advance based on the dominating clay and silt fractions in this interval. Between 500 ± 150 and 300 ± 100 cal BP, a sharp rise in χbulk and a higher concentration of sand and terrestrial plant fragments, suggests higher runoff. Peaking around 400 ± 150 cal BP, we interpret the increase in χbulk to represent the LIA. This glacier advance is also evident from the bog record, where a 4.5 cm thick silt layer with high χbulk values was deposited at 450 ± 200 cal BP. The LIA is also captured by Van Der Bilt et al. (2017) at 300 cal BP (**Figures 11A,B**) and widely documented on the island (e.g., Clapperton et al., 1978, 1989a,b; Gordon and Timmis, 1992; Roberts et al., 2010). Clapperton et al. (1989b) estimated an average ELA depression of 50 m during the LIA on the northeast coast of South Georgia. If applied to our site, the LIA moraine of Diamond glacier is likely represented by moraine M2 (Section Reconstructing the Glacier; **Figure 2B**). Coinciding with the LIA advance, Foster et al. (2016) document rapid cooling from c. 400 cal BP on Annenkov Island, and a millennial-scale cooling trend on James Ross Island culminates at c. 400 cal BP (Mulvaney et al., 2012) (**Figures 11E,F**).

#### Twentieth Century Glacier Advance

Increased χbulk suggests a recent glacier advance during the twentieth century. In the lake record, χbulk increases from AD 1920 ± 10 years, while a silt layer in the bog is dated to AD 1950 ± 10 years. Being a cirque glacier, we argue that Diamond glacier likely advanced in the 1930s, which is the case for cirque and small valley glaciers across the northeastern coast of South Georgia (Hayward, 1983; Gordon and Timmis, 1992; Gordon et al., 2008). The advance concurs with a drop in temperature to 3.3◦C (AD 1931) at Grytviken compared to the mean summer temperature (Dec-Feb) of 4.9 ± 0.8◦C (1905– 2016), and is also recorded in Middle Hamberg lake (Van Der Bilt et al., 2017). According to Landsat imagery, Neumayer Glacier retreated 5600 m between 1999 and 2016 (Pelto, 2016), and we know from our field measurements that Diamond glacier front retreated approximately 150 m between 2008 and 2012. The glacier retreat is reflected in the lake sediments by decreasing χbulk (**Figure 11A**) in the topmost sediments, although not as clearly as in Van Der Bilt et al. (2017)'s record from Middle Hamberg lake (**Figure 11C**). Remobilization of paraglacial sediments in front of the rapidly retreating Diamond glacier, might mask the glacial signal in the lake and bog sediments during the paraglacial period, which typically lasts 10–50 years after a glacier retreat (Bakke et al., 2009). The ongoing glacier recession on South Georgia is well documented by observational data (Gordon et al., 2008; Cook et al., 2010) and follows a pattern of warming and glacier retreat observed from Patagonia to the Antarctic Peninsula (Koppes et al., 2015) as well as the rests of the world (Zemp et al., 2015).

### Glacier Variability on South Georgia vs. Patagonia

The glacier reconstruction from Olsen Valley suggests that Diamond glacier was close to its maximum Holocene extent between 7200 ± 400 and 4800 ± 200 cal BP, 2700 ± 150 and 2000 ± 200 cal BP, 500 ± 150–300 ± 100 cal BP, and in the 1930s. Bertrand et al. (2012) provide a continuous glacier reconstruction from the Northern Patagonian Icefield based on a multiproxy analysis of a fjord sediment core, revealing repeated episodes of centennial scaled glacier advance over the past 5000 years (**Figure 11D**). The two records thus vary on a similar frequency and both show increased glacier activity during the last millennia. The glacier fluctuations documented by Bertrand et al. (2012) seem to be alternately in- and out of phase with those of Diamond glacier, although we find it too early to conclude about timing due to differing glacier response times and dating uncertainties. In general, the two records correspond best during periods of

major glacier advance such as prior to 2000 cal BP, the LIA and the twentieth century advance, while they are out of phase on a millennial time scale during the mid-Holocene (**Figures 11A,D**).

from James Ross Ice Cap, Antarctic Peninsula (F) (Mulvaney et al., 2012).

We speculate that the latitudinal expansion of the westerlies during cold phases causes both sites to be influenced by the same wind circulation system (Lamy et al., 2010), improving the correspondence of the glacier records during these periods. There are also other records from Patagonia to compare with. It turns out that the timing of South Georgia glacier advances concur with several of the Holocene glaciations of the Patagonian Icefield suggested by Anyia (2013) and Strelin et al. (2014), particularly since c. 2700 cal BP. Anyia (2013) proposes five Neoglacial episodes: at 4500–4000 cal BP, 3600–3300 cal BP, 2700–2000 cal BP, 1600–900 cal BP and in the seventeenth to nineteenth centuries century. Based on surface exposure dating, Strelin et al. (2014) suggest six Holocene glacier advances of the Southern Patagonian Icefield at 7700–7200 cal BP, 6000–5000 cal BP, 2500–2000 cal BP, 1500–1100 cal BP >400 cal BP and <300 cal BP. The covariance between South Georgia and Patagonian glacier fluctuations on centennial time scales points to shared regional climate forcing, possibly related to persistent changes in the Southern Hemisphere westerlies which arguably can have a large impact on local variations in precipitation and temperature (cf. Moreno et al., 2014).

### CONCLUSIONS

A 7000-year long continuous cirque glacier reconstruction, based on analyses of glacier-fed lake and bog sediments from the same catchment, suggest that glaciers on South Georgia have varied on centennial time scales for the time interval in question. Evidence for this pattern of variability is obtained by employing the same combination of magnetic, physical, and geochemical sedimentary analyses on the two independently dated sediment archives from the bog and lake, which in sum have enabled us to reproduce the history of the cirque glacier and precisely constrain the timing of the glacier fluctuations. A set of marginal moraines in the catchment outline the past extent of the Neumayer tidewater glacier, which probably filled the lower valley during the Late-Pleistocene, and subsequently advanced to a position just west of Olsen Valley several times, including the past century. In front of the cirque glacier, a mapped moraine sequence provide spatial constraints to the continuous glacier reconstruction. Reconstructed ELAs demonstrate that the recent glacier retreat stands out in a Holocene context although changes are within 100 m altitude. Based on the combination of two different glacial sedimentary archives, a consistent picture emerges of the Holocene glacier evolution in the studied catchment; characterized by deglaciation some time prior to 9900 ± 250 cal BP, and glacier advance close to its maximum Holocene extent between 7200 ± 400 and 4800 ± 200 cal BP, 2700 ± 150 and 2000 ± 200 cal BP, 500 ± 150–300 ± 100 cal BP, and in the 1930s. Changes in glacier size is constrained by a sequence of moraines in front of the

#### REFERENCES

cirque glacier. The reconstructed glacier fluctuations agrees with existing moraine chronologies, two previous continuous glacier reconstructions and various paleoclimate records on South Georgia, highlighting the ability of these sedimentary records to capture concordant shifts. Moreover, the centennial scaled glacier variability documented here are largely in-phase with patterns of retreat and advance of Patagonian glaciers, especially during the past three millennia, although caution is warrant due to chronological uncertainties. This, nevertheless, points to large-scale climate shifts potentially linked to the Southern Hemisphere Westerlies, although this is something that requires further studies before any solid conclusions can be drawn.

### AUTHOR CONTRIBUTIONS

JB, ØP, and LO were responsible for the conceptual design of the work. JB and ØP mapped and the study area in 2008, and JB and colleagues further mapped the catchment and cored the lake and bog in 2012. LO performed the laboratory analyses, and lead the writing and compilation of figures and tables. All authors have given substantial contributions to the interpretations of the data and the writing of the article. In cooperation with LO and JW performed the age-depth modeling and statistical cross-matching (**Figures 7**–**8**), and WvdB designed the maps (**Figures 1**–**3**). All authors approve of the submitted version of the manuscript and thereby agree to be accountable for all aspects of the work, ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

#### ACKNOWLEDGMENTS

This is a contribution from the project "Shifting Climate States of the Polar Regions," funded by the Norwegian Research Council (Grant 210004). We would like to thank Anne Bjune, Bjørn Christian Kvisvik, Åsmund Bakke, Sunniva Solheim Vatle, and Thom Whitfield for their effort during fieldwork, and British Antarctic Survey for their field support. We would also like to thank Sven Lukas and Gunhild Rosqvist for contributing to the manuscript, and Lubna S. J. Al-Saadi and Eivind W. N. Støren for help with grain size measurements. Finally, thank you to reviewers Olga N. Solomina and Bethan J. Davies, whose comments helped improve and clearify the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00002/full#supplementary-material

Okstindbreen, northern Norway: an integrated approach. Quat. Sci. Rev. 29, 1246–1262. doi: 10.1016/j.quascirev.2010.02.012

during the Holocene, northern Folgefonna, western Norway Holocene 15, 161–176. doi: /10.1191/0959683605hl797rp

Antarctica, eds M. R. A. Thomson, L. A. Crame, and J. W. Thomson (Cambridge: Cambridge University Press), 353–359.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Oppedal, Bakke, Paasche, Werner and van der Bilt. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Imprints of Climate Signals in a 204 Year δ 18O Tree-Ring Record of Nothofagus pumilio From Perito Moreno Glacier, Southern Patagonia (50◦S)

Jussi Grießinger 1 \*, Lukas Langhamer 2 , Christoph Schneider 3 , Björn-Lukas Saß1 , David Steger 3 , Pedro Skvarca4 , Matthias H. Braun1 , Wolfgang J.-H. Meier 1 , Ana M. Srur 5 and Philipp Hochreuther 1

#### Edited by:

Davide Tiranti, Agenzia Regionale per la Protezione Ambientale (ARPA), Italy

#### Reviewed by:

Bao Yang, Cold and Arid Regions Environmental and Engineering Research Institute (CAS), China Li Wu, Anhui Normal University, China

> \*Correspondence: Jussi Grießinger<EMAIL_ADDRESS>
#### Specialty section: