The first Australian gravimetric quasigeoid model with location-specific uncertainty estimates

We describe the computation of the first Australian quasigeoid model to include error estimates as a function of location that have been propagated from uncertainties in the EGM2008 global model, land and altimeter-derived gravity anomalies and terrain corrections. The model has been extended to include Australia’s offshore territories and maritime boundaries using newer datasets comprising an additional ∼\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}280,000 land gravity observations, a newer altimeter-derived marine gravity anomaly grid, and terrain corrections at 1″×1″\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1^{\prime \prime }\times 1^{\prime \prime }$$\end{document} resolution. The error propagation uses a remove–restore approach, where the EGM2008 quasigeoid and gravity anomaly error grids are augmented by errors propagated through a modified Stokes integral from the errors in the altimeter gravity anomalies, land gravity observations and terrain corrections. The gravimetric quasigeoid errors (one sigma) are 50–60 mm across most of the Australian landmass, increasing to ∼100\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }100$$\end{document} mm in regions of steep horizontal gravity gradients or the mountains, and are commensurate with external estimates.


Introduction
We describe the computation of the gravimetric quasigeoid component of AUSGeoid2020, herein called Australian gravimetric quasigeoid 2017 (AGQG2017), as well as the computation of error estimates as a function of location (also called geographic specificity by Pavlis and Saleh 2005). These location-specific errors have been propagated from a combination of uncertainties in EGM2008 (Pavlis et al. 2012(Pavlis et al. , 2013, land and altimeter-derived gravity anomalies, and terrain corrections (McCubbine et al. 2017). In a companion paper (in prep), we describe the geometric component of AUSGeoid2020, where AGQG2017 has been distorted by least squares prediction to fit the Australian Height Datum (AHD; Roelse et al. 1971) on land.
AUSGeoid09 ) included a geometric component where the AGQG2009 gravimetric quasigeoid model ) was distorted to fit the AHD using cross-validated least squares prediction (Featherstone and Sproule 2006). This yielded a surface for the more direct transformation of GNSS-derived ellipsoidal heights to the AHD (cf. Milbert 1995;Featherstone 1998Featherstone , 2008Smith and Roman 2001), allowing Australian land surveyors to realise AHD heights directly using GNSS, rather than having to apply post-survey adjustments as with previous Australian quasi/geoid models.
Geoscience Australia (GA), the national geodetic agency, is in the process of moving from the Geocentric Datum of Australia 1994.0 (GDA94) to GDA2020, which is based on the International Terrestrial Reference Frame 2014 (ITRF2014; Altamimi et al. 2016) extrapolated to epoch 2020.0 using Australian ITRF2014 station velocities. This datum change will cause horizontal geodetic coordinates to move ∼1.8 m in a north-easterly direction and ellipsoidal heights to increase by ∼90 mm. For full three-dimensional (3D) implementation of GDA2020, there remains the need to compute an accompanying model of the separation between the AHD and GRS80, as was done for AUSGeoid09, to be called AUSGeoid2020.
Though the geometric component of AUSGeoid09 could be remodelled for AUSGeoid2020 using AGQG2009 and the co-located GNSS-AHD heights now available over Australia, we would also like to profit from the following new data: (i) an additional ∼280,000 land gravity observations; (ii) retracked satellite altimeter-derived marine gravity anomalies that include Jason-1 and CryoSat-2 data (Sandwell and Smith 2009;Garcia et al. 2014;Sandwell et al. 2013Sandwell et al. , 2014; (iii) the 1 × 1 (∼30 m) resolution DEM-H digital elevation model (Gallant et al. 2011) derived from the Shuttle Radar Topography Mission (STRM; Farr et al. 2007); and (iv) improved numerical integration routines  in the 1D-FFT implementation of Stokes's integral (Haagmans et al. 1993).
GNSS users of AUSGeoid09 have had to rely on a wholeof-region uncertainty estimate of 30-50 mm taken from the average of nationwide residuals versus GNSS-AHD ). This single uncertainty estimate is unrealistic in many regions. This is particularly the case in the denserpopulated Australian coastal zones, where AUSGeoid09 has had to rely on altimeter-derived gravity anomalies offshore, which are poor in the coastal zone (e.g., Claessens 2012), and where the ship-track gravity data are sparse and unreliable (Featherstone 2009). Therefore, we have estimated quasigeoid (height anomaly) errors as a function of location (cf. Pavlis and Saleh 2005;Huang and Véronneau 2013) by propagating uncertainties from EGM2008, land and altimeter-derived gravity anomalies, and terrain corrections in a remove-compute-restore (RCR) approach (Sect. 3).

Gravimetric quasigeoid computations 2.1 Preliminaries
All previous Australian geoid or quasigeoid models (cited in Featherstone et al. 2001Featherstone et al. , 2011 have focussed only on the Australian mainland and Tasmania. However, Australia administers several offshore territories and a 200-nauticalmile maritime boundary around each. Since the 1D-FFT requires a rectangular grid, we extended the computation area to 8 • S, 61 • S, 93 • E, 174 • E, resulting in a grid of approximately 15.5 million 1 × 1 nodes, roughly twice that used for AGQG2009. Figure 1 shows the coverage of land gravity observations used in AUSGeoid09, plus those added until May 2016 (used for AGQG2017). Figure 1 also identifies land gravity surveys that have been coordinated using GNSS. Geoscience Australia's online national gravity database (www.geoscience. gov.au/geophysical-data-delivery) uses the outdated AUS-Geoid98 (Featherstone et al. 2001) to transform GNSS ellipsoidal heights to the AHD. Therefore, we transformed all the GNSS-coordinated gravity surveys to the AHD using AUSGeoid09 for the subsequent computation of gravity anomalies.

Gravity data
We experimented with the effect on the quasigeoid of computing the gravity anomalies with respect to an a priori gravimetric quasigeoid model (cf. Amos and Featherstone 2009). The GNSS-levelling data (Sect. 2.4) were used to calculate a geometric component (cf. Brown et al. 2011) that models the difference between the a priori gravimetric quasigeoid and the AHD. This is dominated by a north-south tilt (cf. Featherstone and Filmer 2012) but with some regional distortions. The geometric component was converted to a free-air gravity term and modified-Stokes-integrated using the parameters in Sect. 2.6 to determine the effect on the quasigeoid. The maximum magnitude is less than 3 mm. This is due to (1) the high-pass filtering properties of a modified Stokes kernel over a small (0.5-degree) integration cap (cf. Vaníček and Featherstone 1998), and (2) the isotropic (azimuth independent) modified Stokes kernel is blind to a tilt in the gravity anomalies over the integration cap. That is, a tilt that manifests as positive values over one half of the cap and equal negative values over the other half cancel out and thus make no contribution to the quasigeoid.
For each record in the GA database, we calculated and applied normal gravity (Somigliana-Pizzetti formula, GRS80 ellipsoid (Moritz 1980)) and second-order free air and Bouguer plate corrections to the point gravity values to obtain simple planar Bouguer gravity anomalies. The tensioned spline (Smith and Wessel 1990) routine in the Generic Mapping Tools (GMT; Wessel et al. 2013) with a tension factor of 0.25 was used for the interpolation of point Bouguer gravity anomalies. We used the "reconstruction" technique (Featherstone and Kirby 2000) to generate Molodensky freeair gravity anomalies on the topography (Fig. 2a), using the 1 × 1 DEM-H data from GA. We added a 1 × 1 grid of planar terrain corrections (Fig. 2b;McCubbine et al. 2017) to these to give Faye anomalies and then block-averaged these onto a 1 × 1 grid in GMT for subsequent quasigeoid computations.
Offshore gravity data comprised marine gravity anomalies derived from multi-mission satellite altimetry (Fig. 2a). We chose the grav.img.23.1 model (http://topex.ucsd.edu/ marine_grav/mar_grav.html) (Sandwell and Smith 2009;Garcia et al. 2014;Sandwell et al. 2013Sandwell et al. , 2014 simply because it is accompanied by error estimates as a function of location, which are readily used in the error propagation (Sect. 3). We chose not to use ship-track gravity observations d the May 2016 GA land gravity database used in AGQG2017. Each observation is mapped as a single black dot because they contain long-wavelength errors that cannot be removed by a poorly conditioned crossover adjustment (Featherstone 2009). The land and marine gravity data were merged using exactly the same techniques as for AGQG2009 (Featherstone et al. 2011, Sect. 2.5).

Synthesis on the topography
For AGQG2009 , we used the har-monic_synth.f FORTRAN software (http://earth-info.nga. mil/GandG/wgs84/gravitymod/egm2008/index.html) to synthesise ellipsoidal block-averaged gravity anomalies and quasigeoid heights for use in the RCR technique. After receiving reports of poorer performance in the Great Dividing Range in south-eastern Australia (e.g., Rexer et al. 2013;Sussanna et al. 2016), we discovered that EGM2008 had been synthesised on the ellipsoid instead of on the topography as demanded by Molodensky theory. However, this poorer performance is more likely to be explained by erroneous Australian land gravity data that have now been corrected in the GA database (Sect. 2.6; Fig. 8).
Synthesising high-degree gravity field functionals on the topography dramatically slows the performance of har-monic_synth.f for a 1 × 1 grid over the extents defined in Sect. 2.1, because the numerically efficient recursion routines (Holmes and Featherstone 2002) cannot be implemented on an irregular surface. Therefore, we used the isGrafLab.m Matlab TM software Janák 2013, 2014) instead, which is capable of computing gravity field functionals from spherical harmonic coefficients on an irregular surface quite efficiently using the gradient approach (Hirt 2012) even for very high spherical harmonic degrees. This requires ellipsoidal heights of the topography, computed as follows.
We block-averaged the 1 × 1 DEM-H data onto a 1 × 1 grid (the same resolution as AGQG2017) to yield the computation point locations for the spherical harmonic synthesis on the topography. These block-averaged heights need to be converted to ellipsoidal heights for isGrafLab, thus requiring an iterative scheme. We first synthesised a 1 ×1 quasigeoid grid on the ellipsoid using isGrafLab and then added the blockaveraged DEM-H heights to this to obtain a first estimate of the ellipsoidal heights of the topography. We then recomputed the quasigeoid at this first estimate of the ellipsoidal height of the topography and added the DEM-H heights to gain a second estimate of the ellipsoidal height of the topography. This iteration converged quickly (<1 mm difference between the first and second iterations). Figure 3 shows the differences between calculating gravity anomalies and quasigeoid heights from EGM2008 at the ellipsoidal height of the topography versus on the surface of the GRS80 ellipsoid. The differences reach maximum magnitudes of 42.162 mGal and 0.363 m, respectively. We ran an experiment over the reasonably mountainous island of Tasmania (maximum elevation of ∼1600 m) to identify the magnitude of error that has resulted in AGQG2009 from the theoretically incorrect synthesising EGM2008 on the ellipsoid (Appendix A). It transpires that the RCR technique is rather forgiving of this consistently made mistake. Nevertheless, it remains theoretically superior to synthesise the EGM on the topography for quasigeoid determination, and this was done for AGQG2017.
The reason that this mistake in the EGM2008 synthesis was not detected earlier is the inability of the GNSS-levelling data to clearly discriminate between quasigeoid heights synthesised on the ellipsoid versus on the topography (Sect. 2.4). Using EGM2008 as an example and the 7224 GNSS-levelling points now available over the Australian mainland, the standard deviation of fit to the GNSS-levelling is ±0.188 m when EGM2008 is synthesised on the ellipsoid and ±0.187 m when it is synthesised on the topography. This small difference is also partly due to the generally low Australian elevations, where levelling traverses are less likely to sample the hills (highway effect), but also GNSS observations have tended to be made at more accessible locations at lower elevations. The mean AHD height of the 7224 points is only ∼214 m (Sect. 2.4; Table 1).

Block-averaged ellipsoidal gravity anomalies
The isGrafLab software only computes the spherical approximation of the gravity anomaly g s from the disturbing potential T , via the fundamental equation of physical geodesy where (r, θ, λ) are the geocentric radius, spherical polar co-latitude and longitude of the computation point on the  STD is standard deviation topography, respectively. A more rigorous definition is the ellipsoidal gravity anomaly, which uses partial derivatives with respect to the ellipsoidal normal instead of the geocentric radial and does not approximate the Somigliana-Pizzetti normal gravity field by a spherical normal gravity field (e.g., Jekeli 1981;Grafarend et al. 1999;Hipkin 2004;. For the ellipsoidal gravity anomaly on the topography g e where γ (h) is normal gravity on the topography (Eq. 11, below). Therefore, Eq.
(2) is reformulated so that other gravity field functionals already computed in isGrafLab on the topographic surface can be combined and thus converted to the ellipsoidal form. These comprise the disturbing potential T (r, θ, λ), the gravity disturbance in spherical approximation and the Helmert north-south vertical deflection, also in spherical approximation The chain rule of differentiation converts the partial derivatives of the disturbing potential with prior knowledge of the functionals that isGrafLab generates, viz Inserting Eqs.
(3), (4) and (5) in Eq. (2) gives The partial derivatives ∂r/dh and ∂θ/dh are derived in terms of geodetic coordinates (φ, λ, h) in (Claessens 2006, Chapter 4), where φ is the geodetic co-latitude. (Hirt and Claessens 2011, Eq. 4) use the derivatives at the surface of the ellipsoid, whereas we require them on the topography. In this more general case, they are (Claessens 2006, Chapter 4) and for normal gravity (noting the retained use of geodetic co-latitude) In Eqs. (6) through (12), a is the semi-major axis length of the normal ellipsoid, e is its first numerical eccentricity, f is its flattening, b its semi-minor axis length, and m is the geodetic parameter (ratio of centrifugal to gravitational accelerations at the equator), γ a is normal gravity at the equator and γ b is normal gravity at the poles. Numerical values for the GRS80 normal ellipsoid (Moritz 1980) were used throughout. All isGrafLab computations also included the degree-zero term and accounted for the difference in mass and semi-major axes between the Earth gravitational model (EGM) and GRS80 (cf. Smith 1998). No degree-one terms were included because they are not provided with the EGM2008 coefficients and we desire a geocentric quasigeoid model so these terms become inadmissible.
The above procedure uses the existing outputs of the isGrafLab software to compute ellipsoidal gravity anomalies via Eq. (6), thus avoiding having to recode isGrafLab. However, they are point values and not block-average values required for the RCR technique (cf. . In contrast to gravity anomalies in spherical approximation, block-averaged values of ellipsoidal gravity anomalies cannot be computed through Paul (1978) recurrence relations for integrals of associated Legendre functions . Therefore, we computed a 1 × 1 grid of point ellipsoidal gravity anomalies and then block-averaged these onto a 1 × 1 grid (Fig. 2a) to be subtracted from the land and altimeter-derived gravity anomalies. Figure 4 shows the difference between the block-averaged ellipsoidal and spherical gravity anomalies and their effect on the quasigeoid, computed using a degree-40 modified FEO kernel ) over a 0.5-degree integration radius. This parameter combination was determined from fits to GNSS-levelling data (Sect. 2.6). The largest effect on the quasigeoid is 11 mm, insignificant in relation to the expected error in AGQG2017 (Sect. 3), but possibly significant in the ongoing search for a "1 cm geoid" (cf. Sansò and Rummel 1997;Tscherning et al. 2001).

GNSS-levelling data
The GNSS-levelling data play three roles in the computation of AUSGeoid2020: they are used to (1) select the EGM to be used in the RCR technique; (2) determine the combination of parameter choices in the AGQG2017 gravimetric quasigeoid computations, namely degree of FEO kernel modification and integration cap radius; and (3) generate the geometric component to provide the more direct transformation of GNSS heights to AHD heights (cf. Brown et al. 2011). A near-nationwide dataset of GNSS-levelling points was provided by Australian State and Territory geodetic agencies (Sect. 2.6; Fig. 6d).
We used a high-resolution shoreline map (Wessel and Smith 1996) to identify GNSS-levelling sites on islands, which technically are separate vertical datums (e.g., Filmer and Featherstone 2012;Amos and Featherstone 2009): Tasmania (71 points), Lord Howe Island (1 point) and nearcoastal islands (185 points). These sites were excluded from the tests in roles 1 and 2 above, but will be included in role 3 because the AHD is defined as local mean sea level on these islands, despite the likely presence of vertical offsets in the sea level connections to the mainland (Lord Howe and Tasmania). We detected around 20 consistently identified outliers from comparisons with multiple quasigeoid models (Tables 3, 4), leaving 7224 GNSS-levelling sites for use in roles 1 and 2.
The AHD (Mainland) was fixed to mean sea level (MSL = zero;epoch 1966-1968 at 30 mainland tide gauges in a two-stage least squares adjustment (LSA) in 1971 (Roelse et al. 1971) and is known to contain a north-south tilt (Featherstone and Filmer 2012) and regional distortions. Therefore, we produced heights from three readjustments of the Australian National Levelling Network (ANLN; Filmer et al. 2014), which is an updated version of the levelling data used in the 1971 AHD adjustment: • LSA 1 (MC_NOC) is a minimally constrained (MC) LSA holding only the Johnston origin station in central Australia fixed arbitrarily to its AHD height and using  (1961) normal-orthometric correction (NOC) as used in the AHD. The Australian height systems are reviewed by Featherstone and Kuhn (2006). • LSA 2 (MC_NC) is a minimally constrained LSA holding only the Johnston origin station fixed arbitrarily to its AHD height and using the normal correction (NC). The NC system is theoretically consistent with the quasigeoid. Gravity values at ANLN benchmarks required for the NC were derived from EGM2008, as in Filmer et al. (2010). • Adjustment 3 (CARS) is a constrained LSA, where the 30 mainland tide gauges that were held fixed to MSL in the AHD LSA are now held fixed to the MSL corrected for mean dynamic topography (MDT; aka sea surface topography) from the CSIRO's atlas of regional seas (CARS2009) climatology Ridgway et al. 2002). The MSL was observed for three years in the 1960s and the MDT for the past 50 years with bias to more recent years . This extra constraint to MDT-corrected MSL, which is theoretically aligned with the quasi/geoid, alleviates large (spatial) scale distortions in the ANLN that propagate into the MC LSAs. Normal corrections were applied as per LSA 2.
The idea behind these readjustments is to avoid the tilt in the AHD contaminating the parameter choices (Sects. 2.5, 2.6). However, they cannot account for regional distortions and other errors in the ANLN (cf. Filmer et al. 2014). Table 1 (columns 6-8) shows that there are some quite substantial differences among levelled heights from these different LSAs, due mainly to the MSL constraints used in the AHD at mul-tiple tide gauges . Table 1 (column 9) also shows that the NOC is a reasonable approximation of the NC (cf. Filmer et al. 2010) in relation to the expected precision of the levelling data (Table 2). Another consideration is the inherent precision of the GNSS-levelling data and thus its ability to discriminate among different quasigeoid model solutions. The GNSSlevelling data file provided by Australian State and Territory geodetic agencies is accompanied by error estimates for the ellipsoidal heights (from the GNSS data processing) and levelled heights (from the original 1971 LSA of the AHD). As these are independent error estimates, the precision of the GNSS-levelling-derived (geometric) quasigeoid can be calculated by taking the square root of the sum of the squares (Table 2, column 3). Taking the mean of the STD of all 7224 observations as a proxy for their overall precision, the ability of the GNSS-levelling data to discriminate among quasigeoid models (Sects. 2.5, 2.6) is only ±40 mm, considerably larger than that needed to properly identify a "1 cm geoid".

Choice of EGM for the RCR technique
Two degree-2190 EGMs have been released since EGM2008: EIGEN-6C4 (Förste et al. 2015) and GECO (Gilardoni et al. 2015), both of which include GRACE (Tapley et al. 2004) and GOCE (Drinkwater et al. 2003) data. [EGM2008 only includes GRACE data]. The standard deviations in Table 3 show that there is little difference among these EGMs when compared to the 7224 GNSS-levelling data (all EGM values were synthesised on the topography (cf. Sect. 2.1)), with EGM2008 providing marginally lower STDs. However, bearing in mind the inability of the GNSS-levelling data to discriminate these millimetre-scale differences among quasigeoid models (Sect. 2.4), the choice of EGM2008 is somewhat arbitrary.

Modified Stokes parameter sweeps
The 7224 GPS-levelling sites were also used to choose the optimal combination of Featherstone et al. (1998) FEO kernel modification degree (M) and spherical integration cap radius (ψ 0 ) in the 1D-FFT modified Stokes integration. To add some robustness to the parameter choice, these parameter sweeps were compared with all four variants of the levelling data (AHD, MC_NOC, MC_MC, CARS), without and with solving for a bias and tilted plane (Fig. 5). The tilted plane was included because the minimally constrained adjustments still exhibit a (smaller than AHD) tilt due to systematic levelling errors, and the bias accounts for a combination of the poorly determined zero-degree term and mean offset of the AHD from the 'true' geoid potential, W 0 . From Fig. 5, and as was found for AGQG2009 ) and AUSGeoid98 (Featherstone et al. 2001), the unmodified spherical Stokes kernel is inappropriate because it does not filter out long-wavelength errors from the land and altimeter data (cf. Vaníček and Featherstone 1998). The optimal integration radius is clearly ψ 0 = 0.5 degrees from all comparisons in Fig. 5. The degree of modifi-cation ("FEO-mod" in the Fig. 5 legends) is indistinguishable among degrees ranging from M = 40 to M = 140, so we chose M = 40 as there is no empirical evidence to choose a higher degree of modification. Also, recall that the ability of the GNSS-levelling data to discriminate among quasigeoid solutions is only ±40 mm (Sect. 2.4), so there is quite some leeway in these parameter choices. Figure 5 (bottom panel) demonstrates the relative quality of the different levelling LSAs for quasigeoid testing. Without solving for a bias and tilted plane, CARS has a STD <±0.10 m at 0.5 degree cap radius, whereas the MC_NC and MC_NOC are both >±0.15 m and AHD >±0.18 m. When solving for a bias and tilted plane, the AHD STD decreases to ±0.10 m, while CARS decreases to ±0.09 m, indicating that CARS has already removed most of the AHD's MDTinduced tilt (cf. Featherstone and Filmer 2012). MC_NOC and MC_NC STD decrease to ±0.12 m, but reflect that regional distortions rather than a systematic north-south tilt contribute to errors in these LSAs. This indicates that the CARS-constrained ANLN LSA (cf. Filmer et al. 2014) is the most suitable dataset currently available for testing geoid models in Australia.
We also attempted to determine the optimal parameter combination using a set of 1080 historical vertical deflections observed in the 1950s and 1960s. The precision of these observations is unknown, but they are probably only ±1 . Therefore, they were unable to neither confirm nor   Figure 6a shows the residual terrain-corrected free-air gravity anomalies after the removal of the block-averaged EGM2008 ellipsoidal gravity anomalies synthesised on the topography (Sect. 2.3). Figure 6b shows residual quasigeoid undulations from the M = 40 FEO-modified Stokes integration over ψ 0 = 0.5 degrees. Figure 6c shows the AGQG2017 gravimetric quasigeoid model after restoration of the EGM2008 quasigeoid heights synthesised on the topography. Figure 6d shows the quasigeoid height differences at the 7224 GNSS-AHD stations. The north-south tilt and regional distortions in the AHD necessitate the geometric component (cf. Brown et al. 2011) for the transformation of GNSS ellipsoidal heights to the AHD, until the time that the AHD is redefined. Table 4 shows that, when assessing the gravimetric quasigeoid models against the Australian GNSS-levelling data, there is only a very marginal improvement upon EGM2008. The first qualifier is the uncertainties in the GNSS-levelling data, so a millimetre-scale change among gravimetric models is insignificant when compared to the ±40 mm mean STD of the GNSS-AHD data ( Table 2). The second qualifier is the spatial distribution of the GNSS-levelling data, with relatively few stations at high elevations (only 141 of 7224 greater than 1 km and 806 greater than 500 m) and sparse GNSS-levelling coverage in central and northern Australia, where most new land gravity data have been collected (cf. Figs. 6d, 1b). We considered thinning the GNSS-levelling data, but this became subjective because the standard errors from the LSA do not describe the ∼0.5 m systematic errors in the AHD. Nevertheless, AGQG2017 has used newer and higher-resolution (notably the DEM) data sources than its predecessors. Figure 7 shows the differences between the gravity anomalies used in AGQG2017 and AGQG2009 and the difference between AGQG2017 and AGQG2009 quasigeoid heights. Notable features are: (1) the differences in gravity anomalies offshore between grav.img.23.1 (Sandwell and Smith 2009;Garcia et al. 2014;Sandwell et al. 2013Sandwell et al. , 2014 and DTU08GRAV (Andersen et al. 2010), particularly in the coastal zone (Fig. 7a); and (2) differences in quasigeoid heights in the mountains along the Great Dividing Range in south-eastern Australia (Fig. 7b).
The differences in Fig. 7b are due to a combination of the synthesis of EGM2008 (Sect. 2.3), and the use of the 1 × 1 DEM-H model for both the reconstruction technique (Featherstone et al. 2001) and computation of terrain corrections (McCubbine et al. 2017). In addition, Claessens et al. (2009) identified the possibility of erroneous Australian gravity data in this region, which have since been corrected in the GA gravity database (Fig. 8). The poorer performance of AUSGeoid09 in this region (cf. Rexer et al. 2013;Sussanna et al. 2016) is therefore attributed to a superposition of all the above effects, but the erroneous land gravity data appears to be the major contributor.
In order to quantify the improvement from these corrected land gravity data, we extracted a subset of 173 GNSSlevelling stations in the area bound by 36 .031 m. This shows that the corrected land gravity data have removed a substantial bias from the regional quasigeoid model.

Preliminaries
No previous Australian geoid or quasigeoid model has been accompanied by an error grid, instead relying on an overall error estimate taken from the mean of residuals to GNSS-levelling data (Sect. 1). Notwithstanding least squares collocation (LSC), few previous studies have propagated uncertainties through Stokes's integral (e.g., Sideris and She 1995;Pavlis and Saleh 2005;Huang et al. 2007). LSC was not used for the Australian gravimetric quasigeoid computations simply because of the sheer size and extent of the Australian datasets (Sect. 2.1).
In the most general sense, the terms in Stokes's integral are squared, the gravity anomaly error variances are used instead of gravity anomalies in the integration, and then the square root of the result is taken to give the quasigeoid error standard deviation. This comes from the standard formulas for propagation of independent error variances (i.e., var(x + y) = var(x)+var(y) and var(a.x) = a 2 var(x) for a constant a) applied to the discretised Stokes integral where κ = r/4πγ e for the geocentric radius r of the computation point, ζ i is the quasigeoid height at the computation point i, g j is the mean gravity anomaly at the roving point j, S(ψ i, j ) is the Stokes kernel function dependent on the angular distance ψ i, j between the computation and roving points, and d is the area of the roving point's cell size.
Applying the standard formulas for the propagation of the gravity anomaly error variances, σ 2 ( g j ), into quasigeoid error variances σ 2 (ζ i ) Equation (14) is given in analytical form by Pavlis and Saleh (2005), but without a derivation and noting their omission of the d 2 term, as which was used to produce the quasigeoid error grid that accompanies EGM2008 (Pavlis et al. 2012(Pavlis et al. , 2013. It is important to note that Eqs. (14) and (15) assume that gravity anomaly errors in distinct grid cells are independent, since no covariance terms are considered. Sideris and She (1995) provide error propagation formulas for the 1D-FFT (Haagmans et al. 1993) when using the RCR technique over a region. Huang et al. (2007) provide a spectral form when using the RCR technique and a modified integration kernel over a spherical integration cap. Common to both these approaches is the need to use a satellite-only EGM so that errors in the terrestrial gravity data are independent and thus combined without the need to consider correlations. However, when using a combined EGM such as EGM2008, this assumption breaks down. Therefore, we have devised an alternative approach to error propagation (Sect. 3.3).

Terrestrial gravity error grid
The GA gravity database also contains estimates of the precision of the gravity observations (σ g ), their horizontal locations (σ H ) and gravimeter heights (σ φ ). All are assumed to be standard deviations (i.e., one sigma). A corresponding variance for each point Bouguer gravity anomaly was estimated by propagating the effect of the positioning errors into the gravimetric corrections and combining these with the square of gravity value error estimate where 0.3086 mGal/m is the spherical approximation of the free-air gravity gradient, 0.1119 mGal/m is the Bouguer plate gravity gradient for a topographic bulk density of 2670 kg/m 3 , and 0.000812 sin 2φ is a truncated Chebychev approximation of normal gravity (Moritz 1980). The scattered σ 2 ( g SB ) data were block-averaged and then gridded using GMT at the same 1 × 1 resolution as the land gravity data used in AGQG2017.
In the reconstruction technique (Featherstone and Kirby 2000), the mean free-air gravity anomaly on the topography is recovered by block-averaging the DEM, multiplying by the Bouguer plate term (0.1119 mGal/m), and adding it to the simple Bouguer gravity anomaly grid. The standard deviation of the DEM-H model has been estimated to be ±7.66 m (Gallant et al. 2011;McCubbine et al. 2017). The variance of the 1 ×1 DEM-H heights when block-averaged to 1 ×1 is where N is the number of sums run over all points in each 1 × 1 cell. The DEM-H errors are assumed to be correlated out to ∼100 m (Rodríguez et al. 2006;Gallant et al. 2011). Therefore, the covariance term in Eq. (17) has been modelled with an exponential function of the form where l is the distance between H i and H j , and r = 100 m is the radius beyond which the correlation between the DEM-H errors drops below 95%. This is the same covariance model as used for the removal of correlated DEM-H errors in the gravimetric terrain corrections (McCubbine et al. 2017). Finally, the 1 × 1 terrain correction error variances (Fig. 9a) were added to give the land gravity anomaly error variances Formally, Eq. (19) should also include a covariance term of the terrain correction errors and Bouguer slab 'reconstruction' errors since both are computed from the same DEM. This term is difficult to calculate, so the terrain correction errors and Bouguer slab restoration errors have been assumed to be independent for simplicity. In offshore regions, the error grid that accompanies the marine gravity anomalies from the grav_error.img.23.1 model was concatenated with the land gravity anomaly errors (Fig. 9b) after masking out Australian land areas using the high-resolution shoreline (Wessel and Smith 1996). The large errors on non-Australian land (e.g., New Zealand, Indonesia and Papua New Guinea) are values provided with the grav_error.img.23.1 grid.

Quasigeoid error propagation
In the RCR technique with an unmodified Stokes kernel, the quasigeoid error propagation is where σ 2 (ζ EGM ), σ 2 ( g EGM ) and σ 2 ( g) are the error variances of the EGM quasigeoid height, EGM gravity anomaly and gridded terrestrial gravity data errors, respectively. Importantly, they are not to be confused with the degree or error degree variances calculated from the EGM spher-ical harmonic coefficients or the variance of the anomalous gravity field itself. Haagmans and Gelderen (1991) describe methods for the σ 2 (ζ EGM ) error propagation when the full variance-covariance (VCV) matrix of the EGM geopotential coefficients is available, but this is rarely the case. A practical limitation is the sheer size of the arrays needed to store the full VCV, particularly for a degree-2190 model. Instead, the error degree variances (diagonal of the VCV matrix) are often used to approximate σ 2 (ζ EGM ). In our error propagation, however, we use the error grids of σ 2 (ζ EGM ) and σ 2 ( g EGM ) that accompany EGM2008, and not the error degree variances of a satellite-only EGM (as used by, e.g., Sideris and She 1995;Huang et al. 2007;Huang and Véronneau 2013) and assume independence between error values at grid nodes. For a global integration with an unmodified Stokes kernel, Eq. (20) degenerates to Eq. (15) because the first, third, fourth, fifth and sixth terms on the right-hand side cancel, and the error propagation is solely from the terrestrial gravity data variances σ 2 ( g) (cf. Pavlis and Saleh 2005). However, with a limited integration domain 0 (in our case, a spherical cap) and a modified integration kernel (in our case, Featherstone et al. 1998), the high-pass filtering properties of the integral terms are lessened (cf. Vaníček and Featherstone 1998) and these terms no longer cancel. We therefore make the following two theoretical assumptions for simplicity: (1) the geopotential coefficients of different degrees are independent, and (2) the integration cap behaves as a perfect high-pass filter. Under these assumptions For any modified kernel S and integration truncated to a spherical cap 0 beyond the spherical cap radius ψ 0 and the above two assumptions, Eq. (20) becomes a RCR-like approach to the error propagation The first term on the right-hand side of Eq. (23) is the square of quasigeoid errors taken directly from the EGM2008 website (Fig. 10b). The second and third terms on the right-hand side were evaluated using the 1D-FFT with the same parameter choices as were used to compute AGQG2017 (Sect. 2.6), but now using the squared terms. Adapted from Sideris and She (1995), in 1D-FFT form, the second and third terms in Eq. (23) are, respectively Equation (24) uses the errors propagated from the Australian gravity observations, DEM and terrain correction errors (Sect. 3.2; Fig. 9b). Equation (25) uses the square of the gravity anomaly errors taken directly from the EGM2008 website (Fig. 10a). The EGM2008 quasigeoid height and gravity anomaly errors are provided on 5 ×5 grids, and AGQG2017 is on a 1 ×1 grid. Therefore, the EGM2008 quasigeoid errors were re-sampled to 1 × 1 using GMT's tensioned splines routine (Smith and Wessel 1990) with a tension factor of 0.25. In order to avoid underestimating the error propagation from a smaller integration step, we propagated the 5 × 5 values through Eq. (25) and then re-sampled to 1 × 1 using GMT exactly as above. Figure 11 shows the results of the error propagation for EGM2008 and the terrestrial data, with the latter replacing the former to produce Fig. 12. This novel approach to quasigeoid error propagation allows us not only to avoid the need to use a satellite-only EGM to avoid correlations, but also to avoid the need to use Fig. 11 1D-FFT evaluation of quasigeoid errors from a the EGM2008 gravity anomaly error grid re-sampled to 1 × 1 (Fig. 10a), and b the 1 × 1 grid of land and altimeter gravity anomaly errors (Fig. 9b). Units in m Fig. 12 a The Australian gravimetric quasigeoid model 2017 AGQG2017, and b its associated location-specific error map/chart. Units in m the error degree variances of the geopotential coefficients that only give a single global estimate of the EGM errors. We essentially implement a RCR on the errors by replacing the gravity anomaly errors used in EGM2008 by our own regional estimates on a denser grid, thus reducing the omission and the commission errors in EGM2008.

External validations of the error grid
Our error propagation suggests quasigeoid errors of 50-60 mm across most of mainland Australia, increasing to ∼100 mm in regions of steep horizontal gravity gradients or in the mountains (Fig. 12b). To gauge whether our error estimates are realistic, we compare our values to those reported for Canada by Huang and Véronneau (2013). The maximum AHD elevation is ∼2230 m and the formally propagated quasigeoid error is ∼100 mm in the vicinity of this mountainous region (south-eastern Australia). The maximum elevation in Canada is 5959 m and Huang and Véronneau (2013, p. 785) report geoid error values of 300 mm in the Canadian Rocky mountains. Making the gross assumption that the errors scale linearly with height, these error estimates appear to be reasonably consistent (i.e., 5959 × 100/2230 = 267 mm).
Another check comes from a comparison of the broadscale quasigeoid error estimates deduced from variance component estimation (Filmer et al. 2014). Their error estimate was ±79 mm for ellipsoidal height minus AGQG2009 using 277 GNSS points. The average internal error estimate of the GNSS ellipsoidal heights from the Bernese processing for the dataset used for Filmer et al. (2014) was ±26 mm (scaled by a factor of 10 to provide a more realistic error estimate, e.g., Rothacher 2002), so the quasigeoid component of this error is ∼±75 mm based on variance propagation.
The final check comes from a comparison of our quasigeoid error estimates with the formal errors in the GNSSlevelling data (Sect. 2.4) on a point-by-point basis. We used bicubic interpolation of the error grid to the locations of the GNSS-levelling stations and compared them to the formal error at each of the 7224 points. The mean absolute difference is ∼±25 mm, with the gravimetrically propagated errors being larger than the formal errors in the GNSS-levelling data. The apparently better GNSS-levelling errors are partly due to the disproportionately large number of GNSS points in central-eastern and south-western Australia (Fig. 6d), where the formal errors are smaller. Figure 10a shows a band of relatively low (<1 mGal) marine gravity anomaly error estimates around the Australian coast for EGM2008, which is unique to Australia. When these are propagated through Eq. (25) for a 0.5-degree cap and M = 40 FEO modification, the corresponding quasigeoid error estimates are less than ±10 mm (Fig. 11a). This is contrary to the knowledge that altimeter-derived gravity anomalies are of lower accuracy in the coastal zone (e.g., Deng et al. 2002;Volkov et al. 2007); also see Fig. 11b. Therefore, we have propagated the grav_error.img.23.1 grid of altimeter-derived gravity anomaly errors so as to increase the AGQG2017 error estimates (∼±50 mm) above those of EGM2008 in the coastal zone, expecting this to be more realistic. In other areas, the EGM2008 errors are reduced because the higher-resolution grid reduces the omission error (and part of the commission error) when using the formally propagated terrestrial gravity data errors (cf. Figs. 10b, 12b).

Concluding remarks
We have described the computation of the gravimetric quasigeoid AGQG2017, as well as the computation of location-dependent error estimates that have been propagated from the combined uncertainties in EGM2008, land and altimeter-derived gravity anomalies, and terrain corrections. AGQG2017 has included an additional ∼280,000 land gravity observations over AGQG2009 (1,764,351 in total), re-tracked satellite altimeter-derived marine gravity anomalies that include Jason-1 and CryoSat-2 data, the 1 × 1 (∼30 m) resolution DEM-H digital elevation model for reconstruction of mean free-air gravity anomalies on the topography, computation of terrain corrections that include formally propagated error estimates, and improved numerical integration routines.
Computational refinements were made to the synthesis of EGM gravity field functionals on the topography as is demanded by Molodensky theory (Sect. 2.3). Three separate outputs of the isGrafLab software were combined to generate block-averaged ellipsoidal gravity anomalies, also on the topography. Four least squares readjustments of the Australian levelling data (Sect. 2.4) were used to more robustly select the EGM2008 model for the RCR technique (Sect. 2.5) and parameter sweeps for the degree of kernel modification and integration cap radius (Sect. 2.6). These indicated that a CARS2009 MDT-constrained ANLN is the best currently available levelling adjustment for testing quasigeoid models in Australia. However, the mean error of the GNSS-levelling data is ∼40 mm, making it a less useful tool to discriminate among gravimetric quasigeoid models. Historical deflections of the vertical were not able to assist during the parameter sweeps because of their low precision.
We present a new and alternative technique for the propagation of location-dependent errors without the need to rely on global estimates of the error degree variance of a satelliteonly EGM. Instead, the error grids of quasigeoid heights and gravity anomalies were taken from EGM2008 and used in a RCR-type procedure, where errors in the land gravity anomalies, altimeter gravity anomalies and terrain corrections were propagated through the same modified Stokes integral used to compute AGQG2017 (Sect. 3). The gravimetric quasigeoid errors (one sigma) are 50-60 mm across most of the Australian landmass, increasing to ∼100 mm in regions of steep horizontal gravity gradients or the mountains (Fig. 12b). Our error estimates were validated externally using three separate approaches, indicating them to be realistic.
Finally, we question the band of low EGM2008 gravity anomaly error estimates around the Australian coast, instead increasing the AGQG2017 quasigeoid error estimate based on propagation of errors from the grav_error.img.23.1 grid. On land away from the coasts, the EGM2008 error is reduced through the use of Australian land gravity and DEM data on a high-resolution 1 × 1 grid.
Acknowledgements The current work has been supported financially by the Cooperative Research Centre for Spatial Information, whose activities are funded by the Business Cooperative Research Centres Programme, and by Geoscience Australia. Previous works (e.g., software and algorithm development) that have been used here have been supported financially by the Australian Research Council (ARC) over the past two decades through Grant Numbers A39938040, A00001127, DP0211827 and DP0663020. We thank Scripps Institution of Oceanography (University of California), the US National Oceanographic and Atmospheric Administration and the US National Geospatial-Intelligence Agency for permission to use the marine gravity anomalies from Sandwell et al. (2014). All maps and charts in this paper were produced using GMT  in the Lambert conformal conic projection with two standard parallels or the Mercator projection ( Fig. 8; Appendix A). Nicholas Brown publishes this paper with the permission of the CEO, Geoscience Australia. Finally, we thank the three anonymous reviewers for their constructive and complimentary critiques on the first-submitted version of this manuscript Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: The effect of synthesising an EGM on the topography and on the ellipsoid for use in remove-compute-restore quasigeoid computations
To quantify the effect of incorrectly synthesising an EGM on the surface of the ellipsoid versus correctly synthesising it at the ellipsoidal height of the topography for regional quasigeoid computations in the remove compute-restore (RCR) scheme, an experiment has been conducted over Tasmania (140 • E to 150 • E and 39 • S to 45 • S) with the following stages Remove 1. Synthesise two 1 ×1 grids of EGM2008 gravity anomalies, one on the surface of the ellipsoid and the other at the ellipsoidal height of the topography; 2. Compute two 1 × 1 grids of residual gravity anomalies, the first being a Faye anomaly grid minus the EGM2008 gravity anomaly synthesised at the surface of the ellipsoid, and the second being a Faye anomaly grid minus the EGM2008 gravity anomaly synthesised at the ellipsoidal height of the topography; Compute 3. Compute two 1 × 1 grids of residual quasigeoid heights for both the above residual gravity anomaly grids. So as to replicate the approach taken in Featherstone et al. (2011), the computation used the Featherstone et al. (1998) kernel with modification degree M = 40 and a spherical cap radius of 1 degree; Restore 4. Synthesise two 1 × 1 grids of EGM2008 quasigeoid heights, one on the surface of the ellipsoid and the other at the ellipsoidal height of the topography. The descriptive statistics of the differences are min: −0.004 m, max: 0.180 m, mean: 0.002 m, STD: ±0.010 m; 5. Add these EGM2008 quasigeoid heights to the corresponding residual quasigeoid heights in a consistent manner. That is, the residual quasigeoid from residual gravity anomalies on the ellipsoid are restored to the EGM2008 quasigeoid heights synthesised on the ellipsoid, likewise for the ellipsoidal height of the topography. The difference between these two regional RCR quasigeoid models is shown in Fig. 13. The descriptive statistics of the differences are min: −0.060 m, max: 0.040 m, mean: −0.001 m, STD: ±0.005 m. This is on average 50% less than the difference between the two EGM2008 quasigeoid syntheses (stage 4 above). As such, the consistently incorrectly synthesised gravity anomalies and quasigeoid heights have a marginal impact on the quasigeoid heights from the RCR scheme.