Prospects for Imaging Terrestrial Water Storage in South America Using Daily GPS Observations

Few studies have used crustal displacements sensed by the Global Positioning System (GPS) to assess the terrestrial water storage (TWS), which causes loadings. Furthermore, no study has investigated the feasibility of using GPS to image TWS over South America (SA), which contains the world’s driest (Atacama Desert) and wettest (Amazon Basin) regions. This work presents a resolution analysis of an inversion of GPS data over SA. Firstly, synthetic experiments were used to verify the spatial resolutions of GPS-imaged TWS and examine the resolving accuracies of the inversion based on checkerboard tests and closed-loop simulations using “TWS” from the Noah-driven Global Land Data Assimilation System (GLDAS-Noah). Secondly, observed radial displacements were used to image daily TWS. The inverted results of TWS at a resolution of 300 km present negligible errors, as shown by synthetic experiments involving 397 GPS stations across SA. However, as a result of missing daily observations, the actual daily number of available stations varied from 60–353, and only 6% of the daily GPS-imaged TWS agree with GLDAS-Noah TWS, which indicates a root-mean-squared error (RMSE) of less than 100 kg/m2. Nevertheless, the inversion shows agreement that is better than 0.50 and 61.58 kg/m2 in terms of the correlation coefficient (Pearson) and RMSE, respectively, albeit at each GPS site.


Introduction
The virtual sum of the water stored as natural or artificial lakes, rivers, soil moisture, groundwater, snow and ice, and canopy water compartments is called terrestrial water storage [1].Terrestrial water storage (TWS) provides a good indicator of the total water availability in a region, and knowledge of its spatial and temporal variations is of paramount importance for all nations worldwide since water is essentially required for economic growth.Therefore, monitoring the changes in TWS is crucial because it provides the basis for proper water governance and, consequently, mitigates crises in a variety of water-dependent sectors.However, there is no such dedicated device able to quantify the variations in TWS, considering its definition, and to date, only its separate components have been individually estimated using, for example, piezometers for groundwater assessments.
Moreover, the temporal variations in Earth's gravity field are mainly due to the changes in TWS that are driven by the dynamics of the water cycle [2].Hence, monitoring the former is a viable way to assess the latter.For instance, ground-based measurements using absolute gravimetries have been used to investigate the subsurface water storage dynamics [3][4][5][6].Furthermore, relative gravimeters have been applied to monitor the variations in TWS [3,7,8].Quantum sensors have been developed using a cold-atom interferometry technique, so it might be feasible to monitor TWS with these in the future [9].Conversely, only the Gravity Recovery and Climate Experiment (GRACE) mission has achieved great success in mapping TWS changes from space at global to regional scales at a spatial resolution of approximately 400 km [10].The GRACE mission provided information for the period between April 2002 and June 2017, and its follow-on (GRACE-FO), launched on 22 May 2018, is continuing its legacy [11].However, despite some improvements (e.g., range-rate system), the spatial resolution of inverted TWS fields has still not improved (i.e., 400 km).Nevertheless, several hydrological applications need higher spatial (approximately 25 km) and temporal (daily) resolutions.Although hydrological models meet such requirements, they generally rely on ground data for inputs, which might lead to unsatisfactory results in the South American region.Thus, there is a need for further improvements in the estimations of TWS.
Another geodetic technique used to infer TWS is the Global Positioning System (GPS), as investigated by Ouellette et al. [12], who used GPS observations of radial displacements to predict the water mass loadings primarily driven by snowpack variations after accounting for soil moisture contributions.On the other hand, Argus et al. [13] proposed the inverse technique for quantifying TWS using observations of radial crustal deformation derived from GPS observations.It was shown that, for a GPS dataset used over California, TWS could be inverted with a horizontal scale of 50 km.At the same time, Borsa et al. [14] used a GPS array to investigate the droughts over the western United States by inverting the induced radial deformation to TWS.The inverted TWS depicted the overall water reduction and agreed with the observed decline in precipitation and streamflow.Borsa et al. [14] achieved a spatial resolution of about 200-300 km for the inverted TWS on the basis of the checkerboard test.Likewise, Fu et al. [15] have estimated the TWS series using the GPS-derived radial deformation, and it was revealed that the overall patterns of TWS correlated well with the physiographic properties of the region and also reflected the spatial variabilities of rainfall.Investigations based on a synthetic case study for the western United States were carried out to find the accuracy of snow water storage derived from the knowledge of its induced radial displacements at each GPS station [16].Fu et al. [15] further reported that one of the advantages of GPS-derived TWS is the possibility of increasing the spatial resolution of TWS fields over GRACE observations, and the results can be used as an independent dataset to validate GRACE-derived TWS.
Furthermore, Jin and Zhang [17] also used GPS-observed radial displacements to invert TWS fields over the southwestern United States.Zhang et al. [18] attempted to estimate TWS at a spatial resolution of approximately 100 km using 29 GPS stations over Yunnan Province (approximate area of 394 × 10 3 km 2 ), China.The authors reported an overestimation (of about 18.0%) of the TWS annual amplitudes that were inferred from GRACE monthly solutions.Nevertheless, in regions with limited GPS stations, it is unlikely that an accurate estimation of inverted TWS fields can be achieved.In such cases, GPS-observed radial displacements could be used to assess the timing and duration of TWS anomalies [19,20] and/or properly converted to a drought index [21].For example, Ferreira et al. [21] used the sparse GPS stations over Brazil with the aim of applying the observed radial displacements as indicators to assess hydrological droughts.Such a direct application of observed radial deformation series [19,21] does not require inverse modeling.However, one issue with such a use of GPS series is the time period of the observations, which, in most cases, covers only a few years, though it might be possible to overcome this problem by combining different datasets using a multi-fusion drought index approach [22].For those interested in using individual GPS stations to assess TWS, the studies carried out by Ouellette et al. [12] and Birhanu and Bendick [23] are recommended.
There are relatively few published papers on the topic of GPS-imaged TWS for the United States [13][14][15]17], and there is only one such study in China [18].Thus, this investigation is considered as a second case study outside of the United States.Geometrical observations of the Earth's shape induced by hydrologic loadings are used for imaging changes in TWS fields over South America.Notably, the region experiences one of the largest amplitudes of radial displacements in the world due to its water cycle [24].Furthermore, this study highlights the main challenges facing the inversion technique of the available GPS stations over the region.Hence, this study is the first large-scale attempt to invert TWS from GPS-observed radial displacements.Therefore, the specific objectives of this study are to (i) verify the potential spatial resolution of the imaged TWS considering the available GPS stations using a synthetic checkerboard test, (ii) examine the resolving accuracy of the GPS network and the inversion method using a closed-loop simulation considering a given geophysical model, and (iii) produce an image of TWS using GPS-observed radial displacements over South America.The results of the GPS-imaged TWS are encouraging and show that it is possible to use the method as an independent technique to assess water resources over South America.Notwithstanding the favorable results, there are limitations, as confirmed by the findings, regarding the number of daily GPS stations available for the inversion.
The subsequent sections of the study are structured as follows: Section 2 presents the study area, followed by a description of the datasets and methods in Section 3. The results are presented in Section 4, while the discussion of the results is in Section 5.The conclusions of the study are provided in Section 6.

Hydrology
The South American continent covers approximately 18 × 10 6 km 2 , which is equivalent to nearly 12% of Earth's present land area [25]; hence, it contains many diverse regions and climates.The drainage patterns in South America present an asymmetry due to the tectonic uplift of the Andes, in which the slopes of rivers draining the Pacific and Caribbean are mostly short and steep regarding those draining to the Atlantic, which after steep initial descents from the Andes, runout across sedimentary basins [26].The continent accounts for approximately 20% of global freshwater discharge [27], thus playing a crucial role in the global biogeochemical cycle.Not only is the hydrology/climatology of South America globally significant, physiographically, it also has several interesting natural features, including the Andes (the longest mountain chain), a central lowland (Orinoco, Amazon, and the Paraguay-Parana lowlands), and the Brazilian and Guiana highlands.Consequently, this physiographical diversity results in extreme gradients in water availability and quality [28].The continent also contains the largest river basin in the world-the Amazon river basin-and drier regions, such as eastern Patagonia and the extremely arid Atacama Desert.Furthermore, it has another two important river basins, La Plata and the Orinoco, which are the two largest and highly-productive watersheds on the continent.As a consequence of its abundant river water resources that are suitable for hydropower development, the continent has a large number of constructed dams, although they are less intense and widespread than those in North America, Europe, and Asia ( [29], p. 76).Nevertheless, this figure might change in the near future considering the dams under construction and planned (cf., [30]).However, impounding river waters affect surface and groundwater hydrology, for example by waterlogging lands adjacent to and upstream from reservoirs and by changing streamflow characteristics [31].
On the basis of the major geological, tectonic, and climate features, sixteen hydrogeological provinces have been distinguished in South America (cf., [32]).Consequently, South America has a wide variety of aquifer types, ranging from small aquifers to huge continental aquifers [33].The use of groundwater in the continent is concentrated mainly in the agricultural sector [28].For example, the percentage of area equipped for irrigation using groundwater is approximately 15% of the total area irrigated [34]: this is relatively small compared with that of the other continents.Moreover, as reported by Buytaert and Breuer [28], of the available water resources in South America, agriculture consumes about 70%, which is in overall agreement with the global average.Furthermore, the supply of drinking water derived from groundwater reaches 29% in the urban areas of South America, whereas for smaller towns and rural areas, it reaches more than 50% in Peru and Chile, while this figure is between 25 and 50% for most of the other countries [35].
Between 2002 and 2017, the central lowland region of South America showed a wetting trend in terrestrial water storage (Figure 1) according to monthly GRACE mass concentration solutions (mascons) (see Section 3.1.4for details).This gain in water storage is a recovery from a previous dry period observed over the Amazon and La Plata basins and is consistent with recent reports [36,37].It has been reported that the exceptional drought event in 2005 [38] was mostly caused by an increase in the sea surface temperature in the tropical North Atlantic [39].An observed negative trend in the Brazilian Highlands (Eastern Brazil) (Figure 1) was also reported by Getirana [40] and could suggest that the water resources of the Sao Francisco River basin have been compromised [41].Water depletion was also observed over the Guiana Highlands, which seems to be progressing from a wet period to a dry [42] over this portion of the Orinoco basin.The negative trend over the Patagonian ice field is also notable, confirming the impacts of global warming on the continent [43].The decrease in mass changes over central Argentina is associated with drought periods to some extent [36].The GRACE results over this region, however, could also be masked by the Maule earthquake in Chile, which occurred on 27 February 2010 [44].The superimposition of the possible influence of geodetic perturbations and other physical processes (e.g., subsidence) on the perceived role of the tropical oceans in driving inter-annual variations in land water storage over the region (e.g., [27]) warrants new research directions to improve the temporal assessment of land water storage.The next section highlights South America's geodynamics.

Geodynamics
The main characteristics of South American physiographic are driven by tectonic forcing, which provided the continent with its relief features.On the South American plate, four major domains are recognized, i.e., the continental interior, the western, the Andes, and the eastern side, in which the continental interior (South American Platform) is composed mainly of metamorphic and igneous complexes of Archean age.At the eastern edge of the South American plate, crustal deformation is rare, whereas crustal flexuring, faulting, and volcanism occur as the plate moves away from spreading centers [25].At the northern and southern coasts of the continent, active shearing or transforming margins' seismicity is high, and earthquakes are frequent.Along the present western edge of South America, the Nazca plate collides with and is subducted beneath the South American plate (Figure 1) and often accompanied by intense magmatic and seismic activity.The Nazca oceanic plate is converging eastward at approximately 60-70 mm/year relative to western South America, as shown in Figure 1, and consequently, great earthquakes occur somewhere along the South America's western edge every few years [25]; for example, the M w 8.8 Maule Earthquake (Chile, 27 February 2010) with a slip of about 10 m.GPS data show that the Maule event, maximum horizontal and vertical displacements of about 5.1 m and of 1.8 m (up lift), respectively, was evident nearby the epicenter in central Chile [45].GPS stations located at distances as far as 1000 km from the epicenter also depicted the signature of the Maule event [46].
Volcanism is also a feature of the western edge of the South American plate in the Andes, which is related to the steep subduction of oceanic crust beneath the continental plate [26], and four distinct zones are recognized.These are: the northern (with a latitudinal range of approximately 5 • N-3 • S), the central (approximately 16 • S-27 • S), the southern (approximately 33 • S-46 • S), and the austral (south of approximately 47 • S) volcanic zones.They are generally linked to subduction of the Nazca plate at a steep angle of 30 • -40 • .Some volcanic eruptions succeed some disturbances such as increased seismic activity, ground deformation, gas emission or fumarolic activity, and thermal anomalies [47,48].However, migration of fluids due to the co-seismic displacements caused by earthquakes like the Maule event have been associated with volcanoes' deformation in the southern volcanic zone [49].Conversely, Pritchard and Simons [50] have reported that over the central Andes, only a small fraction of all volcanoes are presently deforming (considered a given threshold).Earthquakes are also induced by surface reservoirs due to the interannual or short-term fluctuations in water levels, which translate to changes in the stress field or the groundwater poro pressure.Water levels in the wells due to the co-seismic and post-seismic changes have been reported, which are attributed to near-surface permeability enhancement due to the ground motion; the magnitude of the crustal strain produced by slip on the fault; and the amplification of a small crustal strain signal [53].For instance, Mohr et al. [54] examined the regional pattern of streamflow response to the Maule event through the diverse gradients of topographic and hydro-climatic regions in Chile.They have analyzed more than 85 streams that have changed their flow, in which 78 mostly increased in discharge shortly after the event due to the groundwater release.Liquefaction has also been associated with the Maule event over several locations in Chile, although the groundwater table was low (Austral summer), and only sites located near a saturation source liquefied [55].
Nonetheless, the nucleus of the South American plate, specifically over the continental cratons exposed in the Brazilian and Guiana shields, are relatively stable areas.Most of river basins with significant intra-annual variability in TWS are located in the South America Platform, and the GPS stations therein (Figure 1) are unlikely to be affected by the post-seismic relaxation movements.However, caution should be taken while using those GPS station placed at the western margin of South America (Section 3.1.1).Furthermore, glacial isostatic adjustment (GIA) effects on the Earth's crust seems to be minimal [56].

Datasets
In this section, the datasets and the pre-processing performed before using them in the experiments are introduced.Basically, the datasets comprise observed radial displacements (Section 3.1.1),a hydrological model to simulate terrestrial water storage (Section 3.1.2),remotely-sensed rainfall (Section 3.1.3),and terrestrial water storage (Section 3.1.4).

GPS Time Series
Daily solutions in Precise Point Positioning (PPP) mode provided by the Nevada Geodetic Laboratory (NGL) and processed with GIPSY/OASIS-II were used for this investigation.The final solutions refer to the geodetic reference frame IGS08 of the International GNSS Service (IGS).IGS08 is a GPS realization of the International Terrestrial Reference Frame 2008 (ITRF2008).The ocean loading was calculated using the coefficients provided by the ocean tide loading for the FES2004 tidal model and reduced from the GPS data.Of the 397 GPS stations used here, some of them cover the maximum period in the dataset of 3 years, 5 months, and 3 days, starting on 15 March 2015 and ending on 18 August 2018.However, some stations present the minimum period of observations of 1 year, 4 months, and 4 days.Only stations with solutions based on an entire 24-h period of observations were considered.Figure 1 shows the distribution of the stations and their respective time-spans (colors of the circles associated with the color bar located within the map frame).
Figure 2a shows the time series of radial displacements at each station (sorted by latitude).The displacement are driven mainly by loading signals, but also by other components, such as transient tectonic signals, GPS technical errors, and site-specific conditions.Among the loading signals, the GPS time series clearly reflects the seasonality of the hydrological mass variations.Furthermore, the period of observations is interesting since it covers the end of the previous GRACE mission and the beginning of the current GRACE mission (GRACE-FO).Hence, any potential inversion using these GPS data could help fill the gap in the TWS series over South America.
The overall variations in radial displacements for the 397 stations depict the seasonality of the Earth's crust induced by surface loadings (hydrology, atmosphere, and ocean), and peaks larger than ±30 mm are observed.Note that the highest variabilities are depicted for those stations between latitudes of 0 • (Equator) and 5.49 • S. The stations therein are mainly spread over the Amazon River Basin, with a cluster of stations in the Andes portion of Ecuador (refer to Figure 1).The stations at the northern part of the Equator are mainly concentrated over the Andean region of Colombia, and they experienced uplifts from March 2015-February 2017 that might be mainly associated with an extreme dry period over that region.Figure 2b shows the number of stations that are available at each epoch for estimating TWS using inverse modeling.As mentioned above, there are 397 stations in total (Figure 1), although some are unavailable for use in this study because of gaps in continuous observations.For example, there were no solutions (Figure 2b) for the following days: 12, 17, 20, 22, 25, 28, and 30 December 2015; 2 January 2016; and 23 March 2017.Therefore, it is impossible to produce an image of TWS for these days.Furthermore, for 2018, 8 and 11 January, 25 and 31 March, and 10 April have only 85, 86, 72, 66, and 60 available solutions, respectively (Figure 2b).Despite this challenge, the data were still considered in the inversion to investigate the overall performance with a reduced number of stations.For the other days, there were 314 stations on average, with a minimum and a maximum of 111 and 353 stations, respectively.
Nevertheless, for the subsequent analysis, the results shown in Figure 2a were reduced by contributions from the non-tidal deformation loading from the atmosphere and ocean to each GPS station, resulting in loading displacements due to hydrological mass variations.Thus, the surface loading products used are those provided by the Global Geophysical Fluids Center (GGFC) of the German Research Center for Geosciences (GFZ), as further described in [57].These products are provided in the center of figure (CF) frame, considering the elastic Earth model "ak135" [58].The loading fields were given on a 0.5-arc-degree regular global grid for time intervals of 3 h for non-tidal loadings (caused by both atmospheric and oceanic mass variations).The global gridded fields were interpolated at the station location by using bi-cubic interpolation.Next, the extracted series were filtered using a sliding window with a length of eight (centered on the current and previous elements) to reduce the three-hourly atmospheric and oceanic effects from each time series.Finally, linear interpolation was conducted for each GPS site since their epochs differ from that for atmospheric and oceanic loads.
As discussed in Section 2.2, South America is relatively stable, whereas its western part is characterized by seismic activities.The GPS time series were de-trended to reduce the long-term tectonic effect.Consequently, traces of those contributions were not reflected in the GPS-imaged TWS.As a result, the inversion reflects only the seasonal time scale variations.To do this, it is assumed that the temporal variations in the radial displacements as observed at each GPS station are composed by a trend, seasonal loading effects, jumps (e.g., antenna changes), and co/post-seismic.Thus, the "trajectory model" proposed by Bevis and Brown [59]: was applied to remove the contributions due to the jumps, co/post-seismic, and linear trends to the GPS series.In Equation ( 1), h(t 0 ) is the ellipsoidal height in the reference epoch t 0 , ḣ is the velocity, t is the time in years, H is the Heaviside function, h i H characterizes the jump that occurs at time t i as an instantaneous displacement, S k and C k are the coefficients for computing the amplitude for each of the two frequencies (here, annual, k = 1, and semi-annual, k = 2, waves), a i is the magnitude of the exponential decay due of the post-seismic change in the period ∆t j in the years since the j earthquake event occurred, and the decay time T was configured with the default value T = 1 year, reaching an excellent fit to the data.Note that this model is insensitive to the effect of assigning a moderately-erroneous value to the T parameter [59].
The parameters of the trajectory model can be estimated by least squares adjustment.It should be noted that the full Equation ( 1) was applied only to the stations most affected by the Maule earthquake (27 February 2010, Chile).For the other stations, only the first four terms of Equation (1) were considered.

Model for Continental Hydrology
The Noah-driven Global Land Data Assimilation System (GLDAS-Noah) is generating a series of land surface state (e.g., soil moisture and surface temperature) and flux (e.g., evaporation and sensible heat flux) products [60], and it was used here to simulate TWS over the continent at a daily scale.The GLDAS-Noah Version 2.1 at a spatial resolution of 0.25 arc-degree and a temporal resolution of 3 h [61] is available from NASA's Goddard Earth Sciences Data and Information Services Center (GES-DISC).GLDAS-Noah TWS is virtually the sum of the water stored as soil moisture in all layers, accumulated snow, and plant canopy surface water.Soil moisture is a significant component of the loading [62], and missing elements of groundwater and surface water (e.g., lakes) might impact the results depending on the location of the computation point.
The 3-hourly GLDAS-Noah TWS series were filtered with a sliding window with a length of eight to describe the series on a daily basis as those of GPS.Furthermore, monthly solutions were de-trended to better agree with GPS for comparing them.Nonetheless, the TWS fields from GLDAS-Noah were considered only in a closed-loop simulation to assess the quality of the inversion algorithm considering the spatial distribution and number of GPS data as presented in Figures 1 and 2b, respectively.Furthermore, it was used in order to check if the GPS-imaged TWS was within an acceptable range.

TRMM Rainfall Fields
The Tropical Rainfall Measuring Mission (TRMM), a joint U.S.-Japan satellite mission, was dedicated to monitoring precipitation and to estimating its associated latent heating over tropical and subtropical regions of the globe [63].The 3-hourly 3B42 TRMM Multi-satellite Precipitation Analysis (TMPA) Version 7 product was used to generate daily fields by simply accumulating the eight 3-hourly values so that the resulting units were mm/day.The rainfall products were used to assess the overall quality of the GPS-observed radial displacements and GLDAS-Noah TWS fields over South America.

GRACE TWS Solutions
Apparently, among the GRACE Level 3 products, that is, the monthly grids with water masses, the mascon solution provided by the Center for Space Research (CSR) seems to have a higher reliability because of its regularized solution without external models [51].The GRACE mascon solutions simplify the use of GRACE-TWS observations for several applications, such as oceanography, land surface hydrology, cryosphere, etc.For example, there is no need for empirical filters and/or scale factors.The temporal resolution is monthly, and the spatial resolution is 0.5 arc-degree; however, this resolution is actually constrained to the band-limited nature of GRACE data, which is likely to be around 300-400 km.More details regarding the derivation of the mascon solutions and their performance metrics were documented by Scanlon et al. [64].Figure 1 shows the overall linear trend of TWS over South America based on the time span between April 2004 and June 2017, although some regions with water depletion and increase were identified (Section 2.1).

Methods
This section introduces the methodology used to invert TWS from observed radial displacements.Since it is necessary to validate the approach, a simulation using synthetic data is described.

Forward Modeling: From Mass Loadings to Radial Displacements
The radial deformation ∆r can be written in terms of the Love number as inferred from Lambeck ([65], p. 99): where V n is the gravitational potential V up to a degree n, h n is the Love number, and g is the gravity.The gravitational potential is caused by the mass of the water, expressed in terms of equivalent water height ∆h.Considering a spherical Earth of radius R ∆h, that is, ∆h/R ≈ 4.7 × 10 −8 , the gravitational potential is given by ([66], p. 3): with ρ w denoting the density of water (approximately 1000 kg/m 3 ) and G representing the gravitational constant.is the distance between the computation point, given by the latitude ϕ and longitude λ, and the running point, given by the latitude ϕ and longitude λ ; dΩ is the surface element of the unit sphere given by: dΩ = cos ϕ dϕ dλ . (4) The reciprocal distance 1/ in Equation ( 3) is expanded in a series as ( [65], p. 99): where the P n (cos ψ) terms represent the Legendre polynomials of the n th degree in cos(ψ).
The hydrological loading now can be expressed by considering Equations ( 3)-( 5), which gives: where G is Green's function for loading.For radial deformation, the function is: h n P n (cos ψ), (7) in which g ≈ GM/R , M is the mass of the Earth and ∆σ = ρ w • ∆h is the surface mass in kg/m 2 .
For each compartment i, the surface masses (∆σ) are replaced by their average value ∆σ.Hence: since the hydrological masses are constant over each compartment, for example 1 • ×1 • .The integrand, G r (ψ), which is constant over the compartment i, may be replaced by its value at the center of the compartment, which gives G r (ψ i ).The final integral is simply the area of the compartment.Finally, Equation ( 8) becomes: Green's function given by ( 7) is pre-computed and available for ψ ranging from 0.0001 • -180 • [67] for different Earth models.Here, we considered Green's functions calculated in the CF frame provided by Wang et al. [68] for the elastic Earth model "ak135".
Equation ( 9) relates the hydrologic surface water loadings ∆σ at a location with ∆r, the crustal radial displacements induced by the loadings.However, what exactly does hydrologic water loading mean here?The terrestrial water storage includes water stored on Earth's surface in the form of snow, ice, and vegetation, as well as water stored within the ground (groundwater and soil moisture), which exert a significant load on Earth's surface, reaching up to 700 kg/m 2 in regions such as the Amazon Basin [21].In this hydrologic context, the Earth's surface changes due to elastic and poroelastic responses.The elastic responses as given by Equation ( 9) reflect the gravitational attraction of all water storages on and below the Earth's surface as indicated by Equations ( 2) and (3), that is mainly due to the variations in surface and sub-surface hydrology.This is somewhat like the GRACE principle, which considers the variations in the Earth's gravity field to estimate the masses driving those variations.However, GPS radial displacements (the observations) are also sensitive to ground subsidence due to poroelastic effects (Figure 3a-c), as well as the irreversible subsidence due to the aquitards' compaction (Figure 3d).The loads (the sum of snow, ice, vegetation, and soil moisture) contribute to the seasonal displacements observed in GPS series, as shown in Figure 2a.Considering the long-term changes in TWS depicted in Figure 1, it is possible to see regions with gain (positive trend) and loss (negative trend) of loadings, which consequently causes a net uplift/subsidence.Thus, monitoring the radial displacements (∆r) using GPS makes it possible to estimate the terrestrial water storages (∆σ).The knowledge of ∆σ provides a key measure of the continental water cycle and available water resources in a given region or river basin [13].
3.2.2.Inverse Modeling: From Radial Displacements to Mass Loadings Notably, Equation ( 9) takes the general form of a forward problem (or direct problem).For instance, given a cause m (mass loadings) and a model G (matrix describing the system response, i.e., Green's function), find the effect d (radial displacements), where d = Gm.Consequently, the goal of the inversion is to find the unknowns m (mass loadings) given the known data d (observed displacements) and the system response matrix G.If perfect and continuous measurements d are available, it is desirable that the solution m is such that Gm = d holds [70].The generalized inverse G † provides such a solution.Nevertheless, in the presence of measurement errors, the generalized inverse cannot be used as such, so regularization is necessary.
Alternatively, the linear system can be solved using the damped least squares problem as ( [71], p. 90): where α is a regularization parameter.The minimization is performed in the usual way by solving the augmented normal equation, which reduces to: where I is the identity matrix characterizing the zeroth-order Tikhonov regularization in (11).
The matrix G is given as: that is a matrix of partial derivatives with respect to each unknown (∆σ i , where m is the number of observations and n the number of parameters.Depending on the discretization (n, number of grid cells, i.e., parameters) of Equation ( 9) and the number of observations (m, e.g., the number of GPS stations shown in Figure 2b), G becomes poorly conditioned and generally rank(G) < min(m, n).Spectral decomposition of G shows that singular values decay gradually toward zero and do not show an obvious jump between nonzero and zero singular values (this is the case here).These problems are generally referred to a "discrete ill-posed problems".The regularization parameter α can be selected using, among other techniques, generalized cross-validation (GCV) [72,73].The premise of GCV is that if an arbitrary element d i of d is removed, the corresponding regularized solution should predict this observation well [74].Hansen [74] illustrated that the GCV method seeks to balance the perturbation and regularization errors; thus, it is related to the corner of the L-curve (another strategy for selecting the regularization parameter (see, e.g., [71], p. 91)).
The model resolution matrix, which maps the true solution to the inversion result [75], for the Tikhonov regularization method is given by ( [71], p. 95): In Equation ( 13), V is an n-by-n orthogonal matrix with columns that are basis vectors spanning the model space.It can be obtained by the singular value decomposition of m G n that provides m U m and n V n orthogonal, and m S n is a diagonal with elements s ii , which are singular values (the diagonal of S is generally arranged in decreasing size), while some elements may be zero.F is an n-by-n diagonal matrix with elements given by filter factors ( f i ) ( [71], p. 92) The trace of R m,α , Tr(R m,α ), is often used as a simple quantitative measure of the resolution.If Tr(R m,α ) is close to n, then R m,α is relatively close to the identity matrix I [71].Note that R m,α is dependent on the particular value of α used in F.

Checkerboard Test: Model Resolution
In the field of inverse problems, especially in tomographic seismological studies, the checkerboard (also spike) test is commonly used to illustrate the quality of inversion results and assess the resolution of inversion [76].Resolution tests with the checkerboard model as in Figure 4 are very commonly used in practice [71].For example, Borsa et al. [14] used the approach to infer the resolution of the inversion of TWS using a dense GPS network.Likewise, Fu et al. [15] conducted a synthetic test using the checkerboard in order to investigate the sensitivity and robustness of their inversion method.In both studies (i.e., [14] and [15]), the authors confirmed the suitability of the checkerboard test to identify the regions in which the inversion could recover the input model well.Such regions are generally related to the proximity of the GPS stations.A checkerboard test in which one applies a forward model to the synthetic loadings with patterns of alternating high (+250 kg/m 2 ) and low (−250 kg/m 2 ) attenuation elementary squares within the South America domain is shown in Figure 4.In order to account for the edge effects while using Equation ( 9), a buffer zone of 5 arc-degree was adopted.A large portion of the buffer zone lies over the oceans, and zeros were assigned to these "squares".The test implemented three different resolutions: 1, 2, and 3 arc-degree.The forward model of the displacements for each GPS site were then inverted for the available stations considering their resolutions (more details are provided in Section 4.1).This is a useful way to visualize the limitations of the observations (d) and the functional model given by Equation (9).

Results
This section presents the main results, which are divided according to the three main experiments shown in the flowchart in Figure 5.The results of the first experiment (Section 4.1) are based on a checkerboard (Section 3.2.3)for three different spatial resolutions, that is 1 × 1 arc-degree, 2 × 2 arc-degree, and 3 × 3 arc-degree.This experiment involved forward (Section 3.2.1)and inverse (Section 3.2.2) models.Next, the second experiment (Section 4.2) used the root-mean-square (RMS) of GLDAS-TWS in two different resolutions, 1 × 1 arc-degree and 3 × 3 arc-degree, and then, forward and inverse models were carried out.Since the number of stations was varied, the impact of the varied number of stations on the inversion, considering both the RMS of TWS and TWS daily fields, was investigated.Finally, the third experiment (Section 4.3) considered the observed GPS radial displacements in an attempt to produce images of TWS; this involved only inverse modeling.Cross-correlation analyses were also performed in the third experiment, and a summary is presented in Appendix A.  The checkerboard test was applied to identify the potential resolution of the inverted TWS products from radial displacements "observed" at each GPS site.Furthermore, it was important for the test to verify possible problems with the distribution of the stations since the network geometry impacts the coefficient matrix G in the linear system in Equation (11).Firstly, synthetic checkerboards were created with different resolutions covering a region beyond 5 arc-degree from the boundary domain of South America.Figure 6a,d,g shows the arrangement of the "squares" for spatial resolutions of 1.0, 2.0, and 3.0 arc-degree, respectively.Next, these loading patterns were used to compute their respective contributions to the displacements at each of the 397 GPS stations using Equation (9).

RMS
Figure 6b shows the results of the inversion of the forwarded displacements at each GPS site, and Figure 6c shows the residuals based on the results of the original synthetic loadings (Figure 6a) and the inverted ones (Figure 6b) using a resolution of 1.0 arc-degree (approximately 111.2 km at the Equator).With this resolution (1.0 • -by-1.0• ), it was not possible to recover the spatial patterns of the input checkerboard given the number of observations and distributions of the 397 GPS stations (Figure 1).The configuration of this inverse problem using a 1.0-arc-degree resolution generated 2862 parameters to be estimated over South America (using a buffer zone of 5.0 arc-degree).The normal matrix in the linear system in this configuration has a condition number with a magnitude of the order of 25, and the rank of the normal matrix is 396, which characterizes an ill-posed problem.
From Figure 6b, it is possible to see deficiencies in the data over the Amazon rainforest, Eastern Patagonia (Argentinian side), a large portion of Peru (Andes Mountain), and Eastern Brazil.A low density of GPS stations is observed over the Amazon rainforest region (see Figure 1), and the inversion results based on the current density of GPS stations failed to deliver meaningful results at a resolution of 1.0 arc-degree.Furthermore, there were no GPS stations available for this study over Peru, and only two stations were available in Eastern Patagonia.At this resolution, the inversion of GPS data seems to yield overall "squares" patterns at the station locations (Figure 6b).
Moreover, Figure 6c depicts the magnitude of the residuals of the inversion.As mentioned above, it was found that near the stations, the inversion was able to reconstruct the overall spatial patterns of the "squares".However, the amplitudes of the residuals were large for almost the whole of South America.This implies that imaging TWS using the actual distribution and quantity of GPS sites might be compromised at a 1.0-arc-degree resolution.Results based on 2.0 arc-degree also failed to produce excellent results from this experiment, as shown in Figure 6d-f, for the condition number with a magnitude of order 22, and the rank of the normal matrix was 386.Thus, the resolution of 2.0 arc-degree will not be further addressed here.Instead, results obtained with a resolution of 3.0 arc-degree were adopted and are presented hereinafter.
Figure 6g depicts the distribution of the "squares" with alternating low and high values for TWS at a 3.0-arc-degree resolution.This configuration presented 317 unknowns with the condition number of the coefficient matrix reaching a magnitude of order 13 and rank equal to 313, characterizing an ill-posed problem.The displacements were computed based on the basis of the loading patterns in Figure 6g, and the inverted results using these displacements are shown in Figure 6h.Overall, it seems to be possible to recover the spatial patterns of TWS using a spatial resolution of 3.0 arc-degree with the given distributions of the 397 stations.The patterns of the "squares" containing alternating TWS are well recovered within the region domain.The residuals (computed as the difference between the input model and the inverted one) are relatively low with magnitudes close to zero (Figure 6i).
Nevertheless, the real scenario is slightly different from this synthetic experiment.Thus, like the results for the checkerboard test, GLDAS-Noah was used to simulate the spatial distribution of TWS at spatial resolutions of 1.0 • -by-1.0• and 3.0 • -by-3.0• .The root-mean-square (RMS) of the TWS series was computed at each grid point to characterize the strength of the water mass loadings over South America.The radial displacements were estimated at each GPS site using Equation ( 9) from the inputs from the RMS of GLDAS-Noah with two resolutions.These radial displacements were inverted in order to recover the input mass loading (RMS of the TWS fields) considering the two resolutions (Figure 7).   Figure 7a shows the inputs used in the forward model with a resolution of 1.0 arc-degree, and the inversion results are shown in Figure 7b.Overall, the RMS of the GLDAS-Noah-derived TWS showed amplitudes of about 180 kg/m 2 over the region of the Amazon River Basin.However, the inverted values at the resolution of 1.0 arc-degree could not depict an accurate image of the water mass patterns, as shown in Figure 7b.The residuals between the input model and the inverted ones showed a correlation with the GPS site distribution (black dots in Figure 7c).The RMSE for those residuals was about 38.34 kg/m 2 with a mean value of 12.86 kg/m 2 (the buffer zone due to the edge effect was removed in this depiction).
The same experiment was carried out, however, considering the RMS of TWS fields given in a spatial resolution of 3.0 arc-degree (Figure 7d).As depicted by the checkerboard test using a 3.0-arc-degree resolution (Figure 6h), the results of the inversion recovered the spatial patterns and amplitudes of the input model, as shown in Figure 7e.Indeed, the results for the inverted mass using GLDAS-Noah inputs, which provide a natural pattern of TWS distribution, are slightly better than those for the synthetic model, as can be seen over the buffer zone by comparing Figures 6h and 7e.However, the buffer zone was not considered since the cells therein were used only because of the edge effects in the integration of Green's function.As can be seen in Figure 7f, the amplitude of the residuals of the differences between the inputs (Figure 7d) and inverted (Figure 7e) model was at 1.10 kg/m 2 , which is almost negligible relative to all of South America.However, the northeast region presented relatively high uncertainties, and the density of the GPS stations there was low.

Experiment 2: How About the Impact of the Varied Number of Daily GPS Data?
From Figure 2b, it is known that the actual number of stations with daily GPS solutions over South America varied, and a few days lacked solutions.Further, despite the total number of stations being 397, only about 300-350 were indeed available for each daily imaging throughout the period study.However, in some cases, the number of days was about 50-150, such as in the second half of 2017 and the first half of 2018, as shown in Figure 2b.Thus, in order to know what the quality of the inversion in this scenario could be, forwarded radial displacement fields using the RMS of the daily TWS derived from GLDAS-Noah were used to assess the impact of the varied number of stations.At this stage, only the 3.0 arc-degree resolution for the inversion was used, and the number of columns in the coefficient matrix was n = 317; that is, the number of parameters was considered, while the number of rows m changed according to the solutions available for each day (Figure 2b).
First, the grid with the RMS of TWS (as given in Figure 7d) was used to forward the daily radial displacements from each available GPS station with a solution.Then, these daily solutions with a varied number of stations covering the period ranging from 15 March 2015 to 18 August 2018 were inverted in an attempt to reconstruct the original input model (the one displayed in Figure 7d).Figure 8a shows the results of RMSE versus the number of stations.Overall, there was a correlation between the number of stations with daily solutions and the RMSEs (Figure 8a).However, for a daily number of solutions equal or superior to 317, the RMSEs were on the order of 1.0 kg/m 2 ), without a correlation between them.Below 317 to the minimal number of daily solution (equal to 60), there was a linear correlation between both (Figure 8a).For a given number of stations equal to 60, the RMSE reached approximately 57.1 kg/m 2 , and generally, increasing the number of stations reduced the RMSE.However, even when using 317 stations, the RMSEs presented some discrepancies, with the maximum value reaching 4.93 kg/m 2 .The variations in the RMSEs for a fixed number of daily solutions were probably due to the arrangements of the observations.This affected the model space in the coefficient matrix, hence providing different values for the inversion despite the same number of observations.Moreover, the rank p of the coefficient matrix m G n was also affected and, consequently, so were the inverted values.
Secondly, an experiment with the daily TWS fields from GLDAS-Noah compatible with each GPS daily solution was carried out.This could depict the real results more realistically.The daily inputs from GLDAS were used to forward the daily radial displacements from each available GPS station with a daily solution.Next, they were inverted to the corresponding respective daily TWS grids.The RMSEs were computed for the differences between daily observed TWS (GLDAS-Noah) and those inverted (GPS-imaged TWS), and the results are summarized in Figure 8b.Overall, the results depicted in Figure 8b show almost the same patterns as those of Figure 8a.
The results depicted in Figure 8a,b suggest that it is possible to produce images of daily TWS using the daily variations in radial displacements.The RMSEs correlate with the number of available daily solutions as observations to the inverse problem, and the optimal number was close to 317, which is the number of unknown parameters to be estimated.However, a significant conclusion from Figure 8a,b is that the variability of TWS and, consequently, the radial displacements, through the wet and dry seasons in the different regions of South America, seemed not to impact the inversion results.That is, one could produce images of TWS using GPS with almost the same accuracy in different seasons.
Nevertheless, for a minimum number of stations with available daily solutions, the RMSE was approximately 100.89 kg/m 2 .For the available daily solutions at 317 stations, the maximum RMSE was 0.95 kg/m 2 .However, about 300 stations could also provide reasonable RMSEs depending on their distribution; the values reached up to 8.89 kg/m 2 .Furthermore, there was a bias between the model input and the imaged model; for example, for the 60 stations with a daily solution, the bias reached 15.46 kg/m 2 .

Experiment 3: What Is the Result of Considering the Actual GPS Data?
From Figure 8b, it is clearly seen that, for some stations with 300 or more daily solutions, producing an image of TWS over South America looks feasible, with an overall RMSE value that was better than 11.53 kg/m 2 .However, actual radial displacements as those depicted in Figure 2b had high uncertainties, and other signals were difficult to model (see, for example, [18]).There are two ways in which variations in TWS may affect observed radial displacements (Section 3.2.1).First, water mass variation may induce crustal deformation, as it is the basis of this study.Second, groundwater pumping may cause soil compaction and thus subsiding of the observation point [7], which generally occurs within a region of 10 km from the cone of the depression of the water table.The GLDAS-Noah TWS series were extracted from each GPS site location, and both radial displacements and TWS correlated (anti-correlation) well (Figure A1).This correlation supports the idea that radial displacements might be mainly driven by the seasonal mass loadings and unloadings at the GPS sites, as shown in Figure 1, without evidence of a poroelastic component.
Therefore, in comparison with the radial displacements predicted by a hydrological model like GLDAS-Noah, a closed-loop simulation introduces certain errors.Thus, the actual observations of the radial displacements after appropriately accounting for the atmospheric and ocean loadings, as already discussed in Section 3.1.1,were used to produce images of TWS.Exactly 1244 daily GPS-imaged TWS solutions were inverted from the available GPS observed radial displacements.Figure 9a summarizes the inversion results of those predicted by GLDAS-Noah TWS in terms of RMSEs.It is noteworthy that GLDAS-Noah was used here only to check whether the GPS-imaged results were within an acceptable range since there were missing water compartments in this hydrological model, such as groundwater and surface waters (e.g., rivers and reservoirs).
The minimal RMSE value (40.07 kg/m 2 ) was found using the GPS solutions available on 12 February 2017 and was the result of 314 GPS stations (Figure 9a). Figure 9d shows the distribution of these 314 stations.Figure 9b shows the GPS-imaged TWS, and Figure 9c shows GLDAS-Noah TWS.
Overall, the GPS-imaged TWS (Figure 9b) depicted the most important patterns of TWS (Figure 9c) variations as the amplitudes were lower than those of GLDAS-Noah (Figure 9d).There were significant discrepancies between the two, particularly over the Amazon and Tocantins basins.The mean of the differences was approximately 16.79 kg/m 2 (i.e., GPS-imaged TWS underestimated the amplitudes determined with GLDAS-Noah TWS) with a standard deviation of 36.50 kg/m 2 .
However, only 74 (approximately 6%) GPS-imaged TWS fields presented RMSE values below 100 kg/m 2 , as shown in Figure 9a.Nevertheless, from the results depicted in Figures 6c and 7c, the assessment of the inversions showed acceptable residuals near the sites used for forward modeling (Section 4.1).Thus, it is imperative to check whether the GPS-imaged TWS series is comparable to that simulated by GLDAS-Noah at a site.For each GPS location, linear interpolation was used to extract the series for both the inversion results and GLDAS-Noah at each daily solution.However, this could be carried out for only the stations lying within at least two grid points in each direction (north-south and east-west).Furthermore, residuals larger than 100 kg/m 2 and less than −100 kg/m 2 were masked out.Only 207 out of 397 sites fulfilled these criteria.Figure 10a depicts the results of the differences between the two series at each location of the 207 GPS sites.
Generally, the residuals of the differences presented an RMSE for each site in the range of 28.17-61.58kg/m 2 , as summarized in Figure 10b.About 8.7% of the presented RMSE values were less than 50 kg/m 2 .
Moreover, Figure 10b also depicts the correlation coefficient (Spearman) between the GPS-imaged TWS and GLDAS-Noah TWS.Approximately 82.6% of the stations presented a correlation coefficient equal to or larger than 0.50, which indicates relatively good agreement between the two series.A negative correlation coefficient was found for a few stations (approximately 3.4%), one of which is depicted in Figure 10b with the lowest RMSE (28.17 kg/m 2 ).This shows that GPS-imaged TWS captured the main characteristics of the hydrological cycle, as simulated by GLDAS-Noah.

Discussion
The assessment of TWS has been possible only with the observations acquired by satellite gravimetries (e.g., GRACE mission [10,11]) and, to some extent, by terrestrial gravimetries [3][4][5][6][7][8].Recently, the GPS-observed radial displacements driven by the hydrological cycle have been applied in the United States [13][14][15]17] and, at a small scale, in China [18], to produce images of TWS fields.Therefore, for GPS observations, this study is the second example outside the United States that aims to produce images of TWS using observations of radial displacements driven by TWS over South America.Implementing this technique in South America makes it unique since some GPS stations in the region experienced one of the largest amplitudes in radial displacements in the world [24].Hence, the region could benefit from the inversion of TWS from GPS observations since it can provide high temporal and spatial resolutions for those values derived from GRACE measurements [13].Furthermore, GPS-imaged TWS could be used to fill in the temporal gap between the GRACE and GRACE-FO missions, as well as validate their estimations.However, this might be possible only for regions with a dense GPS network [16], although reasonable results have been reported for a relatively small number of stations [18].Here, 397 GPS stations were considered from which only (at most) 353 were available, with an average of 314 stations having daily solutions.That is, there was an variations in the number of stations due to issues such as the lack of a complete 24-h data record, repair, and antenna changes (cf., Section 3.1.1).
Considering these factors, should the 397 stations be available for solving the inverse problem-that is, estimating TWS from the inversion of the radial displacements-then experiments based on checkerboard tests (Section 4.1) show that a resolution of 3.0 arc-degree is feasible (see the summary in Table 1).At this resolution, the spatial patterns and amplitudes of the TWS could be almost exactly recovered on the basis of the closed-loop simulation (Figure 6e).However, when using resolutions higher than 3.0 arc-degree, for example 1.0 arc-degree, the inversion failed to recover the input model, as shown in Figure 6b.Of course, changing the resolution implies increasing/decreasing the number of unknown parameters in the inverse problem (11).In the present study, the number of parameters was n = 2862 and n = 317 for 1.0-and 3.0-arc-degree resolutions, respectively.The respective ranks for the normal matrices were p = 396 and p = 313.For the case of 1.0-arc-degree resolution, despite p < n (rank-deficient problem), the trace of the model resolution matrix (R m,α ) was approximately 396; that is, R m,α was not close to an identity matrix using this arrangement.Furthermore, the off-diagonal terms of R m,α cannot be negligible, and they contributed to biases in the inverted model, which could be avoided by changing G to improve the solution [76].

RMS of TWS
Daily TWS * * Considering daily solutions in Figure 9a with RMSE less than 100 kg/m 2 .
Nevertheless, the resolution of 3.0 arc-degree was considered since this choice affected R m,α , and its trace was exactly 317, given an identity matrix with almost zero for the off-diagonal elements.In such a case, the estimated model can be perfectly resolved in synthetic experiments [71].Indeed, optimal results are possible using the spatial resolution of 3.0 arc-degree to recover the spatial patterns and amplitudes of the input signal.This is depicted in Figure 6h,i, respectively, where the number of parameters estimated is n = 317 and the rank of the normal matrix is equal to 313, characterizing an ill-posed problem (there is no obvious jump between nonzero and zero singular values).The same investigation was performed using the RMS of the GLDAS-Noah TWS series as the input model with the resolutions of 1.0 and 3.0 arc-degree (Table 1).The actual water patterns described by the RMS provide almost the same conclusion as that obtained from the alternating patterns in the checkerboard test.That is, Figures 6h and 7e present the resolving power of the GPS dataset available over South America.
However, the real scenario regarding the availability of daily GPS solutions is entirely different from the number used in the experiments in Section 4.1, i.e., 397 observations.Thus, the GLDAS-Noah TWS was considered in two different ways (Section 4.2).Firstly, only the RMS (as depicted in Figure 7d) was used in the inversion for the actual number of stations available on each day, as shown by Figure 2b.Overall, the RMSEs present a correlation with the number of GPS stations available per day (Figure 8a).This is observed only when the case that the number of stations was less than approximately 317.As expected, the largest RMSE (57.6 kg/m 2 ) was calculated for the minimum number of observations, i.e., 60 observations on 10 April 2018.On the other hand, RMSEs better than 2.1 kg/m 2 were found for several stations (equal to or larger than 317).Actually, for the same observations (e.g., 317), different RMSEs were obtained, which implies that the arrangement of the stations, that is the network geometry, could deliver different results that are no worse than 2.1 kg/m 2 .The results depicted in Figure 8b present the RMSEs considering each daily solution provided by the GLDAS-Noah TWS series.Further, in the presence of noise-free GPS data, the same conclusion might be drawn for real GPS-imaged TWS, and the intra-annual variability in inverted TWS can be feasibly recovered.
Nonetheless, the results for the real GPS data (Section 4.3) are not optimistic like the simulations described in Sections 4.1 and 4.2.It was found that only approximately 6% of the daily GPS-imaged TWS presented RMSE values below 100 kg/m 2 (Figure 9a).The inversion results using 314 GPS sites with daily solutions on 12 February 2017 were in good agreement with GLDAS-Noah TWS, although differences between the two are expected due to the missing water compartments in the hydrological model [60].The RMSE for this epoch was 40.07 kg/m 2 , which agrees with previous results for Yunnan Province, China, published by Zhang et al. [18], of about 48-52 kg/m 2 at annual time scales based on monthly values, which are smoother than daily solutions used in the present study.Likewise, RMSEs in the range of 25-77 kg/m 2 across the mountains of the western United States were reported by Enzminger et al. [16], after using synthetic data to assess the inversion of radial displacements into TWS.Moreover, on the basis of the results of the checkerboard test and the closed-loop simulation, it was found that generally, the results were relatively good near the GPS sites (Figures 6c and 7c).Thus, the extracted values of GPS-imaged TWS series for each GPS location and those for GLDAS-Noah presented a good correlation coefficient (≥0.5) for about 82.6% of the sites with RMSE values lower than 61.58 kg/m 2 (Figure 10).This is somewhat equivalent to the approach used by Ouellette et al. [12], where six GPS stations were used to estimate snow water equivalent at each site location.
In this study, only a limited number of GPS stations were considered over South America (see Figure 2b) to produce images of daily TWS fields.A more significant number of sites could lead to a better generalization of the results found here, and it could be directly compared with the findings of Borsa et al. [14] and Argus et al. [13].The contributions of GPS uncertainties to the inversion was not investigated, although this knowledge could perhaps improve the results since they could work as weights.All GPS raw data are being re-processed in order to improve the estimations of the coordinate time series, as well as their uncertainties.More GPS data are being collected, especially over Brazil, in order to expand the distribution of the sites.
Moreover, the investigation of the impact of groundwater variations on the results was not considered since there are no stations placed close to excessive pumping sites (near-field).Regions of South America have experienced several drought episodes [21,36,40,41] and, consequently, reduced groundwater levels.Thus, the drought's contribution to the heads of aquifers can result in land subsidence, and although a good correlation between modeled land water storage (GLDAS-Noah) and GPS-observed radial displacements (Figure A1) was observed, this deserves further investigations.For example, de Luna et al. [77] has reported vertical displacements with a rate of approximately 0.68 mm/year based on the spirit leveling in Recife, northeast of Brazil, as a consequence of the excessive exploitation of groundwater.Recife lies over a multi-layered aquifer system and has experienced anthropogenic pressure on groundwater resources due to the intense water pumping to alleviate drought impacts, and there is a large number of illegal wells [78].
Although the need for further studies within this field is paramount, the results obtained were optimal and confirmed previous findings that GPS can be used as an independent tool to assess water resources.The results can be extended to South America, even though they are based on site locations.From a hydrological point of view, radial displacements sensed by GPS offer a new way to monitor the hydrological process at horizontal scales of approximately 300 km.

Conclusions
The advancement in space-borne sensors has benefited South America, a data-poor region, by providing relevant information about its water resources.One example of such a contribution is the results from the past GRACE mission, which provided monthly fields of TWS at a spatial scale of approximately 400 km.Likewise, GPS also has provided the opportunity to map TWS by exploring the elastic response of the Earth's crust induced by the dynamics of the water cycle.Thus, for the first time, an image of the TWS over the continent was produced using the radial displacements sensed by GPS.Firstly, experiments based on simulations showed that it is feasible to produce images of TWS at a spatial scale of 300 km and with a temporal resolution of one day, which presents advantages in comparison with the results of GRACE.Secondly, while considering the actual GPS-measured radial displacements, which contain errors, the findings indicate that it is possible to produce an image of TWS only for a few days with an RMSE better than 100 kg/m 2 in comparison with a hydrological model (GLDAS-Noah).Considering the values of the daily GPS-imaged TWS, the results showed that 82.6% of the stations presented correlation higher than 0.50 and RMSE values lower than 61.58 kg/m 2 .Therefore, GPS shows potentiality for monitoring TWS in different physiographic regions of South America, which contains one of the driest (Atacama Desert) and wettest (Amazon Basin) regions of the world.Hence, GPS can currently extend the capabilities of remotely-sensed products such as TWS derived from GRACE in the sense that they could be combined to take advantage of both temporal resolution and spatial coverage, to enhance the estimations of TWS over South America.Future work

Figure 1 .
Figure 1.Distribution of GPS stations (circles) over South America.The colors of the circles depict the number of years of observations, as shown in the color bar (within the map's frame) for each GPS station.The colored surface shows the long-term changes in terrestrial water storage (TWS) for the period between April 2002 and June 2017 using GRACE data[51].The plate boundaries were provided by Bird[52].The arrows indicate horizontal plate motions[46].The scale refers to the parallel 20 • S in Albers projection.

Figure 2 .
Figure 2. (a) The time series of radial displacements of each GPS station sorted by latitude from north to south.The right y-axis shows the approximated values of the latitude ranges.(b) The number of stations available at each epoch (day) available for the inversion.

Figure 3 .
Figure 3. (a) Recharged groundwater and the Earth's surface; (b) subsidence induced by groundwater drop (e.g., pumpage); (c) shows the uplift due to the groundwater recharge; (d) permanent subsidence due to the irreversible deformation caused by the compaction of aquitards.Source: Adapted from Galloway et al. [69].

Figure 4 .
Figure 4. Schematic of the experimental design using a closed-loop simulation.The top shows the arbitrary model used for the computation of synthetic loading data with a pattern of alternate high (+250 kg/m 2 ) and low (−250 kg/m 2 ) attenuation elementary squares, which is the pattern used in the so-called "checkerboard test".The middle panel shows the displacements due to the loading induced by the synthetic loading patterns at each GPS station.The bottom panel shows the inverted model obtained from the synthetic data corresponding to the model on top, which is identical to the initial model.

Figure 5 .
Figure 5. Flowchart presenting the main steps of the experiments based on the methods and datasets described in Sections 3.1 and 3.2, respectively.Comparisons are presented in the following sections accordingly.4.1.Experiment 1: What Is the Potential Spatial Resolution of GPS-Imaged TWS over South America?

Figure 6 .
Figure 6.Top panels show (a) synthetic checkerboard patterns with positive and negative mass loadings, (b) inverted checkerboard patterns, and (c) residuals of the differences between "observed" and inverted mass loadings at a spatial resolution of 1.0 • × 1.0 • .Middle Panels (d,e,f) show the same results as Panels (a-c), respectively, but for a spatial resolution of 2.0 • × 2.0 • .Bottom Panels (g-i) show the same results as Panels (a-c), respectively, but for a spatial resolution of 3.0 • × 3.0 • .The dashed lines indicate a buffer zone of 5.0 • w.r.t. the boundary of South America.

Figure 7 .
Figure 7. Top panels show (a) gridded mass loadings expressed by the RMS of the TWS series from GLDAS-Noah, (b) inverted TWS recovered from the 397 forward displacements, and (c) errors quantified as the difference between the input (a) and inverted loads (b) at a spatial resolution of 1.0 • × 1.0 • .Bottom Panels (d,e,f) show the same results as Panels (a-c), respectively, but for a spatial resolution of 3.0 • × 3.0 • .The dashed lines in Panels (a,b,d,e) indicate a buffer zone of 5.0 • with respect to (w.r.t.) the boundary of South America, and the dot-marks in Panels (c) and (f) show the distribution of GPS sites.

Figure 8 .
Figure 8.(a) Daily errors in terms of the RMSE for the GPS-imaged TWS w.r.t. the input model of the RMSs of the GLDAS-Noah TWS series.(b) The same as (a), but the input model is the GLDAS-Noah TWS series for each day corresponding to the available GPS solutions.

Figure 9 .
Figure 9. (a) Daily errors in terms of the RMSE for the GPS-imaged TWS w.r.t.GLDAS-Noah "TWS".(b) The GPS-imaged TWS over South America for 12 February 2017, equivalent to the inversion with the minimum RMSE as shown in (a).(c) GLDAS-Noah TWS for 12 February 2017.(d) Residuals as the differences between GPS-imaged TWS and GLDAS-Noah TWS.The black dots indicate the distribution of the GPS sites used for the inversion.

Figure 10 .
Figure 10.(a) Daily differences between GPS-imaged TWS and GLDAS-Noah TWS at each GPS site.(b) The coefficient correlations (circles) and the RMSEs (colored circles with amplitudes depicted by the color bar) between the GPS-imaged TWS and GLDAS-Noah TWS at each GPS site.

Figure A1 .
Figure A1.(a)The maximum correlations (circles) between GPS-observed radial displacements and GLDAS-TWS; the color bar depicts their respective lags (circles' colors).Positive lag values mean GLDAS-TWS leads GPS radial displacements, and conversely, negative lag values mean GLDAS-TWS lags GPS radial displacements.(b) The maximum correlations (circles) between GPS-observed radial displacements and precipitation (TRMM), and the respective lags (circles' colors) in days.(c) The maximum correlations (circles) between GLDAS-TWS and precipitation (TRMM) and the respective lags (circles' colors) in days.

Table 1 .
Summary of the results presented in Sections 4.1-4.3.The check marks (ticks) and stand for "yes" and "no", respectively.