GPS Imaging of Global Vertical Land Motion for Studies of Sea Level Rise

We estimate the rates and patterns of vertical land motion (VLM) on all locations on Earth's land surface using GPS Imaging. The solution is based on a large database of uniformly processed GPS data from solutions that are aligned to the International Terrestrial Reference Frame. We provide global maps and estimates of VLM at all tide gauges of the Permanent Service for Mean Sea Level to better constrain the difference between geocentric and relative sea level rise. To enable critical assessment of the VLM estimate, the temporal and spatial contributions to rate uncertainty and variability are generated and included for every gauge. Seasonality and trends of uplift are assessed and found to be strongly correlated with observations from gravity data suggesting that loading from the terrestrial hydrosphere is a dominant driver of non‐glacial isostatic adjustment (non‐GIA) VLM. Although stations are dominantly concentrated at subsiding parts of continents, GPS Imaging geographically balances VLM signals, correcting for bias associated with network distribution. This allows us to make a global assessment of the budget of uplift and subsidence attributable to GIA and non‐GIA sources. We show that the surface motion of the continents is on average upward, implying that the unobserved areas (composed of the ocean basins and ice‐covered areas) move on average downward with respect to Earth center. However, after correcting for the GIA the reverse is true, and observed areas subside on average implying that the unobserved areas undergo net non‐GIA‐related uplift.

Many processes contribute to vertical motion of Earth's surface (Pfeffer et al., 2017). The one with the greatest impact on global sea level is the steady motion arising from glacial isostatic adjustment (GIA) following the loss of the continental ice sheets, which has horizontal wavelengths of thousands of kilometers (e.g., Caron et al., 2018;Milne & Mitrovica, 1998;Peltier et al., 2015). However, many other processes contribute to VLM and these have a variety of rates, patterns, and wavelengths. These include tectonic uplift (e.g., Beavan et al., 2010;Bürgmann et al., 2006;Marshall et al., 2017;Serpelloni et al., 2013), isostatic adjustment following erosion (England & Molnar, 1990;Small & Anderson, 1995), sediment loading (Ivins et al., 2007), and changes in mantle flow or dynamic topography (Arnould et al., 2018;Faccenna & Becker, 2010;Hoggard et al., 2016;Kreemer et al., 2020). Shorter term and unsteady VLM changes can arise from elastic response to contemporary ice loss owing to climate change (e.g., Argus et al., 2011;Jiang et al., 2010;Larsen et al., 2005), trends in storage within aquifers and hydrocarbon reservoirs (Faunt, 2009;White & Morton, 1997), surface loading owing to changes in the terrestrial hydrosphere (Amos et al., 2014;Borsa et al., 2014;van Dam et al., 2001), sediment compaction (Carbognin et al., 2005;Dixon et al., 2006;Johnson et al., 2020), and earthquake cycle deformation (e.g., Burgette et al., 2009;Mazzotti et al., 2007;Smith-Konter et al., 2014), to name a few. Some of these processes, such as aquifer depletion, can be caused or exacerbated by human activity (Sneed et al., 2013). Thus, while GIA represents the single largest factor driving contemporary VLM, many other processes contribute, making direct measurement to capture its full complexity and time variability necessary. Here, we evaluate the budget of VLM on a global scale to assess how much of it is attributable to non-GIA processes.
Increasingly, global data sets and models are being brought together to improve understanding and project regional sea level rise into the future (Hamlington et al., 2020). To best support and accelerate these efforts, VLM products that are available with global scope are needed, with the greatest possible coverage and precision. Recent improvements have resulted from the continued proliferation of high-precision GPS networks, which are trending toward global coverage with stations having longer records (up to 25 years) at a greater number of locations. Network improvements are coupled with developments in analysis, software, data products, and shifts toward open data policies which continue to improve the estimates of trends in the data for multidisciplinary science applications . These improvements allow imaging of VLM with GPS at resolutions that have utility for regional assessment of sea level rise in a large fraction of Earth's coastlines. However, owing to the heterogeneous distribution of stations and complexity of the signals questions remain about how precise are the VLM estimates at any given location.
Here, we present an analysis of a globally extensive database of GPS data, from which we derived vertical rates using the MIDAS trend estimator . From these rates, we estimate a VLM field using GPS Imaging and a weighted median interpolation method . These methods capture VLM signals by relying on robust nonparametric (nonleast squares) methods which are insensitive to seasonality, undocumented steps in the time series, outlier positions, or velocities. The estimates provide VLM at any location, the precision of which depends on several factors, including distance to nearest GPS station, number of nearby stations, duration of observation at those stations, and spatial and temporal variability of land motion. We pay special attention to quantifying the uncertainties in VLM and provide metrics that allow users to identify sources of uncertainty. These results can be used to correct trends seen at tide gauges, identify gauges not effected by VLM, or in analyses where trends from tide gauges, and satellite altimetry are compared to VLM (Cazenave et al., 1999;Nerem & Mitchum, 2002;Ray et al., 2010). They can also be combined with constraints from other techniques (e.g., InSAR and leveling) that work at shorter geographic scales to align constraints to the global reference frame (e.g., Bekaert et al., 2018;Hammond et al., 2018;Shirzaei & Bürgmann, 2018). and other products from the Jet Propulsion Laboratory. All the products are based on, and the solutions are aligned to, the IGS14 realization  of the ITRF14 reference frame . To correct for tropospheric effects, we used the VMF1/ECMWF mapping function with elevation-dependent data weight, nominal dry and wet zenith delays from the numerical weather model (Boehm et al., 2006). Antenna adjustments were made with the IGS14 ANTEX calibration table that has a new algorithm that corrects for antenna and radome fields available in RINEX headers (Schmid et al., 2007). The models include adjustments for solid Earth tides, ocean and atmospheric tidal loading, but not nontidal ocean or atmospheric loading. Additional details about treatment of metadata, data editing, ambiguity resolution, antenna phase center calibrations, and estimation strategy including the atmospheric parameters are provided in Kreemer et al. (2020) and are not repeated here.
To estimate the rate of vertical movement, we use the nonparametric MIDAS algorithm  on the IGS14-aligned time series. MIDAS calculates the median trend from the difference in pairs of vertical position data with intersample duration of ∼1 year. This makes the MIDAS rates insensitive to seasonality and steps in the time series even if they are undocumented, as long as the steps are sparse enough in time. The method is robust in the sense that the final result is insensitive to outliers up to a breakdown point (up to ∼25% incorrect data). The resulting velocity field is characterized by fewer stations with outlier velocities, making them more suitable for imaging, geophysical modeling, and interpretation (e.g., Hammond et al., 2016). We base the analysis presented here on the MIDAS file downloaded from our own GPS data products system on December 22, 2020. This file has 17,172 stations with sufficient data to obtain a MIDAS velocity ( Figure 1) and forms the basis for estimation of a VLM field.
The updated models and products used in our analysis have, compared to the previous generation of solutions in IGS08, resulted in significant reductions in RMS residual scatter of vertical position time series. The improvement is in part due to more effective partitioning of signals between atmospheric effects and real Earth motion (Martens et al., 2020), which translate directly to velocity precision and reduction in bias associated with sources not related to solid Earth motion. All the derived solutions produced by NGL are available at http://geodesy.unr.edu, including time series, time series plots, station maps, tables of site coordinates and durations of observation, and MIDAS velocities.

Site Selection and Coverage
We omit from the data set stations that have time series duration shorter than 3.5 years and mean completeness (number of observations days per year divide by 365 days) of less than 30%. We omit these stations because they offer less constraint on VLM and keeping them in our database of station locations may give an overly favorable picture of the spatial distribution of stations which we analyze below. After these exclusions, 90% of the stations have completeness greater than 68%. The median/maximum time series duration is 10.2/26.5 years, 90% have their first position on or after the year 2002.0, and the median midtime of the HAMMOND ET AL.  time series is 2014.7 ( Figure 2). Thus, the data constrain VLM that occurred mostly after 2002.0. We also omit stations whose motion is known to us to not be representative of the solid Earth. For example, stations on flowing ice that have been deployed to study or serve as reference on glaciers in Antarctica or Greenland, and stations on oil platforms, but these collectively affect only 29 stations.
GPS Imaging relies on Delaunay triangulation of the network, which performs best when stations are not collocated. Therefore, prior to imaging we decluster the stations, creating single representative metastations with a vertical rate that is the weighted median value in the cluster when stations are within 0.1 km of one another. The uncertainty of the cluster rate is determined as if it were a weighed mean of HAMMOND ET AL.  the contributing rates, where the weights are one over the uncertainties, making large uncertainty rates unlikely to contribute. The median difference between the maximum and minimum rate in a cluster (excluding cases with only one station in the cluster) is 0.7 mm year −1 and 90% are less than 2.9 mm year −1 . Though a few maximum differences are large they are usually filtered in this declustering stage or later if they are rate outliers. Most clusters have only two stations so the one with the lower uncertainty will be selected via the weighted median. This and the selection criteria described above results in a reduction in the effective number of stations to 13,006, where 3.2% of stations represent clusters. In this set, 98% of the vertical rate uncertainties are less than 2 mm year −1 , and 68% are less than 1 mm year −1 (Figure 2).
Most of the VLM signal lies within a narrow range, with 97% between −10 and 10 mm year −1 . A histogram of the velocities is shown in Figure 2. The median vertical rate at GPS stations is −0.60 mm year −1 . The slightly negative median is significant with respect to the uncertainty in the estimate of the sample median. However, it is negative owing to the bias introduced by station concentrations in subsiding areas such as the central United States which is widely affected by GIA fore-bulge collapse over a large segment of the country. Europe subsides in the southern half of the continent for similar reasons (Lambeck et al., 1998;Peltier et al., 2015). GPS Imaging estimates VLM on a grid which mitigates the bias associated with station concentrations. This is important because, as we show below, the median VLM of the continent is strongly positive and would be misrepresented if the heterogeneous station distribution were not taken into account.

How Much of the Land Area and Coasts Are Covered?
Owing to the heterogeneous station distribution, the quality of VLM constraint varies greatly among and within the continents. Figure 1 shows that the density of station coverage is much greater in the Northern Hemisphere, especially in North America and Europe. The continental interiors have large areas that are poorly sampled, especially in Asia and Africa. To objectively portray the variability in geographic coverage, we show in Figure 3 the distance from every point on land to the nearest GPS station that passed our exclusion criteria. These distances range from about zero to a maximum of 13.9° (∼1,540 km) in the Saharan Desert of northern Africa near the border between Chad and Libya.
The percentage of land area that is within a given distance to a GPS station in shown in Figure 3b. This chart is created by summing the area of Figure 3a that has values less than or equal to the value on the horizontal axis of Figure 3b, excluding oceanic areas. For example, on the major land masses, 52.7% of land area has no GPS station within a distance of 2° (∼222 km), and 84.5% of land area has a GPS station within 5° (blue bars in Figure 3b). However, the perimeter of continents tends to be better covered, because of the concentration of populations (and GPS stations) near coasts. So, we repeat the analysis except by taking steps of median size 28 km along the coastlines, summing the length of linear coasts that are within the given distances on the horizontal axis of Figure 3b, divided by the total length of all perimeters (red bars in Figure 3b). This indicates that across the globe 73.9% percent of the linear coastlines have a GPS station with 2°, and 93.0% have a GPS station within 5°.
Despite the presence of hundreds of GPS stations on islands, the oceans are generally poorly covered. In the ocean basins, the maximum nearest distance from a permanent GPS station is 26.7° (∼2,970 km), about twice the maximum distance on land. That location is in the southeast Pacific Ocean (at −49.75° latitude, −137.75° longitude). For this reason, we omit the ocean basins from the imaging. While sea-floor geodesy is becoming a more viable mode for observing Earth deformation (e.g., Bürgmann & Chadwell, 2014), the GPS networks are unlikely to proliferate onto the sea floor at the rate they have on land. Thus, we expect the ocean/land sampling bias to be a long-lived feature in the analysis of global VLM. For similar reasons, we block out areas covered with ice, including the interiors of Greenland and Antarctica, since there are no bedrock observations there. However, in Antarctica, we preserve the area of the Transantarctic and Ellsworth mountains which have some bedrock GPS sites that can be used to measure uplift (Argus et al., 2011;Bevis et al., 2009).

GPS Imaging
A benefit of GPS Imaging is that stations in the continental interiors are used to obtain insights about variability of signals and processes that effect VLM at the coast. The method has been documented elsewhere , but we briefly describe the essential aspects here since they may impact the interpretation of the images. There are two major steps. In the first step, we apply weighted median spatial filtering to the vertical MIDAS GPS rates, and in the second step, we perform the imaging, where the same algorithm is applied to points on a grid.
In both steps, we require a spatial structure function (ssf) which represents the degree of similarity in VLM rates as a function of distance between locations. The ssf itself is a function of distance Δ ij between GPS station i and evaluation point j, which varies from 0° to 180°. The form of the ssf is based on the expectation that two GPS stations in the same location will have the same rate, and stations separated by great distance are expected to have weak similarity. Stations in the same location will be similar to the degree expected from their aleatory uncertainties. Stations separated by a great distance will have less similarity but are sampled by the global GPS network, so there is still some degree of expectation of similarity (e.g., described by the histogram in Figure 2e). In our construction of the ssf, it is forced to have a form such that ssf(0°) = 1 and ssf(180°) = 0. We determine the values of the ssf between Δ ij = 0° and 180° empirically from the data. It can have any shape, but we require that it decreases monotonically. In the supplementary materials, we describe in detail how we generate the ssf for each continent and additionally one for the entire data set to use at locations not on major land masses. Once the ssf is determined, we use it to estimate the weighted median velocity v j  at each evaluation location j: where v i are the set of vertical MIDAS velocities at GPS stations i that are connected to the evaluation point j via Delaunay triangulation. For each evaluation location, the weights are computed using Equation 2 using the MIDAS rate uncertainty  i . To perform the filtering in the first step, we replace each velocity v i with the weighted median of the original velocities at neighboring locations.
In the second step, we estimate the VLM rate at grid point locations using the filtered data from the neighboring GPS stations as connected via an ad hoc Delaunay triangulation with the grid point location temporarily added to the list of station locations. We calculate the weighted median of the filtered MIDAS station velocities adjacent to the evaluation point. This estimate is free from prior assumption about the functional form of the signal since it is made independently at each grid point. This is desirable when imaging complex fields that may contain signals from multiple processes that operate at different amplitudes and wavelengths, possibly by orders of magnitude (e.g., GIA vs. aquifer deformation). It also makes the estimate insensitive to velocity outliers and geographically local, sensitive to the nearest GPS stations, and resolution dependent more on station density than on grid parameterization.

Spatial Variability of VLM Signals
In practice, the importance of the ssf is low when the GPS network is dense and velocity uncertainties are similar, so the weights will vary little among the stations connected to the evaluation point. However, the importance of the ssf is larger where the GPS station spacing is greater than the horizontal length scale of the VLM signal, and/or the station density is variable. We can quantify this spatial part of the VLM rate uncertainty as which is unitless since it is the median ssf value (all between 0 and 1) at the set of distances Δ ij from stations i connected to the evaluation point j. In Figure 4, we show at every grid point on land the value of ˆj. This map illustrates a factor unrelated to the uncertainties in the GPS velocities that influences how well the VLM is resolved at each location. Maps of ˆj can be thought of as a replacement for a checkerboard test of the kind commonly used in geophysical imaging. It is more direct than a checkerboard test because it does not presume any specific length scale or pattern of VLM signal (such as checkers) against which to compare resolution. Rather it compares the network to the dominant wavelength of actual signal observed in the VLM field. We note that the nearest GPS station may have a larger ssf value and contribute to the local VLM estimate, so ˆj may be a somewhat pessimistic measure in some locations. When we create an alternative map (not shown) keyed to the minimum ssf instead of the median ssf, it has higher values in the sparsely covered regions, suggesting better resolution of VLM but based only on the nearest station. We construct Equation 4 to represent something close to the center and body of constraint on VLM at each location, instead of to its single nearest station.
In Figure 4, blue domains represent places where the GPS stations are generally dense enough to resolve horizontal wavelengths of VLM signal observed in the data. These areas include North America (excluding the Arctic Archipelago and Greenland), most of Africa near the coast, southern Europe, Britain, Japan, and northeast New Zealand. Some areas including most of Asia (aside from South Korea) and parts of Antarctica are not well covered. Some areas with many stations have only medium quality of coverage according to the ˆj parameter, these include South America, Australia, and parts of Antarctica.

Seasonal VLM Variability
Seasonality in vertical motion of the Earth surface has its source in multiple processes, including loading from the terrestrial hydrosphere, atmosphere, thermal changes, and aquifer pumping (e.g., Blewitt & Lavalleé, 2002;Bos et al., 2010;Davis et al., 2012;Prawirodirdjo et al., 2006;van Dam et al., 2001). Here, we estimate the amplitude and phase of annual vertical motion from the residual time series for each station after removal of other signals discussed above such as trends, steps, and exponential transients from large earthquakes. We use the relation: where x(t) is the vertical position as a function of time t. We solve for sine S 1 and cosine C 1 terms then convert them into annual amplitude A u and phase ϕ u of the vertical signal: where atan2 is the four-quadrant inverse tangent, taken modulo 2π. For convenience in interpretation of the phase in terms of seasons, we multiply by the factor  365/2 to convert from radians to the day of calendar year ϕ u . To simplify comparison to gravity data, which exhibit seasonal maximum equivalent water heights when the seasonal GPS station height attributable to loading is at minimum, we add π so that the resulting day represents the time where GPS station elevation is minimum. Because S 1 and C 1 are inferred from multiyear time series, these are the mean annual amplitudes and phases for the station across the period of observation.
We then use GPS Imaging to perform robust filtering and interpolation to map ϕ u and A u on the continents ( Figure 5). The amplitudes are large (8 mm and greater) in locations where seasonal precipitation is large such as in the Amazon Basin, southern Alaska, southeast Asia, and sub-Saharan West Africa. The HAMMOND ET AL.  amplitudes are small (2 mm or less) in relatively dry regions such as the Saharan Desert, southwest United States, and southern Europe.
These signals show a strong correspondence to the annual terms resolved from Gravity Recovery and Climate Experiment (GRACE) satellite data (Luthcke et al., 2013;Watkins et al., 2015). We used time series from the mascon solutions (Loomis et al., 2019) to solve for the annual amplitude and phase using the same method as for the GPS time series. We interpolate the result from mascon centers to a regular grid to map the seasonal amplitude and phase (the day of year the equivalent water height is at a maximum- Figure 5). For the amplitudes, we see that there is a very strong spatial correlation between GPS and GRACE seasonal amplitudes, suggesting the dominance of seasonal mass loading changes driving the GPS seasonal vertical motion.
For the phase signals, there is a strong dependence on latitude, and for large sectors of the Earth's terrestrial surface, the timing is strongly associated with timing of seasonal precipitation. There is a very strong correlation between the timing of GPS and GRACE changes, again consistent with seasonal hydrological loading being the driving mechanism. Areas with a large proportion of their precipitation occurring in Northern Hemisphere winter such as northern North America and Europe have their lowest positions in late spring when the most water has accumulated and driest periods in the late autumn. Areas with a large proportion of their precipitation budgets coming from summer monsoon precipitation such as southeast Asia, the southwest United States, Mexico, and northern Africa have their driest periods in the spring, out of phase with other areas. There is a lack of symmetry between the hemispheres. For example, Amazonian South America (between −30° and 0°) and northern North America (between 40° and 70°) are both at minimum elevation near day of year 150. This is the time of Northern Hemisphere spring and Southern Hemisphere autumn, but in both locations is the time of year that ends the season with heaviest precipitation.
The seasonal linking of VLM and hydrological loading indicates that they are also linked during multiannual changes in hydrological loading, such as drought periods. Thus, trend changes in terrestrial hydrology HAMMOND ET AL.  (Loomis et al., 2019). There is a very strong global correspondence between the amplitude and phase of seasonal signals in GPS and GRACE. Note that factor of 10 change in scale between GPS and GRACE and that phase of GPS refers to day of year the station is at its minimum height in its cycle. See text for discussion.
can impart trend changes in VLM, and the duration of those changes depends on the duration of change in precipitation. These inflections can be observed with GPS as nonlinear signals (Adusumilli et al., 2019;Borsa et al., 2014;Hammond et al., 2016;Silverii et al., 2016) which we examine on a global scale in the next section.

Nonseasonal VLM Variability
Owing to various factors including hydrologic loading, there is often uncertainty about the extent to which trends in GPS data reflect long-term steady geologic movement. An indication that GPS data may not represent long-term geodynamic processes is when its velocity is not constant, suggesting that long-term motion is overprinted by transient processes. We can quantify the steadiness of GPS velocities by computing for each time series an index of rate variability. For each time series, we take 2-year long subintervals and compute the MIDAS rate for each interval. We start with the first 2 years of the time series, then select a new interval by shifting the window forward in time by 1 month and repeat the process until the window reaches the end of the time series. The 2-year window size is selected because it is long enough to include enough data pairs to compute the MIDAS rate, yet short enough so that many rate samples can be generated for the time series in our data set. The result is a time series of velocity for the station which indicates how the rate varied over the observation period. Because of the properties of the MIDAS trend estimator, this variability index is insensitive to outliers and annual periodicity. Hence, this index is a robust measure of nonseasonal velocity variability (υ i ) in velocity.
The index is computed as the median absolute deviation (MAD) of the values in the velocity time series:  Figure 7 shows a histogram of υ i for all stations, indicating 50% are less than 1.4 mm year −1 . A histogram of υ i normalized by the uncertainty in the MIDAS rate for the entire time series and multiplied by 1.4826 to scale the MAD so it represents a standard deviation in a Gaussian distribution indicates that 43.2% have values greater than 2. This indicates that the temporal variability of vertical rate is significant for a large fraction of the stations. However, only ∼5% have variability higher than station QUIN. To see where on Earth variability most affects the rates we perform the imaging on υ i to obtain gridded values ˆj (Figure 7). The robust imaging prevents individual stations with outlying υ i dominating, so only variability present in multiple adjacent stations and representative of large-scale land motion are seen.
We see high variability in locations that are known to be influenced by hydrological mass changes (Figure 7). For example, previous work based on GRACE and GPS saw that hydrological loading changes drive seasonal and multiannual motions in the Amazon basin of South America (Ferreira et al., 2020;Fu et al., 2013). Since ˆj represents nonseasonal changes in motion the Amazon variability indicates multiyear loading changes that affect rates on a near-continental scale. There is also high variability in other areas affected by hydrological mass changes, including the Greenland coast, Antarctica, and Central California. In some cases, variability can be high in regions where VLM is low, and vice versa. For example, we see low variability in central Canada, where uplift is caused by GIA, which is steady over time.

Spatial Scatter in VLM Rate
An independent indicator of uncertainty in VLM estimates is the spatial variability of vertical rates within clusters that are used to image rates at evaluation points. Because our estimates are based on medians of these clusters, this nearest neighbor variability is a measure of noise that is mitigated by the spatial filtering.
Here, we quantify this spatial variability and show it has very weak dependence on the temporal variability, which indicates that it must be separately accounted for when assessing the quality of VLM estimate.
We define the nearest neighbor spatial variability ζ i at station i as the MAD of the rates at neighboring stations that are connected to i by Delaunay triangulation: This kind of spatial variability is strongly affected by interstation separation distance and is expected to be noise that is very different in structure compared to the signals encoded in the ssf which is estimated HAMMOND ET AL.
10.1029/2021JB022355 11 of 26 from rates from stations separated by large distances. Figure 8 shows the relationship between ζ i and υ i in a scatterplot of the values in a logarithmic scale. The plot shows that there is a relatively weak correlation (r = 0.240) between the log 10 (ζ i ) and log 10 (υ i ). The correlation between ζ i and υ i is even weaker (r = 0.028). Thus, the spatial and temporal variabilities are generally distinct and separate indicators of noise contributions, and both should be considered when evaluating the reliability of GPS-based estimates of VLM trends. Below, we show an example of how to use the variabilities as a tool for diagnosing the source of VLM uncertainty, showing how reliability of rates can be evaluated.

Summary of Uncertainty Factors
The uncertainty components have symmetries (summarized in Table 1) that we note in an attempt to make them easier to understand and remember. We distinguish the sources of uncertainty as temporal (σ, υ) or spatial (ζ, η) in origins, and also whether the uncertainty is aleatory, that is, estimated from the scatter around the estimated value (σ, ζ) or from signal variability in the (spatial or temporal) neighborhood in which the estimate is made (η, υ). We distinguish signals at GPS stations (e.g., v i ) from signals imaged from multiple stations with a carat (e.g., ˆj v ). Having these defined clarifies the sources of uncertainty in VLM at specific locations. As an example, we consider the tide gauge in San Francisco, CA (longitude −122.465°, latitude 37.807°). The constraints on VLM at this location are given in Table 2. The table lists the GPS stations that contribute to the estimate of ˆj v (imaged GPS rate at the tide gauge location), including the important values that characterize the rate and uncertainty. In this case, the five stations (CADY, CCSF, MDHL, PBL1, and UCSF) have rate uncertainties that vary from 0.61 to 0.95 mm year −1 . The stations are all relatively close to the tide gauge compared to the wavelength of the signal, as indicated by η i = 0.913 for all stations, with the closest station (MHDL) having η i = 0.938. The stations contribute with slightly different weights (varying from 0.172 to 0.268) to the median value of ˆj v = −0.64 mm year −1 . Nearest neighbor variability ˆj and temporal variability ˆj are in the middle of the ranges for these parameters (Figures 7  and 8). Thus, the overall estimate of VLM is typically well constrained at the San Francisco tide gauge, since the contributing uncertainties are near the middle of their ranges. Below, we illustrate with an example where the VLM estimates at the tide gauges are ranked in quality.

Results
We show the result of GPS Imaging the vertical rates with the oceanic areas masked out because coverage in the ocean basins is generally poor (Figure 9). We do not exclude stations on islands in the ocean basins from the analysis but mask out the resulting estimate of the vertical rate field. Thus, stations on islands near the coast can contribute if they are close enough to have a significant weight. In this section, we discuss some features of the signals, recognizing that some have been described previously, others not. While a complete discussion of all aspects of VLM on Earth is beyond the scope of what is possible to cover in this manuscript, we point out some important features of the rate field. In Discussion, we discuss the global balance of vertical motion across the Earth's surface, and how much of this motion is attributable to non-GIA sources.

GIA
Immediately clear and of greatest impact are the GIA uplift signals of Fennoscandia in Northern Europe, northern North America, and Antarctica that have been documented in many previous studies (reviewed in Whitehouse (2018)). These signals have been measured extensively using multiple techniques, including GPS, leveling, gravity, raised shorelines, and tide gauges, and dominate the uplift field with maximum values over 10 mm year −1 that span thousands of kilometers. VLM from GIA also includes large areas undergo-HAMMOND ET AL.  7) consistent with the processes being a long-term steady feature of mantle flow.

Subduction Zones
Coastal areas at convergent tectonic plate boundaries experience vertical motions associated with interseismic strain accumulation at subduction zones. The islands of Japan and Sumatra are heavily affected and have high degrees of velocity variability (Figure 7) indicating effects from their high rate of recent seismic activity. While ˆ is insensitive to the coseismic steps themselves, these events initiate postseismic mantle relaxation which cause vertical velocities to evolve over months to decades (e.g., Pollitz et al., 2006;Yamagiwa et al., 2014). The Nazca subduction zone in South America exhibits the Earth's longest section of onshore interseismic uplift that is nearly continuous from southern Chile to southern Peru (Figure 9). At the northern end of this zone, the imaging of interseismic uplift becomes equivocal in Peru where there are fewer coastal GPS stations available and also exhibits a medium degrees of velocity variability in the stations that do exist (Figure 7). At this latitude, horizontal GPS measurements indicate that the degree of coupling of the plate interface undergoes a transition from locked to creeping (Chlieh et al., 2011), so onshore uplift may not extend further north. In the Cascadia subduction zone of the Pacific Northwest of the United States, we also see interseismic uplift. However, it is difficult to see in a map of the scale of Figure 9, especially since the uplift occurs only near the southern and northern ends of Cascadia, where the trench is close enough to Figure 9. GPS Imaging of vertical rate on land areas, blue is downward, red upward.   (Burgette et al., 2009;Mazzotti et al., 2007). Velocity variability of Cascadia coastal GPS stations is generally low, consistent with steady uplift and the lack of recent slip of the megathrust.
Coastal subduction-related uplift or subsidence may, according to our GPS measurement of VLM variability, appear steady. However, it is transient when considering time scales longer than the seismic cycles, for example, ∼500 years in Cascadia (Atwater & Hemphill-Haley, 1997), when earthquakes reset the cycle so that areas with net uplift will subside and vice versa. Thus in these cases, projections of sea level rise over time scales longer that the period of GPS observation would be misled by the GPS rates, which should be considered as valid until the time of the next large earthquake.

The Amazon Basin
The Amazon Basin has been shown with GPS and GRACE data to experience very large seasonal fluctuations in vertical and horizontal motions (Ferreira et al., 2020;Fu et al., 2013). We find, in addition to the large seasonality of vertical motion, a very high degree of nonseasonal uplift rate variability, the largest and most variable on Earth in the imaging with the greatest variability located at the Amazon River ( Figure 6). This indicates that changes in rates in the Amazon not only are seasonal but also have trends that inflect upward and downward, changing sign on multiannual time scales between drought and wet periods. This may explain why when considering overall trends Knowles et al. (2020) found an uplift trend in the Amazon, while our results show a combination of uplift and subsidence. This is consistent with the imaged temporal variability and indicates that the trends in the Amazon Basin are likely not permanent features of the vertical GPS rate field, and the sign of uplift depends on the time period of observation.

Australia
In Australia, uplift rates change in sign depending on the time interval observed. Han (2017), for example, found that the Australian continent was depressed by the La Niña water load in 2010-2011 and afterward rebounded as the water load evaporated. When considering the overall trends, however, we find mostly subsidence in Australia over our period of observation, where the median value of the image is −0.8 mm year −1 , and 97.8% of the pixels in the image have negative vertical rate. There is low variability in the central and western part of the country and moderate variability on the eastern part of the country, consistent with some impact from the climate-related loading sources. Nearly, all Australia subsides (with the exception of a patch of the Northern Territory) but the trend is likely not a long-term effect. Also, these rates have been shown to be inconsistent with the effects of GIA (Riddell et al., 2020), and the moderate variability suggests a time-variable source, which is also inconsistent with GIA. Riddell et al. (2020) also discussed numerous potential noise sources that could result in a downward trend for all of Australia, including biases in GPS from draconitic period signals, reference frame origin and geocenter drift, colored noise, but did not settle on single likely source. If the land surface in Australia continues to subside at the imaged rate most of Australia, which is the lowest continental land mass in the world, would be submerged in half a million years. This is not likely to happen since similar subsidence has not occurred in the recent geological past (Sandiford, 2007). However, subsidence over the time scale of the glacial cycles would not continue as long and may be feasible. If VLM oscillated with periods longer than observed with GPS, then the GPS data would not be sensitive to this variability. GRACE data indicate that the western half of Australia has experienced steady loss in equivalent water height in 2003-2016, while the eastern half experienced mass gain (Luthcke et al., 2013;Watkins et al., 2015) suggesting the continent-wide subsidence is not a single process related to hydrological loading or poroelastic deformation. Given that GRACE and vertical rate variability differ, it may be that the imaged subsidence is driven by different factors in the west and east halves of Australia.

Southwest Africa
The GPS network in Angola has 16 stations, almost all of which provided data for the same interval, from 2010 to 2015. During that time nearly all of the stations experienced subsidence at rates between 2 and 6 mm year −1 and are responsible for the subsiding area in southwest Africa (Figure 9). Hydrogeology data for Angola indicate the presence of a group of large moderate-to-high productivity aquifer systems in unconsolidated and Tertiary-Quaternary sediments (Africa Groundwater Atlas, 2020;MacDonald et al., 2012). While the GPS data consistently indicate subsidence, they are geographically sparse and have near the minimum of completeness and duration we allow in our analysis. Their median completeness is 31%, that is, on average 116 days of observation per year, which should be enough to prevent bias in the estimate of the rate and its variability. The temporal variability is moderate in this part of Africa (Figure 7) further suggesting that the aquifer affects vertical motions. Observations of trends in GRACE data from 2010 to 2015 in mascons in Angola indicate a complex time series of water equivalent height changes that increase and decrease over time but have a general decrease over these years. Decrease in overall water mass associated with surface subsidence and high temporal GPS rate variability is all consistent with the effects of porous aquifer depletion. Similar signals of subsidence and mass loss at other locations have been attributed to aquifer depletion and compaction, for example, in the California Central Valley (Famiglietti et al., 2011;Faunt, 2009).

Southern Africa
In Southern Africa, there is a spatially coherent anomaly with peak uplift rate of about 1 mm year −1 in the northern part of the country of South Africa, southern Zimbabwe, and eastern Botswana. The amplitude of uplift decreases toward the coastline on both east and west coasts of South Africa. The anomaly is well imaged by the network, with 79 stations having a median observation duration of 13.4 years. The neighboring countries of Namibia, Botswana, Zimbabwe, and Mozambique have fewer stations in the NGL database, with between 0 and 10 stations in each country, so constraints on VLM are not as strong. The East Africa Rift (EAR) has many GPS stations but does not show an uplift anomaly. The transition between upward and downward moving domains lies along the EAR plate boundary where it trends along Lakes Tanganyika and Malawi toward the Indian Ocean. This zone is the extensional boundary between the Rovuma microplate and the Nubian plate (Saria et al., 2013;Stamps et al., 2018). At this boundary, the extensional EAR gives way southward to the Kalihari Craton, where the lithosphere and sublithospheric mantle have faster seismic velocities (Paysanos & Nyblade, 2007) indicating it is a major transition in lithospheric structure across the Rift.
South of the EAR, within the uplift domain, Figure 7 shows low to medium levels of temporal variability, but to the north of the EAR variability is higher. This suggests that the EAR may be affected by climate-related uplift signals that are not sufficient to drive uplift trends. Trends in equivalent water thickness from GRACE in South Africa and the EAR are mixed in sign, including signals of both increasing and decreasing anomalies in equivalent water height. Thus, spatially the VLM trends are not consistent with the mass change trends, suggesting that the uplift anomaly is not a consequence of mass unloading from climate-related or hydrological loading effects. While the uplift domain has proximity to the aquifer-driven subsidence in Angola (discussed in Section 4.2.4.), it does not surround the zone of subsidence and therefore does not have the spatial characteristics of an unloading feature from water loss centered on the Angola aquifers. However, GPS constraints to the north and northwest of Angola are poor to nonexistent, so a radially symmetric uplift zone could be present but is not detected.
The large size, regional coherence, and low variability of the uplift anomaly in southern Africa indicate that it may be a geodynamic consequence of mantle flow, dynamic topography, and/or lithosphere-scale processes. The location of this uplift domain south of the EAR plate boundary corresponds to a part of the broad region of high elevations and shallow bathymetry termed the south African Superswell (Nyblade & Robinson, 1994). Support of the Superswell has been attributed to buoyancy and consequent vertical flow in the mantle as derived from seismic velocity anomalies (Lithgow-Bertelloni & Silver, 1998). While the exact locus in the mantle of where the support originates, that is, whether primarily from the deep lower mantle low seismic velocity zone or shallower upper mantle features, has not be conclusively determined (Gilfillan et al., 2019), it is roughly collocated with the current location of the Quathlamba hotspot and plume (Hartnady, 1985) suggesting a connection. However, the coincidence of the boundary between uplift and subsidence with the southern EAR suggests a component that could be influenced by shallower lithospheric processes, and thus the entire effect may result from an interaction between asthenospheric mantle flow and surface plate kinematics. But in any case, this class of explanation implies that the uplift is a persistent feature and that the VLM will be a steady contribution that should be accounted for in sea level rise projections for southern Africa.

The Himalaya and Asia
The high mountain region of the Himalayan Mountain chain was built in response to collision between the Indian subcontinent and Asia, and this is evident in its active uplift (Avouac, 2003). This region is not near coasts where sea level rise is a concern, so we do not discuss these signals at length. However, the smearing of these signals is among the most extreme on Earth so bears a brief discussion since it illustrates how GPS Imaging behaves when station density is highly variable. In general, GPS Imaging exhibits no preference for zero rates, that is, the image is in no way biased toward zero, so areas with poor coverage tend to take on signals from nearby measured areas. In China, for example, there are very few GPS stations in the NGL database but there are many on its periphery. Thus, the imaged uplift is largely a function of stations outside China. There are enough stations in Nepal (which has a higher η j , Figure 4) so that the uplift is adequately separated from the signal of subsidence in India (Fu & Freymueller, 2012). Thus, smearing from Nepal is mostly to the northeast into China and south into India.

VLM at Tide Gauges
Although changes in global mean sea surface height are important for understanding the impacts of sea level rise at coasts (Church & White, 2011;Milne et al., 2009), it is the change in sea level relative to the land (RSL) that is, most relevant. Networks of tide gauges provide the most widely used way to measure RSL (e.g., Church et al., 2004;Davis & Mitrovica, 1996;Douglas, 2001;Pugh & Woodworth, 2014). However, tide gauge measurements are affected by movement of the land and of the sea surface, leaving the driving factor behind RSL ambiguous. Sea surface altimetry, which has global coverage, can be used to assess the relationship between RSL and VLM (e.g., Cazenave et al., 1999;Nerem & Mitchum, 2002;Ray et al., 2010) but connecting open ocean sea height measurements to the coast can be challenging (Wise et al., 2018). Regional VLM measurements can be combined with tide gauge and satellite altimetry to better reconcile the budget related to RSL changes by combining all three techniques (e.g., Hawkins et al., 2019). In any case, improving knowledge of VLM at the tide gauges is critical for reconciling the various constraints on RSL and its impact on coastal communities.
Here, we estimate the rate of VLM at the 2,357 tide gauge stations of the Permanent Service for Global Mean Sea Level (PSMSL, Holgate et al., 2013) using GPS Imaging. The estimation points are the tide gauge coordinates instead of the regular grid points used in the GPS Imaging and include locations on islands in the ocean basins ( Figure 10). The patterns of VLM at the tide gauges follow those seen in the individual GPS stations and in the imaging (Figures 1 and 9). Subtracting the rates of VLM from the rates of sea level rise estimated from tide gauges tends to increase the mean rate of sea level rise by 0.2-0.3 mm year −1 (Figure 11) because of the concentration of tide gauges in areas with positive VLM. This remains true whether considering the entire set of the tide gauge data, only data from 1993 to the present, or only data from 2000 to the present. The increase in rate of sea level rise owing to the VLM correction is less than the increase of rate that has occurred in recent decades owing to climate change (Nerem et al., 2018).
The uncertainty of the GPS VLM rate at the tide gauges is a function of the data at the nearest GPS stations, including their number, distance, uncertainty, and variability of their trend estimates. These are given by its uncertainty vector (ˆj, ˆj, ˆj, ˆj) as described in Table 1. Supplementary Table S1 presents a data record Figure 10. Vertical rates at tide gauges of the Permanent Service for Global Mean Sea Level (Holgate et al., 2013). See Supplementary Table S1 for rates at each tide gauge along with list of GPS stations that contributed to the rate and associated temporal and spatial uncertainty parameters.
for each tide gauge with the same information as in Table 2, so that the uncertainty may be assessed. The completeness of this table allows it to be used in subsequent analysis to objectively rank order tide gauges by VLM quality constraints and provides a means to trace back to which sources of uncertainty are most important for a given gauge. Exactly which uncertainty components are most important depends on the goals of a subsequent analysis. For example, in some cases, low temporal variability stations may be required or in others low spatial variability.
We provide an example approach where the tide gauges are rank ordered according to quality constraints, that is, which have good, medium, or poor quality VLM estimates. We use a balanced set of criteria based on the distributions of the four kinds of uncertainty. Figure 12 shows a histogram of each uncertainty parameter, with vertical bars that indicate how each parameter is divided into its tercile ranges (lowest, middle, and highest thirds of the data) and the 5th and 95th percentile ranges. We consider a VLM constraint at a tide gauge "good" if it none of its uncertainty measures lies in the worst tercile. This corresponds to the condition ˆ < 0.42 mm year −1 , ˆ < 0.86 mm year −1 , ˆ < 1.54 mm year −1 , and ˆ > 0.57. In this case, the 33.1% of the tide gauges fall into this "good" category. We call a VLM constraint "poor" if any of its uncertainty parameters ˆ, ˆ, ˆ, ˆ are in the 95% percentile, that is, among the worst 5% of the values. In this case, the 15.5% of the tide gauges fall into this "poor" category. The remaining tide gauges have medium quality, and these comprise 51.4% of the gauges. Figure 13 shows the locations of the gauges with good, medium, and poor constraints. VLM tends to be well constrained (mostly "good") at gauges where GPS networks are dense, but this is not always the case. For example, in California, there are 26 tide gauges in the PMSL, 13 of which have "good" VLM constraint and 13 have "medium" constraint. In Japan, despite the excellent coverage with GPS, about 33% of gauges have "poor" VLM constraint. None of the tide gauges with "poor" VLM constraints are classified as such because of a low ˆ, consistent with the excellent geographic coverage of GPS networks in Japan. About two thirds of these are classified as "poor" because they exceed the threshold HAMMOND ET AL.

10.1029/2021JB022355
18 of 26 for temporal rate variability (ˆ), consistent with the seismic cycle in Japan causing large changes in uplift rates over the observed time.

Discussion: Integrating Global VLM Constraints
A benefit of GPS Imaging is that it balances the contributions from stations that are very heterogeneously distributed. For example, the median vertical rate of the GPS stations is −0.60 mm year −1 , a value that is, significantly different than zero given the large number of data used to estimate it. However, this signal reflects the number of stations that are concentrated in areas that experience subsidence, for example, the large fraction of the central United States that undergoes fore-bulge collapse (Figure 1). While uncertainties in the imaged field ( Figure 9) are not the same everywhere, it does represent a field that is, geographically balanced and can be used to evaluate properties of Earth's global VLM field. Some recent studies have provided global-scale interpolated VLM based on the NGL database (e.g., Husson et al., 2018; Schumacher Vertical bars indicate terciles separating lower, middle, and upper third of data (green dashed lines) and limits of 5% and 95% percentiles of the data (red dashed lines). These are used to divide tide gauge stations into sets with good, medium, and poor GPS estimates of VLM (see Figure 12 and text for discussion). et al., 2018) but tend to focus their concern on the GIA process. Detailed comparison of results from the GPS Imaging and Bayesian surface reconstruction methods has been conducted elsewhere (Husson et al., 2018) so we do not repeat such comparisons here.
We wish to assess the impact of non-GIA so we ask the questions: how much of the global VLM is non-GIA-related? and how significant is it on a global scale? To do this, we first quantify the impact of the VLM field by integrating to estimate vertical volume flux rate: where  r denotes position and ∂S can be any part or all of Earth's surface. This kinematic formula is valid regardless of the driving process behind the VLM, be it tectonics, GIA, elastic loading, viscoelastic flow, poroelastic rebound, etc. If the volume encapsulated by the entire Earth's surface is constant over time, then  F = 0 when ∂S covers all of Earth. Several of the most important processes driving VLM are conserving in this way. For example, mantle flow and surface loading are  F = 0 since the loads on Earth surface move from place to place and never leave Earth. Integrating the ICE-6G_C GIA model (Argus et al., 2014;Peltier et al., 2015) results in total  F ∼ 0. This means that while GIA currently drives uplift in some places, elsewhere it is matched by subsidence to keep the Earth volume constant. Similarly, hydrological loading changes are net  F = 0, since water moving of the surface changes a load location, but the net deflection is the same. The same is true for some other processes but not all. For example, water that is, removed from an aquifer will cause a change in loading location, but also causes a separate contribution to subsidence owing to compaction and/or poroelasticity. The poroelastic effect is not volume conserving because the water does not necessarily return to another aquifer, resulting in net subsidence. Thus, the total Earth surface is not strictly volume conserving but is approximately so since the largest mechanisms are GIA and hydrological loading, and others have relatively minor effects.
When we perform the integration in Equation 10 over Earth's land areas, we find that the net vertical volume flux of the observed land areas is 22.1 km 3 year −1 ( Table 3). The land areas are in net uplift, implying that the unobserved areas (e.g., ocean basins and masked areas under existing ice sheets) experience net subsidence. Averaged over the observed land area the net uplift is 0.16 mm year −1 , while the subsidence in the unobserved areas needed to compensate is −0.06 mm year −1 , a slower absolute rate because the unobserved areas have a larger area.
But how much of this is attributable to non-GIA processes? We can begin to address this by correcting the VLM field using a GIA model. The GIA model we consider (ICE-6G_C- Argus et al., 2014;Peltier et al., 2015) has a net upward vertical flux in the imaged land areas of 87.5 km 3 year −1 and downward flux of −13.7 km 3 year −1 for a net upward flux of 73.8 km 3 year −1 ( Table 3). The net flux of the GIA model is far greater than the net Figure 13. Map indicating by color whether vertical land motion estimate at tide gauge location is good (green), medium (yellow), or poor (red) according to criteria discussed in text.
22.1 km 3 year −1 found in the imaging of the same land area, suggesting that non-GIA processes drive net subsidence on the continents. Subtracting the predictions of the GIA model from the VLM field ( Figure 9) removes most of the uplift and fore-bulge subsidence in North America and Europe and changes signals in Antarctica ( Figure 14). The residual volume flux of −51.7 km 3 year −1 is the strongly negative (downward) net signal that non-GIA vertical motion has on the land areas. That these trends are largely caused by elastic response of the Earth to changes in the location of hydrological loading is evident when comparing the VLM trends to trends in GRACE data. Figure 14 shows the trends from GRACE, with the scale flipped so that drying is portrayed as red so it may be compared to locations where GPS-measured motion is upward, and blue can be compared to locations where GPS-measured motion is downward, consistent with loading. While differences in GPS and GRACE resolution exist, in many locations there are similarities between GRACE and GPS signals which are consistent with the effects of elastic hydrological loading, just as they were similar in their seasonality ( Figure 5). In many areas, the GPS and GRACE trends do not agree, suggesting processes other than GIA or hydrological loading (e.g., tectonics, mantle flow, aquifer dynamics) drive trends in VLM. Also, disagreements between GPS and GRACE seem to be in the areas that are not well covered by GPS. Clearly, having a global model for time-variable hydrological elastic loading to correct the GPS and GRACE data is desirable, so we could evaluate these trends. Generating or evaluating such models is beyond the scope of this study, but progress is being made by several services that predict loading displacements globally using models of the atmosphere, oceans, and terrestrial hydrology, for example, those of the EOST/IPGS Loading Service (Mémin et al., 2020) and GFZ Earth System Modeling group (Dill & Dobslaw, 2013).
If non-GIA-related VLM observed on land is on balance negative, then this implies that much of the unobserved areas rise owing to non-GIA sources to maintain a near net-zero Earth volume (subject to the caveats mentioned above). This unobserved area includes the oceans and the currently ice-covered areas of Greenland and Antarctica, which over the period of observation with GPS have been losing mass at high rates owing to warming climate conditions with corresponding elastic unloading response that is not accounted for in the GIA model. Thus, it is likely that much of the uplift needed to balance the observed subsidence is under these HAMMOND ET AL.   ice sheets. It could also possible that the sea floor experiences net uplift owing to non-GIA-related processes. However, groundwater systems around the Earth are experiencing trends of depletion (e.g., Rodell et al., 2018) so a portion of the subsidence may be associated with aquifer compaction. Quantifying the proportion attributable to aquifer compaction will require more detailed study that includes comparisons with other data sets.
Regardless of the details in these arguments regarding budgets, the net volume flux of non-GIA VLM of −51.7 km 3 year −1 is nearly as great in absolute value as the GIA contribution of 73.8 km 3 year −1 . We conclude that while non-GIA VLM does not have as large an overall impact as GIA, it is still a globally substantial component of the total VLM, and it should be considered in analyses of RSL. This applies when correcting tide gauge records or closing the budgets between other techniques that measure trends in sea surface height such as satellite altimetry.

Conclusions
We use a database of vertical component position time series from globally distributed GPS stations to obtain a geographically balanced estimate of global VLM. The signals show a diversity of processes driving uplift and subsidence, including GIA, tectonic and seismic cycle deformation, seasonal and multiannual variations in hydrological loading, and many others. The imaged field is the observed VLM from all process-HAMMOND ET AL.  Figure 9 except with predictions of glacial isostatic adjustment (GIA) model (Argus et al., 2014;Peltier et al., 2015) subtracted. (b) Trends in equivalent water height from Gravity Recovery and Climate Experiment (GRACE) with the scale flipped so that water loss is reckoned positive. This results in a color scale that is red when GPS-measured motion is upward and GRACE indicates drying and is blue for GPS motion downward and increased GRACE-measured water load. Note that the GPS in is units of mm year −1 and GRACE in cm year −1 . This GRACE product is provided by Loomis et al. (2019) with the contribution from GIA removed. es combined. Analysis of the net vertical flux from GPS Imaging averaged over the continents indicates that the non-GIA sources have about 70% of the impact of GIA's contribution to VLM.
Characteristics of the GPS signals help partition the sources of uncertainties in VLM, which are both spatial and temporal and come from aleatory uncertainty and signal variability. Considering all GPS stations in the database allows us to quantify these factors and summarize them for each tide gauge, providing the ability to identify sources of VLM uncertainty at each gauge. These results can be used to improve projections of present and future RSL rise from tide gauge and satellite altimetry data, which constrain RSL.
Comparison between GPS and GRACE data shows that most of the trends, seasonal changes, and multiannual variability in VLM come from interaction of Earth's surface with its fluid envelopes. Hydrological loading affects nearly all of Earth's surface.
VLM of the continents is on average upward, implying that the unobserved areas (composed of the ocean basins and areas currently under the large ice sheets) move on average downward with respect to Earth center. However, after correcting for GIA the reverse is true, and continents subside on average suggesting that the ocean basins and ice-sheet-covered areas undergo net non-GIA-related uplift. A large portion of this is likely attributable to elastic unloading from contemporary ice loss that is not accounted for in GIA models.