Comparison of state-of-the-art GNSS-observed and predicted ocean tide loading displacements across Australia

We seek to quantify and understand the residual signal in GPS and GLONASS estimates of ocean tide loading displacements (OTLDs) after removing state-of-the-art model estimates. To consider contributions over a broad spatial scale, we estimate OTLD over the Australian continent using ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}5.5 years of continuous GPS and GLONASS data from 360 sites. We compare these with modelled estimates, with a focus on the lunar semidiurnal M2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document} and diurnal O1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_1$$\end{document} constituents. We observe spatially coherent patterns of residual OTLD in each of the east, north, and up coordinate components after the removal of tidal loading using elastic models. We subsequently assess the impact of including anelastic dispersion in the model and show a 0.2 mm reduction in range of the up component residuals at coastal sites. A similar reduction at all sites is observed in the east and north components. Of the seven ocean tide models used, we find that three recent models, FES2014b, GOT4.10c and TPXO9.v1, perform similarly, noting these comparisons are made in the CE frame. However, we show that the latter contains centre of mass (CoM) biases in amplitude up to 0.2 mm and 0.5 mm for M2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document} and O1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_1$$\end{document}, respectively, due to the assimilated altimetry data having not been corrected for geocentre motion. We find OTLD estimates are sensitive to the chosen orbit and clock products used in our analysis, with differences of up to 0.5 mm in the east component between solutions using the JPL and either of ESA or CODE products (GPS-only). Our analysis shows that current GNSS estimates of OTLD over Australia are typically accurate to ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}0.2 mm at which point we are unable to explain spatially coherent residuals when compared to modelled quantities. These may depend on the appropriate treatment of the CoM variation, anelasticity and/or three-dimensional Earth structure. As such, we recommend great care is taken when interpreting OTLD at the level of 0.1 or 0.2 mm, even if it is regionally coherent.


Introduction
Data from continuous GPS sites have been used to estimate the ocean tide loading displacements (OTLDs) of the solid Earth since the early 2000s (Allinson 2004;Khan and Tscherning 2001;Schenewerk et al. 2001). These have been used to yield insights into ocean tide models (King et al. 2005;Penna et al. 2008) and solid Earth rheology, notably asthenosphere anelasticity Martens et al. 2016;Martens and Simons 2020;Wang et al. 2020;Ito and Simons 2011). The use of satellite geodesy to infer OTLD has mainly been limited to the use of the GPS constellation as it was the only complete Global Navigation Satellite System (GNSS) from the late 1990s until 2010 when the GLONASS constellation was fully restored. Since that time, many GNSS networks started migrating to multi-GNSS receivers and antennas, deploying new sites with GPS+GLONASS or multi-GNSS capable receivers (e.g. Villiger and Dach 2020, Sect. 3.1). As such, there is now often one decade of both GPS and GLONASS data available at hundreds of sites globally (e.g. Prange et al. 2017). Besides providing additional data to assist in observing millimetre or sub-millimetre level displacements (e.g. Geng et al. 2017), GLONASS has some advantages over GPS in terms of the absence or reduction of systematic errors at tidal frequencies (Abbaszadeh et al. 2020;Matviichuk et al. 2020).
Most previous studies related to OTLD have focused on regions over the scale of a few hundred km which makes it difficult to separate widespread GNSS systematic errors from more localised OTLD model error. Some of the exceptions include Bos et al. (2015) and Martens et al. (2016) though these studies used solely GPS data. We focus here on the tidal deformation of Australia, making use of GPS and GLONASS data from a continental array of multi-GNSS receivers. Assessing displacements over a large continental region such as Australia provides insights not only into the ocean tides and solid Earth rheology but also into potential systematic errors in the space geodetic technique used to compute the estimates themselves. Separating this systematic error component can be particularly challenging over smaller spatial scales such as coastal study areas. As such, we consider site displacement timeseries based on a range of satellite orbit and clock products, the choice of which was demonstrated by Matviichuk et al. (2020) to affect estimated OTLD in western Europe. The authors demonstrated that the magnitudes of differences between M 2 and O 1 solutions computed using different orbit and clock products reached 0.5 mm in that region. Following various previous studies, we focus on the major lunar semidiurnal and diurnal constituents (M 2 and O 1 , respectively) which are thought to be mostly free from major GNSS systematic errors (e.g. Yuan et al. 2013). We do so using GPS+GLONASS data acquired from sites across the Australian continent. In our study, we remove from the sites' observed displacements the modelled solid Earth body tides and OTLD based on purely elastic Earth models and focus our investigation on the determination of the source of the residual signal. The body tides, while not being perfectly modelled, are often stated to have accuracy of around 1 % (e.g. Yuan and Chao 2012) and are not analysed in this study.

GNSS data and processing
After a preliminary quality assessment, data from 360 Australian continuous GNSS sites were selected for processing, with data spanning from the beginning of 2015 to mid-2020. Most sites are known to be well connected to the local bedrock, although some sites are located on buildings or are on less stable geology. The site locations are shown in Fig. 1 (see Table S1 for coordinates of all sites). There were two selection criteria: a GLONASS-capable receiver with both GPS and GLONASS data availability and near-continuous recording since DOY 001 in 2015. The total data span is close to 2000 days which is in excess of that used by Penna et al. (2015), Martens et al. (2016) and Matviichuk et al. (2020) who demonstrated ∼1000 days was sufficient for the reliable estimation of OTLD. Average data availability is 95 % over all sites with most experiencing a 1-month gap in December 2017, likely due to a system issue in the Geoscience Australia data server. This gap should not impact the results (see Penna et al. 2015).
We use the kinematic precise point positioning (PPP) technique (Zumberge et al. 1997) using NASA JPL's GipsyX software suite v1.3 (Bertiger et al. 2020) to process daily 30 s RINEX files following the approach of Matviichuk et al. (2020). The analysis involves estimating site coordinates and Zenith Wet Delay (ZWD) every 5 min as random walk processes, with process noise settings discussed below. Receiver clock terms are estimated as a white noise process. ZWDs were mapped to the elevation angle of the satellites using the VMF1 mapping function and using a priori zenith wet and hydrostatic delays derived from the European Centre for Medium-range Weather Forecasts (ECMWF) datasets (Boehm et al. 2006). We modelled solid Earth and pole tides according to the IERS2010 Conventions (Petit and Luzum 2010). Our approach to estimate OTLD is discussed below.
A range of final fiducial satellite orbit and clock products from three separate analysis centres (CODE, ESA and JPL) were used. All products were designated 'final' and included CODE MGEX (Dach et al. 2020), ESA Operational (in the sense of being aligned to the most recent homogeneous reanalysis) (Villiger and Dach 2020), and final JPL repro3.0 (NASA JPL 2020) (a non-IGS release, distributed through http://sideshow.jpl.nasa.gov). The products are each provided with different origins: ESA and CODE products (as submitted to the IGS) are provided in the centre of mass of the solid Earth (CE) frame, and JPL products are in the centre of mass of the solid Earth + fluid masses (CM) frame (Blewitt 1989;Petit and Luzum 2010). Both CODE and ESA products were used for processing GPS, GLONASS and GPS+GLONASS solutions in which ambiguities were kept as non-integer estimates, while JPL products were used to provide GPS solutions with ambiguities both floating and resolved to integers (hereon termed AR). The ambiguities were fixed using GPS wide lane and phase bias estimates that are disseminated as part of JPL's native product format (Bertiger et al. 2010). The JPL products do not yet include orbits, clock corrections or wide lane and phase bias estimates for GLONASS; hence, we were unable to produce any solutions including AR with JPL products for GLONASS or Red dots represent site locations. Note the different colour bar ranges between panels. The phase shown with black contour lines is relative to Greenwich. The site 'BRO1', used in process noise optimisation tests, is marked with a red cross GPS+GLONASS. The elevation cutoff angle was set to 7 • for both GPS and GLONASS constellations. Abbaszadeh et al. (2020) reported that solutions using the GLONASS constellation may benefit from decreased satellite cut-off angle for sites at latitudes within ±50 • but our additional tests using these sites with elevation cutoff angles of 3 • and 5 • demonstrated negligible difference with GLONASS-derived OTLD estimates. We did not closely examine the amount of data at these low elevation angles.

OTLD estimation from GNSS
The resulting kinematic timeseries were then resampled to 30 min via averaging (thus reducing data volumes without loss of information at the frequencies required). The resultant time series were then analysed with ETERNA v3.30 as per Penna et al. (2015) with settings recommended by Wenzel (1996) for timeseries of over 1 year in length: M 2 and O 1 were analysed within 1.914129−1.950419 and 0.911391−0.947991 cycles per day, respectively (see ETERNA config file, Listing S1, in the Supplementary Information). ETERNA makes use of the Tamura tidal potential which has 1200 tidal harmonics, some of which mitigate the 18.6-year modulation. This tidal potential also contains harmonics that correspond to the smaller/nonlinear tides such as M 3 and M 4 . The uncertainties on each constituent amplitude and phase are derived from the spectrum of the residuals of the least squares adjustments, using average spectral amplitudes over individual tidal bands (e.g. Wenzel 1996). We used a remove-restore procedure of a priori OTLD values computed using the FES2004 ocean tide model and the Gutenberg-Bullen Earth model as per Penna et al. (2015). These modelled values, though computed using outdated models, are sufficient to reduce the tidal displacements at the computation stage. Figure S16 provides a 7-day snapshot of the JPL GPS AR kinematic solution and the modelled OTLD timeseries at the BRO1 site, similar to (Martens et al. 2016, Fig. 10). The modelled timeseries were generated using the HARDISP program, part of the IERS software collection (Petit and Luzum 2010), using the constituent parameters determined by ETERNA via the procedure outlined above.
We focus on M 2 and O 1 tidal constituents, avoiding the complexities of solar-related constituents, S 2 , K 2 , K 1 , P 1 that possess additional unmodelled effects. Satellite orbits have biases at S 2 , producing biases in both GPS (King 2006;Thomas et al. 2007) and GLONASS OTLD estimates (Matviichuk et al. 2020). Also, the S 2 constituent has an atmospheric component which can be modelled (Tregoning and van Dam 2005) but models are unable to fully explain the residuals when compared to OTLD models. The K 2 and K 1 constituents are known to be problematic with GPS due to the proximity of the repeat periods of individual satellites and the overall constellation geometry. This therefore exposes the constituents to bias from multipath signals (Urschl et al. 2005). GLONASS-only has been found to perform better than GPS for K 2 and K 1 (up to 1 mm smaller residuals in the up component than with GPS or GPS+GLONASS) but the recovered residuals were far larger than the values at other major constituents (Matviichuk et al. 2020). The displacements at lunar N 2 and Q 1 constituents are relatively small which may limit the accuracy of the recovered phase values (Matviichuk et al. 2020). Assessment of the M 2 and O 1 combination has been done previously through gravimetric (Bos and Baker 2005) and GPS observations (Martens et al. 2016;Martens and Simons 2020).
The residual OTLD is computed via vector differencing the observed and modelled OTLD, i.e. Z res = Z obs -Z otl , according to the naming conventions of Yuan and Chao (2012), where Z obs denotes the observed estimates of amplitude and phase with the previously removed FES2004-based OTLD restored prior to estimation.
We consider Z res to consist mainly of the residual OTLD since the body tide that is applied at the PPP processing stage has an accuracy of better than 1 %  away from regions of substantial lateral variation in Earth's structure such as subduction zones (Zürn et al. 1976). With an M 2 body tide amplitude of 80-150 mm over Australia from south to north, a 1 % error can lead to ∼1 mm residuals, which may be significant for this application. Bos et al. (2015) obtained very small GNSS-derived Z res in the middle of Europe, far from the coast indicating that the modelling of the solid Earth tide seems to be quite accurate. We will show in Sect. 3.4.1 that we also consider the body tides to have an accuracy exceeding 1 % over this study region in Australia.

OTLD prediction
A set of seven global ocean tide models was used in the analysis, excluding FES2004 that was used in the initial processing: FES2012 (Carrere et al. 2013), FES2014b (Lyard et al. 2021), GOT4.8 and GOT4.10c (Ray 2013), TPXO8, TPXO9.v1 and TPXO9.v2a (Egbert and Erofeeva 2002). The older models were included for consistency with previous studies and also to demonstrate their deficiencies as they are still being actively used by the geodetic community. During its creation, the GOT4.10c model was explicitly corrected for the effect of altimeter observations being in the CM frame (hence the 'c') (Desai and Ray 2014). The FES2014b loading tide used GOT4.8ac data in the data assimilation process for the preliminary version of the FES2014b ocean tide model (Lyard et al. 2021), with GOT4.8ac also using the corrected altimeter data. By contrast, the TPXO9 models assimilated altimetry datasets that were not corrected for geocentre motion (pers. comm., Svetlana Erofeeva, 02 Feb 2022).
OTLD predictions were computed using the ocean tide models listed above with a set of Green's functions (see e.g. Bos and Scherneck 2013) computed with each of the PREM (Dziewonski and Anderson 1981), STW105 and S362ANI (Kustowski et al. 2008) Earth models (Farrell 1972) using the CARGA software (Bos and Baker 2005) within the free ocean tide loading provider (http://holt.oso.chalmers.se/ loading). We replaced the water layer in PREM with the density and elastic properties from the underlying rock layer since our GPS stations are on land. S362ANI is a transversely isotropic seismic tomographic model for the upper mantle. We computed the mean shear velocity of S362ANI at various depths underneath Australia and used these values together with the density and compressional velocities to construct a loading Green's function. We used a land-sea mask based on the GMT high resolution coastline (Wessel and Smith 1996), converted to a regular longitude latitude grid of 0.00137×0.00137 • .
In addition to purely elastic Earth models, we also consider the effects of mantle and asthenosphere anelasticity. We do this following the approach of Bos et al. (2015) by accounting for the frequency-dependence of the shear modulus as the period of the stress cycle is increased at tidal frequencies relative to the 1 Hz reference frequency of seismic Earth models. The load Love numbers are computed using complex-value elastic constants, resulting in real and imaginary-value Green's functions. However, since the latter only changes the phase-lag by less than 0.2 • , only the real-valued Green's function is used. The single M 2 constituent was selected to represent the tidal frequency band as the shear modulus difference between the diurnal and semidiurnal frequency bands was found to be negligible. The Q quality factor provided with the reference Earth models was assumed constant, i.e. independent of frequency (Lambeck 1988;Bos et al. 2015). The depth-dependent values of Q used per model are listed in Table S4. In all cases we use radially symmetrical Earth models.
We also considered the effects of spatially varying water density and water compressibility (Ray 2013). Spatial water density and compressibility corrections were demonstrated to decrease the M 2 Z res in New Zealand at the same level as the dispersion correction which was ∼0.2 mm in the up component (Matviichuk et al. 2021). The spatial water density values were extracted from a 0.25×0.25 • grid of the World Ocean Atlas 2013 (Zweng et al. 2013) and included using the method of Ray (2013). We add 'd' and 'c' suffixes to the Green's function name (e.g. STW105dc) to reflect when anelastic dispersion ('d') and spatial water density and compressibility corrections ('c'), were applied in CARGA. The OTLD is modelled in the centre-of-mass frame that is consistent with each of the orbit and clock products (Fu et al. 2012;Blewitt 2003): CM for JPL and CE for ESA and CODE. We note that individual analysis centre choices in generating orbit and clock products may affect the representation of CM and CE in a complex way such as due to imbalanced station geometry or orbit mismodelling, and the representation may vary over time such as due to changing station geometry at the orbit and clock estimation stage. The testing of the kind undertaken in this study sheds some light on these respective differences.

Results and discussion
We first discuss the results based on each of ESA, CODE and JPL products in terms of determining optimal process noise and intercomparison relative to the FES2014b and STW105dc OTLD modelled values. Second, we compare the results with the OTLD based on the different ocean tide models and Green's functions, and assess the residuals computed with various products. This is followed by the detailed assessment of the OTLD residuals computed using the three constellation mode solutions: GPS-only, GLONASS-only and GPS+GLONASS. Finally, we analyse the observed biases related to the sub-daily centre of mass (CoM) tidal variations and assess inter-model differences, and the effects of spatial water density and water compressibility corrections.

Process noise sensitivity tests
Our first step was to establish the optimum process noise settings for the coordinate and ZWD parameters, repeating similar work done for the UK Matviichuk et al. 2020) who found optimal settings for coordinate and ZWD process noise values of 3.2 and 0.1 mm/sqrt(s), respectively. Martens et al. (2016) focused on a South American dataset and used an alternative approach for determining optimum process noise through the analysis of the OTLD RMS misfits and found optimum values of 0.5 and 0.05 mm/sqrt(s) for coordinate and ZWD process noise, respectively.
We tested optimal process noise settings using BRO1, a site with the maximum OTLD amplitude experienced in the chosen Australian network. Its location is marked with a red '+' in Fig. 1. We followed the approach of Penna et al. (2015) and Matviichuk et al. (2020) and chose the optimum coordinate process noise value (3.2 mm/sqrt(s)) as the value immediately after the magnitude of the residuals stabilised for all eight major tidal constituents while also minimising the standard deviation of the coordinate timeseries. The choice of the ZWD process noise value (0.1 mm/sqrt(s)) was governed by the minimum standard deviation of the differences between the ZWDs computed with kinematic and static PPP ( Fig. S1.1). These process noise values are the same as found for the UK by Penna et al. (2015). The tests were repeated using JPL, ESA and CODE products, and the conclusions remained the same giving confidence that the process noise settings were robust. These values were used for all further tests.
Variations to the optimal settings may occur given the very different ZWD conditions in other parts of Australia. Nevertheless, the rate of ZWD variation is most critical in setting process noise parameters and these rates do not vary substantially across the planet. This is further supported by our identification of similar optimal settings to those found using BRO1 at other sites located in a very different climatological setting in western Europe and, specifically, Great Britain (Matviichuk et al. 2020;Penna et al. 2015). We do note, however, a difference to those settings found in South America (Martens et al. 2016) which may warrant further investigation in the future.

OTLD absorption by tropospheric parameters
We expanded the test above to assess the presence of tidal signals within the estimated ZWD timeseries. Such signal could indicate leakage of OTLD due to overly constrained (i.e. too tight) coordinate process noise. We did this by vector differencing signals at tidal frequencies extracted from a priori ECMWF-derived ZWD timeseries and those derived from the estimated ZWD at site BRO1 ( Fig. S1.2). We found that tidal frequency signals in estimated ZWD are insensitive to variations of the coordinate process noise. Iterating over various ZWD process noise values showed very low levels of propagation of tidal signal into resulting ZWD when using process noise values below the default value of 0.1 mm/sqrt(s). Tidal signals in parameter estimates did not change with higher process noise values and this held for all assessed constituents. From this analysis, we consider our optimal process noise values to be robust. Figure 1 shows the amplitudes of M 2 and O 1 OTLD based on the mean of seven of the most recent global ocean tide models and S362ANI Green's function corrected for dispersion at M 2 and their respective standard deviation (SD) maps. The OTLD amplitudes range from 0.1 to 31.0 mm for M 2 and from 0.15 to 8.5 mm for O 1 (Fig. 1a, b), with maximum OTLD for both constituents in the north west and west of Australia. The SD between OTLD grids at M 2 and O 1 is negligible for most of the continent but reaches over 3 mm in localised areas. These variations in OTLD and OTLD SD are reflected in differences in these regions between ocean tide models. Following Bos et al. (2015), we also computed the mean ocean tide models and SD maps for both M 2 and O 1 constituents after resampling all ocean tide models to a consistent 0.05×0.05 • grid (Fig. S2). Areas with high ocean tide model SDs may be regions where errors exist in the ocean tide models that may impact the estimated OTLD, though this does not imply that all models have large errors in these regions. We regard agreement between models as an indication the models are robust. In areas of disagreement, we interpret each model having potential errors in the relevant region.

Comparison of ocean tide models around Australia
While the spatial-mean SD of the ocean tide models over the total area is less than 0.15 m for M 2 and 0.1 m for O 1 , some shallow shelf areas demonstrate localised SD anomalies. Regarding M 2 , high SD values (over 0.2 m) are observed in the north and north-east of Australia (Fig. S2c) where the M 2 amplitude exceeds 2.5 m. Modelled O 1 has large, highly-localised SDs of up to 0.2 m in the north of Australia ( Fig. S2d) but is generally otherwise small. Considering the areas of localised SD between ocean tide models at M 2 and O 1 , these correspond to the areas of high SD between OTLD grids (comparing Fig. 1c, d and Fig. S2c, d).

Observed OTLD
OTLD estimated from the JPL GPS AR solutions, after restoration of the signal removed during processing is shown for the 360 sites across Australia in Fig. 2 for M 2 (left) and O 1 (right) for each of the east, north and up components. Each panel shows strong spatial coherence giving confidence in the robustness of the estimates. This is unsurprising given that the amplitudes regularly exceed 10 mm in all coordinate components for both M 2 and O 1 . As expected, signals are generally largest nearest the coast and decay inland, although substantial signal exists in the interior for various constituents and components (e.g. O 1 north, M 2 east). Particularly, striking spatial patterns are evident along the densely observed east coast with M 2 OTLD in the up component reaching around 15 mm. As with the OTLD predictions, the largest observed M 2 OTLD is found in the Broome region of north-west Australia where site BRO1 shows the maximum M 2 OTLD amplitude in the up component of 30 mm. Given negligible spatial variation in the formal errors, the OTLD confidence ellipsoid shown in Fig. 2 is the mean ellipsoid from all sites. The same approach is used in Figs. 3, 4, 6 and 7.
For M 2 up, the phase of the signal is shown to rotate along the east coast and for several hundreds of km inland where the amplitude reduces to ∼3 mm in central Australia. The decay in amplitude towards central Australia is also seen in the east and north components of M 2 though with smaller amplitudes.
The observed O 1 amplitudes in the up component decay towards the east coast, while in the east and north components all sites demonstrate OTLD of ∼4 mm in amplitude with no particular change across the continent. The phase of O 1 OTLD in the east and north components is each quite spatially uniform (0 • and −90 • , respectively, Fig. 2 While the OTLD field in Fig. 2 is limited to GPS AR solutions, results from any constellation are broadly similar when plotted at this scale. Thus, we proceed with analysis of the GNSS OTLD estimates differences with modelled OTLD and discuss potential systematic errors.

Comparison of solutions based on different orbit and clock products
In Fig. 3a Fig. 3 show the up component Z res but for GPS+GLONASS solutions using ESA (Fig. 3c, d) and CODE (Fig. 3e, f) products. M 2 panels (Fig. 3a, c, e) demonstrate a high degree of similarity in Z res between the different solutions. The M 2 ESA Z res values are rotated along the east coast and are much larger in northern Australia than with JPL or CODE solutions, while CODE M 2 Z res is larger in western Australia. The M 2 Z res in the horizontal components is also close between solutions (Fig. S4.1).
The mean of the M 2 residual magnitudes, Z res , is ∼0.22, ∼0.23 and ∼0.40 mm in the east, north and up components, respectively, independent of which orbit and clock products were used in generating GPS+GLONASS solutions. Comparable residual magnitudes were observed when using JPL GPS AR solutions. The M 2 up value is slightly smaller for both sets of GPS+GLONASS solutions compared to GPS AR adding further to evidence that adding multiple constellations improves estimates of tidal displacements. The M 2 Z res values from solutions based on non-integer ambiguity GPS-only and GLONASS-only solutions are mostly larger than the above GPS AR or GPS+GLONASS values independent of the products and component (see Table 1). Despite the inter-solution similarity, the estimates based on JPL products (independent of AR) in the east component (Fig. S4.1a) show a lateral change of Z res ("pivoting" relative to 130 • E longi- tude) between eastern and western Australia, relative to both ESA (Fig. S4.1d) and CODE (Fig. S4.1g) GPS+GLONASS Z res solutions. In these cases, the difference in residuals, Z res , between JPL and either of ESA or CODE reaches 0.4 mm. We do not understand the origin of this signal, although it must originate in the orbit and clock products.
We also directly compared Z res estimates derived from GPS+GLONASS and GPS AR by computing the mean of the respective Z res vector difference magnitudes, Z res . The CODE GPS+GLONASS Z res in the north and up components are close to that for JPL GPS AR solutions with Z res of ∼0.1 mm while ESA GPS+GLONASS Z res is close to the JPL GPS AR value only in the north component ( Z res <0.02 mm). This is different from the results reported in Matviichuk et al. (2020) where the Z res based on ESA GPS+GLONASS solutions were found to be closer, over all coordinate components, to those derived from JPL AR solutions, compared to those from CODE solutions. This change may be related to the differences in products used - Matviichuk et al. (2020) used purely IGS repro2 products, whereas in this study we use operational products post-repro2 and these may include analysis developments since then.
Regarding the O 1 constituent, the Z res phasors show rotations and magnitude variations between solutions computed using JPL, ESA and CODE products (Fig. 3b, d, f). In particular, JPL AR and ESA residuals appear rotated by about 90 • along the east coast (comparing Fig. 4d and e). Visually, estimates based on CODE products have the smallest O 1 Z res in the up component (Figs. 3f and 4f). As shown in Table 1, solutions based on CODE products demonstrate the smallest Z res in the up component, while ESA-based solutions have the smallest Z res in the horizontal com- Our tests in this section show that there are small but important differences in GNSS OTLD estimates that depend on product selection and choice of satellite constellations and ambiguity resolution that are not yet fully understood. Future products that allow multi-GNSS with ambiguity fixing in GipsyX may help resolve which is the more important of these choices.

Sensitivity of modelled OTLD to different ocean tide models
Now considering different ocean tide models used in the modelling, OTLD computed relative to GOT4.10c are generally comparable to those computed relative to FES2014b, as described above in Sect. 3.2. There are, however, localised differences between modelled OTLD values in the up component. These are concentrated in M 2 and mostly in Bass Strait (Fig. 1c) and several areas in the north and north west of Australia with Z res reaching 1.5 mm in the very north, at the Thursday Island site (TITG) in the Torres Strait (see Fig. S6.1a-c). In terms of O 1 , the differences in modelled displacements between the two tide models do not exceed 0.1 mm and are usually half of that, with Z res only 0.04 mm (Fig. S6.2a-c).
The M 2 up results using TPXO9.v1 and TPXO9.v2 models show similar results to those using FES2014b and GOT4.10c, with reduction of the observation residuals in Bass Strait and the east coast (comparing GPS+GLONASS Z res in Fig. S4.1 and S5.1). The residuals are still larger in Bass Strait than in other locations, pointing to residual tide model errors in the two TPXO9 variants. No noticeable difference is present in the horizontal components of M 2 (Fig. S6.1d-f). The O 1 TPXO9.v1 OTLD values are very close to those computed with FES2014b (mean difference of 0.03 mm) with marginal increase in Z res for the up component (0.08 mm) as demonstrated in Fig. S6.2d-f. In terms of vertical OTLD over Australia, the difference between TPXO9.v1 and TPXO9.v2a is negligible, thus we limit our discussion to TPXO9.v1 (noting that by extension, we expect the conclusions to apply for both versions).

Sensitivity of modelled OTLD to Green's functions and corrections
The three Green's functions used in the analysis, computed with PREM, STW105 and S362ANI, do not produce a noticeable difference in terms of residual magnitude distribution for any given ocean tide model (see Fig. S9.1−9.3 and Fig. S10.1−10.3 for M 2 and O 1 , respectively). One of the explanations for this is the high number of inland stations, away from the coast, where OTLD is small, thus the effect from changing the Green's function is insignificant. We now assess, in turn, the effects of considering the corrections to the previously defined Green's function due to anelastic dispersion correction at M 2 . In addition, instead of using a constant mean sea water density value, we use a spatially varying one that is a function of the depth of the ocean. The deeper the ocean, the higher the mean density in each column of water. We also take the compressibility effect into account which increases the water density value even further, see Ray (2013). Figure 5 shows the combined effect from such corrections to M 2 (a, c, e) and O 1 (b, d, f) in the up component. These tests were conducted using FES2014b and STW105 and have very similar impacts with the other tide models and one-dimensional Earth models.
Large differences between OTLD values computed using a Green's function based on a seismic Earth model valid at 1 Hz and a Green's function that has been adjusted for anelastic dispersion will only occur if the OTLD value itself is large. Since many sites are located well inland where the OTLD values are small, this difference for the M 2 constituent in the up component is generally negligible (Fig. S8.1a, S9.1−9.3). The dispersion correction is only significant in the up component at the 50 coastal sites located along the east and north-west coasts of Australia (marked with * in Table S1). The impact of the dispersion correction is less uniform over Australia than in the previous studies of western Europe Matviichuk et al. 2020). Results from the coastal sites demonstrate a 0.2 mm improvement of the mean Z res at M 2 due to the anelastic dispersion correction as previously reported by Matviichuk et al. 2020Matviichuk et al. , 2021. This also reduces the variance regardless of the underlying Earth or ocean tide model (Fig. S9.4−9.6). The effect at M 2 from the dispersion correction in the east component is largest in the south-east and north of Australia (Fig. 8.1a). The correction phasors in the north component rotate in south-east Australia, close to Bass Strait, and are aligned, with a magnitude of up to 0.3 mm in the west of Australia (Fig. S8.1b).
Altering the elastic constants of the Earth model from 1 Hz to the M 2 frequency due to dispersion also affects other tidal constituents including O 1 . The dispersion correction at O 1 is usually negligible (∼0.03 mm) over all components, reaching 0.2 mm in the up component in the very north of Australia, close to Torres Strait (Fig. S8.2a-c, S10.1−10.3). We also tested altering the elastic constant from 1 Hz to the O 1 constituent frequency, but the impact from the correction on O 1 and other constituents was found to be negligible (<0.02 mm).
The spatial water density and compressibility corrections demonstrate a uniform effect over both inland and coastal sites, but the average effect is small. For M 2 , the average magnitude of the correction effect is 0.1 mm in the up and east, and negligible (∼0.04 mm) in the north component (see Fig. S8.1d-f). The correction reaches a maximum of 0.2 mm in the up component at the sites that experience maximum OTLD. At O 1 , the average correction magnitude is half of that of M 2 but shows a 0.1 mm magnitude in the up component at the sites in western Australia (Fig. S8.2d-f).

Comparison of solutions using GPS, GLONASS or GPS+GLONASS
We now provide an assessment of the impact of the choice of constellation configuration on the derived OTLD at M 2 and O 1 . For this we fixed the OTLD model predictions to the combination of the FES2014b tidal model and STW105dc Green's function. Previous analysis (e.g. Abbaszadeh et al. 2020;Matviichuk et al. 2020) has shown that GPS+GLONASS performs close to GPS AR at lunar constituents, especially at M 2 , while GPS AR has greatest impact in the east component. Figure 6 shows residual phasors derived from GPSonly (blue), GLONASS-only (orange) and GPS+GLONASS (green) solutions using ESA products. In addition, direct differences between GLONASS-based and GPS-only constellation modes were computed. The results are illustrated in Fig. 7, which provides a detailed map of differences between GLONASS-based (GLONASS-only and then also Fig. 5 The modelled changes in OTLD values when the effects of dispersion in the asthenosphere and a spatially varying water density are taken into account for the east (top), north (middle) and up (bottom) components at M 2 (left panels) and O 1 (right panels) GPS+GLONASS) derived OTLD and GPS-only derived OTLD. All constellation modes perform similarly in the north component at both M 2 and O 1 (Fig. 6c, d) with the most visible difference in the east component (Fig. 6a, b).
Z res , a Z res vector difference computed at each site between different constellation mode solutions, reveals that GPSonly and GPS+GLONASS OTLD estimates are close: in the up component, the M 2 mean Z res between GPS and GPS+GLONASS is ∼0.13 mm for either ESA and CODE products, while for O 1 it is ∼0.15 mm for ESA and ∼0. ESA and CODE GPS+GLONASS Z res estimates have the smallest uncertainties (SD) at both M 2 and O 1 , producing the smallest mean confidence ellipsoid over all components (Fig. 6). The confidence ellipsoids were estimated from amplitude and phase uncertainties derived from  (Figs  S4.1, S4.2) presumably due to the impact of GPS ambiguity fixing on the east component (Blewitt 1989). The GLONASS mean uncertainty is the largest, independent of the component. Note that uncertainties here are derived from ETERNA, reflecting timeseries noise, and are not reliant on the formal uncertainties of the coordinate timeseries. The average uncertainty of the estimated amplitude values for both constituents is less than 0.01 mm, while the average uncertainty of the phase values is on the level of 2-5 • . Table 2 thus ignores the mean amplitude uncertainty values and presents only the uncertainty of the phase values.
We further illustrate the variations in the OTLD residuals as the constellations and products are varied by using a mod-ified empirical cumulative distribution functions (ECDFs) as per Martens and Simons (2020) (Fig. 8). Curves that rise most quickly to 1.0 reflect a relatively better overall fit of the modelled values with the observations. Note that each ECDF curve has had its respective mean Z res value subtracted (normalised). This removes possible common mode errors/biases and shows the pure distribution of fit (see also Martens et al. 2016). Unmodified ECDF curves are demonstrated in Fig. S11. Figure 8 shows that the horizontal components of the M 2 and O 1 residuals show close agreement with the models, with the ECDFs converging to 1.0 at 0.3 mm with the only exception being GLONASS-based solutions that reach 1.0 at 0.5 mm in the cases of M 2 (ESA) and O 1 (both ESA and CODE). Note that, with the mean Z res removed  (normalised), no obvious difference between GPS AR and GPS+GLONASS fit in the east component is present (as reported in Sect. 3.5), independent of the products used. Without the normalisation, the ECDF for GPS+GLONASS CODE solutions diverge in the 0.25 to 0.5 mm window at O 1 (Fig. S7f). This suggests the presence of a common signal in the east component of the CODE GPS+GLONASS solutions at O 1 which is affecting a significant number of sites. The steepening of the ECDF slope after the normalisation suggests that similar effects are present to a lesser degree in all the solutions, for both constituents. For the up component, the normalised Z res magnitude using JPL GPS AR solutions shows better fit (ECDF is 1.0 at 0.8 mm) than CODE or ESA GPS+GLONASS solutions (ECDF is 1.0 at 1mm) for M 2 , with similar performance for O 1 (see Fig. 8). GPS-only (non-AR, solid blue lines) M 2 ECDFs are similar across ESA and JPL products (Fig. 8a,  b), as reported previously in Matviichuk et al. (2020) for western Europe. However, there is a slightly reduced agreement with CODE that demonstrates GPS-only performance similar to that of GLONASS solutions. The GPS-only O 1 solutions demonstrate better convergence (sharper rise of the ECDF curve) than GLONASS solutions in the up component, independent of the satellite products. GLONASS solutions demonstrate the largest residuals at both M 2 and O 1 , while performing similarly to GPS-only in the east component at M 2 and O 1 and in the up component at O 1 .

Large-scale bias in TPXO9
Next, we fix the Green's function to STW105dc while iterating through the ocean tide models. Using the TPXO9.v1 ocean tide model, residuals at M 2 are reduced relative to those computed with FES2014b and GOT4.10c along the east coast of Australia and around Bass Strait by up to 0.2 mm for sites at the coast (Fig. S5.1 and S5.2). This, however, is only observed with ESA and CODE product solutions. JPL product solution shows consistently higher magnitude Z res (bias, for simplicity) when using the TPXO9.v1 ocean tide model rather than FES2014b or GOT4.10c at M 2 and O 1 by ∼0.2 and ∼0.3 mm, respectively, at all sites (Figs S5.1a-c and S5.2a-c). To illustrate the issue, we computed a vector difference between FES2014b and TPXO9.v1 modelled values (CM) at M 2 and O 1 (Fig. 9a, d). We next explore potential contributors to the biases evident with TPXO9.v1-based models.

Assessment of CoM values
We now investigate the absence of the large-scale residuals bias when applying TPXO9.v1 in the CE frame. With ESA and CODE solutions, its similarity between sites suggest that the difference is related to the CoM correction associated with far-field differences between TPXO9.v1 and both FES2014b and GOT4.10c. This is also confirmed by the absence of significant ocean tide models anomalies in the region (Fig. S2). First, we computed the vector difference between FES2014b and TPXO9.v1 OTLD values each given in the CM frame (Fig. 9a, d). Next, we computed CoM values by subtracting CE OTLD values from their respective CM for M 2 and O 1 (Fig. 9b, e). The alignment of both phasor maps strongly suggests that the observed difference found between the modelled values is in the CoM correction. Note that both CM and CE OTLD are computed through convolution of the modelled tidal mass with respective CM/CE Green's function. As such, Subtracting the TPXO9 CoM bias from the OTLD difference makes the resulting TPXO9/FES2014b vector difference close to that when using solutions computed with CE satellite products, with localised differences in the up component related to the Bass Strait, Timor Sea and the Great Barrier Reef at M 2 and to the east coast at O 1 (Fig. 9c, f). Exactly the same effect is observed for the horizontal components but with less localised differences for both M 2 and O 1 (Fig. S12.1 and S12.2, respectively). The observed CoMrelated bias reflects global-scale mass changes as opposed to the more regional or local-scale variations that dominate the loading deformations.
Fundamentally, the tidal variations in the CoM of all water masses on Earth are counter balanced by tidal variations in the position of the CoM of the solid Earth. The tidal displacements of the CoM (an integral part of the CM frame) were computed for the M 2 and O 1 constituents and the TPXO9.v1, FES2014 and GOT4.10c models to assess the difference between the ocean tide models in terms of the resulting CoM (Table 3). Next, for each position in our area of interest, we converted this tidal motion of the CoM into local east, north and up displacement components. The largest components of CoM displacement are shown in Fig. 10 as amplitudes and phase-lags for M 2 up and O 1 north. We then rotated each of the in-phase and out-of-phase CoM components to the XYZ coordinate system, which is effectively the same set of XYZ values at each site over the spatial scale of the Australian network. The observed noise associated with numerical precision at the level of <0.05 mm was removed with averaging of values over sites. We also observed small CoM biases created by introducing spatial water density and compressibility corrections (∼0.05 mm for M 2 and ∼0.02 mm for O 1 ) in each of X, Y and Z components (Fig. S14.1). We consider these to be negligible. We also explored the effect of water mass conservation but we found it to be negligible as well (see Fig. S14.2).

Regional sources of the CoM differences
We now explore the potential of large-scale tide model differences to produce the observed TPXO9.v1 CoM differences. Such differences may occur due to subtle differences in handling altimeter data assimilated into many of the models. For instance, Martens et al. (2016) have demonstrated that predicted OTLD is sensitive in important ways to how assimilated altimeter data was treated in the generation of the tide  model. In particular, there is a small but important effect due to the correction (or not) of the altimeter for CoM variations.
To explore this further, we assess the effect of distant tides (i.e. distant from Australia) by focusing on five ocean tide models of which three are the most recent: FES2014b and GOT4.10c, with their earlier versions, FES2012 and GOT4.8, respectively, and TPXO9.v1. The two older models were included to highlight the effect of the CoM correction applied in the most recent models. We divided the global tide models into seven polygons covering the main ocean areas, replicating the conventional geographical division (Fig. S15). The resulting CM and CE values, in ENU, were vector differenced to obtain CoM values and rotated to XYZ. CoM values from each global water region were then used to construct complete CoM phasors (Fig. S13). Overall, the Pacific region has the largest impact on the CoM, however the Atlantic and Indian oceans have comparable impact in the Y and Z components of M 2 (Fig. S13b, c). The effect from the Pacific is also the largest in O 1 with the exception of the X component where the Indian and Atlantic water regions are comparable to that from the Pacific (Fig. S13d).
All reconstructed CoM phasors are similar at M 2 while showing minor differences at O 1 (Fig. S13). The most significant differences in O 1 were in the Z component where TPXO9.v1 diverges from FES2014b and GOT4.10c (Fig. S13f).
We next compare TPXO9.v1 with the other four selected ocean tide models by vector differencing with each regional phasor to better study the region-specific differences. If the CoM bias is concentrated in a single area, we should see it as a single phasor in the vector difference. Figure 11 clearly demonstrates that the TPXO9.v1 CoM deficiency relative to the two global ocean tide models CoM (FES2014b and GOT4.10c) is distributed over multiple regions of the globe. FES2012 CoM shows the smallest variations relative to TPXO9.v1 (Fig. 11). In addition, one can observe how FES2014b and GOT4.10c are similar at M 2 and O 1 in all components with the exception of M 2 in the Z component, presumably associated with the Antarctic, Arctic, and Atlantic regions. In all other cases, FES2014b and GOT4.10c are in closer agreement than TPXO9.v1 and any other model.
The absence of the geocentre motion correction in the assimilated altimetry datasets used in the computation of the TPXO9 models may be the reason for the observed difference, especially considering the similarity between TPXO9.v1 and the preceding models. The exact source of the CoM bias in TPXO9.v1 is not certain. We note the bias is also present in TPXO9.v2a and we have not examined earlier versions of the TPXO tide model series. While both TPXO9 variants can safely be used for OTLD modelling in the CE frame, the CM correction computed with these models should be ignored in favour of the correction computed from either of FES2014b or GOT4.10c. An example of such a use-case is processing a PPP solution using orbits and clocks in the CM frame as in the case of native JPL products. In such cases, the CM correction is favoured to ensure OTLD are applied in a comparable frame.

Discussion and conclusions
This paper investigates GNSS-derived ocean tide loading displacements (OTLDs) of semidiurnal M 2 and diurnal O 1 lunar constituents across 360 sites distributed across continental Australia. By comparing GNSS-derived OTLD against those from various models, we sought to gain insight into model performance and possible impact of GNSS-specific systematic errors. We assessed time series computed using GPSonly, GLONASS-only and GPS+GLONASS using ESA, CODE and JPL orbit and clock products. The processed dataset reveals large scale biases, in particular in the east coordinate component across western Australia (specific to JPL GPS-only solutions). We also discovered high-resolution coherent residual OTLD variations over a dense network of sites adjacent to Bass Strait (southern Australia).
Our study fits between global studies that investigate a broad range of tide and Earth structures (e.g. Yuan et al. 2013) and more specific smaller-scale studies that investigate the response to complex tides (e.g. Wang et al. 2020;Penna et al. 2015;King et al. 2005;Yuan et al. 2009;Martens and Simons 2020) or complex Earth structures (Matviichuk et al. 2021;Ito and Simons 2011). Our work fills a knowledge gap in OTLD interpretation over regional scales in a similar way to the work of Yuan and Chao (2012) over the United States, Martens et al. (2016) over South America and Bos et al. (2015) over western Europe.
The deficiencies of modern ocean tide models are highly localised, generally limited to shallow and polar water regions (Stammer et al. 2014). While most ocean tide models and Green's functions resulted in similar modelled displacement values, the TPXO9.v1 model resulted in a noticeable bias at all sites when using products in the CM frame. This bias was found to originate in the centre-of-mass (CoM) motion. We assume that the bias is caused by the fact that satellite altimetry data used in generation of this tide model were not corrected for the CoM motion. Note that TPXO9.v1 can outperform other tide models in some areas (even with the CoM bias present) due to better bathymetry, higher spatial resolution of the grid and tide gauge data used during the assimilation process.
While the residual OTLD magnitudes for both M 2 and O 1 are well within the reported range from the previous OTLD studies in western Europe , South America (Martens et al. 2016) and Alaska (Martens and Simons 2020) independent of the orbit and clock products used, no single set of products consistently provides the best solution (see Table 1). The differences between solutions vary through constituents and constellation modes, producing rotations in residual phasors specific to an analysis centre and constellation mode, and dependent on the coordinate component-negligible differences in the north component, increasing for the east and up components, respectively. We did, however, notice a regional change in M 2 residual amplitude of up to 0.5 mm in the east component from east to west of Australia between solutions using JPL and ESA or CODE products (Fig. S9.1-S10.3). These conclusions do not align with those previously reported by Matviichuk et al. (2020), which suggests these findings are regionally specific (noting GLONASS-only solutions were confirmed to perform poorly, independent of the products provider). Similarly, while not being recommended on its own, GLONASS clearly adds information content to GPS alone, as solutions using the GPS+GLONASS constellation mode demonstrated comparable results to those of GPS with ambiguities resolved (AR).
The study also provides strong evidence for the inclusion of dispersion and spatial water density and compressibility corrections into modelled OTLD, specifically for the coastal sites. The addition of a dispersion correction to the OTLD models reduced the residual OTLD amplitude of M 2 by 0.2 mm in the up component at the coastal sites and in the east component at inland and coastal sites, independent of the GNSS products used. Applying spatial water density and compressibility corrections produced a further 0.2 mm reduction of residuals at the coastal sites in the up component. These findings at the coastal sites align with those reported previously in Matviichuk et al. (2021Matviichuk et al. ( , 2020, however, demonstrate better behaviour of the corrections for inland sites in every coordinate component (as per Bos et al. (2015), however, their study only considered the up component) even more so than a study of South America (Martens et al. 2016), due to, for example, complexities in modelling OTLD in the Amazon River delta and on the Patagonian shelf.
Accounting for the effect of anelasticity, in the form of a dispersion correction, and for the effect of the spatial water density and compressibility, can explain 0.2−0.5 mm of unmodelled signal depending on a site's distance from the coast. If these effects are not accurately treated they will propagate into conventional 24-hr position solutions (Penna et al. 2007). Correcting for these effects should contribute to the enhancement of satellite products, prompt further studies of Earth's anelastic response, and result in increased inclu-sion of GNSS-derived tidal estimates into the computation and assessment of tide models.
Overall, our results show that while GNSS-derived estimates of ocean tide loading displacements show generally close agreement with models in Australia (∼0.2−0.4 mm, see Table 1), important systematic differences (up to 0.5 mm) remain that depend on both the GNSS products and the ocean tide models, in addition to contributions from uncertain elastic and anelastic Earth properties and body tides. These systematic differences are evident in continental scale observations that may not be otherwise clear in smaller-scale networks most commonly studied to date (e.g. Wang et al. 2020;King et al. 2005;Yuan et al. 2009;Martens and Simons 2020;Matviichuk et al. 2021). This motivates future research to understand the exact origin of these errors. The accuracy of the body tide also requires further investigation (e.g. Lau et al. 2017). Such a task requires a global set of stations, ideally far from the coast to avoid errors in OTL, to eliminate the degree-1 motion (CoM) and see the degree-2 (body tide) deformation more clearly. The small magnitude of residual OTLD observed at sites in central Australia (well distant from the dominant effects of OTL) is consistent with modelling of the body tide with an accuracy of better than 1 % -likely closer to 0.5 % in this region.
Residual errors driven by systematic differences in clock and orbit products are important to understand given they have similar magnitudes to OTLD residuals themselves. If GNSS-estimated OTLDs are used to investigate Earth rheological structures, they will be affected by orbit and clock product-related discrepancies at some level. Examples of datasets that could largely benefit from the reduced discrepancies due to orbit and clock product inconsistencies include those that have complex regional Earth structure, especially if this includes proximity to the active continental margins, e.g. South America (Martens et al. 2016), Alaska (Martens and Simons 2020) and New Zealand (Matviichuk et al. 2021). The availability of a reprocessed multi-GNSS set of final orbit and clock products from JPL and support of IGS phase bias products (Schaer 2018) could help resolve the high variability of the residual OTLD values. This could decrease the differences between solutions computed with different ACs' orbit and clock products and, potentially, reduce the effect of the regionality of the product-related conclusions and enable detailed geophysical interpretation of the observed residuals.