Present‐Day Vertical Land Motion of Australia From GPS Observations and Geophysical Models

The secular rate of Australia's vertical surface deformation due to past ice‐ocean loading changes is not consistent with present vertical velocities observed by a previously sparse network of Global Positioning System (GPS) sites. Current understanding of the Earth's rheology suggests that the expected vertical motion of the crust should be close to zero given that Australia is located in the far field of past ice sheet loading. Recent GPS measurements suggest that the vertical motion of the Australian continent at permanent sites is between 0 and −2 mm/year. Here we investigate if vertical deformation due to previous ice sheet loading can be recovered in the time series of Australian GPS sites through enlarging the number of sites compared to previous studies from ~20 to more than 100 and through the application of improved data filtering. We apply forward geophysical models of elastic surface displacement induced by atmospheric, hydrologic, nontidal ocean, and ice loading and use independent component analysis as a spatiotemporal filter that includes multivariate regression to consider temporally correlated noise in GPS. Using this approach, the common mode error is identified, and subsequent multivariate regression leads to an average reduction in trend uncertainty of ~35%. The average vertical subsidence of the Australian continent is substantially different to vertical motion predicted by glacial isostatic adjustment and surface mass transport models.


Introduction
Accurate knowledge of vertical land motion is key to understanding sea level change and variability in both regional and global contexts (Wöppelmann & Marcos, 2016). To be able to determine regional and global patterns of absolute sea level change, corrections for vertical land motion need to be made to relative sealevel estimates from tide gauges (Bouin & Wöppelmann, 2010;King et al., 2012;Pfeffer & Allemand, 2016;Wöppelmann & Marcos, 2016). Vertical land motion originates from both geophysical and anthropogenic sources and is an important factor to consider when inferring global changes from local observations Frederikse et al., 2018) and for predicting future sea-level changes (Kopp et al., 2014). In the current Intergovernmental Panel for Climate Change Fifth Assessment Report , only vertical land motion due to glacial isostatic adjustment (GIA), the solid Earth's response to ice-ocean loading changes (Mitrovica & Milne, 2002;Tamisiea & Mitrovica, 2011), is accounted for in the regional and coastal mean sea level rise projections. Neglecting non-GIA-related sources of vertical land motion is likely to bias interpretation of projections in some areas Han et al., 2015).
The main geophysical process driving spatially coherent vertical velocities above 0.1 mm/year globally is GIA. Current geodetic estimates of Australian vertical land motion do not agree with the motion predicted by earlier models of GIA, such as ICE-5G (VM2) (Peltier, 2004), with a discrepancy at the level of a few millimeters per year (e.g., Ostanciaux et al., 2012). Across Australia, recent model predictions of GIA are ambiguous as to the sign of vertical land motion (Figure 1), although predictions have values of <±0.5 mm/year over the Australian landmass (Table 1). The differences between GIA models are partly due to their process of construction and theoretical differences in their computation (Whitehouse, 2018). Recent sparsely distributed Global Positioning System (GPS) measurements suggest that the vertical motion of the Australian continent at a dozen well-distributed permanent sites is between 0 and −1 mm/year (Altamimi et al., 2016;Burgette et al., 2013;King et al., 2012;Santamaría-Gómez et al., 2012). Here we investigate vertical land motion at more than 100 sites across the Australian continent.
The present-day vertical velocity field from GPS demonstrates spatial coherency across the Australian continent. If the GPS velocities are correct, this suggests the nonnegligible influence of other large-scale geophysical phenomena, for example, far-field deformation due to present-day ice unloading and postseismic deformation. If the derived quantities of land motion are biased, they will in turn bias GPS-corrected tide gauge estimates of sea level rise or adversely affect other geophysical interpretation (e.g., surface mass transport, ice loading, postseismic deformation, thermal expansion, groundwater, and oil pumping). If biases in vertical land motion estimates are coherent over large spatial scales and sustained over significant time periods, reference frame accuracy is also likely to be affected. Given that GIA models and vertical land motion observations are used to correct tide gauge measurements of relative sea level (White et al., 2014), resolving this discrepancy is important for both regional and global sea level change studies as well as for improving our understanding of crustal deformation and dynamics. The reduction of uncertainty achieved here allows for a clearer understanding of spatial and temporal solid Earth deformation patterns in Australia.
Systematic error or drift in the GPS reference frame could be an important and widespread source of velocity bias (Argus, 2012;Wu et al., 2011). The Earth's center of mass (CM; geocenter) is not stationary, and its secular and seasonal motion is driven primarily by changes in mass on the Earth's surface and within the interior of the Earth's crust and mantle. GPS station height time series are inherently sensitive to the stability and accuracy of a global reference frame, especially the variations of the origin with respect to the Earth's surface (Serpelloni et al., 2013). The secular origin of the International Terrestrial Reference Frame (ITRF) approximates the CM of the total Earth system, including the fluid envelope. However, various ITRFs are not perfect realizations of CM even at secular timescales. Argus (2012) found the uncertainty of the ITRF2008 origin translation parameters to be ±0.4 mm/year for the X and Y components and ±0.9 mm/year for the Z component (95% confidence limit). ITRF2014 shows a modest reduction in uncertainty compared with ITRF2008 where the uncertainty of the X and Y components was found to be ±0.3 mm/year, and a reduction in the Z component uncertainty was found to be ±0.7 mm/year (95% confidence limit) (Riddell et al., 2017).
Velocities from GPS in an ITRF and from forward GIA models are not strictly directly comparable due to differences in frame origin. Forward GIA models are computed in the CM of the solid Earth (CE) frame, which is not accessible by observation (Blewitt, 2003). Over seasonal and shorter time scales the ITRF reflects a center of figure (CF) frame (Blewitt, 2003;Collilieux et al., 2009;Dong et al., 2003). Jiang et al. (2013) found that the majority (90 %) of root-mean-square variations between loading models computed on a grid in CE and CF frames were below 1%. There are multiple reasons for differences in the uncertainty of the velocity of the center of the Earth (i.e., the velocity between CM and CE). These include the implicit sources of error in satellite laser ranging (SLR), which is the only geodetic technique currently used to define the Earth's CM. The SLR Z translation rate uncertainty is ±0.7 mm/year with a 95% confidence interval (Riddell et al., 2017). By minimizing the differences in geodetic observations of vertical land motion rates from GPS and radial land motion rate predictions from a postglacial rebound model, the velocity of CE can be estimated (Argus, 2007;Argus et al., 2010). There is ambiguity as to the magnitude of the absolute motion of the geocenter. Argus et al. (2012)   ±0.2 mm/year in all three translational components (X, Y, and Z) over the period 2005.0-2015.0. We also note that other estimates based on analysis of Gravity Recovery and Climate Experiment (data over the period January 2003 to August 2016 (Sun et al., 2019) suggest a slightly higher value of ±0.5 mm/year and that this is an area of ongoing investigation. A recent study by Métivier et al. (2020) has obtained a geocenter velocity reaching 0.9 ± 0.5 mm/year in 2013 with a Z component of 0.8 ± 0.4 mm/year from spherical harmonic coefficients estimated from ITRF2014 Global Navigation Satellite System (GNSS) velocities, over the period 1994.0-2015.1, which is larger than previous estimates and attributed to accelerations of recent global ice melting as a purely elastic response. We consider this further below and note that accurately determining small rates of change from GPS coordinate time series is complicated by reference frame errors and inconsistencies.
Other geophysical phenomena and measurement noise introduces short-term variability in the GPS time series, which can bias velocity estimates and/or increase velocity uncertainties. Spatially correlated noise is the largest source of error in a regional network, and this is also commonly referred to as "common mode error" (e.g., Dong et al., 2006;Liu et al., 2018). To our knowledge, common mode error across a spatially dense field of Australian GPS sites has not previously been investigated. Common mode error within a network of GPS stations can be associated with nontectonically induced deformation from uncorrected or mismodeled surface mass transport, including atmospheric, hydrologic and nontidal loading; tectonic-related crustal deformations; reference frame realization bias and mismodeled effects of polar motion; satellite errors related to clock, orbit and antenna phase center variations; and systematic errors caused by software and data processing strategies.
In this paper, we seek to identify whether GPS estimates of GIA in Australia are robust and if using decomposition techniques allows the identification or isolation of GIA as a common cause of vertical land motion. To do this, we first remove the short-term variability from GPS time series at Australian sites to avoid biasing trend estimates computed from the vertical component. We then compare our observed velocities of vertical land motion to predicted GIA radial rates. By attempting to obtain present-day velocities reflective of long-term processes, we aim to elucidate a signal with geophysical origin above the inherent noise. In the context of this study, "noise" refers to those quantities that are not of geophysical origin. The methodological advances in this research include the use of multivariate regression with power law noise, a unified approach to vertical velocity and uncertainty estimation at Australian sites (including common mode error), and the quantification of the discrepancy between observed and modeled vertical rates of solid Earth motion. Figure 2 shows the processing flow and structure of this paper. Section 2 presents the data used for the comparison of the GPS vertical velocities with forward surface mass transport models. The models are then subtracted from the GPS to remove short-term variability from the time series. Section 3 describes the use of spatiotemporal filtering methods to further filter common mode error. Finally, we compare the final GPS vertical velocities and GIA, followed by a discussion of the data, models, assumptions, and conclusions.

GPS Data and Geophysical Models
2.1. Data 2.1.1. GPS GPS data from sites in the Australian regional GNSS and AuScope networks (Ruddick & Dawson, 1990), listed in supporting information Table S1, were processed as daily precise point position solutions using GIPSY v6.3 combined with the Jet Propulsion Laboratory clock and orbit products (Bertiger et al., 2010;Zumberge et al., 1997). A more thorough description of the parameters, constraints, and models are detailed in the supporting information (Text S1). Times of offsets in the resultant time series were identified by considering site log file entries, supplemented by visual identification (Gazeaux et al., 2013). The daily nonfiducial position time series were transformed into ITRF2014 (Altamimi et al., 2016).
Some sites in the Australian network have been operational for more than 22 years, with time series beginning in 1995, but the majority were installed in the period 2010-2013. To maximize the number of sites with a common-data period, we only used data from 2013 onward, giving a total of 113 sites. This presents a significant increase in the number of sites used in previous studies, where four sites were used by King et al. (2012) to investigate vertical land motion at tide gauge stations; 12 sites were used in Burgette et al. (2013), with a focus on GPS sites near Australian tide gauges; 13 sites were used by Santamaría-Gómez et al. (2012), again with a focus on GPS sites located near the coast; and 36 Australian sites were used in Altamimi et al. (2016) for the realization of ITRF2014.
The data used in this analysis span from 2013.0 to 2019.0 (6 years). Data spans of 3 years or greater are sufficient to determine a velocity with realistic uncertainty. Bos et al. (2010) took into consideration the effect of colored noise and found that 4.5 years of daily GPS data is sufficient to reduce the trend uncertainty to below 1 mm/year (estimated with an annual signal). Consequently, a 6-year data span is considered sufficient to determine a robust linear trend and uncertainty in the presence of power law and white noise where the spectral index (κ) is also estimated (i.e., not fixed to κ = −1 [flicker noise] or κ = −2 [random walk noise]). To estimate velocities from the GPS data, we use the maximum likelihood estimator in Hector v1.6 software , with a power law plus white noise model, simultaneously estimating annual and semiannual signals as well as offsets with times as predetermined.
Before GPS velocities are compared to GIA, the GPS velocities need to be transformed into reference frames with compatible origins (Argus, 2012). To correct GPS vertical velocities for the CM drift, we use a range of forward models of surface mass transport as an approximation of CE-CM corrections, including those from GRACE with elastic loading due to present-day ice mass loss, ideally achieving vertical land motion values in a secular CF frame. The CF frame can be considered as closely approximating a CE frame (Dong, Dickey, Chao, & Cheng, 1997) and is theoretically realized by an infinitely dense and uniformly distributed network of station positions (and velocities) on the Earth's surface. The difference between CM and CE is predominantly due to present-day mass distributions from ice, terrestrial water storage, and the atmosphere.
It is important that all geophysical corrections made to GPS data are also consistent with the geophysical time scale of interest. Changes in the motion of the Earth's rotation pole will change geodetic vertical velocities (King & Watson, 2014). Conventionally, the "pole tide correction" (Petit & Luzum, 2010) only removes deformation associated with periods of less than a few years. Longer-term deviations in the data are neglected, whereas for comparison to GIA models, for instance, they need to be referred to a reference polar motion compatible with millennial time scales. We discuss our treatment of this bias below.

Forward Surface Mass Transport Models
We use two alternative combinations of geophysical forward models to remove the effects of present-day surface loading changes, with the aim of reconciling GPS observations of vertical land motion and motion predicted by GIA models. The first combined forward model accounts for atmospheric, hydrologic, and nontidal ocean loading. The second approach uses surface mass loading determined from the GRACE mission. We assess both to determine the optimal approach and then adopt the time series corrected using the preferred approach in subsequent analysis.
Subtracting models of surface deformation due to surface mass transport is commonly used to reduce time series variance (evidenced by a reduction in root-mean-square error) (e.g., Li et al., 2017). This technique aims to reduce uncertainty and improve accuracy, but using the root-mean-square error as an indicator of reduced variance is not particularly effective at enabling the refinement of secular velocities (Santamaría- Gómez & Mémin, 2015). With a focus on the vertical component, other studies have shown that the impact of mass surface transport caused by continental water (Fritsche et al., 2012), atmospheric pressure (Dach et al., 2011), and nontidal oceanic loading  can bias geodetic velocity estimates depending on the time series length. Reduction of the temporal variability is beneficial for isolating transient signals in the time series, and so we use this approach to start our investigation of vertical land motion at Australian GPS sites.
Atmospheric loading displacements were estimated at each site across the Australian AuScope network using the International Mass Loading Service (Petrov, 2015) Modern-Era Retrospective analysis for Research and Applications (MERRA-2) products (Gelaro et al., 2017). Terrestrial land water storage was computed using land water storage pressure from the numerical weather model MERRA2 (Reichle et al., 2017). Nontidal ocean loading, based on the spherical harmonic transformation of the ocean bottom pressure, were computed with the OMCT05 model (Dobslaw & Thomas, 2007). The atmospheric, land water storage, and nontidal ocean loading surface displacements were downloaded from http://massloading.net/ using Stokes coefficients of the harmonic expansion of the contribution of the atmosphere to the geopotential and load Love numbers derived on the basis of the Preliminary Reference Earth Model (Dziewonski & Anderson, 1981). All surface loading products are downloaded with reference to the center of solid mass of the Earth (CF), which is appropriate for comparison with GPS time series in an ITRF that contain nonlinear motion (Dong et al., 2003).

GRACE-Determined Surface Mass Loading
As an alternative approach, we considered surface displacements based on GRACE data. The redistribution of surface mass causing vertical displacement can be well predicted by the GRACE mission (Davis et al., 2004;Wang et al., 2017). Surface loading deformations at Australian GPS sites were extracted from monthly GRACE data and background models where the surface displacement solution is Release-06 from Center for Space Research (CSR; at University of Texas, Austin), up to degree 96 using degree-1 values from Swenson et al. (2008). The selection of the CSR solution is arbitrary, with differences between alternative solutions from the Jet Propulsion Laboratory and German Research Centre for Geosciences (GFZ) being small and within the error bounds of the GRACE solution itself; for further information, see Sakumura et al. (2014). The degree-2,0 values were replaced following Cheng and Ries (2017), and spatial filtering following Kusche (2007) (the DDK7 model) was applied. A GIA model is not applied in the CSR spherical harmonic solution.
The surface loading contributions were split up into components from glaciers and ice sheets, and a hydrological module, each of which is accompanied by a self-consistent ocean response to conserve mass at a global scale. The deformations are provided in the CF frame and represent the elastic response of Preliminary Reference Earth Model to the estimated surface loads. The secular linear trend was estimated over 2002.3-2017.5 and removed from the modeled surface deformation.

Glacial Isostatic Adjustment
Three GIA forward models are used for comparison to the GPS vertical velocities. Rates of radial displacement on a 0.2 × 0.2°grid were downloaded from http://www.atmosp.physics.utoronto.ca/~peltier/data. php (last accessed 7 May 2018) based on the ICE-6G_C ice history model with VM5a Earth model (Argus et al., 2014;Peltier et al., 2015). The model of Purcell et al. (2016) provides an alternative realization of this model using the same ice and Earth models but a different computational framework, with results provided on a 1°grid of GIA rates that have been modeled with alternative Stokes coefficients and is hereafter called ICE6G_C_ANU.
We supplement these two conventional forward GIA models with a more recent model determined using a Bayesian framework to generate expected GIA estimates and formal uncertainties globally (Caron et al., 2018). This model is based on 128,000 forward models of varying ice history and 1-D Earth structure constrained by moraine positions, dating history of the ice sheet collapse, tide gauge-relative sea level measurements, and vertical land motion trends derived from GPS. There are minimal data constraints in Australia. The model output is on a 0.4°grid. The uncertainties from this model over the Australian continent are relatively consistent at~0.12 mm/year.

Mass Loading Comparison
The GPS time series are "corrected" with the modeled surface deformation data to account for interannual variability and to obtain a vertical velocity that is not predominantly induced by terrestrial water storage.
10.1029/2019JB018034 Table 2 shows the median values over the Australian AuScope network where the mass transport models have been removed from the GPS time series. We found that truncating the time series to the time span 2013.0-2019.0, to maximize the number of sites with co-temporal data, does not degrade or significantly change the noise characteristics or velocity estimated from the individual site time series as shown in the supporting information (Table S2).
Colored noise is known to yield quasiperiodic patterns and variations within time series, interfering with the real periodic terms, making the detection of signals with small amplitude and power difficult to isolate (Amiri-Simkooei, 2007; Amiri-Simkooei, 2013). Using the maximum likelihood estimator in Hector software, we assess the geophysical models-(i) forward surface mass transport models and (ii) GRACE-derived mass transport-against the GPS time series to attenuate the influence of mass surface loading. The combination of power law and white noise was chosen as the representative noise model for GPS time series (Bos et al., 2007;Santamaría-Gómez et al., 2011;Williams, 2003), based on both the Akaike information criterion and the Bayesian information criterion (He et al., 2019).
When the respective surface mass transport loading models are removed from the GPS series (GPS minus SMT), the median unfiltered GPS trend shows a reduction of the trend by −0.6 mm/year (with a median uncertainty of 0.5 mm/year) ( Table 2). Removal of surface mass transport does not substantially change trend uncertainties, with the exception of the instance when surface displacement from GRACE is removed and the median uncertainty increases from 0.4 to 0.9 mm/year (1σ). The dominant negative trends across the network suggest subsidence is the dominant signal even after consideration of present-day mass loading signals.
The amplitude of power law and white noise is reduced when mass loading series are removed from the GPS time series, an example being the median reduction from 8.5 (with a median uncertainty of 0.2) to 7.0 (0.2) mm/year −κ/4 , when the combined surface mass transport loading displacement product is removed from the GPS series. The largest change in spectral index is an increase from a median value of −0.6 (0.1) for the GPS series to −0.9 (0.1) when the GRACE series are removed from the GPS time series, representing a coloring of the noise. It is also interesting to note that regressing the atmospheric loading displacement increases the spectral index to −0.9 (0.1).
After following the standard approach of removing mass loading from the GPS time series, we note that there is remaining residual signal with some power at annual and semiannual frequencies ( Figure 3). We test an idea that the loading models are mis-scaled, due to, for example, GRACE or forward model resolution issues, meaning the mass change is not sufficiently localized. Consequently, as a variation on the standard removal approach, we choose to regress the mass loading time series against the unfiltered GPS time series within Hector while estimating the other parameters and noise model as usual. The multivariate regression   considers the relationship between the surface mass transport loading and the GPS time series in the presence of noise other than white (as would be the case for simply removing the surface mass transport series from the GPS time series).
Regression of each of the forward mass loading models against the GPS time series results in a larger reduction of the trend uncertainty than that by simply removing the forward model from the GPS time series (  (Table 4), but the level of noise remains too high to identify if GIA is a contributing factor to the distinct pattern of subsidence. The velocity of vertical motion after the respective surface mass transport forward models have been regressed against the GPS height series provides an improvement in the match between GPS and GIA, but the GPS velocities are biased toward the GIA models. This comparison does not account for differences in the reference frame between GPS and GIA. Significant residual noise remains in the time series, which may mask the geophysical signal of interest (Figure 3). To address this potential effect, in the next section, a common mode filtering approach is investigated.

Residual Analysis
Using a spatiotemporal filter, we now assess the residuals remaining after including the combined surface mass transport loading (SMT REG from Table 2), as a dependent variable in the regression, with regression parameters estimated in the presence of power law and white noise. Spatiotemporal filtering is an effective method of assessing the relative error contribution of various geodetic data sets and testing if the scaling of the variates is appropriate. We seek to identify and reduce residual common mode error in the Australian AuScope network.

Spatiotemporal Filtering
Using the surface mass transport loading deformation at each site within the regression provides a reduction of the power law and white noise amplitude and trend uncertainty, but substantial noise remains in the solution (Figure 3). We therefore attempt to further reduce the noise through empirical filtering. In this section we investigate spatially correlated common mode error using principal component analysis (PCA) and independent component analysis (ICA) methods for decomposing GPS time series into temporal components and their associated spatial responses. PCA methods are well known and are described, for instance, in He et al. (2015).
Improving the signal-to-noise ratio of geodetic coordinate time series by separating and removing the source signals that contribute to the reduction of precision is an effective way of identifying small or weak deformation signals. Two commonly used matrix factorization methods for reducing dimensionality in data are PCA and ICA (Hyvärinen & Oja, 2000). Component analysis approaches have been used to identify and remove common mode error in GPS networks for both the spatial and temporal domains (e.g., Dong et al., 2006;Gruszczynski et al., 2016;He et al., 2015;Li et al., 2015;Liu et al., 2015;Liu et al., 2018;Ming et   By applying spatiotemporal filtering using the PCA and Karhunen-Loeve expansion algorithms, the primary objective of Dong et al. (2006) was to identify local site effects within a network of stations. Liu et al. (2015Liu et al. ( , 2017 focused on investigating if ICA could be used to identify signals with geophysical meaning related to surface mass loading. For ICA computations, we use the FastICA fixed-point optimization algorithm (Hyvärinen, 1999) to maximize computational efficiency with the pow3 approach, where the derivative function g = G′ is termed the nonlinearity. Using the classical kurtosis optimization approach gives the nonlinearity g(u) = u 3 (pow3). The implementation of ICA follows the description in Hyvärinen et al. (1999) and Comon (1994): x ¼ As where x is an observed m-dimensional vector, s is an x-dimensional random vector whose components are mutually independent (the sources), and A is a constant m × n matrix to be estimated (mixing matrix). The original sources (s) can be recovered by multiplying the observed signals (x) with the inverse of the mixing matrix (W = A −1 ), also known as the unmixing matrix.
The separation of independent components from multidimensional data is the objective of ICA, but a wellknown issue is that because both the separation matrix and independent components are estimated simultaneously, neither the amplitudes nor the ordering of the components can be established. Hendrikse and Spreeuwers (2007) consider ordering the components based on the power of each component as a contribution to the independent components. Strong sources are those that provide a larger contribution, and weaker sources are those that present a smaller contribution to the overall signal. Selection of the components with the highest power mimics the ordering of eigenvectors in PCA, which provides the best selection of independent components in the realm of least squares. We apply this pseudo-ordering of ICA components to identify the largest contributors. Gaps in the time series of individual sites were initially infilled with a spatial mean and then updated with reconstructed series from the six leading PCA modes.
Each GPS position component (north, east, and up) is assessed individually, considering that spatial cross correlations between components are not likely to be significant (Amiri-Simkooei, 2013). Cross correlations between the position components of the Australian AuScope GPS network were confirmed as insignificant and only the results of the vertical time series are presented here.

Results
We find that the leading six modes of PCA are representative of the common mode error in the GPS time series, noting the common spatial pattern and proportion of variance explained. For ICA the appropriate number of principal components used needs to be resolved. Here the number of principal components used in ICA was determined using a trial-and-error approach (Barnie & Oppenheimer, 2015). Careful consideration needs to be given using this approach as too many principal components retained can result in overfitting sources as isolated   spikes, and too few can result in discarding useful information .
The temporal and spatial representations of the leading modes are shown in Figure 4, following principal and independent analysis of the AuScope GPS vertical time series. We normalize each time series to allow for the comparison of PCA and ICA time series. Each spatial response is normalized by the maximum absolute value to enable comparison and are scaled to lie between −100% and +100%.
For large networks, common mode error is generally represented by smooth transitional patterning of the spatial response (Serpelloni et al., 2013). The common sign of the spatial eigenvector in Figure 4 provides a strong indication that the modes from each analysis demonstrate a coherent and common spatial pattern that can be interpreted as a signal above the noise floor. The temporal series of the leading principal component (PC1) and the first independent component (IC1) show a quasiperiodic signal with similar amplitude that is prevalent in the series even after removing modeled deformation due to atmospheric, land water storage, and nontidal ocean surface mass transport. This suggests short-term quasiperiodic signals that remain in the GPS height series are spatially coherent over length scales more than 4,000 km; whether these are due to surface loading model deficiencies or other effects requires further investigation.
We removed the leading modes from the GPS time series and examined its effect on temporal and spatial variability. We found that this approach did not substantially improve the time series and note that this approach of simply subtracting common modes from time series neglects the temporal correlation of noise and uncertainty propagation. To improve on this, we took the approach of regressing the principal modes against the original GPS time series while also coestimating other parameters and noise models as above.
The leading modes from both PCA and ICA were detrended and then regressed individually and collectively (PC1-2, PC1-3, PC1-4, etc.) to assess the change in noise characteristics and to determine which component analysis method is more effective at reducing noise and trend uncertainty. The regression coefficients are generally less than 1 with a median value of 0.8 for PCA components and 0.7 for the ICA components (Table  S1). The modes can also be used to make downstream determinations of individual site anomalies and regional deformation patterns.
We compared the estimated noise characteristics for the time series after removing regressed leading modes, PC1 and IC1, against those for the unfiltered GPS time series. The results in Table 3 show the larger reduction of power law noise amplitude, spectral index, and trend uncertainty for the PC1 regression, where the power law noise amplitude of the unfiltered GPS series is changed from 8.5 (0.2) to 8.7 (0.2) and 7.8 (0.2) mm/year −κ/4 after the leading modes from ICA and PCA have been removed and regressed, respectively. As more of the leading modes are included in the regression, (e.g., top four components, 1:4), the difference in power law noise amplitude and trend between PCA and ICA is reduced from a 10% difference to a 5% difference.
A general concern with PCA is the orthogonality constraint imposed by the statistical definition of the method. By defining that the primary and secondary axes of the maximum variation must be orthogonal, a constraint is placed on analysis that can result in overinterpretation or underinterpretation of real spatial and temporal variability. We address this by providing results from ICA as a comparison and noting that they are similar. The combination of the top four independent components (Table 3) demonstrated a larger relative reduction of power law amplitude, reducing the median power law and white noise from 8.7 (0.2) mm/year −κ/4 of the first independent component to 7.2 (0.2) mm/year −κ/4 for the leading four independent components combined.
The regression coefficient from the multivariate regression analysis demonstrates a similar scaled magnitude and spatial pattern as the spatial eigenvectors from each PCA and ICA, providing confidence in the regression technique. A check of the independent component ordering also showed that the regression of the first independent component against the individual GPS time series showed the largest reduction in power law noise amplitude, the smallest trend uncertainty, and the lowest spectral index when compared to the other independent components, providing validation that the pseudo-ordering of the independent components is reasonable.
The median value of the GPS vertical velocities with the common mode error removed (the top four modes from ICA) is −0.5 mm/year with a median uncertainty of 0.4 mm/year compared with the same statistic computed from the unfiltered velocity field of −0.8 (0.4) mm/year (Table 3). ICA and PCA provided similar velocities when the residual time series were regressed against the original GPS time series over the period 2013.0-2019.0, and so we compare only the GPS velocities computed after removal of regressed ICA to GIA. By reducing common mode error in GPS time series, we can better isolate the remaining signal with the assumption that this represents underlying vertical land motion and compare the velocities with those predicted from GIA models.

GIA Velocity Comparison
The median uplift velocities from the three GIA models are shown in Table 4, after interpolation of the GIA grids to the GPS site locations. The GPS and GIA models all agree within 1σ, but this assumes that the uncertainties are robust. In fact, the median residual trend after filtering remains as subsidence within the network. The mean trend of GIA radial displacement from the ICE6G_C model, interpolated at each of the GPS sites is −0.12 mm/year (uncertainty is not provided). The remaining signal cannot clearly be attributed to GIA and could be the product of other large-scale geophysical phenomena, reference frame errors or systematic errors within the GPS series (such as effects related to aliasing at draconitic periods).
To mitigate the effects of reference frame errors, we recompare the GPS and model outputs by considering the spatial covariance of the estimated trends. The empirical semivariogram for the GPS ( Figure S1) suggests the correlation length to be around 500 km, while this is larger than 3,000 km for the GIA models ( Figures S1  and S2).

Discussion
The contribution of GIA to vertical land motion velocities in Australia is expected to be relatively small, considering that Australia lies in the far field of the past major ice sheets (i.e., North America, northern Europe, and Greenland). However, as the different contributions of surface mass loadings improve, the discrepancy 10.1029/2019JB018034 Journal of Geophysical Research: Solid Earth between GPS and GIA velocities is worthy of investigation, noting the importance of both in sea level studies. GIA models are commonly used as a correction in both altimetry and tide gauge studies (Tamisiea, 2011), and so the validity of the modeled deformation is paramount for reliable sea level estimates. In Australia, previous studies have noted that GIA is close to zero, whereas the GPS observations show a general rate of subsidence that is currently unexplained.
We have demonstrated that the noise level in residual time series from Australian GPS sites, combined with reference frame uncertainty remains too high to confidently separate competing GIA models even after cleaning and filtering. Mismatches between the decomposed and filtered GPS observations and GIA models could result from other sources of uncertainty in GPS such as draconitic aliasing with seasonal signals (Amiri-Simkooei, 2013;Bogusz & Klos, 2015;Ray et al., 2007); other unconsidered geophysical processes such as thermal bedrock expansion (Wang et al., 2018;Yan et al., 2009); and site-specific effects such as multipath or antenna phase center corrections (King & Watson, 2010) and reference frame alignment (Argus, 2012;Riddell et al., 2017;Santamaría-Gómez et al., 2017).
Due to the shortness of the time series, we are not able to separate the solar seasonal and draconitic signals in a harmonic analysis. We expect that the contribution from draconitic signals would be small. Using an extended set of global sites, yet a comparable GPS analysis strategy to our approach, Amiri-Simkooei et al. (2017) find that the mean amplitude of the first eight draconitic signals have an mean amplitude of less than 1 mm for the horizontal components and 2 mm for the vertical component. Colored noise is known to alias with periodic patterns and variations (Amiri-Simkooei, 2007;Amiri-Simkooei, 2013), mixing with the real periodic terms in time series, making the detection of signals with small amplitude and power difficult to isolate.
We note that there is a limitation associated with the length of the time series that is related to the date of when the majority of the Australian sites were installed. Unfortunately, while there are a small number of sites that have operated for an extended period (eight sites since the mid-1990s), the majority (more than 80%) of the sites were installed after 2012. Klos et al. (2017) demonstrated using a global network of 174 GPS sites processed with GIPSY that 13 years of data was able to attain a trend uncertainty of less than 0.1 mm/year for the horizontal components and less than 0.15 mm/year for the vertical components (simultaneously estimating time-variable seasonal signals). Assuming the noise model is time constant and representative, using a simulation of temporally correlated noise in GPS time series with amplitudes ranging from 1 to 25 mm/year κ/4 , we find it would take >10 years for the trend uncertainty to reach~0.1 mm/year at an individual site. Using the median noise amplitude from our data of 8.5 mm/year κ/4 , we find approximately 15 years is required to reach ±0.1 mm/year. Initial investigations were performed with fewer sites over longer periods of time, but the spatial resolution of the principal and independent components (common mode error) was compromised by the reduced number of sites and length of separation between the available sites. Here we aim to maximize the number of sites while still robustly estimating the trend uncertainty, and so we use the sites with time series available from 2013.0 to 2019.0.
GPS reference frame bias will likely contribute to the mismatch of median values between the GPS and GIA models since an ITRF does not necessarily accurately represent CM. Past ice loss and present-day ice mass loss are not likely to provide a significant difference in the CE and CM (Argus, 2012), with the assumption that plate interiors are not drifting with respect to CE. However, a drift can exist between the reference frame origin and the true Earth CM, which manifests as a bias between measured and modeled vertical velocities (van der Wal et al., 2015). The reference frame uncertainty over Australia is likely larger than the GPS uncertainties, where Argus (2012) estimated the CE velocity from GPS vertical observations to be ±0.4 mm/year for the X and Y components and ±0.8 mm/year at the upper bound for the Z component for ITRF2008. By comparing estimates from very long baseline interferometry and SLR, Wu et al. (2011) found that the uncertainty of the scale drift of ITRF2008 was ±0.16 mm/year, which may also contribute to the observed differences between GPS and GIA. Using the results from Riddell et al. (2017), we calculate the expression of ITRF2014 geocenter drift in the vertical direction at Alice Springs (central Australia) to be~0.5 mm/year. GIA models are computed in a CE system, whereas GPS (referenced to an ITRF) is representative of CM at temporal scales longer than a few seasonal cycles. Although the differences between the reference frames are considered to be small, they may contribute to the discrepancy between GPS trends and GIA velocities.

Journal of Geophysical Research: Solid Earth
Postglacial rebound models have an implicit assumption that the plate interiors are moving vertically relative to the CE (Argus, 2012). Due to the small ongoing differences in ocean water distribution due to GIA, the difference between GIA predictions in CM and CE is minimal (Argus, 2007;Klemann & Martinec, 2011). Station velocities from GPS are reference frame dependent, being affected by both rotation and translation parameters. Uncertainty in the reference frame origin translation parameters will likely contribute to GPS error and are larger than the bias associated with different CE and CM frames.
Variations in surface mass transport loading due to hydrologic forcing and terrestrial water storage is not well modeled primarily due to the absence of a global monitoring network and a lack of understanding of the temporal and spatial unpredictability (Lettenmaier & Famiglietti, 2006). Terrestrial water storage load is known to deform the Earth's surface and influence the annual signal of GPS time series, particularly the vertical component. Although GRACE is able to provide a global picture of mass redistribution, it is not able to recover variations in the total mass of terrestrial water storage (degree-0 coefficients) (Meyrath et al., 2017). Hydrologic models are also limited by the omission of polar ice sheet contributions and poor representation of groundwater storage. Groundwater has not been considered here despite knowing that extraction and recharge signals can significantly contribute to vertical land motion signals, for example, in the Perth Basin (Featherstone et al., 2015).
Interannual variability of the terrestrial water storage deformation signal in Australia is largely related to the El Niño Southern Oscillation and the Indian Ocean Dipole (Forootan et al., 2012;García-García et al., 2011). This signal was observed as a spatial pattern in the second principal component and independent component ( Figure S2), where the spatial eigenvectors of sites toward northern Australia focused around the Gulf of Carpentaria are flipped in sign, compared to other sites. Forootan et al. (2012) were able to demonstrate that climate teleconnections, such as the El Niño Southern Oscillation, and independent components of GRACE in northern Australia were significantly correlated. In Australia, the apparent vertical land motion due to Neogene tilting as a north-down, SSW-up differential reflects contributions from both dynamic topography and the geoid evidenced by paleo-shorelines (Sandiford, 2007).
Each of the GIA models will also contain levels of uncertainty due to discrepancies in approach and assumptions made during computation. Uncertainties from the definition of the spatial and temporal variations in the ice load and uncertainties in the viscoelastic rheology of the solid Earth will propagate into the GIA modeled radial rate (Martín-Español et al., 2016). Unlike the model from Caron et al. (2018), the ICE-6G_C (VM5a) and ICE6G_C_ANU forward GIA models are not provided with uncertainty estimates, which makes it difficult to determine if the radial velocities in Australia are considered statistically significant from that model. Assumptions made about the viscosity of the mantle will propagate as differences between the respective GIA models and each models' respective difference from the median GPS velocity. Uncertainties for GIA models are difficult to compute given the unknown unknowns, and while Caron et al. (2018) are able to provide uncertainty estimates, it is unclear whether these are realistic and representative of the fit between GPS data and the GIA model in regions such as Australia where minimal GPS data is used to constrain the model.
Using a median value as the primary metric of comparison between rates of land motion presents some limitations over an irregularly distributed and sparse network. This assumes that there is minimal spatial autocorrelation that we know is not the case over the Australian continent but is sufficient here as we are conducting a first-order comparison. A more complex comparison accounting for spatial autocorrelation is planned for future research.
The improved techniques (multivariate regression to account for variation and trend estimation in the presence of colored noise) yield improved geophysical interpretation given a reduction in uncertainty in longterm rates of vertical land motion. Our work provides an improved quantification of these land motions, highlighting the finding that the observed land motion cannot be explained by GIA or other elastic modeling over shorter time scales (we contend this is an important incremental advance in understanding land motion in Australia), yet also provides insight in the techniques also applicable elsewhere.
Improved trend estimates of vertical land motion directly improves our understanding and interpretation of the velocity field dynamics in Australia. We were able to achieve a reduction of the variability in the Australia GPS series, providing a refined estimate of trend and uncertainty following spatiotemporal 10.1029/2019JB018034 Journal of Geophysical Research: Solid Earth filtering. A comparison of the traditional practice of removing the common mode error with the novel approach taken in this paper by regressing the common mode error demonstrates that including the leading modes in the regression analysis as dependent variables further reduces the noise. Providing realistic uncertainty estimates of height changes is paramount for downstream geophysical interpretation and allows for an improvement of the signal-to-noise ratio to isolate small deformation signals.
The installation date and location of the Australian AuScope sites has an influence on the application of spatial filtering methods. The AuScope GNSS network was designed and installed with the intent of 200-km nominal spacing between stations (Johnston et al., 2008). This has been achieved in some areas, but there remain large portions of the continent where there is noticeably sparse coverage (i.e., Western Australia). Amiri-Simkooei (2013) demonstrated that spatial correlations of individual components (north, east, and up) between stations are significant over angular distances of 30°(~3,000 km) and that maximum spatial correlations are found between the nearest sites, suggesting that noise in GPS time series has a physical basis. Gruszczynski et al. (2016) looked for the common spatial response of 87 GPS sites in central Europe using PCA where maximum site separation was less than 1,900 km. The maximum distance between any two sites in the AuScope network is <1,000 km (WLAL to WARA is approximately 995 km). There is presently a trade-off between using longer time series to improve the signal-to-noise ratio of the vertical time series and the spatial density of sites given the number available.
Considering that GPS is a point source measurement and only representative of the change at a single location, it would be relevant to investigate the use of other space geodetic techniques in the future that are capable of continental-wide motion, such as interferometric synthetic aperture radar.

Conclusions
This paper investigates the present-day velocity and spatial variability of GPS-derived vertical land motion across Australia and finds distinct patterns that can be partially explained by only GIA. These are compared with models of GIA in an attempt to separate competing GIA models that show either subsidence or uplift across Australia. Time series of 113 Australian GPS sites are analyzed to quantify the residual noise magnitude after removing geophysical models of surface mass transport from atmospheric, water, and nontidal ocean influences, as well as mass loss from past ice sheets and glaciers. Previous studies of Australian vertical land motion have relied on less than 20 GPS sites.
By using multivariate regression of surface mass transport in the presence of power law noise, we reduce the temporal variability of the GPS series and then test the impact of regressing the leading ICA mode against the residuals from the GPS analysis. Simply removing the leading mode from the time series and then comparing the noise properties of the resulting time series do not consider the aliasing of colored noise with periodic patterns. In our analysis, once the common mode error has been identified, the leading modes are regressed against the original GPS time series. Spatiotemporal filtering provides a methodology for reducing correlated errors between stations within a regional network, allowing for the suppression of noise. This results in better estimates of solid Earth deformation in the Australian region and evidence of common mode error.
The technique of multivariate regression is advanced by the use of colored noise to account for variation and trend estimation. We regress the determined component time series against the GPS time series as a method of detecting change in the noise properties in the presence of power law and white noise. We then compare the noise characteristics of the original time series with the multivariate-regressed time series and note a substantial reduction of power law and white noise amplitude from 8.5 (0.2) to 7.0 (0.2) mm/year −κ/4 .
The average vertical subsidence of the Australian continent is substantially different to vertical motion predicted by postglacial rebound models. This is based on the evidence that the median vertical trend of land motion across the Australian continent is −0.8 mm/year with a median uncertainty of 0.4 mm/year and the first and third quartile values of −1.2 and −0.4 mm, respectively. The background rates of motion from postglacial models across Australia range between −0.20 and 0.03 mm/year. We note that it is challenging to have a statistically significant difference given the irregular spatial distribution of sites across the continent, the non-Gaussian distribution of rates and their associated uncertainties, and the unknown uncertainty of the ICE 6G/VM5a vertical rate predictions.
We note that uncertainty in the absolute reference frame is likely at the level of 0.5 mm/year and will likely contribute to the discrepancy between GPS and GIA vertical velocities. Here we have considered that the reference frame errors would express in a spatially coherent way across the spatial scale of the Australian continent. Spatiotemporal filtering is an effective method of assessing the relative error contribution of various geodetic data sets to geophysical interpretation. This leads to a more rigorous determination of vertical land motion, which is particularly important for sea level estimation at tide gauge sites along the Australian coastline.