Next Article in Journal
Impacts of Climate Change and Human Activities on NDVI in the Qinghai-Tibet Plateau
Next Article in Special Issue
Sentinel-1 Response to Canopy Moisture in Mediterranean Forests before and after Fire Events
Previous Article in Journal
Ensemble Machine Learning for Mapping Tree Species Alpha-Diversity Using Multi-Source Satellite Data in an Ecuadorian Seasonally Dry Forest
Previous Article in Special Issue
On the Choice of the Most Suitable Period to Map Hill Lakes via Spectral Separability and Object-Based Image Analyses
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Detection of Crustal Uplift Deformation in Response to Glacier Wastage in Southern Patagonia

1
Laboratorio de Geomática Andina, Instituto Argentino de Nivología, Glaciología y Ciencias Ambientales (IANIGLA), Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET), Mendoza CP 5500, Argentina
2
Departamento de Geografía, Universidad de Chile, Santiago 1234, Chile
3
Facultad de Ciencias Exactas y Tecnología, Universidad Nacional de Tucumán, San Miguel de Tucumán T4000, Argentina
4
Instituto Geográfico Nacional (IGN), AAD Buenos Aires C1426, Argentina
5
Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, La Plata B1900, Argentina
6
International Center of Earth Sciences, Universidad Nacional de Cuyo, Mendoza CP 5500, Argentina
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(3), 584; https://doi.org/10.3390/rs15030584
Submission received: 31 October 2022 / Revised: 16 December 2022 / Accepted: 17 January 2023 / Published: 18 January 2023
(This article belongs to the Special Issue Remote Sensing in Geomatics)

Abstract

:
The Southern Patagonian Icefield (SPI) is the largest continuous ice mass in the Southern Hemisphere outside Antarctica. It has been shrinking since the Little Ice Age (LIA) period, with increasing rates in recent years. An uplift of crustal deformation in response to this deglaciation process has been expected. The goal of this investigation is to analyze the crustal deformation caused by ice retreat using time-series data from continuous GPS stations (2015–2020) in the northern area of the SPI. For this purpose, we installed two continuous GPS stations on rocky nunataks of the SPI (the GRCS near Greve glacier and the GBCS close by Cerro Gorra Blanca). In addition, ice elevation changes (2000–2019) were analyzed by the co-registration of the SRTM digital elevation model and ICESat elevation data points. The results of the vertical components are positive (36.55 ± 2.58 mm a−1), with a maximum at GBCS, indicating the highest rate of crustal uplift ever continuously recorded in Patagonia; in addition, the mean horizontal velocities reached 11.7 mm a−1 with an azimuth of 43°. The negative ice elevation changes detected in the region have also accelerated in the recent two decades, with a median Δ h (elevation change) of −3.36 ± 0.01 m a−1 in the ablation zone. The seasonality of the GPS signals was contrasted with the water levels of the main Patagonian lakes around the SPI, detecting a complex interplay between them. Hence, the study sheds light on the knowledge of the crustal uplift as evidence of the wastage experienced by the SPI glaciers.

1. Introduction

The Southern Patagonian Icefield (SPI) represents the major temperate ice mass of the Southern Hemisphere outside Antarctica. Together with nearby smaller glaciers the SPI cover an area of 14,151 km2 [1] with a maximum ice thickness of 1.5 km [2]. During the last century, there has been an increased tendency for ice mass loss in glaciers and ice sheets around the world as a consequence of climate change [3]. The Patagonian Andean glaciers are not the exception since they have suffered strong retreat from the Little Ice Age (LIA) period at the end of the 19th century, with an accelerated trend toward the end of the 20th century and the beginning of the 21st century ([4,5,6,7,8,9,10,11,12,13], and others). As a result, the SPI experienced an ice mass loss of 10.7 Gt a−1 from 2010 to 2019 [14] with a mean contribution to sea level rise throughout the last 360 years (1650–2010) equivalent to <0.005 ± 0.0001 mm a−1 [15].
As ice sheets and glaciers melt, and water is redistributed at the global ocean scale, the deformation of the Earth’s crust generates a complex pattern of three-dimensional movements on the Earth’s surface [16]. Although the current sea level changes are mainly a consequence of the massive global ice sheet melting, the mean sea level remains in continuous change due to late viscoelastic responses of the Earth by surface mass redistribution [17]. Consequently, this is a delayed effect of continental ice sheets since the Last Glacial Maximum (LGM, 21,000 years ago) [18]. Those changes originate in the solid Earth through an isostatic readjustment or GIA (glacial isostatic adjustment), one of the geophysical processes that can easily be observed without resorting to sophisticated scientific measurement techniques. The currently observed uplift rates are the sum of an immediate elastic response of the Earth’s crust to current ice mass changes and the viscous response to ice load changes in the past [19]. The crustal deformation survey by the Global Positioning System (GPS) geodetic network has been widely used to monitor crustal variations [20] since the 1990s thanks to it being an inexpensive and precise technique [21]. Furthermore, the system allows for the detection of responses to depth and superficial deformations that have been observed over time ([22,23], and others). Currently, geodetic time series are widely used for detecting the Earth’s deformation processes.
In the SPI region, there are evidences of crustal deformation responses to ice mass volume change during the last 150 years [24]. Because of these processes, an exceptionally fast crustal elevation with an interplate tectonic movement has been experienced by the SPI [19]. In recent decades, GPS stations have been placed on the edge, both sides, of the SPI with the objective of monitoring crustal deformation of the South American plate ([25], among others). This allows researchers to establish that GIA is one of the primary components explaining the crustal observed deformation, with another factor being the western interseismic tectonic deformation field related to plate subduction [26]. In general, previous studies have placed the GPS stations (continuous and episodic ones) at the margin of the SPI, with few discrete observations at the higher SPI plateau.
The goal of the present study is to analyze the crustal movements in the northern and higher part of the SPI by GPS continuous measurements during the 2015–2020 period at two locations selected mainly by accessibility (see Figure 1). Due to extreme conditions prevailing in the area, both stations were installed in 2015 on a rigid metal structure fixed to rock outcrops. The first GPS system was called GRCS (Greve Continuous Station), and recorded for a total of 378 days. The second one, named GBCS (Gorra Blanca Continuous Station), registered 911 days in total. Two episodic GPS sites (PUMA and LPAS) located near the front of the Upsala glacier (see station details in Table 1) were also included. For a better validation of the seasonality analyses, two continuous GPS stations called CHLT and VOGH belonging to the International GPS Service (IGS) network were added. Scientific processing of GPS data yielded network station velocities in N, E, and Up, allowing the tendency, seasonality, and other periodic signals of the time series to be examined. The hydrologic system is driven by the loading processes and the transfer of mass from the upper to the lower parts of the SPI basin where the glacial lakes are key natural elements. Consequently, we believe that the trend and seasonality of the observed altimetric GPS signal are related to the water levels of the lakes. For this purpose, we use the gauged lake data of Viedma, O Higgins, and Chico, and the La Leona river.
In addition, to analyze the ice loss in the region, we estimated ice elevation changes by differentiation of the Space Shuttle Radar Topography Mission (SRTM) model from February 2000 and altimetric data profiles from the Ice, Cloud, and Land Elevation Satellite (ICESat) mission-1 (2005–2009), and mission-2 (2018–2019). Consequently, the study contributes to the knowledge of the crustal uplift in the northern part of the SPI by using geomatics/remote sensing techniques.

Study Area

The SPI has an average width of 40 km and north–south length of 350 km (Figure 1) [27], at elevations between 1000–1500 m a.s.l. (above sea level) [28]. The Equilibrium Line Altitudinal (ELA) at the glaciers of the western slope varies between 800 to 1000 m a.s.l., while at the eastern side, the ELA can reach 1300 m a.s.l. [29]. The large hydrographic basin of this region is the Santa Cruz river basin, including the lakes Viedma and Argentino. The lake Argentino is the largest in Argentina, with an east–west extension of approximately 100 km, an area of 1330 km2, and a water depth exceeding 600 m near the Upsala glacier front [30]; while lake Viedma presents an approximate length of 80 km east–west, an area of 1088 km2, and a maximum water depth of 919 m [31]. The region is located on the South American plate bounded to the west by the Nazca and Antarctic plates, forming the triple junction called Chile Triple Junction (CTJ) whose position is at ~−46.5° S (Figure 1).

2. Geodetic and Satellites Input Data

2.1. Continuous GPS Stations REFUGIO GREVE and GORRA BLANCA: Field Activities

The continuous GPS stations are part of a local network located in the high plateau at the northern zone of the SPI. The GRCS station is placed in the accumulation area (see Figure 1) at 1438 m a.s.l., and it was installed on a nunatak. A Trimble receiver model NetR9 and a Trimble Zephyr geodetic antenna (TRM41249.00), without a radome was used. The system was powered by a solar panel installed and mounted on a 1.7 m steel pole attached to the rock. On the other hand, the GBCS station was located on the edge of the high plateau, in a rocky outcrop at approximately 1600 m a.s.l. (see Figure 1). A Trimble receiver model NetRS and a Trimble Zephyr geodetic antenna (TRM41249.00) without a radome was used. A similar solar panel was installed and mounted on a 1.7 m steel pole similar to the GRCS. Both stations were configured with a logging interval of 15 s and elevation mask of 5°.
In total, the GRCS recorded 378 discontinuous days from 25 October 2015 to 18 June 2018 when the station ceased operation due to low battery charge and damage to the solar panel due to strong winds. GBCS recorded in total 911 days from 2015 to November 2020 (see Table 1) with some interruptions due to energy failures. These data collection interruptions were a consequence of the extreme weather conditions prevailing in the area.

2.2. GPS Episodic Stations: Upsala Glacier

Two episodic GPS stations (ES) installed on rocky outcrops located in the vicinity of the Upsala glacier front (LPAS and PUMA) were added to the geodetic network (Figure 1), following the fixed anchor State Survey Mark [32] minimizing the leveling error [33]. The minimum observation period at these stations was 48 h per campaign, corresponding to two 24 h observation sessions. The measurements in all periods were carried out at the same time of the year in order to reduce the influence of the seasonal effects on the estimation of site velocities. In total, the PUMA station recorded data during 13 days between 2015 and 2017. Likewise, LPAS registered 12 days during the same period. Table 1 shows the GPS site occupation time of each of the continuous and episodic stations of the geodetic network.

2.3. GPS Accessory Data at Viedma Glacier

Surface elevation of the glacier lower tongue was also surveyed in April 2015 with a Trimble 5700 dual-frequency GPS receiver using a Differential GPS (DGPS) mode linked to the CHLT station.

2.4. ICESat Data Profiles

The ICESat-1 mission provided multiyear elevation data for ice sheet mass balance determination from 2003 to 2009 [34]. The GLAH14 product of the GLAS/ICESat L2 Global Land Surface Altimetry data set, version 34 in HDF5 format, was used for this study (Table 2). Heights have a centimeter order vertical error [35]. All elevation data are referenced to the WGS84 ellipsoid. We downloaded the available data from 2004 to 2008 collected at the upper zone of the Upsala and Viedma glaciers (see Figure 1). These products were obtained from the official NASA website (https://search.earthdata.nasa.gov/, accessed on 21 December 2021). We estimated ice elevation changes by comparing the data profiles to the SRTM elevation model from 2000 to 2004–2009.
The ICESat-2 mission has been used in several areas of the planet for monitoring the cryosphere since September 2018. Land ice heights measured through the Advanced Topographic Laser Altimeter System (ATLAS) represent the average height of the land ice surface averaged over ~17 m ground segments and spaced 20 m apart. ATL06 processing calculates ground ice surface heights and stores them per ground track so that subsequent measurements on the same reference ground track can be easily compared. Heights are referenced to the WGS84 reference system (ITRF2014 reference frame) [36] and have a vertical precision of centimeter order [37]. For this study we selected data collected from November 2018 to November 2019 (Table 2). HDF5 files were downloaded from the official NASA site (https://search.earthdata.nasa.gov/search, accessed on 15 December 2021), considering data from the SPI’s northern region, including the Upsala, Viedma, O’Higgings, and Chico glaciers. We estimated ice elevation changes by comparing the data profiles to the SRTM elevation model from 2000 to 2018/19.

2.5. SRTM Model

Map products derived from SRTM data were sampled to 1 arc-sec by 1 arc-sec grid (30 m by 30 m approximately), with a linear relative height error of less than 10 m [38]. We used a mosaic based on six SRTM Global Version 3 files, Digital Object Identifier (DOI) number:/10.5066/F7PR7TFT. The SRTM model was used to estimate the elevation change over time ( Δ h ) with respect to the ICESat profiles. To keep consistency between the data sets, elevation values were converted to WGS84 ellipsoid with EGM96 geopotential model [13]. The final product was projected according to UTM system, zone 18S, with a resolution of 30 m.

2.6. Gauging Data of Lake Viedma, O’Higgins and Chico

The water levels data from lake Viedma and La Leona river were obtained from the Argentine National Undersecretariat of Water Resources (https://www.argentina.gob.ar/obras-publicas/infraestructura-y-politica-hidrica, accessed on 21 December 2021). Lake Viedma: daily scale reading values from January 2015 to 2021; the limnometric scale in La Leona river: daily values of scale from 2015 to 2021. La Leona river connects lake Viedma with Argentino lake. Daily gauge data series collected by the Chilean Water Cadaster (DGA) at the lakes O’Higgins and Chico were also used (www.dga.cl, accessed on 21 December 2021).

3. Geospatial Data Processing

3.1. Scientific Processing of the GPS Network

Data daily solutions were obtained at GRCS, GBCS, LPAS, and PUMA stations that were referenced to the International GPS Service 2008 (IGS08) terrestrial reference frame by scientific processing of GPS data with BERNESE v5.0 software [39]. Accurate ephemerides and Earth rotation parameters consistent with the reference frame were used, as well as corrections to the phase center of the antennas. The FES2004 model was applied to perform the oceanic load corrections. For the tropospheric correction we used the Vienna Mapping Function with zenithal delay estimation every 2 h [40], and for the ionospheric corrections we followed Fritsche et al. [41]. Finally, the ambiguity resolution strategy was QIF (quasi-ionosphere-free).
The daily GPS time series allowed linear displacement trends to be estimated for each site after outlier removal. From these series, repeatability analysis of each component of the continuous station coordinates sites were performed to estimate the uncertainty in positioning and velocities. Uncertainty of 10 mm in the antenna position for the three components was added for the episodic stations since the antenna was mounted in each campaign. The total uncertainty obtained for each site was propagated according to the time interval between the first and last observation at each site. In particular, for the vertical component (Up) an uncertainty of 2 mm a−1 was added in agreement with Bevis and Brown [42].

3.2. Geodetic GPS Data: Linkage to the Reference Frame

In order to obtain precise positions of the geodetic network, the data were linked to the terrestrial reference frame IGS08 [43], providing consistency with the IGS products, such as satellite orbits, corrections to antenna phase centers, Earth orientation and rotation parameters, etc. Eight continuous GPS stations of the IGS global network were linked to the proposed geodetic network used in this study. These stations generate the regional network and are BRAZ, FALK, LPGS, OHI2, PLAM, PARC and UNSA (IGS Stations), and CHLT (IGN-Ar Station) (Figure 1). Data from these stations were downloaded in RINEX format from the Scripps Orbit and Permanent Array Center (SOPAC) repository (http://garner.ucsd.edu/, accessed on 21 December 2021), except for the CHLT station, which belongs to the Argentine National Geographic Institute (IGN-Ar). Table 3 shows the stations used, as well as their location, geodetic coordinates, tectonic plate where they are located, and the agency to which they belong.

3.3. GPS Time Series Tendency and Seasonality Analysis

GPS data series derived from the positions of continuous GPS stations allow the determination of variations over time by extracting the seasonality of the signals [44]. In addition, seasonal variations can be manifested as annual and semi-annual oscillations [45]. We use the Empirical Mode Decomposition (EMD) method for seasonality analysis, which decomposes the signal in the time domain. The way it decomposes ensures integrity, where the original signal can be represented as the sum of the new basis functions found. The process is useful for the analysis of natural signals because it is an adaptive method, where the definition of the basis is derived from the data without imposing linearity or stationarity conditions [46]. This method assumes that, for each measurement instant, the time series can have several simple oscillatory modes of different frequencies superimposed on each other. Each of these components is called an intrinsic mode function (IMF) and must achieve the following requirements: (a) throughout the data set, the number of extremes and the number of zeros must be equal or differ by at most one; (b) for each data value, the mean value of the envelope defined using the local maxima and local minima is close to zero.
Unlike other non-parametric methods, this method does not require a zero reference in the original signal. During the screening process, the trend of the original data is obtained as an IMF. The number of functions depends on the original data series and the parameters for the envelope calculation. The obtained IMFs can then be grouped (summarized) as representing the same physical phenomenon.

3.3.1. Robust Time-Series Data: Outlier Removal and Interpolation

Before calculating the IMFs, it is essential to prepare the signals by eliminating outliers and interpolating gaps. The interpolation for missing data is an essential step before processing GNSS time series with other advanced methods [47], as missing data may occur for various reasons [48], e.g., due to weather contingencies among other causes, and could affect the accuracy of the trend analysis. To achieve a robust data series, we decided to combine the signals from the GRCS and GBCS stations. As both stations are less than 50 km apart and above the same icefield, we can assume that the sites have a similar crustal uplift response. The study period considers 2154 days in total, with the GBCS series covering 42% and the GRCS series covering 17% of this period of time. In a first step, the time series of both stations were taken to the same altitudinal scale through a least-squares adjustment method since they had an offset of approximately 200 m. From now on, we will call the new signal GBGR-CS.
For the analysis of the trend and seasonality of the GBGR-CS station, we annexed two continuous stations belonging to the IGS network, CHLT and VOGH. We removed outliers with the three-sigma criterion. Subsequently, an interpolator was used to fill in missing data in each series. Next, the average year without a trend was calculated. This was achieved by comparing and calculating the daily average for each site. The new incorporated data do not come from a simple linear interpolation of the previous and subsequent points. However, they are calculated from the average year and the trend of the site for that time.

3.3.2. Sifting Process

The algorithm allows for finding the IMFs from a time series of data, which we will call x ( t ) [49]. First, we searched for local extrema and classified them into relative maxima or minima. Second, two time series were assembled, one with each of them. To form the upper envelope emax (t), the maxima are connected with a cubic spline., and the minima are connected to form the lower emin (t). Next, the mean m1,1 is calculated between the two envelopes. Finally, the first proto-mode (h) is computed by subtracting the mean from the original series as Equation (1) shows:
h 1 , 1 = x ( t ) m 1 , 1
The procedure is repeated k times starting from Equation (1) until the following equation is obtained:
h i , k = x ( t ) m i , k
The purpose of Equation (2) is to eliminate background signals on which the IMF is built, and then bring it closer to the requirements of its definition. The first IMF component is taken as Equation (3), which should contain the short period variations of the original signal x(t).
c i = h i , k
Finally, the residual is calculated as the difference between the input signal and the found IMF (Equation (4)). This residual is treated as the new input data set and the procedure is repeated until a cut-off criterion is satisfied [50]. Note that the criteria for stopping the iterative process are satisfied when the cut-off parameter is below a set value.
r i = x ( t ) c i
The input data can be thought of as the sum of n intrinsic mode functions plus a residual (Equation (5)):
x ( t ) = i = 1 N c i ( t ) + r N ( t ) .

3.3.3. Seasonality Components

From the IMFs, a frequency analysis of the seasonal components was performed through the Fourier transform. Moreover, in order to assess the link between the altimetric GPS signal and the water level of the lakes, we calculated the correlation coefficients between the IMFs obtained and the gauge data of lakes Argentino, Viedma, and O’Higgins and the La Leona River.

3.4. ICESat Profile Processing

Data filtering was performed in three steps: (1) spatial filtering of profiles with data collected on Upsala, Viedma glaciers, and in the north region of SPI (near GRCS and GBCS stations). The glacier basin outlines were obtained from the 2000 Randolph Glacier Inventory (RGI v6.0) provided by GLIMS (https://www.glims.org/RGI/rgi60_dl.html, accessed on 4 July 2021); (2) in the case of ICESat-1, the parameters used to filter the data were those used by Vacaflor et al. [13]; ICESAT-2 data were used as long as their elevations had the condition atl06_quality_summary = 0, indicating that no elevation quality problems were identified [36]; (3) only elevations with ICESat values differing by less than 50 m from SRTM were accepted since even after applying step (2), anomalous elevations with positive differences against SRTM were observed. The SRTM elevation model was co-registered to ICESat-1 elevation profiles using the methodology and parameter calculations described by Vacaflor et al. [13] to estimate elevation changes.
For a better estimation of elevation changes, two filters were performed: first, it is well known that the error in high mountainous areas increases with a slope >25° [51]; hence, all points with values equal or higher than this threshold were filtered out. Second, to obtain the seasonality elevation change, we separated the results according to the ablation and accumulation area defined, respectively, by the ELA of the main glaciers [29]. We estimated the relative errors by comparing the differencing result of the SRTM and elevation profiles over stable terrain (off-glacier zones) using the median and the normalized median absolute deviation (NMAD) [52]. Note that the median as the central measurement, and NMAD as the dispersion measurement are robust statistical indicators [53]. The relative random error ( σ Δ h r a n d ) of the median difference for each profile ( Δ h ^ (m)) was calculated using the NMAD of Δ h , Equation (6), in off-glacier zones for all the dates:
σ Δ h r a n d = N M A D n
where n is the number of observations.

4. Results

4.1. Velocities and Tendency of the Data Series

The estimated velocities for the GPS network sites along with their uncertainties for a 3- and 5-year observation period are shown in Table 4. The study period was 2154 days with some interruptions; through the temporal resolution given by the continuous data from the stations, it was possible to estimate the displacement for the sites in the three components (North, East, and Up) (Figure 2). The highest value in Up was at the GBCS station located on the periphery of the SPI, reaching a velocity of 37.16 ± 2.10 mm a−1; then, in the north area of SPI, the uplift rates found were high. On the other hand, the lowest value recorded was for the PUMA episodic station with an Up velocity of 20.30 mm ± 5.03 mm a−1 in the southern region of the network (Figure 3a). The results of the data series in the three components N, E, and Up of the GBCS and GRCS stations are shown in Figure 2, the trend of which in the Up component is positive, which would indicate crustal uplift at all the sites.
The estimated horizontal components have a mean velocity of 11.7 mm a−1 with an azimuth in the northeast direction with respect to the IGS08 frame. The continuous stations GRCS and GBCS have azimuth values of 17° and 22°, respectively, and the episodic stations LPAS and PUMA, which are located to the south, have a more easterly orientation with values of 54° and 64°, respectively. This trend in velocities responds to the regional crustal deformation of Patagonia in the South American plate (SOAM) [54]. It should be noted that the area of the SPI where the network is located does not belong to the rigid part of the plate but is affected by a regional-scale deformation [26].
A topographic profile was performed between the stations located in the SPI (GRCS and GBCS), and one in the periphery (CHLT) and another more distant (UNPA) in the Patagonian plateau, the latter belonging to RAMSAC (Red Argentina de Monitoreo Satelital Continuo) [55]. This profile was derived from the Shutter Radar Topography Mission (SRTM) digital model. Figure 3b shows the profile plotted between the four stations, which denote an altitudinal delta of approximately 1600 m. Comparing the trends, we observe that at those stations located at high altitude within SPI and CHLT, the slope has a marked positive upward trend (Figure 3c). Note that GBCS has a slope value of 38 mm a−1, GRCS of 35 mm a−1, and CHLT of 25 mm a−1. In contrast, UNPA station near the Atlantic Ocean has a slightly negative vertical displacement rate with a value of −0.077 mm a−1. These results are consistent with those expected, indicating a crustal uplift of the SPI region and its periphery.

4.2. Tendency and Seasonality Analysis

The GBGR-CS combined series improved the amount of data by 17%, covering a total of 49% of the study period. Figure 4 shows the result for the CHLT, VOGH, and GBGR-CS time series. In red is shown the original series without outliers, and in blue, the interpolation points are added to continue with the processing for seasonality estimation.
To assess the seasonality of the series, the signals were decomposed through the IMFs. Figure 5 shows the results of the GBGR-CS, CHLT and VOGH stations’ time-series without gaps. They were decomposed into three intrinsic mode functions: long period, medium, and high frequency. According to these results, we decided to emphasize the long period component because it showed a positive linear trend, but at the same time it overlaps with an annual periodic component of 10 mm. These few millimeters of signal amplitude with an annual or semi-annual period are partially caused by atmospheric and hydrological loads, the thermal expansion of soil and monumentations, or variations in multiple paths [56]. It should be noted that the high and medium frequency components have not been used for the seasonality analysis. The combined continuous stations located at the SPI plateau (GBGR-CS) and VOGH in the periphery, exhibit a similar tendency and period behavior; however, the CHLT station shows a linear and positive trend with a very slight seasonality. The estimated slopes for the stations according to a linear fit with a 95% confidence interval were 36.3 mm a−1 for GBGR-CS, 24.8 mm a−1 for CHLT, and 21.4 mm a−1 for VOGH.
The results of the Fourier transform and amplitude in the IFM’s long period signal show that in two stations, the frequency is associated with an annual variation (f = 1 [1/a]) (See Figure 6). This signal predominates in the magnitude of the spectrum at stations GBGR-CS and VOGH, which may be linked to the hydrological cycle. Analyzing interannual frequency, only VOGH shows a frequency at f = 2.0 [1/a] that corresponds to a period of six months and the other to a period of 4 months. Likewise, an analysis of the yearly variations shows that the amplitudes of the seasonal signal are 0.0081 m (0.0162 m peak to peak) for GBGR-CS, and 0.085 m (0.017 m peak to peak) for VOGH. Although the adjustment for CHLT is poor, the amplitude reached 0.0019 m (0.0038 m peak to peak). As is shown in Figure 6, the annual signal does not predominate over the interannual signal as in the other stations.

Lake Variability Correlation with GPS Signal

The maximum and minimum epochs of the gauging data series Chico and O’Higgins in Chile, and Viedma and La Leona in Argentina are coincident with some delays. The influence of lake variability on the GPS signal are shown in Table 5 and Figure 7. The Pearson’s correlation coefficients between the seasonal component of the vertical displacement time series and the level variations in the two lakes exceed 67% in GBGR-CS, while VOGH yielded at Chico 48% and at O’Higgins 64%, respectively. For the CHLT station, the correlation of the seasonal component with respect to gauge data does not exceed 22%. The highest correlation was observed at lake Viedma, and Chico achieved the lowest. The long period component of the signal in CHLT is dominated by the linear component and, in the seasonal variation found, the annual signal does not predominate over the interannual ones, as in the other stations.

4.3. Ice Elevation Changes

Altimetric ice changes obtained from the co-registration of the SRTM model and ICESat profiles (missions 1–2), indicate that the study area has experienced negative changes, with a marked loss of ice (Table 5). The altimetric profiles of ICESat-1 are concentrated in the high plateau zone between the Upsala and Viedma glaciers. The profiles cover an altitude range between 545 m and 2473 m a.s.l., including Upsala, Viedma, and smaller tributaries. Ice elevation changes during the first period of study (2000–2004/8) yielded minimum values of −1.53 ± 0.16 m a−1 in the accumulation zone, and maximum values of −5.49 ± 0.14 m a−1 in the ablation zone.
Nevertheless, the ICESat-2 mission from 2018/2019 included the north of the SPI where the continuous GPS stations are located. These ICESat-2 profiles are located near the fronts of the Upsala, Viedma, Chico, and O’Higgings glaciers with altitude ranges between 191 m to 3342 m. According to the obtained results (Table 6), the altimetric changes yielded ~−1 m a−1 in the accumulation zone. However, the ablation zone reached a median value of −3.36 ± 0.01 m a−1. The results from 2000–2018 are higher than the second period, but the 2018 data were only collected in November. Moreover, the profiles cover an area that corresponds to glacier tongues. Likewise, for the year 2019, the data sampling has a better temporal distribution since it comprises several months of the year with more than 110,000 data samples. The profiles were generally located on the eastern slope of the SPI between the Upsala glaciers and the Jorge Mont accumulation zone.
Figure 8 shows the altimetric changes that have occurred in four of the main SPI glaciers, where there is a generalized negative change trend, especially below the ELA. The decrease in elevation changes at the lowest elevation is related with the frontal retreats of calving fronts. Note that the parallel patterns shown in the graphs are associated with the complete disappearance of the glacier terminus, especially at Upsala and Viedma whose fronts have retreated considerably during the last decade. In addition, we analyzed the data from the GPS collected on the glacier surface in 2015 and compared them with the SRTM 2000 model. The results yielded a Δ h equal to 96 m in 15 years, which is a significant change. Note that the glacier disappeared in this area in November 2018 due to a massive calving event.

5. Discussion

The continuous SPI ice losses since the middle of the 20th century have accelerated during the 21st century, particularly with high frontal ablation experienced by freshwater calving glaciers. Although glacial retreat of the SPI is heterogeneous, the main contributor to the total mass losses is calving. Calving glaciers represent 91% of Patagonian Icefields [14]. The spatial pattern of the SPI ice elevation change is complex and non-uniform comparing different basins [11]. Studies by Gomez et al. [57] using geodetic data suggest that the El Niño–Southern Oscillation (ENSO) 2015–2016 had a major impact on SPI ice mass changes that has continued through 2018. In this sense, the results we have achieved in the period 2000–2019, both in the plateau (accumulation) and in the frontal ablation zone of the SPI eastern slope, are coincident with those reported by Minowa et al. [14].
From the analysis of ICESat-1, we estimated a maximum Δ h rate of −2.58 ± 0.26 m a−1 for the period 2000–2006 for the high plateau of the Upsala glacier. However, the highest thinning rate was found at the ablation zone of the Upsala glacier, yielding −5.49 ± 0.14 m a−1 in the 2000–2007 period. By using ICESat-2, the value reached −3.36 ± 0.01 m a−1 for the 2000–2019 period. Despite the difference in time intervals, our results are similar to those reported by Malz et al. [58], Jaber et al. [11], and Minowa et al. [14], showing that a significant retreat and mass losses have taken place in the SPI glaciers during the last decade. On the other hand, our results are higher than those estimated by Vacaflor et al. [13] who estimated the geodetic mass balance of the Santa Cruz River basin 1979–2018 period, where the Upsala and Viedma glaciers reached maximum rates of −2.16 ± 0.11 m a−1 and −1.59 ± 0.11 m a−1, respectively. The discrepancy could be attributed firstly to the longer time span and, secondly, because they used optical sensors (Hexagon KH-9 and ASTER) for the estimation of elevation changes.
The four most important glaciers on the eastern and northern side of the SPI (O’Higgins, Chico, Viedma, and Upsala) are shrinking rapidly. Our results show negative exposed trends are intensified in areas below the ELA (1300 m a.s.l.), and the ice elevation has changed by more than 100 m in almost two decades (Figure 8). During the period 2000–2015/16, these four glaciers have recorded thinning rates of up to 2 ma−1 in the accumulation zones, while in the tongues, the Chico reaches up to −6 m a−1 and Viedma up to −4 m a−1 [58]. Jaber et al. [11], found a high surface elevation loss rate in the ablation areas (lower terminus) of the Upsala and Viedma glaciers in the summer of 2011/2012. Instead, the Viedma glacier has suffered an alarming retreat since 2015, with mean velocities at its terminus of 3.5 m d−1 between April 2014 and April 2016 [59]. In the last seven years, it has lost an area of 2.5 km2 in the frontal zone, equivalent to what it had previously lost in 50 years.
Consequently, the acceleration in ice loss would produce an unequivocal acceleration in the elastic response of the Earth. The volumetric decrease in SPI ice masses in the last 150 years could be linked to the Earth’s crustal uplift [24]. Geodetic studies to detect the crustal deformation have been carried out in the SPI region. Dietrich et al. [25] found that in the area near the SPI plateau, the uplift displacement was 39 mm a−1 during the period 2003–2006. Considering a longer time (1999–2014), Lange et al. [19] and Richter et al. [26] observed crustal uplift in the Nunatak Viedma and Cerro Gorra Blanca of 31 and 41 mm a−1, respectively. It should be noted that the Model B of vertical crustal motion prediction proposed by Lange et al. [19] matches with the results found in our study. In consequence, the collected GPS data over the SPI accumulation zone, the first of its type, represent a novel approach to improving future viscoelastic glacial isostatic adjustment models in the region.
Our GPS stations are placed in a particular geological zone, affected by the subduction of the spreading-ridge system composed of the Nazca and Antarctic plates, from the Miocene to the present, which led to the formation of an asthenospheric window (Figure 1) [60,61], contemporary with the onset of the Patagonian glaciations [62]. This gap in the asthenosphere generated thermal and chemical anomalies and consequent changes in the lithosphere–asthenosphere boundary below Patagonia [63], as well as a thinning of the lithosphere [19,64]. Above this slab window, the viscoelastic response of the solid Earth to ice mass changes is fast [19]. This could be one of the causes of fast crustal uplift, as shown by our results. Hence, this slab window provides key evidence supporting the previously hypothesized connection between post-Little Ice Age anthropogenic ice mass loss and rapid geodetically observed glacial isostatic uplift [65].
The region is in a particular tectonic context since there is a point where three tectonic plates (South American, Nazca, and Antarctic) coexist at the north of the SPI (The triple junction point at Peninsula Taitao, see Figure 1). This interaction generated a window in the lithosphere through which the asthenosphere could flow and combine with surface layers, having geomechanical and geodynamic consequences such as regional deformation and a higher uplift rate than in other regions, e.g., Antarctica or Greenland [25,61,66,67]. The horizontal displacements obtained in our study show a trend in a northeasterly direction with a rate of 14 mm a−1. This indicates a regional-scale deformation, indicating that the SPI area can be considered as a non-South American craton sector that is the stable rigid part of the South American plate [26].

Lake Seasonal Variability in GPS Time-Series

Lake-level variations may help us to understand the interactions between glaciers, lakes, and solid earth in an environment affected by ongoing climate change, rapid ice mass losses, and intense GIA [68]. Until present, there has been a lack of studies in Patagonia related to crustal deformation by the Up-GPS signal and lake seasonality interactions. Pereira et al. [69] obtained a weak correlation for some stations in the Patagonian lakes with GRACE (2002–2017 period).
The hydrological system surrounding the GPS stations is comprised of glacial lakes of the region, some of them the largest in Patagonia. However, the area is indeed under the influence of a complex hydrological system that refers to various sources of water, including glaciers, snow, canopy water, surface water (lakes, rivers, and reservoirs), soil moisture, and groundwater. In addition, the natural processes that we are investigating are dominated by an annual period, even the GPS signal. Distinguishing those signals in geodetic measurements is critical for monitoring spatiotemporal variations in continental water storage and tectonic signals [70]. GPS records are related to the integration of all tectonic and non-tectonic sources, including the mass loading from the atmosphere, ocean, snow, soil moisture, and groundwater [71]. Therefore, isolating the seasonal signals that demonstrate the interactions between solid Earth and the hydrological cycle remains a difficult task, especially in Patagonia [72].
Missing data and gross errors always exist in GPS time series due to inevitable factors of the environment and human activities, such as the replacement of the receiver antenna, bad observation conditions, signal interruptions between the satellite and receiver, among others [47]. These interruptions increase signal instability; moreover, GPS velocity solutions tend to be unstable until a 2-year data span is exceeded [73]. In spite of the interruptions, our two continuous GPS stations continuously recorded crustal deformations for more than two years providing novel insights into the complex natural factors’ interplay taking place in the region. The possible signal instabilities and errors due to the harsh environment conditions affecting GPS receivers and antennas were overcome by the signal combination and data gap filling applied method. The results of the seasonality analysis derived from gauge data and long period signals (IMFs) show dissimilar correlation values between those GPS signals, where the mean correlations for GBGR-CS, CHLT, and VOGH are 0.43, 0.09, and 0.61, respectively. However, thanks to this approach, the results are consistent with theory and aligned with previous studies [71]. In summary, our results deserve further analysis by combining more field data with models that can address the complexities in the interconnected hydrological–glaciological–tectonic Patagonian natural system.

6. Conclusions

In this paper, we presented the first continuous GPS data obtained at two stations drilled into rock outcrops within the Southern Patagonian Icefield. These stations operated between 2015 and 2020 yielding a vertical uplift of 36.55 ± 2.58 mm a−1, the highest rate of crustal deformation continuously recorded in Patagonia. This crustal deformation is a response to glacial ice adjustment that has taken place in the region due to strong and fast ice wastage as was detected in this work by the comparison of ICESat-2 and SRTM data in the ablation area between 2000 and 2019, resulting in a median annual ice thinning of −3.36 ± 0.01 m a−1 in the ablation zone, and −0.94 ± 0.01 m a−1 in the accumulation zone. The main results of this study are similar in magnitude, trend, and signs compared to previously reported uplift rates mainly obtained by episodic surveys carried out in the region. The results of lake seasonal components with respect to GPS stations showed annual variations that may be linked to the hydrological cycle of the region, resulting in correlation coefficients exceeding 62% between GBGR-CS and O’Higgins/Chico lake. The correlations were smaller when comparing combined GPS data with lake Viedma data. This study illustrates the magnitude of the glacier changes taking place in Patagonia using geomatics techniques, and deserves further studies for a more detailed understanding of the complex interplay between climate change, glacier dynamics, the hydrological cycle, and crustal deformation.

Author Contributions

M.G.L. carried out the conceptualization, design, writing—original manuscript preparation, methodology; A.R. installed the GPS stations, collected data, writing and editing; M.D.: GPS station design, field data measurements, GPS processing, analysis; P.V.: ICESat data processing; M.C.: processing of seasonal data analysis; E.L.: GPS station design, field data measurements, analysis of the manuscript; M.G.: GPS and seasonal methodology, evaluation, and formal analysis of the manuscript; L.L. reviewed the article, and acquired funding. All authors have read and agreed to the published version of the manuscript.

Funding

This research has been funded by the Argentina’s National Science and Technology Agency (PICT 1995–2013 and PICTO 050-2016). Andrés Rivera was funded by FONDECYT 1171832 ANIDA, Chile.

Data Availability Statement

The study did not report any data.

Acknowledgments

We appreciate the collaboration in the field of the Andean Geomatics Lab members. We thank the Los Glaciares National Park for permits to explore the study area. We thank the collaboration of those who helped install and operate the GPS stations, including Carlos Fouilloux, Jonathan Garcés, Alejandro Silva, Sebastián Pulgar, and Jorge Hernández among many others. The collaboration of DGA and CECs is also appreciated.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Rivera, A.; Bown, F.; Napoleoni, F.; Muñoz, C. Balance de Masa Glaciar; Ediciones CECs: Valdivia, Chile, 2016; 203p. [Google Scholar]
  2. Millan, E.; Rignot, A.; Rivera, V.; Martineau, J.; Mouginot, R.; Zamora, J.; Uribe, G.; Lenzano, B.; De Fleurian, X.; Li, Y.; et al. Ice thickness and bed elevation of the Northern and Southern Patagonian Icefields. Geophys. Res. Lett. 2019, 46, 6626–6635. [Google Scholar] [CrossRef] [Green Version]
  3. IPCC. Summary for Policymakers. In Climate Change; The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, MA, USA, 2021; in press. [Google Scholar]
  4. Aniya, M.; Skvarca, P. Characteristics and variation of Upsala and Moreno glaciers, southern Patagonia. Bull. Glaciers 1992, 39–53. [Google Scholar]
  5. Leiva, J.C.; Lenzano, L.; Cabrera, G.A. Geodetic and glaciological work at Glaciar Chico, South Patagonian Ice Field. Glob. Planet. Change 2004, 59, 10–16. [Google Scholar] [CrossRef]
  6. Rivera, A.; Bown, F.; Acuña, C.; Ordenes, F. Chilean glaciers as indicators of climate change. Terra Glacialis 2008, 11, 193–207. [Google Scholar]
  7. Rivera, A.; Corripio, J.; Bravo, C.; Cisternas, S. Glaciar Jorge Montt (Chilean Patagonia) dynamics derived from photos obtained by fixed cameras and satellite image feature tracking. Ann. Glaciol. 2012, 53, 152. [Google Scholar] [CrossRef] [Green Version]
  8. Willis, M.; Melkonian, A.; Pritchard, M.; River, A. Ice loss from the southern Patagonia Ice Field, South America, between 2000 and 2012. Geophys. Res. Lett. 2012, 39, L17501. [Google Scholar] [CrossRef] [Green Version]
  9. Mouginot, J.; Rignot, E. Ice motion of the Patagonian Icefields of South America: 1984–2014. Geophisical Res. Lett. 2015, 42, 2661. [Google Scholar] [CrossRef] [Green Version]
  10. Moragues, S.M.; Lenzano, M.G.; Lo Vecchio, A.; Falaschi, D.; Lenzano, L.E. Surface velocities of Upsala Glacier, Southern Patagonian Andes using cross correlation satellite imagery: 2013–2014 Period. Andean Geol. 2018, 45, 87–103. [Google Scholar] [CrossRef] [Green Version]
  11. Jaber, W.; Rott, H.; Floricioiu, D.; Wuite, J.; Miranda, N. Heterogeneous spatial and temporal pattern of surface elevation change and mass balance of the Patagonian ice fields between 2000 and 2016. Cryosphere 2019, 13, 2511–2535. [Google Scholar] [CrossRef] [Green Version]
  12. Richter, A.; Groh, A.; Horwath, M.; Ivins, E.; Marderwald, E.; Hormaechea, J.L.; Dietrich, R. The rapid and steady mass loss of the patagonian icefields throughout the GRACE era: 2002–2017. Remote Sens. 2019, 11, 909. [Google Scholar] [CrossRef] [Green Version]
  13. Vacaflor, P.; Lenzano, M.G.; Vich, A.; Lenzano, L. Co-Registration Methods and Error Analysis for Four Decades (1979–2018) of Glacier Elevation Changes in the Southern Patagonian Icefield. Remote Sens. 2022, 14, 820. [Google Scholar] [CrossRef]
  14. Minowa, M.; Schaefer, M.; Sugiyama, S.; Sakakibara, D.; Skvarca, P. Frontal ablation and mass loss of the Patagonian icefields. Earth Planet. Sci. Lett. 2021, 561, 116811. [Google Scholar] [CrossRef]
  15. Glasser, N.F.; Harrison, S.; Jansson, K.N.; Anderson, K.; Cowley, A. Global sea-level contribution from the Patagonian Icefields since the Little Ice Age maximum. Nat. Geosci. 2011, 4, 303–307. [Google Scholar] [CrossRef]
  16. Coulson, S.; Lubeck, M.; Mitrovica, J.X.; Powell, E.; Davis, J.L.; Hoggard, M.J. The global fingerprint of modern ice-mass loss on 3-D crustal motion. Geophys. Res. Lett. 2021, 48, e2021GL095477. [Google Scholar] [CrossRef]
  17. Peltier, W.R. Global Sea level rise and glacial isostatic adjustment. Glob. Planet. Change 1999, 20, 93–123. [Google Scholar] [CrossRef]
  18. Peltier, W.R. Ice age paleotopography. Science 1994, 265, 195–201. [Google Scholar] [CrossRef]
  19. Lange, H.; Casassa, G.; Ivins, E.; Schroeder, L.; Fritsche, M.; Richter, A.; Groh, A.; Dietrich, R. Observed crustal uplift near the Southern Patagonian Icefield constrains improved viscoelastic Earth models. Geophys. Res. Lett. 2014, 41, 805–812. [Google Scholar] [CrossRef]
  20. Even-Tzur, G. GPS vector configuration design for monitoring deformation networks. J. Geod. 2002, 76, 455–461. [Google Scholar] [CrossRef]
  21. Segall, P.; Davis, J. GPS applications for geodynamics and earthquake. Annu. Rev. Earth Planet. Sci. 1997, 25, 301–336. [Google Scholar] [CrossRef] [Green Version]
  22. Johnson, H.; Agnew, D.C. Correlated noise in the geodetic time series. U.S. Geol. Surv. 1995, 102, 591–603. [Google Scholar]
  23. Khan, S.A. Surface Deformations Analyzed Using GPS Time Series. PhD Thesis, Danish National Space Center, Lyngby, Denmark, 2005. [Google Scholar]
  24. Davies, B.J.; Glasser, N.F. Accelerating shrinkage of Patagonian glaciers from the Little Ice Age (~AD 1870) to 2011. J. Glaciol. 2012, 58, 1063–1084. [Google Scholar] [CrossRef] [Green Version]
  25. Dietrich, R.E.R.; Casassa, G.; Lange, H.; Wendt, J.; Fritsche, M. Rapid crustal uplift in Patagonia due to enhanced ice loss. Earth Planet. Sci. Lett. 2010, 289, 22–29. [Google Scholar] [CrossRef]
  26. Richter, A.; Ivins, E.; Lange, H.; Mendoza, L.; Schröder, L.; Hormaechea, J.L.; Horwath, M. Crustal deformation across the Southern Patagonian Icefield observed by GNSS. Earth Planet. Sci. Lett. 2016, 452, 206–215. [Google Scholar] [CrossRef]
  27. Casassa, G.; Rivera, A.; Aniya, M.; Naruse, R. Current knowledge of the southern Patagonia icefield. In The Patagonian Icefields; Springer: Boston, MA, USA, 2002; pp. 67–83. [Google Scholar] [CrossRef]
  28. Aniya, M. Holocene glaciations of Hielo Patagónico (Patagonia Icefield), South America: A brief review. Geochem. J. 2013, 47, 97–105. [Google Scholar] [CrossRef] [Green Version]
  29. De Angeli, H. Hypsometry and sensitivity of the mass balance to changes in equilibrium-line altitude: The case of the Southern Patagonia Icefield. J. Glaciol. 2014, 60, 14–28. [Google Scholar] [CrossRef] [Green Version]
  30. Moragues, S.N.; Lenzano, M.G.; Rivera, S.A.; Oberreuter, J.G.; Vich, A. Characterization and reconstruction of the Agassiz land- slide using geospatial data. Southern Patagonia, Argentina. Andgeo 2021, 48, 557. [Google Scholar] [CrossRef]
  31. Rivera, A.; Lenzano, M.; Lanutti, E.; Moragues, S.; Lenzano, L.; Lenz, A.J. Vich Estudio de la profundidad del lago Viedma, Parque Nacional Los Glaciares, Argentina. GEOACTA 2022, 43, 4–6. [Google Scholar]
  32. Normandeau, J.; Meertens, C.; Bartel, B. GPS antenna monuments and mounts supported by UNAVCO: Options and Effectiveness. AGU Fall Meet. Abstr. 2008, 2008, G41B-0627. [Google Scholar]
  33. Durand, M.; Lannutti, E.; Lenzano, M.G.; Lenzano, L.E. Evaluación del uso de diferentes antenas y la influencia de la no verticalidad en mediciones GPS. GEOACTA 2017, 42, 1–12. [Google Scholar]
  34. Lee, J. GLAS_HDF Standard Data Product Specification. 2012. Available online: https://icesat.gsfc.nasa.gov/icesat/hdf5_products/docs/GLAS_HDF_SDP.pdf (accessed on 5 October 2019).
  35. Abshire, J.B.; Sun, X.; Riris, H.; Sirota, J.M.; McGarry, J.F.; Palm, S.; Yi, D.; Liiva, P. Geoscience Laser Altimeter System (GLAS) on the ICESat Mission: On-Orbit Measurement Performance. Geophys. Res. Lett. 2005, 32, 24028. [Google Scholar] [CrossRef] [Green Version]
  36. Smith, B.; Fricker, H.A.; Gardner, A.; Siegfried, M.; Adusumilli, S.; Csathó, B.; Holschuh, N.; Nilsson, J.; Paolo, F.; the ICESat-2 Science Team. ATLAS/ICESat-2 L3A Land Ice Height Version 3; NASA National Snow and Ice Data Center Distributed Active Archive Center: Boulder, CO, USA, 2020.
  37. Neumann, T.A.; Martino, A.J.; Markus, T.; Bae, S.; Bock, M.R.; Brenner, A.C.; Brunt, K.M.; Cavanaugh, J.; Fernandes, S.T.; Hancock, D.W.; et al. The Ice, Cloud, and Land Elevation Satellite–2 Mission: A Global Geolocated Photon Product Derived from the Advanced Topographic Laser Altimeter System. Remote Sens. Environ. 2019, 233, 111325. [Google Scholar] [CrossRef] [PubMed]
  38. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The Shuttle Radar Topography Mission. Rev. Geophys 2007, 45, 183. [Google Scholar] [CrossRef] [Green Version]
  39. Dach, R.; Hugentobler, U.; Fridez, P.; Meindl, M. Bernese GPS software version 5.0. Astron. Inst. Univ. Bern 2007, 640, 114. [Google Scholar]
  40. Kouba, J. Implementation and testing of the gridded Vienna Mapping Function 1 (VMF1). J. Geod. 2008, 82, 193–205. [Google Scholar] [CrossRef]
  41. Fritsche, M.; Dietrich, R.; Knöfel, C.; Rülke, A.; Vey, S.; Rothacher, M.; Steigenberger, P. Impact of higher-order ionospheric terms on GPS estimates. Geophys. Res. Lett. 2005, 32, L23311. [Google Scholar] [CrossRef]
  42. Bevis, M.; Brown, A. Trajectory models and reference frames for crustal motion geodesy. J. Geod. 2014, 88, 283–311. [Google Scholar] [CrossRef]
  43. Rebischung, P.; Griffiths, J.; Ray, J.; Schmid, R.; Collilieux, X.; Garayt, B. IGS08: The IGS realization of ITRF2008. GPS Solut. 2012, 16, 483–494. [Google Scholar] [CrossRef]
  44. Chen, B.; Bian, J.; Ding, K.; Wu, H.; Li, H. Extracting seasonal signals in GNSS coordinate time series via weighted nuclear norm minimization. Remote Sens. 2020, 12, 2027. [Google Scholar] [CrossRef]
  45. Wang, K.; Jiang, W.; Chen, H.; An, X.; Zhou, X.; Yuan, P.; Chen, Q. Analysis of seasonal signal in GPS short-baseline time series. Pure Appl. Geophys. 2018, 175, 3485–3509. [Google Scholar] [CrossRef]
  46. Huang, N.E.; Shen, Z.; Long, S.R.; Wu, M.C.; Shih, H.H.; Zheng, Q.; Yen, N.-C.; Tung, C.C.; Liu, H.H. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis. R. Soc. Proc. Math Phys. Eng. Sci. 1998, 454, 903–995. [Google Scholar] [CrossRef]
  47. Bao, Z.; Chang, G.; Zhang, L.; Chen, G.; Zhang, S. Filling missing values of multi-station GNSS coordinate time series based on matrix completion. Measurement 2021, 183, 109862. [Google Scholar] [CrossRef]
  48. Shen, Y.; Li, W.; Xu, G.; Li, B. Spatiotemporal filtering of regional GNSS network’s position time series with missing data using principle component analysis. J. Geod. 2014, 88, 1–12. [Google Scholar] [CrossRef]
  49. Yang, X.; Cheng, G.; Liu, H. Improved empirical mode decomposition algorithm of processing complex signal for IoT application. Int. J. Distrib. Sens. Netw. 2015, 11, 862807. [Google Scholar] [CrossRef]
  50. Huang, N.E.; Wu, Z. A review on Hilbert-Huang transform: Method and its applications to geophysical studies. Rev. Geophys. 2008, 46, RG2006. [Google Scholar] [CrossRef] [Green Version]
  51. Liu, Z.; Zhu, J.; Fu, H.; Zhou, C.; Zuo, T. Evaluation of the vertical accuracy of open global DEMs over steep terrain regions using icesat data: A case study over hunan province, china. Sensors 2020, 20, 4865. [Google Scholar] [CrossRef] [PubMed]
  52. Höhle, J.; Höhle, M. Accuracy Assessment of Digital Elevation Models by Means of Robust Statistical Methods. ISPRS J. Photogramm. Remote Sens. 2009, 64, 398–406. [Google Scholar] [CrossRef]
  53. Huber, P.J.; Ronchetti, E. Robust Statistics, 2nd ed.; Wiley Series in Probability and Statistics; Wiley: Hoboken, NJ, USA, 2009; ISBN 978-0-470-12990-6. [Google Scholar]
  54. Sánchez, L.; Drewes, H. Crustal deformation and surface kinematics after the 2010 earthquakes in Latin America. J. Geodyn. 2016, 102, 1–23. [Google Scholar] [CrossRef]
  55. Piñón, D.A.; Gómez, D.D.; Smalley, R.; Cimbaro, S.R.; Lauría, E.A.; Bevis, M.G. The History, State, and Future of the Argentine Continuous Satellite Monitoring Network and Its Contributions to Geodesy in Latin America. Seismol. Res. Lett. 2018, 89, 475–482. [Google Scholar] [CrossRef] [Green Version]
  56. Klos, A.; Bos, M.S.; Bogusz, J. Detecting time-varying seasonal signal in GPS position time series with different noise levels. GPS Solut. 2018, 22, 21. [Google Scholar] [CrossRef] [Green Version]
  57. Gómez, D.D.; Bevis, M.G.; Smalley, R.; Durand, M.; Willis, M.J.; Caccamise, D.J.; Casassa, G. Transient ice loss in the Patagonia Icefields during the 2015–2016 El Niño event. Sci. Rep. 2022, 12, 1–8. [Google Scholar]
  58. Malz, P.; Meier, W.; Casassa, G.; Jaña, R.; Skvarca, P.; Braun, M.H. Elevation and mass changes of the Southern Patagonia Icefield derived from TanDEM-X and SRTM data. Remote Sens. 2019, 10, 188. [Google Scholar] [CrossRef] [Green Version]
  59. Lenzano, M.G.; Lannutti, E.; Toth, C.; Rivera, A.; Lenzano, L. Detecting glacier surface motion by optical flow. Photogramm. Eng. Remote Sens. 2018, 84, 33–42. [Google Scholar] [CrossRef]
  60. Breitsprecher, K.; Thorkelson, D.J. Neogene kinematic history of Nazca–Antarctic–Phoenix slab windows beneath Patagonia and the Antarctic Peninsula. Tectonophysics 2009, 464, 10–20. [Google Scholar] [CrossRef]
  61. Russo, R.M.; VanDecar, J.C.; Comte, D.; Mocanu, V.I.; Gallego, A.; Murdie, R.E. Subduction of the Chile Ridge: Upper mantle structure and flow. Gsa Today 2010, 20, 4–10. [Google Scholar] [CrossRef]
  62. Rabassa, J.; Coronato, A.; Martinez, O. Late Cenozoic glaciations in Patagonia and Tierra del Fuego: An updated review. Biol. J. Linn. Soc. 2011, 103, 316–335. [Google Scholar] [CrossRef] [Green Version]
  63. Ávila, P.; Dávila, F.M. Heat flow and lithospheric thickness analysis in the Patagonian asthenospheric windows, southern South America. Tectonophysics 2018, 747, 99–107. [Google Scholar] [CrossRef]
  64. Ávila, P.; Dávila, F.M. Lithospheric thinning and dynamic uplift effects during slab window formation, southern Patagonia (45°–55° S). J. Geodyn. 2020, 133, 101689. [Google Scholar] [CrossRef]
  65. Mark, H.F.; Wiens, D.A.; Ivins, E.R.; Richter, A.; Ben Mansour, W.; Magnani, M.B.; Marderwald, E.; Adaros, R.; Barrientos, S. Lithospheric erosion in the Patagonian slab window, and implications for glacial isostasy. Geophys. Res. Lett. 2022, 49, e2021GL096863. [Google Scholar] [CrossRef]
  66. Ivins, E.R.; James, T.S. Bedrock response to Llanquihue Holocene and present-day glaciation in southernmost South America. Geophys. Res. Lett. 2004, 31, L24613. [Google Scholar] [CrossRef]
  67. Folguera, A.; Gianni, G.; Sagripanti, L.; Vera, E.R.; Colavitto, B.; Orts, D.; Ramos, V.A. Active Deformation, Uplift and Subsidence in Southern South America Throughout the Quaternary: A General Review About Their Development and Mechanisms. In Marine Isotope Stage 3 in Southern South America; 60 KA BP-30 KA BP; Springer: Cham, Switzerland, 2016; pp. 107–127. [Google Scholar]
  68. Richter, A.J.; Marderwald, E.R.; Hormaechea, J.L.; Mendoza, L.P.O.; Perdomo, R.A.; Connon, G.C.; Dietrich, R. Lake-level variations and tides in Lago Argentino, Patagonia: Insights from pressure tide gauge records. J. Limnol. 2016, 75, 62–77. [Google Scholar] [CrossRef] [Green Version]
  69. Pereira, A.; Cornero, C.; Matos, A.C.; Pacino, M.C.; Blitzkow, D. Detection of total water mass changes in the Patagonian glaciers area by satellite gravimetry. Geofísica Int. 2021, 60, 161–174. [Google Scholar] [CrossRef]
  70. Nicolas, J.; Verdun, J.; Boy, J.-P.; Bonhomme, L.; Asri, A.; Corbeau, A.; Berthier, A.; Durand, F.; Clarke, P. Improved Hydrological Loading Models in South America: Analysis of GPS Displacements Using M-SSA. Remote Sens. 2021, 13, 1605. [Google Scholar] [CrossRef]
  71. Yan, J.; Dong, D.; Bürgmann, R.; Materna, K.; Tan, W.; Peng, Y.; Chen, J. Separation of sources of seasonal uplift in China using independent component analysis of GNSS time series. J. Geophys. Res. Solid Earth 2019, 124, 11951–11971. [Google Scholar] [CrossRef]
  72. Carabajal, C.; Boy, J.P. Lake and reservoir volume variations in South America from radar altimetry, ICESat laser altimetry, and GRACE time-variable gravity. Adv. Space Res. 2021, 68, 652–671. [Google Scholar] [CrossRef]
  73. Blewitt, G.; Lavallée, D. Effect of annual signals on geodetic velocity. J. Geophys. Res. Solid Earth 2002, 107, ETG-9. [Google Scholar] [CrossRef]
Figure 1. Location of GPS stations used in the study in the SPI. The black line in the South American map at Patagonia denotes the asthenospheric window. In red square is represented the study area. The continuous reference GPS stations are indicated in yellow circles. On the right, the map shows the employed ICESat profiles belonging to the SPI subregion considered for spatial filtering (see Section 3.4). The hillshade on which the ICESat profile are mapped were made from the ALOS Global Digital Surface Model (AW3D30, https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm, accessed on 12 July 2020).
Figure 1. Location of GPS stations used in the study in the SPI. The black line in the South American map at Patagonia denotes the asthenospheric window. In red square is represented the study area. The continuous reference GPS stations are indicated in yellow circles. On the right, the map shows the employed ICESat profiles belonging to the SPI subregion considered for spatial filtering (see Section 3.4). The hillshade on which the ICESat profile are mapped were made from the ALOS Global Digital Surface Model (AW3D30, https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm, accessed on 12 July 2020).
Remotesensing 15 00584 g001
Figure 2. Estimated speed from GPS stations GRCS and GBCS.
Figure 2. Estimated speed from GPS stations GRCS and GBCS.
Remotesensing 15 00584 g002
Figure 3. (a) Velocity vectors of the continuous stations GBCS and GRCS, as well as the episodic stations PUMA and LPAS. (b) Altimetric profile of the GBCS, GRCS, and we added CHLT and UNPA located in the Patagonian plateau. (c) Up time series of the four stations denoting the slope of each of the GPS series.
Figure 3. (a) Velocity vectors of the continuous stations GBCS and GRCS, as well as the episodic stations PUMA and LPAS. (b) Altimetric profile of the GBCS, GRCS, and we added CHLT and UNPA located in the Patagonian plateau. (c) Up time series of the four stations denoting the slope of each of the GPS series.
Remotesensing 15 00584 g003aRemotesensing 15 00584 g003b
Figure 4. Time series GPS data with the application of interpolation to fill out the gaps.
Figure 4. Time series GPS data with the application of interpolation to fill out the gaps.
Remotesensing 15 00584 g004
Figure 5. Results of IMFs for seasonality analysis of GBGR-CS, CHLT, and VOGH stations. Note that the units on the abscissa axes are meters.
Figure 5. Results of IMFs for seasonality analysis of GBGR-CS, CHLT, and VOGH stations. Note that the units on the abscissa axes are meters.
Remotesensing 15 00584 g005
Figure 6. Results of frequency analysis of GBGR-CS, CHLT and VOGH stations. Upper: seasonal components in the GPS stations’ height time series derived from IMFs. Bottom: Fourier transform.
Figure 6. Results of frequency analysis of GBGR-CS, CHLT and VOGH stations. Upper: seasonal components in the GPS stations’ height time series derived from IMFs. Bottom: Fourier transform.
Remotesensing 15 00584 g006
Figure 7. Analysis of the seasonality between the temporal sequences of gauge data and the Up time series of the GPS stations under study.
Figure 7. Analysis of the seasonality between the temporal sequences of gauge data and the Up time series of the GPS stations under study.
Remotesensing 15 00584 g007
Figure 8. Elevation changes for O Higgins, Chico, Viedma, and Upsala glaciers in the 2000–2019 period.
Figure 8. Elevation changes for O Higgins, Chico, Viedma, and Upsala glaciers in the 2000–2019 period.
Remotesensing 15 00584 g008
Table 1. GPS occupation time.
Table 1. GPS occupation time.
SiteTypeHeight
(m a.s.l.)
YearStart DateEnd DateTotal Days
GBCSContinuous163020158 October30 November911
20164 January5 March
20179 October31 December
20181 January26 June
9 October31 December
20191 January28 May
28 July31 December
20201 January6 June
7 October23 November
GRCSContinuous1436201524 October31 December378
20161 January11 July
20179 October31 December
20181 January18 June
LPASEpisodic754201518 April20 April12
201619 April21 April
23 October25 October
201717 April19 April
PUMAEpisodic807201518 April20 April12
201618 April20 April
22 October25 October
201717 April19 April
Table 2. ICESat profiles database.
Table 2. ICESat profiles database.
File NameAcquisition
Date
Laser
Campaign
(I(ICESat-1))
/G GTN * Used
(I(ICESat-2))
Footprint Size
Major Axis Mean
± St. Dev.
Eccentricity Mean
± St. Dev.
GLAH14_634_2107_002_0043_0_01_000125 February 2004L2B89.52 ± 4.930.822 ± 0.045
         2107_003_0043_0_01_000126 May 2004L2C88.37 ± 19.120.892 ± 0.044
         2109_002_0043_0_01_000112 October 2004L3A55.79 ± 0.430.567 ± 0.043
         2111_002_0043_0_01_000127 February 2005L3B79.53 ± 11.550.753 ± 0.051
         2113_002_0043_0_01_000130 October 2005L3D52.04 ± 1.060.523 ± 0.010
         2115_002_0043_0_01_00012 March 2006L3E52.31 ± 1.600.483 ± 0.040
         2119_002_0043_0_01_000120 March 2007L3H55.61 ± 0.480.521 ± 0.019
         2121_002_0043_0_01_000111 October 2007L3I57.28 ± 0.570.590 ± 0.013
         2123_002_0043_0_01_000125 February 2008L3J58.66 ± 1.520.575 ± 0.036
         2125_002_0043_0_01_000112 October 2008L3K51.99 ± 1.120.611 ± 0.036
ATL06_20181126060903_08950109_005_0126 November 20181-2-3-4-5-617 m
       20190206150633_06130213_005_016 February 20191-2-3-5-6
       20190225014905_08950209_005_0125 February 20193-4-5-6
       20190508104610_06130313_005_018 May 20191-2-3-4-5-6
       20190526212840_08950309_005_0126 May 20191-2-3-4-5-6
       20190807062551_06130413_005_017 August 20191-2-3-4-5-6
       20190825170829_08950409_005_0125 August 20191-2-3-4-5-6
       20190905050156_10550413_005_015 September 20191-2-3-4-5-6
        20190923154435_13370409_005_0123 September 20196
       20191124124821_08950509_005_0124 November 20193-4-5-6
       20191205004146_10550513_005_015 December 20191-2-3-4-5-6
* GTN, Ground Track Number. This table was constructed using the metadata included in the HDF5 files and the official metadata table provided by the ICESat Science Investigator-led Processing System (I-SIPS), 2014.
Table 3. Regional GPS network. GPS continuous stations involved in the processing. The CHLT and VOGH stations were used for the seasonality analysis.
Table 3. Regional GPS network. GPS continuous stations involved in the processing. The CHLT and VOGH stations were used for the seasonality analysis.
StationCity—CountryHeight
(m a.s.l.)
PlateAgency
BRAZBrasilia, Brazil1106South American (SOAM)IGS
UNSASalta, Argentine1257South American (SOAM)IGS
FALKIslas Malvinas50South American (SOAM)IGS
LPGSLa Plata, Argentine29South American (SOAM)IGS
OHI2O’Higgins, Antarctica33Antarctica (ANTA)IGS
PALMEstación Palmer, Antarctica31Antarctica (ANTA)IGS
PARCPunta Arenas, Chile22South American (SOAM)IGS
CHLTSanta Cruz, Argentine485South American (SOAM)IGN-Ar
VOGH *Chile301South American (SOAM)IGS
* Note that the VOGH station was only used by validation of seasonality gauge.
Table 4. GPS continuous and episodic station network velocities.
Table 4. GPS continuous and episodic station network velocities.
SiteVel N
(mm a−1)
σN
(mm a−1)
Vel E
(mm a−1)
σE
(mm a−1)
Vel Up
(mm a−1)
σUp
(mm a−1)
GBCS7.870.503.332.1137.162.10
GRCS10.200.633.151.1233.032.14
LPAS8.701.0912.301.2334.005.03
PUMA5.401.0911.301.2320.305.03
Table 5. Results of the Pearson’s correlation values between the GPS signal from the stations and water levels.
Table 5. Results of the Pearson’s correlation values between the GPS signal from the stations and water levels.
Lake GaugesGBGR-CSCHLTVOGH
Chico0.67360.06030.4858
O’Higgins0.6269−0.07960.6431
La Leona0.37970.14290.6718
Viedma0.02640.22720.6365
Table 6. Results of ice thickness changes ( Δ h ) between the SRTM2000 and ICESat-1 and -2 data profile analyzed by accumulation and ablation zone.
Table 6. Results of ice thickness changes ( Δ h ) between the SRTM2000 and ICESat-1 and -2 data profile analyzed by accumulation and ablation zone.
ICESat DateICESat Data Quantity Δ h ^
(m)
Δ h   ^ Annual (m) Time
Elapsed
Time Elapsed
(Years)
Δ h ^
(m a−1)
Δ h ^   Annual
(m a−1)
Accumulation Zone
SRTM-ICESat-1 (2000–2004/8)
February 200436−6.23 ± 1.11−6.61 ± 0.694.04.31−1.56 ± 0.28−1.53 ± 0.16
May 200438−7.94 ± 1.944.25−1.87 ± 0.46
October 200431−5.19 ± 1.334.67−1.11 ± 0.28
February 200575−10.60 ± 0.64−10.55 ± 0.485.05.33−2.12 ± 0.13−1.98 ± 0.09
October 200599−10.47 ± 0.455.67−1.85 ± 0.08
March 200692−15.67 ± 1.58−15.67 ± 1.586.086.08−2.58 ± 0.26−2.58 ± 0.26
March 200751−13.68 ± 0.34−13.13 ± 0.667.087.38−1.93 ± 0.05−1.78 ± 0.09
October 2007117−13.07 ± 0.627.67−1.70 ± 0.08
February 200851−17.57 ± 1.07−16.78 ± 0.638.08.34−2.20 ± 0.13−2.01 ± 0.08
October 200868−16.19 ± 0.888.67−2.87 ± 0.10
SRTM-ICESat-2 (2000–2019)
26 November 20183771−24.47 ± 0.13−24.47 ± 0.1318.7418.74−1.31 ± 0.01−1.31 ± 0.01
6 February 2019671−18.09 ± 0.31−18.13 ± 0.0418.9319.31−0.96 ± 0.02−0.94 ± 0.01
25 February 20192187−10.85 ± 0.1718.99−0.57 ± 0.01
8 May 201923,036−19.61 ± 0.0619.18−1.02 ± 0.01
26 May 2019675−14.85 ± 0.3319.23−0.77 ± 0.02
7 August 20198898−22.61 ± 0.0819.43−1.16 ± 0.01
25 August 20193147−31.28 ± 0.1819.48−1.61 ± 0.01
5 September 20197361−1.80 ± 0.1419.51−0.09 ± 0.01
5 December 20192872−15.32 ± 0.2419.76−0.78 ± 0.01
Ablation zone
SRTM-ICESat-1 (2000–2004/8)
February 2004127−10.07 ± 0.59−10.54 ± 0.374.04.31−2.52 ± 0.15−2.45 ± 0.08
May 2004131−10.90 ± 1.054.25−2.56 ± 0.25
October 2004111−10.42 ± 0.704.67−2.23 ± 0.15
February 200547−22.26 ± 0.80−18.14 ± 0.585.05.33−4.45 ± 0.16−3.39 ± 0.11
October 200572−17.25 ± 0.535.67−3.04 ± 0.09
March 200632−30.94 ± 2.68−30.94 ± 2.686.086.08−5.09 ± 0.44−5.09 ± 0.44
March 200744−27.21 ± 0.37−40.49 ± 1.017.087.38−3.84 ± 0.05−5.49 ± 0.14
October 200727−43.54 ± 1.307.67−5.68 ± 0.17
February 200830−39.08 ± 1.04−33.75 ± 0.788.08.34−4.08 ± 0.17−4.05 ± 0.09
October 200848−16.19 ± 0.888.67−4.51 ± 0.12
October 200938−33.61 ± 1.73−33.61 ± 1.739.67 −3.48 ± 0.18−4.03 ± 0.18
SRTM-ICESat-2 (2000−2019)
26 November 201810,901−76.27 ± 0.07−76.27 ± 0.0718.7418.74−4.07 ± 0.01−4.07 ± 0.01
6 February 20191841−46.96 ± 0.19−64.97 ± 0.0418.9319.31−2.48 ± 0.01−3.36 ± 0.01
25 February 20193359−57.86 ± 0.1418.99−3.05 ± 0.01
8 May 201911,842−50.51 ± 0.0619.18−2.63 ± 0.01
26 May 20196859−93.08 ± 0.1019.23−4.84 ± 0.02
7 August 201910,754−74.33 ± 0.0719.43−3.83 ± 0.01
25 August 20197824−62.21 ± 0.1119.48−3.19 ± 0.01
5 September 20198549−53.83 ± 0.1319.51−2.76 ± 0.01
23 September 2019256−120.07 ± 0.4219.56−6.14 ± 0.02
6 November 2019897−118.63 ± 0.2419.73−6.01 ± 0.01
24 November 20195833−44.22 ± 0.1719.76−2.24 ± 0.01
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Lenzano, M.G.; Rivera, A.; Durand, M.; Vacaflor, P.; Carbonetti, M.; Lannutti, E.; Gende, M.; Lenzano, L. Detection of Crustal Uplift Deformation in Response to Glacier Wastage in Southern Patagonia. Remote Sens. 2023, 15, 584. https://doi.org/10.3390/rs15030584

AMA Style

Lenzano MG, Rivera A, Durand M, Vacaflor P, Carbonetti M, Lannutti E, Gende M, Lenzano L. Detection of Crustal Uplift Deformation in Response to Glacier Wastage in Southern Patagonia. Remote Sensing. 2023; 15(3):584. https://doi.org/10.3390/rs15030584

Chicago/Turabian Style

Lenzano, María Gabriela, Andrés Rivera, Marcelo Durand, Paulina Vacaflor, Micaela Carbonetti, Esteban Lannutti, Mauricio Gende, and Luis Lenzano. 2023. "Detection of Crustal Uplift Deformation in Response to Glacier Wastage in Southern Patagonia" Remote Sensing 15, no. 3: 584. https://doi.org/10.3390/rs15030584

APA Style

Lenzano, M. G., Rivera, A., Durand, M., Vacaflor, P., Carbonetti, M., Lannutti, E., Gende, M., & Lenzano, L. (2023). Detection of Crustal Uplift Deformation in Response to Glacier Wastage in Southern Patagonia. Remote Sensing, 15(3), 584. https://doi.org/10.3390/rs15030584

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop