GRACE Data Explore Moho Change Characteristics Beneath the South America Continent near the Chile Triple Junction

: The internal and external mass migration and redistribution of the Earth system are usually accompanied by changes in the gravity ﬁeld, and the Gravity Recovery and Climate Experiment (GRACE) has been proven to be able to effectively monitor and evaluate such changes. The Chile Triple Junction (CTJ) is the convergence point of the Nazca plate, the Antarctic plate and the South American plate. Subductions of different forms and rates in the north and south of the CTJ have varying degrees of impact on the surface and underground material changes of the South American plate. In this study, GRACE data are used in the estimation of the comprehensive mass changes in the South America Continent (SAC) Near the CTJ (~15 ◦ range). In addition, surface movement changes constrained by GNSS data cannot fully explain the GRACE results after deducting hydrological information, which indicates that residual signals might be attributed to mass changes beneath the crust, that is, the Moho interface deformation. After eliminating surface movement and hydrological signals from the comprehensive mass changes of GRACE, this study obtains the deep structural information and calculates the Moho changes of three signiﬁcant regions with rates of − 2.12 ± 0.67 cm/yr, 0.18 ± 0.19 cm/yr and − 6.46 ± 1.31 cm/yr, respectively. Results have demonstrated that the subductions of the Nazca plate and the Antarctica plate have an effect on the uneven deformation of the Moho interface beneath the SAC. The Moho beneath the SAC mainly shows a deepening trend, but it is uplifted in some areas north of CTJ. On the whole, the rate of Moho changes is greater in the south than in the north. The relationship between Moho changes and surface changes also indicates that a longer timescale may be needed for maintaining isostatic balance.


Introduction
The Chile Triple Junction (CTJ) is formed by the Subduction of the Chile ridge (Nazca-Antarctic spreading center) beneath the South American plate [1][2][3]. Since 17 Ma, a series of South Chile Ridge (SCR) segments parallel to the trench have subducted in the southern Patagonian Andes, leading the CTJ to migrate northward as a whole [2,[4][5][6]. Since about 6 Ma, three relatively short ridges of 6 Ma, 3 Ma, and 0.3 Ma collided closely in the Golfo de Penas area, the CTJ has come to its present position of about 46 • 30 west of the Taito Peninsula [7]. The Antarctic plate (AN), to the south of the CTJ, subducts along the E-W direction almost vertically (with the ridge) at a rate of~2 cm/yr, while the Nazca plate (NZ), to the north of the CTJ, subducts obliquely in the trench N-80 • direction at a rate of~7 cm/yr [6][7][8][9][10][11][12]. Subduction in the northern part of the CTJ is associated with the formation of lower South American topography (~2300 m elevation), a small number of Cenozoic molasse sedimentary rocks in the east, and results in the lack of a fold and thrust. Meanwhile, subduction in the southern part of the CTJ is associated with higher topography (~4000 m elevation) with exhumed Patagonian batholith and pre-Jurassic rocks, the late Cenozoic molasse deposits and the development of a fold and thrust belt which has been shortened by 25-45 km [13,14]. The continuous subduction, with apparently different forms and rates, results in different degrees of surface changes due to internal mass redistribution in the north and south of the South American continent (SAC). At present, many studies have explored the external and internal geoscience issues of the SAC, such as the impact of the subduction of the NZ on the lithospheric erosion of the SAC, the changes of horizontal and vertical surface displacement, the distribution of stress and strain fields in Patagonia and the crustal thickness and Moho depth of the SAC, all of which have obtained satisfactory results. [15][16][17][18][19][20][21][22][23].
In addition to the various structural information of mass changes, the SAC also yields complex hydrological information, particularly the glacier changes in the Patagonia region. The Patagonian ice sheet is the second largest ice body in the Southern Hemisphere and is covered by numerous glaciers [24,25]. In recent years, many scholars have focused on the study of the mass changes of Patagonia glaciers, and their results show that the glaciers have undergone an alarmingly rapid retreat, which has an important impact on the flow of lakes and rivers on the Earth's surface [25][26][27][28][29]. Besides, the rapid global climate change (drought, flood and extreme weather, etc.) will also affect the hydrological balance of the SAC. Therefore, in order to properly understand the Earth system mass changes in this region, scientifically accurate monitoring is much needed for both the underground tectonic activities and the surface hydrological changes.
With the continuous development of modern geodetic and remote sensing technologies (such as Global Satellite Navigation System (GNSS), Interferometric Synthetic Aperture Radar (InSAR), Gravity Recovery and Climate Experiment (GRACE) and Remote Sensing Satellites, etc.), these technologies provide more and more observation data are available for the Earth science research. At present, these data have been widely used to explain the Earth surface deformation and the underground mass migration as well as other problems and important research results have been achieved. In particular, the GRACE gravity satellite, jointly developed and launched by the National Aeronautics and Space Administration (NASA) and the German Aerospace Center, has been observing the surface mass distribution and subsurface mass flow of the entire Earth system at a global regional scale for more than 10 years, and these observations provide explanations for processes such as polar ice, soil moisture (SM), terrestrial water storage (TWS), large-scale underground mass migration, and ocean mass dynamics [25,[30][31][32][33][34][35][36][37][38][39][40][41]. Mikhailov et al. [42] inquired whether it was possible to identify the temporal variation of gravity field caused by tectonic processes (TECs) (geodynamic signals) from satellite gravity data (GRACE and future gravity satellites). By studying subduction-driven mechanics, they found that satellite gravity data could not only be used to detect and distinguish the geodynamic signals generated by subduction-driven mechanics but could also be used to study the gravity of other geodynamic targets over time. It has been confirmed that some TECs can be detected by GRACE at present, such as the coseismic and post-seismic deformation generated by the 2004 Sumatra-Andaman earthquake, the 2010 Maule earthquake, the 2011 Tohoku-Oki earthquake, and the crustal subsidence caused by the 2013 Okhotsk deep-focus earthquake [43][44][45][46] as well as the mass changes in glacial isostatic adjustment (GIA) regions like Fennoscandia, Antarctica and Greenland [31,[47][48][49][50]. Yi et al. [51] combined GRACE observations and numerical models to analyze and explain the changes of Moho and crustal thickness in eastern Tibet. Jiao et al. [52] extracted tectonic signals of the Tibetan Plateau by using GRACE data, GPS data and hydrological model and further interpreted the tectonic signals as the result of the Moho changes. Rao and Sun [53] considered the impact of erosion and sedimentation on the Tibetan Plateau documented by the previous works and used GRACE data to conduct a further research on the changes of Moho interface and crustal thickness beneath the Tibetan Plateau. Previous studies have manifested that GRACE observations, combined with other Earth science data, can provide effective estimates and explanations for Earth surface information changes caused by subsurface mass migration.

Motivation and the Study Objective
Although many studies have examined crustal deformation, lithospheric erosion, crustal thickness and other issues beneath the SAC, we have not yet fathomed the rate of Moho uplift (subsidence), nor other subsurface processes that may influence the uplift (subsidence). Considering the impacts on the tectonic changes beneath the SAC caused by subduction of different rates and forms between northern and southern CTJ, this study focuses on the Moho changes beneath the SAC near the CTJ, and the factors affecting Moho changes also are discussed in the study area. By dint of existing research techniques, GRACE observations, hydrology data, remote sensing data, precipitation data and GNSS data are combined to reasonably explore the Moho changes caused by subsurface mass migration in South America. Firstly, this study calculates the comprehensive mass changes (including tectonic and hydrological changes) in the study area using GRACE data and then figures out the effects of other mass sources changes, such as the loading effect, crustal movement, GIA, glaciers, SM and groundwater (GW) from different data sources. After removing these mass changes from GRACE results, the residual mass changes caused by interior mass migration of the Earth can be obtained. Finally, the Moho depth changes can be estimated by using the corresponding formula transformation.
The scope of the study area (the SAC (yellow)) locates roughly in 80 • W-60 • W and 54 • S-39 • S, as shown in Figure 1. The main content of this paper is as follows. Section 2 explains why we conduct this study. Section 3 introduces data sets and data processing methods. Section 4 shows the results. Section 5 shows the discussion. Finally, conclusions are given in Section 6.

GRACE Data Processing
The study uses the GRACE Level 2 (Release 06) monthly gravity solutions, which are provided by the Center for Space Research (CSR) at the University of Texas, Geo-ForschungsZentrum (GFZ) in Potsdam and the Jet Propulsion Laboratory (JPL) (all data can be obtained from http://icgem.gfz-potsdam.de/home, accessed on 8 November 2020). Gravity solutions are spherical harmonic coefficients (SHs) solutions of order 60, and the effects of non-solid tidal atmosphere, high-frequency ocean signals, various tides, Earth tides, solid polar tides and other planetary gravity disturbances are deducted [54]. The average value from January 2003 to June 2017 is deleted from the monthly SHs to obtain the monthly gravity change value.
The GRACE gravity solutions do not provide the Degree-1 SHs of the geocentric motion (due to the use of geocenter as the reference origin in gravity field definition [55]), so we add back the geocentric correction term calculated by Swenson et al. [56]. The Degree-2 zonal terms (C 20 ) of GRACE data can be replaced by satellite laser ranging observations to improve the accuracy [57,58]. Due to the limitations of satellite design and observation data quality, the obtained data have large noise in the high-degree SHs, and the correlations between different degree SHs lead to "strips" errors in GRACE data [59]. Therefore, the 300 km Gaussian smoothing filter is used to reduce the noise of higher degree SHs, and the sliding window algorithm proposed by Swenson and Wahr [59] is adopted to suppress the strips. The numerical model proposed by Geruo et al. [60] is used to adjust the influence of GIA. Finally, the SHs are transformed into the monthly changes in the equivalent water height (EWH) anomaly (Equation (1)) [30], and the trends are obtained by using the least squares method.
2n + 1 1 + k n (∆C nm cos mλ + ∆S nm sin mλ)P nm (cos θ) (1) where h w is EWH; (θ, λ) is the latitude and longitude of ground points; (∆C nm , ∆S nm ) are the changes of SHs with degree n and order m; P nm (·) is the fully normalized associated Legendre function with degree n and order m; k n is the Load Love number of degree n; ρ e is the average density of solid Earth; ρ w is the density of water; R is the average radius of the Earth. The above processing steps can reduce the noise and strips errors in the GRACE data, but they may also attenuate the real signal. Therefore, signal leakage and signal attenuation mainly caused by filtering and truncation should be recovered [37,61]. The iterative recovery (ITER) method proposed by Chen et al. [62] is used to recover GRACE observed signals, that is, to restore the original model by iteratively adjusting the simulation model to approximate the observed values. The specific steps are as follows: (1) The EWH observation value f 0 with an initial resolution of 1 • after smoothing in the study area is regarded as the original mass rate m 0 ; (2) The original mass rate m 0 is expanded into SHs, truncated at degree 60, and then the 300 km Gaussian smoothing is applied. In this way, the forward model apparent mass rate m 1 can be acquired; (3) The difference between the apparent rate observed by GRACE and the forward model apparent mass rate (m 0 − m 1 ) is added to the original model m 0 ; (4) Repeat steps (2) and (3) until there is an agreement between the apparent rate observed by GRACE and the forward model apparent mass rate, or a certain number of iterations are reached. Usually, to speed up the convergence rate, a factor of 1.2 is added in Step (3), that is, m 0 = m 0 + 1.2 × (m 0 − m 1 ).

GLDAS Hydrological Model
Global Land Data Assimilation System (GLDAS) is a global high-resolution hydrological model jointly established by NASA and the National Centers of Environmental Prediction (NCEP). This hydrological model includes four version models (NOAH, VIC, CLSM and MOSAIC). For data acquisition reasons, the data sets provided by NOAH, VIC, and CLSM are used in this study. Here, the average snow water equivalent (SWE) of these three versions of models are used as the final SWE estimation. The period analyzed is from January 2003 to June 2017, the range of longitude and latitude covered is 90 • N-60 • S and 180 • E-180 • S, the spatial resolution is 1 • × 1 • grid, and the time resolution is one value per month [63].

CPC Hydrological Model
The Climate Prediction Center (CPC) of the National Oceanic and Atmospheric Administration (NOAA) provides a dataset of CPC hydrological models [64]. Compared with the GLDAS model, the CPC model uses more measurement data so that it greatly enhances the usability of real-time data [53]. Therefore, the SM data provided by CPC data is adopted in data processing. The spatial resolution of the data is 0.5 • × 0.5 • , the temporal resolution is one value per month, and the period is from January 2003 to June 2017.

Water GAP Hydrological Model
The Water GAP Hydrological Model (WGHM) is a global hydrological model, which is designed to assess water resources changes on a global scale [65,66]. The latest version of the model currently released is Water GAP2.2d, which is established based on the Water GAP2.2 version released in 2014 to achieve model calibration, data update and algorithm improvement [67]. GW data, lake water storage (LWS) data and SWE data are provided by it in this study. The spatial resolution of the data is 0.5 • × 0.5 • , the time resolution is one value per month, and the period is from January 2003 to November 2016.

Precipitation Data
As the only input variable in the water cycle, precipitation plays an important role in the global and regional changes of surface water, SWE, GW, etc. The precipitation data in this study are drawn from the PRECipitation REConstruction over Land (PREC/L) of the NOAA Physical Sciences Laboratory (PSL) [68]. The spatial resolution of the data is 1 • × 1 • , the time resolution is one value per month, and the period is from January 2003 to June 2017.

Remote Sensing Data Results for Glaciers
Glaciers are an important component of the water cycle and one of the hydrological signals that should be considered in this study. Glacier change information is mainly manifested as mass changes of the Patagonia ice sheet in the study area. Dietrich et al. [18] pointed out that the mass changes signal of the Patagonia ice sheet calculated by GRACE must be divided into (negative) mass changes caused by glaciers decline and (positive) mass change caused by crustal uplift. Therefore, accurate monitoring of glaciers mass changes is crucial to the data processing results. Foresta et al. [28] calculated the glaciers mass changes of the Patagonia ice sheet from April 2011 to April 2017 by using cryosat-2 strip radar altimetry data. The Patagonia ice sheet was divided into nine sub-blocks and the glacier mass change rate of each sub-block was obtained. The total glacier mass budget of the Patagonian ice field was −21.29 ± 1.98 Gt/yr based on the mass changes of nine glacier sub-areas in their study. Here, the glaciers mass changes rates of the nine sub-regions are evenly distributed in a regular grid of 1 • × 1 • , respectively.

Vertical Velocities Caused by Surface Movement
The Global Position System (GPS) has been widely used in the field of geodynamics since the mid-1980s. It is mainly used to monitor plate movement and crustal movement deformations at global or regional scales. At present, the accuracy of GPS monitoring plate movement and crustal deformation can reach 1~2 mm/yr in the horizontal direction and 2~4 mm/yr in the vertical direction. The distribution of GPS stations in the SAC is irregular, as most of them locate in the northern part of South America, and the number of available stations in the southern part is limited. Richter et al. [20] provided some useful GPS site data in the study area, and Nevada Geodesy Laboratory (NGL, http: //geodesy.unr.edu/, accessed on 8 May 2021) also set up a corresponding number of GPS sites in the southern region of South America. So, the GPS vertical velocities data are provided by them. The Kriging interpolation method is used to obtain the regular grid vertical velocity field data in the study area, which is usually suitable for data interpolation with uneven distribution [51,52,69]. The deformation characteristics of GPS observations are usually the superposition of GIA, tectonic activity, and surface loading deformation [20,53]. We need to separate or correct surface loading information and GIA changes to obtain mass changes caused by tectonic activities. Surface loading signals are usually caused by changes in hydrological information, and the effect of loading deformation usually reflects changes in gravity or displacement. Here, the data provided by GRACE is used to simulate vertical crustal displacement due to redistribution of Earth's surface mass [70][71][72][73], as shown in Equation (2). The vertical displacement change caused by GIA is estimated by using the model established by Geruo et al. [60].
where def. is elastic deformation; h n and k n are the Load Love numbers [74,75]. The GPS velocities used are the results of removing the surface loading and GIA's contribution, so they reflect the plastic deformation of the surface, that is, when assuming that the crustal density is constant, the mass change caused by the surface deformation may be related to its volume change [51,52,72,76].

Hydrologic and Remote Sensing Data Processing Methods
In order to compare with GRACE results, GLDAS data, CPC data, Water GAP data and glacier data are processed by spherical harmonic expansion, truncated at order 60, and then the 300 km Gaussian smoothing filter is applied. The ITER method is used to restore the signal, and the least squares method is used to obtain the mass trend.
It should also be noted that due to the difference in time periods of each observation model, the lengths of data sets results obtained are not consistent. But after calculation, we find that all of these data show a trend of stable change. Therefore, the final results should reflect the tectonic trends of the study area even if the time periods are different.

Tectonic Information Acquisition
The mass changes observed by GRACE usually involve mass redistribution in the surface and interior of the Earth, which can be divided into hydrologic mass changes and TECs. Hydrological factors mainly include changes in SM, SWE, GW, LWS, glacier changes and other signals, while TECs usually lead to the vertical movement of the crust, including surface uplift, surface erosion and sedimentation and Moho depth changes, etc. They generally satisfy the following relationships [52,53]: Equation (3) indicates that the gravity changes observed by GRACE (∆GRACE) mainly include hydrologic (∆Hydrological) and tectonic (∆Tectonic) changes. Equation (4) gives the main components of hydrological signals, where ∆SM is the change of soil moisture; ∆SWE is the change of snow water equivalent; ∆GW is the change of groundwater; ∆LWS is the change of lake water storage and ∆Glacier is the change of glacier mass. This study is interested in the changes of tectonic information. Hence, combining the Equation (3) with the Equation (4), we can obtain this by deducting hydrological information from the results of GRACE calculations.

Mass Changes Caused by Surface Movement
Through the Kriging interpolation method mentioned in Section 3.5, the interpolated vertical displacement change of the study area can be obtained based on the existing GPS data. Jiao et al. [52] showed that the mass change (EWH) caused by the surface deformation of the grid could be calculated by the following Equation (5): where ∆h water represents the EWH of grid mass variation; ∆h gps is the average vertical displacement change of grid calculated by GPS; ρ crust is the crustal density of grid, data sets come from the CRUST1.0 global crustal model (https://igppweb.ucsd.edu/~gabi/ crust1.html, accessed on 30 May 2021) [77], and we weight it according to the thickness to obtain the average density; ρ water is the density of water, we take 1 g/m 3 . Through the above Equation (5), the vertical mass change observed by GPS can be converted into EWH.

Estimation of Moho Interface Changes
The tectonic information obtained by the calculation method provided in Section 3.7 mainly includes surface uplift (subsidence) changes, sedimentation and erosion changes, GIA effect and Moho changes [53]. The method of Section 3.8 is used to estimate the mass change caused by surface deformation (EWH), and then the surface deformation information will be deducted from the tectonic information to get the mass changes caused by Moho changes (EWH). The Moho fluctuation changes (EWH) can be considered to be caused by the density difference between the crust and the mantle. Under the Bouguer plate assumption, the following Equation (6) can be used to establish the relationship between mass changes ∆h water (EWH) and Moho interface changes ∆h Moho [52]: where ρ water is the density of water; ∆ρ is the difference of density between the mantle and the crust.

The Impact of Filtering and Truncation on Data
In order to maintain consistency with the GRACE data processing results, truncation and filtering are performed on the hydrological data, glaciers data, GPS data, etc. Usually, truncation and filtering will cause signal leakage and signal attenuation, so we need to understand how big the impact is-whether the impact can be ignored, or that the signal must be restored. Here, this study makes a simulation to show the effect of truncation and filtering on the original data. Assume that the 5 • × 5 • range of 50 • S-46 • S and 73 • W-69 • W produces a 1 cm water layer fluctuation in the study area, as shown in Figure 2a. The signal of this water layer is expanded to 60-order SHs, as is shown in Figure 2b. On the basis of Figure 2b, 300 km Gaussian smoothing filter is performed, and the result is shown in Figure 2c. The results after truncation and filtering are expanded along the black dotted line shown in the figure, and the cross-sectional value is shown in Figure 2d. It is obvious that after truncation and filtering, the signal strength becomes weaker and the spatial influence range compared with the original state appears enlarged, which means that the signal does attenuate and leak. The signal attenuation caused by truncation is about 2~10%, and the attenuation of filtering is about 50~60%. Therefore, in order to obtain the original model data, as shown in Figure 2a, the inversion method is needed to perform signal recovery.

GRACE Observations
There are multiple signal sources (such as lakes, glaciers, soil water, structures, etc.) in the study area, and the mass change information usually observed by GRACE is the comprehensive result of these signal sources. The GRACE SHs solutions provided by CSR, GFZ and JPL are used to calculate the long-term mass trends, and the results are shown in Figure 3. The difference between the results in the study area is within the range of ±0.7 mm/yr, which shows that the accuracy of the solutions provided by different institutions is similar.
Considering the effects of signal attenuation and leakage in the GRACE data results, the unconstrained ITER method is used to restore the GRACE signal in the study area. Jiao et al. [52] showed that the unconstrained ITER method required neither signal prior information nor consideration of the spatial distribution of each independent signal, which could be used to realize the recovery processing of GRACE signals in the study area. In order to reduce the difference of GRACE data provided by different institutions, the average value of mass rate change calculated by the products of the three institutions is regarded as the initial value of the ITER method, and the recovery results of GRACE signals are shown in Figure 4.
The ITER method is used to recover GRACE observations, and the specific process is as follows: Firstly, the average GRACE observations ( f 0 ) provided by CSR, GFZ and JPL are regarded as the initial mass rates (m 0 ), as shown in Figure 4a. Secondly, the initial mass rates (m 0 ) are expanded into SHs, truncated at degree 60, and then the 300 km Gaussian smoothing is applied. The forward model apparent mass rates (m 1 ) can be obtained in this way. Thirdly, the differences between the initial mass rates and the forward model apparent mass rates (m 0 − m 1 ) are added to the initial mass rates m 0 . Fourthly, repeat steps (2) and (3) until the observations mass rates ( f 0 ) and the forward model apparent mass rates (m 1 ) (Figure 4c) are agreement or a certain number of iterations are reached. Finally, the true mass rates can be obtained, as shown in Figure 4b. Figure 4d represents the differences between GRACE observations ( f 0 ) and the forward model apparent mass rates (m 1 ). By calculating the differences between the observed and modelled apparent mass rates, we find that the largest positive difference in the study area is less than 0.07 cm/yr, and the smallest negative difference is greater than −0.09 cm/yr, which indicates that the signal has been better recovered in the study area. Comparing Figures 4a and 4b, it can be seen that the recovered signal is significantly enhanced, and the signal appears to converge to the inside. This also shows that the ITER method can effectively suppress signal leakage while recovering the signal.  The "true mass rates" recovered by the ITER method; (c) The apparent mass rates calculated from the "true mass rates" recovered from (b); (d) The differences between the observed apparent mass rates and the modeled apparent mass rates (i.e., (a-c)). Note the different color scale.
Based on the GRACE observations, this study selects point A and B as research representatives because of the significant signal changes they indicate and then calculates the time series of the mass changes from 2003 to 2017, respectively. Among them, point A is located on the Patagonian ice sheet, and point B is located in south-central Chile. Figure 5 shows the time series calculated by GRACE for these two points (including seasonal, nonseasonal and long-term trends) and its relationship with precipitation data, hydrological data, GIA data and GPS data. In order to distinguish different change information, seasonal signal and long-term signal are fitted with a mathematical model containing annual, semiannual and linear trend term, and the specific expression is as follows [78]: where M(t) is the mass change; t is time; a is the constant; b is the long-term change; c k , ω k and ϕ k are amplitude, frequency and phase, respectively; k = 1 represents annual change, k = 2 represents semi-annual change; ε(t) represents the modeled residual term. In particular, a 161-day Tidal Aliasing term is added to the above Equation (7) at point A after referring to Chen et al. [25] analyzed glacier mass change in Patagonia.   Figure 5a. The Patagonian Ice Sheet, where point A is located, is the second-largest ice body in the southern hemisphere. This area is mainly a gathering area of glaciers. Therefore, the mass change calculated by GRACE mainly reflects the glacier change, and the long-term trend of point A shows a negative increase, which is consistent with previous reports that the Patagonia ice sheet had experienced a general glacial decline over the past decade [25,27,[79][80][81][82]. At the same time, point B is located in central-southern Chile. This area contains a lot of information changes (such as precipitation, evapotranspiration, and population factors, etc.), and these factors have been observed by GRACE. The time series analysis of precipitation data and GRACE results is shown in Figure 5b. It can be seen that the annual precipitation in the area of point A is uniform, but there is no obvious interaction between the precipitation data and the GRACE results, and precipitation change has a weak direct influence on GRACE results. Chen et al. [25] calculated that the total mass loss rate of the Patagonian ice sheet was about −24.3 ± 4.3 km 3 /yr. Among them, the change in TWS was about −5.4 km 3 /yr and the postglacial rebound (The contribution of PGR) of the Earth mantle was about 9 ± 8 km 3 /yr, and the glacier melting rate −27.9 ± 11 km 3 /yr (all measured by EWH). This also shows that the GRACE observation of the area, where point A is located, reflects the comprehensive changes of multiple different mass sources. The time series of mass change calculated by GRACE in the area, where point B is located, is in good agreement with the precipitation data in most cases. The GRACE observation also shows an increasing trend during the period when the precipitation is increasing. At the same time, we also observe a feature that the winter precipitation is greater than that in summer. Cobain climate classification indicates that central Chile usually has a Mediterranean climate. Under the influence of this climate, there is usually an obvious precipitation peak in winter and a scarcity of precipitation in summer, which also explains the abnormal changes of precipitation and GRACE results in this area. Compared with GRACE mass change rate, the comprehensive mass change rate of hydrology, GIA and GPS is smaller than that of GRACE, as shown in Figure 5c.

Hydrological Results
The acquired hydrological data are shown in Figure 6, where Figure 6a is the change data of SM content obtained from the CPC data, Figure 6b is the change result of LWS from WGHM, and Figure 6c is the result of the change in GW provided by WGHM. The glacier mass changes in Figure 6d is the result of using CryoSat-2 strip radar altimetry. SWE result from WGHM is shown in Figure 6e, and the SWE result in Figure 6f is provided by the GLDAS model. It is found that the glacier mass shows an obvious change, while the other data (except the SWE provided by GLDAS) show weak changes. Meanwhile, we also find that SWE data provided by GLDAS show strong signals in Patagonia and even reach the mass changes rates of glaciers, but the rates seem unreasonable and unreliable. Jiao et al. [52] also observed a similar phenomenon in the SWE data obtained by using the GLDAS hydrological model in the Tibetan Plateau-the amplitudes of SWE signals exceeded those observed by GRACE and they thought it unreliable. Therefore, this study abandons the SWE data provided by GLDAS and only use the SWE results provided by WGHM.  According to Equation (2), TECs can be obtained by subtracting the mass changes of hydrological information sources from the mass changes observed by GRACE. Figure 7 shows the process of subtracting hydrological signals from GRACE observations. Figure 7a shows the superposition result of all the hydrological signals shown in Figure 6 (this study considers all the hydrological signals that may exist), Figure 7b is the mass changes rates observed by GRACE, and Figure 7c shows the GRACE residual signal after subtracting the hydrological signal. There are four distinct signal extreme regions in the study area, and the changes from north to south are negative, positive, negative and positive, respectively. The redistribution of surface mass is mainly affected by various forms of terrestrial water change and changes in atmospheric pressure and ocean transports at different time scales [70,83]. The GRACE satellite data can be used to detect changes in gravity caused by mass migration on the Earth's surface or the interior. When only considering the regional gravity changes due to surface mass migration, we subtract the surface mass changes from the results of GRACE observations, and the remaining signal size should approach zero. However, after deducting the influence of surface mass changes, there are four abnormal signal areas in the study area, and some of the signal values even exceed those observed by GRACE. Therefore, the residual signals of GRACE are considered to come from the changes generated by TECs.

Mass Changes Caused by Surface Movement
The surface movement can be obtained with the help of GPS observations, and the accuracy is usually in millimeters or even sub-millimeters, which can meet the accuracy requirements of the research results. Meanwhile, it can be converted into the corresponding mass changes according to the introduction in Section 3.8.
Firstly, the collected GPS data are used to obtain the vertical velocities in the study area according to the interpolation method in Section 3.5, as shown in Figure 8a (blue dots represent the GPS point distribution). Secondly, GRACE data is used to calculate the loading deformation effect caused by hydrological factors on the Earth's surface, and the results are shown in Figure 8b. At the same time, this study considers that the vertical velocities results calculated by GPS include the GIA influence and Figure 8c shows the uplift changes caused by the GIA model established by Geruo et al. [60]. After deducting the loading deformation information and the GIA deformation effect from GPS, the free-loading-GIA uplift results are obtained, as shown in Figure 8d. CRUST1.0 Global Crust model provides crustal density data for the study area, which is obtained by using a weighted average method, as shown in Figure 8e. Finally, according to the free-loading-GIA uplift results in Figure 8d and the crustal density data in Figure 8e, Equation (5) is used to obtain the spatial distribution results of free-loading-GIA uplift mass changes information (EWH).

Moho Changes Beneath the Study Area
The residual result after deducting the hydrological signal from the GRACE observation is shown in Figure 9. Three regions with significant changes in the residual signal are selected to obtain the mass changes rates, as shown in Figure 9a  Comparing the results of these two sets of data, it can be seen that the surface mass change rate obtained by GPS is only 10% the result of GRACE residual signal in region 1; The change of surface mass accounts for about 67% of the residual signal of GRACE in region 2; The change of surface mass obtained by GPS exhibits a converse trend to the GRACE residual signal in region 3. Therefore, it also means that the surface movement cannot fully explain the tectonic changes in the SAC. Jiao et al. [52], Sun et al. [76], Braitenberg and Shum [84] and Chen et al. [85] had studied that TECs usually consisted of three parts: (1) surface movement; (2) crustal density changes; (3) Moho changes. It is assumed that the crust is in a quasi-equilibrium state after millions of years of development. Therefore, considering the characteristics of the large density difference of the Moho interface, the Moho changes through gravity observation data can be studied under the assumption that the Earth's density remains unchanged [52,76]. Figure 9c is the deformation (EWH) estimated by the GIA model. Therefore, Moho changes (EWH) in the study area can be obtained by deducting the mass changes rates generated by surface movement and the uplift information estimated by the GIA model from the GRACE results after deducting hydrological signals, as shown in Figure 9d. Hydrologic signals and mass change signals caused by surface movement are deducted from GRACE observations, and the final residual signal results (EWH) are converted by Equation (6) to obtain the Moho interface depth changes in the study area, as shown in Figure 10. Figure 10a shows the remaining signals of GRACE in the study area. During data processing, we truncate and filter the data when necessary and then deduct it from GRACE results so that the final estimation of the Moho changes (EWH) is also the result that has been truncated and filtered. Consequently, the recovery of the real signal is needed. The ITER method is still used, and Figure 10b is the result after recovery. Figure 10c is the spatial distribution of the density difference between the mantle and the crust in the study area, where the density data of the crust and mantle are both from the CRUST1.0 Global crust model. Figure 10d is the spatial distribution of final Moho depth changes in the study area. Among them, Moho depth presents significant rates of change in the northeast area of the study area, central Chile region and Patagonia region, showing a decrease, an increase and a decrease trend, respectively. The other regions also show varying degrees of Moho depth changes rates, which indicate that different degrees of tectonic activities are taking place beneath the SAC. The Moho depth changes rates of the three significant signal regions are calculated, and the results are shown in Table 1.

Interpretation of Moho Deformation
The uneven deformation of the Moho indicates that there are different TECs taking place in different parts of the SAC, and the factors affecting Moho changes are also different. The factors contributing to the Moho depth changes of the three significant signal regions are analyzed in the study area.
By analyzing the results in region 1, we notice that the hydrological signal superimposed in the area shows an increasing trend as a whole, however, the GRACE observation presents a negative trend in the area as a whole. We suppose that the surface loading caused by the increase of hydrological mass in this area has a certain impact on the deformation of the crust and the activities of underground tectonics.
Moho uplift may indicate a thinning trend in the crust. After analyzing the possible reasons for the uplifting of the Moho, Rodriguez and Russo [23] and Yi et al. [51] discussed potential unloading processes as follows: (1) surface erosion, (2) glacial isostatic adjustment (GIA) due to glacier shrinkage after the Last Glacial Maximum (LGM), and (3) subduction erosion. Firstly, region 2 is mainly located in Chile. Carretier et al. [86] had shown that the effective erosion rate in the Chile area is less than 0.5 mm/yr, and the erosion in the Chile area does not counterbalance the continuous rapid rock uplift, so it's not very possible. Secondly, this study investigates the glacier situation in this region by using the Randolph Glacier Inventory (RGI 6.0 version, available at http://www.glims.org/RGI/randolph60. html, accessed on 27 June 2021) [87]. The results show that there are a few glaciers along the latitude in this region, and the areas of most glaciers are less than 1 km 2 . If it is caused by glacier decline, the impacts of range and value are much smaller than current observation results. Therefore, the GIA is not a possible reason. Region 2 is near the Chile trench, where the lithosphere of the NZ subducts into the mantle beneath the South American plate. In general, subduction erosion is more likely to be the reason than surface erosion or GIA.
What is the cause of the deepening of the Moho? Jiao et al. [52] studied the Moho changes beneath the Tibet plateau and found that the Moho changes in the middle of the Tibet plateau showed a deepening trend and GRACE results showed a negative trend after deducting hydrological signals, which is similar in region 3. They pointed out three potential reasons: (1) the increase of surface loading; (2) the compression and shortening of the crust, which may lead to the crustal deformation and thickening; (3) the underthrusting of subducting plate. Region 3 is located in the Patagonian ice field. If the surface loading increases, it is more likely to be caused by the increase in the mass of glaciers, which leads to the subsidence of the Earth's crust in the short term and the deepening of the Moho surface. However, Chen et al. [25], Aniya et al. [26] and Rignot et al. [27] indicated that Patagonian ice sheet glaciers were in a period of decline and the rate of glacial melting was increasing, so the first hypothesis fail to explain the results. Richter et al. [20] used GNSS observation data to study the crustal deformation in the southern Patagonia region and their strain rate results showed that the crustal deformation showed a gradual transition from the compression of E-W direction in the west to the extension of E-W direction in the east, while the Moho changes mainly show a deepening trend in region 3. Therefore, the deformation of the crust itself cannot directly and fully explain the Moho deformation in this region, and the second hypothesis does not possess sufficient explanation. Breitsprecher et al. [88] indicated a ductile lower crust and the low viscosity of the mantle favored shortening and uplift in the south of CTJ. Chen et al. [25] pointed out that Patagonia Icefield was located in a unique tectonic setting, where the upper mantle viscosity is much lower. Thus, it may be a result of many causes-the underthrusting process of AN, the shortening of the crust and low mantle viscosity.
In addition to these three regions with significant changes, other regions also show varying degrees of Moho change characteristics. In the north of CTJ, the Moho change rate is smaller than that in the south.

Erosion and Sedimentation Effects
Bagherbandi et al. [21] discussed the influence of the distribution of sedimentation on the gravity field and even the geometry of the Moho in the SAC. From their research results, we find that marine sedimentation is mainly distributed in the eastern and northern continental margins of the SAC, and continental sedimentation and river sedimentation are also mainly located in the north. In the study area, it is shown that the vast majority of sedimentation have little effect on the thickness of the crust. Rao and Sun [53] pointed out that the spatial distribution of sedimentation thickness was considered to be spatially correlated with sedimentation rates. Considering the distribution of sedimentation thickness in the study area, the impact of sedimentation rates on the results is ignored in the study. The factors affecting the erosion rates are complex. Carretier et al. [86] studied the erosion rates of 485 areas in Chile, and the results showed that the main erosion rates were between 0.001 mm/yr and 0.5 mm/yr, and the Chilean Andes was much less than that of Taiwan, Japan and the Himalayas. At present, there are little erosion data in South America, so the influence of erosion is also not considered in the previous calculation process.

Isostatic Status of the Study Area
The Airy hypothesis states that in order to maintain isostatic balance, the mass outside the geoid is compensated internally in the form of mountain roots (anti-mountain roots) [89], so the uplift of the surface also leads to changes in the Moho interface [51,90]. If it is in an isostatic balance state, for every 1 cm rise of the surface in the study area, the corresponding Moho should deepen to about 7 cm (The mantle density is about 3.3 g/cm 3 , and the crustal density is about 2.85 g/cm 3 .). For example, the surface is uplifting at a rate of 0.97 ± 0.65 cm/yr so that in order to maintain 100% isostatic balance in region 3, the Moho would have to deepen at a rate of about −6.79 cm/yr. In reality, Moho is actually deepening −6.46 ± 1.31 cm/yr, which approaches the Moho change required to maintain isostatic balance as a whole, but the regional characteristics are quite different. This also indicates that the factors affecting Moho changes are more complex underground in southern Patagonia. This study also selects the measurement results from the northern region of Patagonia (72 • W-68 • W, 42 • S-46 • S), in which the surface uplift rate is 0.21 ± 0.19 cm/yr, and the Moho deepening rate should be about −1.47 cm/yr to maintain isostatic status in theory. However, the actual deepening rate is −1.32 ± 1.68 cm/yr, which is about 90% of the Moho change required to maintain isostatic status. Thus, the agreement cannot be reached between the theoretical Moho changes calculated by the equilibrium formula from the measured surface uplift data and the results from GRACE data in most regions, with only a few regions showing isostatic status.

Conclusions
The SAC is subject to the rapid eastward subduction of the Nazca plate in northern CTJ, whereas in the south, the Antarctic plate is subducted at a low speed with an expanding ridge. These two types of subduction result in different deformations beneath the SAC. Using satellite geodetic data, hydrological models, remote sensing data, and GIA models to constrain the movement of the study area, this study captures the signal changes caused by TECs, and estimates the Moho depth changes beneath the SAC near the CTJ.
The GRACE observations provided by CSR, GFZ and JPL are put in comparison, and the results show that the solutions of the three institutions are similarly precise, so that their average values are used in data process. This study also analyzes the influence of signal attenuation and leakage caused by truncation and filter in data processing, and the ITER method is used to recover the mass rates changes in the study area. The residual results of GRACE deducting hydrological signals cannot be fully explained by the vertical displacement changes observed by GPS, and the variation characteristics are different in the SAC. Therefore, the study concludes that signals other than GPS may be caused by TECs. After deducting the influence of surface mass movement and GIA from residual GRACE observations, the study estimates the Moho changes rates in three regions with significant changes characteristics, which are −2.12 ± 0.67 cm/yr, 0.18 ± 0.19 cm/yr and −6.46 ± 1.31 cm/yr, respectively and Moho changes also show different rates in other regions. The Moho mainly deepens in the southern SAC of the CTJ with large rates of change, while the local uplift of the Moho appears in the north and the overall deepening rates are relatively small. At the same time, in the northern part of Patagonia, the Moho depth change constrained by space geodesy accounts for 90% of the theoretical value which could be approached in the southern part. Either the SAC is moving away from isostatic balance at present or the timescale for maintaining isostatic balance is longer than the length of geodetic measurements.
The uncertainty of hydrological models and GIA model, the interpolation accuracy of GPS vertical displacement change, the accuracy of GRACE inversion of loading deformation and the accuracy of remote sensing data have an important influence on the final mass changes rates. The effects of surface erosion and sedimentation should be considered more carefully in subsequent studies.