Australian quasigeoid modelling: review, current status and future plans

Key Points: a) We review the history of quasi/geoid modelling in Australia, from the first model produced in 1967 to the most recent AGQG2017. b) The results of more than 10,000 trial quasigeoid solutions are presented, computed to assess whether methods used to produce AGQG2017 could be improved. c) Details of a new Australian vertical reference system, and future planned efforts to improve the underlying quasigeoid model are given.

In recent years, a few tens of new and more precise vertical deflections have been collected using contemporary digital equipment (CCD star cameras integrated with GPS) and used to assess gravimetric-only quasigeoid models in Perth, Western Australia, where the quasigeoid gradient is very steep (Schack et al., 2018). The historical vertical deflection data, used by Fischer-Slutksy and Fryer, are also still of some relevance. They have been used to evaluate geoid and quasigeoid models determined entirely from gravity data (e.g., Featherstone and Morgan, 2007;Featherstone et al. 2018a) and experiments have been conducted into the a posteriori fitting of a gravimetriconly geoid model to these historical vertical deflections (Featherstone and Lichti, 2009). The estimated accuracy of these historic vertical deflections is variable and not available on a perstation basis, with estimates by different authors ranging from ±0.22 arc-sec to ±0.76 arc-sec, which are cited in .
In 1971, the Australian Height Datum (AHD, Roelse et al., 1971) was realised as the first and only vertical datum for the Australian mainland from a least squares adjustment of a continent-wide set of differential levelling observations, constrained to zero height at 30 tide gauge observations of mean sea level distributed around the coast. In 1983, differential levelling in Tasmania was least squares adjusted by holding two tide gauges fixed. Though technically two separate vertical datums, both are called AHD. Several attempts have been made to unify the AHD between the Australian mainland and Tasmania (Rizos et al., 1991;Featherstone, 2000;, but none of these have led to a revision of AHD heights in Tasmania or the mainland. Gravity observations were not collected along all the levelling traverses at that time, so GRS67 normal gravity was used instead via truncated versions of the equations in Rapp (1961) (Roelse et al., 1971). This provided a normal-orthometric height (cf. Featherstone and Kuhn, 2016) of H=566.3 m at the Johnston Origin, which was coupled with the ANS ellipsoidal height of H=571.2 m (legally < A c c e p t e d M a n u s c r i p t > fixed at that time) to provide a new reference value for the geoid at the Johnston Origin of +4.9 m.
This was subsequently used to re-align Fryer's astrogeodetic geoid model with a shift of +10.9 m.
At around the same time, Mather (1969) produced a gravimetric co-geoid model which, unlike the models produced by Fischer-Slutsky and Fryer, was referenced to the geocentric GRS67 ellipsoid. Rather than using astrogeodetic deflections, Mather used free-air gravity anomalies determined from a nation-wide set of gravity observations that had been collected primarily for resource exploration (Fraser et al., 1976;cf. Gilliland, 1987). A remove-compute-restore technique was employed using a combined global geopotential model, determined from both satellite and terrestrial gravity data. Residual free air gravity anomalies were then transformed in to geoid undulations using Stokes's integration with an unmodified kernel. After a translation for the geocentre between the ANS and GRS67 (Mather, 1970), Mather's model was compared to Fischer-Slutsky's model and showed differences exceeded 10 m and had an RMS of ±3 m (Kearsley and Govind, 1992). The differences between the models were attributed to inadequacies (long wavelength control and gravity data density) in the gravity data used by Mather and the limited number of vertical deflections used by Fischer-Slutsky (Kearsley and Govind, 1992). Mather's model was technically a free-air co-geoid because direct and indirect topographic effects (e.g., Heiskanen and Moritz, 1967) were not considered.
Shortly after, Grushinsky et al. (1971) produced a gravimetric quasigeoid model by evaluating Stokes integral, using direct numerical integration, in a remove compute restore scheme. As this work was done in Russia, where the Molodensky theory was prevalent, the gravity anomaly data used in the computation were a zero order approximation to the Molodensky gravity anomaly (free air anomalies) computed from gridded Bouguer gravity anomalies (supplied by the Bureau of Mineral Resources) and a topographic map of Australia. They estimated the accuracy of the gridded free air gravity anomalies to be between 5 and 10 mGal. Their model was (i) compared to that < A c c e p t e d M a n u s c r i p t > produced by Mather (1961) which they report had differences with an RMS of ±5 m and (ii) the astrogeodetic model by Fischer and Slutsky (1967) -and they report differences with an RMS ±3.7 at a 51 points.
In the mid-to-late-1980s, Global Positioning System (GPS) technology started to become prevalent for geodetic surveying through the use of carrier-phase measurements, something that GPS was never designed to do. Therefore, in contrast to the earlier need for a regional model of the geoid with respect to the ANS, the GPS user community needed a model to determine physically meaningful heights from GPS-derived ellipsoidal heights that refer to a geocentric ellipsoid, nominally WGS84. In recognition of this, academics at the University of New South Wales and the University of South Australia independently computed gravimetric co-geoid models. Co-geoid models were produced by Kearsley (1988aKearsley ( , 1988b using the 1980 release of the Australian Geological Survey Organisation's national gravity database, which comprised ~360,000 terrestrial gravity measurements the Australian mainland and ~75,000 measurements offshore. This and subsequent co-geoid models used "ring integration" (Kearsley, 1985(Kearsley, , 1986Kearsley et al., 1998) to evaluate Stokes's integral where the integration radius from the computation point was limited to 0.5 degrees and with an unmodified kernel. Kirby et al. (1997) showed that the ring integration approach gives results commensurate with quadrature integration of rectangular blocks. Gilliland (1989) used same release of the gravity database to calculate the free air anomalies on a 0.1° by 0.1° grid. The Kearsley and Gilliand models employed a remove-compute-restore technique.
Residual free air gravity anomalies were then transformed to residual geoid heights on a 1° by 1° grid using Stokes's integration, evaluated by direct numerical integration over 3° by 3° blocks, with an unmodified kernel. Compared to Mather's model, Gilliland's model had significant differences, including a large N/S trend, which were attributed to the improved gravity data quality and < A c c e p t e d M a n u s c r i p t > coverage (Kearsley and Govind, 1992). Both Kearsley's and Gilliland's models were also free-air cogeoids.
In 1991, the then Australian Surveying and Land Information Group (AUSLIG) formally recognised the need for a national geocentric gravimetric geoid model for use in conjunction with GPS. They subsequently used a gravimetric geoid model called AUSGeoid91 (Kearsley and Govind, 1992). This model was referenced to the WGS84 ellipsoid and was again computed with the 1980 release of the national gravity database using an unmodified Stokes integral (evaluated by bruteforce numerical ring integration) with the remove-compute-restore technique. For this calculation, global model OSU89A (Rapp and Pavlis, 1990) was used which contained spherical harmonic coefficients up to degree and order 360 and the integration was performed within a 0.5 degree radius about each computation point -this was so that the terrestrial gravity data only contributed to wavelengths shorter than those in the global model. Two years later, AUSGeoid91 was superseded by AUSGeoid93 (Kearsley and Steed, 1995), which was computed using the same computational technique but with a refined global model OSU91A to 360 (Rapp et al., 1991).
AUSGeoid93 was then used as the national standard for converting GPS measurements to AHD heights for the following five years (Steed and Holtznagel, 1994). AUSGeoid91 and AUSGeoid93 were also free-air co-geoids.
Throughout the mid-1990s and early 2000s, a number of experimental studies were undertaken around refining the computational techniques in the Australian setting. In a PhD dissertation, Zhang (1997) looked at the merits of evaluating Stokes's integral in the Fourier domain. Featherstone and Kirby (2000) detailed a method to utilise a digital elevation model to reduce spatial aliasing when interpolating pointwise gravity anomalies onto a regular grid. Various gridding/interpolation methods and their impact on geoid computations was assessed by Goos et al. (2003) and Zhang and Featherstone (2004). Some local geoid models were also produced for < A c c e p t e d M a n u s c r i p t > example by Freund et al. (1997) just for the Australian Capital Territory, a series of geoid models over Western Australia were produced by  to explore the impact of modifying Stokes kernel, and Vella and Featherstone (1999) computed a quasigeoid over Tasmania using a modified Stokes kernel with the 1D discrete Fourier technique (Haagmans et al. 1993). Alternative computation methods were also experimented with, particularly in the Perth region. Claessens et al. (2000) investigated point mass modelling of the geoid, but the results were contaminated by erroneous marine gravity data (cf. Featherstone, 2009). Darbeheshti andFeatherstone (2009, 2010) investigated non-stationary least squares collocation (LSC), borrowing from spatial statistics to form covariance functions in a region where the gravity field is not a stationary stochastic variable. Soltanopur et al. (2006) investigated non-orthogonal second generation wavelets as an alternative to LSC for creating a hybrid geoid by combining gravimetric and geometric geoid models. However, none of these local models were formally released or officially adopted in the AUSGeoid models, but the investigations informed the computation processes.

Contemporary quasigeoid modelling
In 1998, AUSLIG in collaboration with academic sector computed a new geoid model, AUSGeoid98 . It also brought together the three most prominent gravimetric geoid modellers in Australia at that time: Kearsley, Gilliland and Featherstone, who jointly attracted funding from the Australian Research Council (ARC) for the project. The methods were provided, free-of-charge, to AUSLIG who ran the final computations and released and administered the model. This practice has continued ever since, and allows a single point of contact for the model(s).
Compared to previous Australian geoid models, the computation of AUSGeoid98 was (i) more theoretically rigorous, (ii) contained substantially more gravity data, gridded at a comparatively higher resolution, and (iii) used the 1D Fourier method to evaluate Stokes's integral in conjunction < A c c e p t e d M a n u s c r i p t > with a deterministically modified kernel  and parameters specifically chosen to fit GPS-levelling ground truth data. Moreover, it included gravimetric terrain corrections and their indirect effects, so was a geoid and not a co-geoid. AUSGeoid98 geoid undulations were specifically referenced to the GRS80 ellipsoid (Moritz, 1980) to be directly compatible with the Geocentric Datum of Australia 1994 (GDA94; Featherstone, 1994) and also GPS.
Prior to AUSGeoid98, previous Australian geoid models had been computed by evaluating Stokes's Integral with free air gravity anomalies only. This approach does not properly account for the effect of topographic masses external to the geoid (Heiskanen and Moritz, 1967). For AUSGeoid98, topographic effects were accounted for more precisely by first adding terrain corrections to the free air anomalies  evaluating Stokes integral and then applying the primary indirect effect given by the quadratic term of Wichiencharoen's (1982) formula, to produce the final geoid model.
(ii) Onshore, over 690,000 terrestrial gravity observations (from the 1996 release of the national gravity database) were reduced to Faye gravity anomalies (Heiskanen and Moritz, 1967) and combined with, 134,539 shipborne gravity data and satellite altimetry-derived gravity anomalies offshore and gridded on a 2 by 2 arc-minute grid using tensioned splines (Smith and Wessel, 1990). This gravity anomaly grid was then used in a remove-computerestore scheme to compute the residual geoid undulations with full 360 degree and order spherical harmonic expansion of the Earth Gravity Model 1996 (EGM96; Lemoine et al., 1998) acting as the reference model.
(iii) Stokes's integral is a convolution integral used to transform gravity anomalies (Δ ) into geoid undulations or N-values (Eq. 1.1.1). Prior to AUSGeoid98, Australian geoid models computed from gravity data evaluated Stokes's integral by brute-force numerical integration using an unmodified Stokes kernel ( (Ψ)). For AUSGeoid98, the convolution

< A c c e p t e d M a n u s c r i p t >
integral was evaluated in the Fourier domain, where the convolution of two functions (the gravity anomaly and modified Stokes's kernel) becomes their pointwise products -which is substantially more efficient to evaluate (Haagmans et al., 1993;Zhang, 1997). To further reduce the computational load, and to accommodate the latitude dependency of the computation point in spherical distance (Ψ) / Stokes's kernel, the two dimensional convolution was evaluated as the sum of one dimensional convolutions (Eq. 1.1.2).
In addition to this, the geoid undulations were partially high-pass filtered by only evaluating the integral over a restricted radius from the computation point ( <̂) and by adapting Stokes's kernel to produce the modified kernel, ̂( Ψ) described by . Vaníček and Featherstone (1998) discuss the partial filtering properties of limited integration caps and deterministically modified kernels. A range of geoid models were produced using a number of capped integration radii ( ̂) with the modification degree set to M=20. The intention of the [partial] high pass filtering was to alleviate the contribution of any long-wavelength errors present in the national gravity database.
The range of potential AUSGeoid98 models were finally compared to a dataset of 900 geometric geoid undulations calculated from differentially levelled AHD heights H and GPS derived ellipsoidal heights h. The comparison was used to identify the optimal choice of integration cap radii and to gauge the precision of the geoid model for its intended main use in GPS-levelling. With an RMS in the GPS-levelling/geoid differences of ± 0.364 m, the integration cap radius was determined to < A c c e p t e d M a n u s c r i p t > be 1 ∘ . To align the AUSgeoid98 model with the AHD a zero degree term was determined from the 93.7 cm mean of the GPS-levelling/geoid differences and subtracted from the model, which also absorbs any constant offset between the AHD and the geoid.
Whilst AUSGeoid98 was a substantial improvement over its predecessors, the ±0.364 m RMS of the comparison of the model to GPS minus AHD heights was much larger than its international counterparts . The differences were determined to be largely attributable to the discrepancies in the way AHD was realised (e.g., Kearsley, 1988;Morgan, 1992;Featherstone, and Stewart, 1998;Featherstone and Filmer, 2008;Filmer and Featherstone, 2009).
One prominent feature of the differences was a large N/S tilt , which was identified as being due to mean sea surface topography (aka the ocean's time-mean dynamic topography) embedded in AHD due to referencing the levelling observations to 30 estimates of mean sea level  around the Australian mainland. There also a number of local distortions (> 150 mm) present in the levelling network, which is to be anticipated when differential levelling is conducted over a vast area such as mainland Australia to third order tolerances (12 root k).
The primary intention of the AUSGeoid98 model was to enable GPS users to determine AHD heights, but due to the significant differences this was not possible to an acceptable degree of accuracy for most uses . Incorporating the difference between the AHD reference level and the geoid model, by interpolating the GPS-AHD and geoid differences onto a grid and including this as second "layer" to the model, has been proposed in a number of studies (e.g. Featherstone 1988, 2000, Featherstone and Sproule, 2006. It was formally implemented in 2009 with the computation of a new, hybrid quasigeoid model, AUSGeoid09. AUSGeoid09 was computed in two stages (Featherstone et al. 2011, Brown et al., 2011. Before this nationwide implementation of hybrid models, a small hybrid patch in Perth, Western Australia, had been < A c c e p t e d M a n u s c r i p t > applied to the released version of AUSGeoid98 because of the steep geoid gradients in this area (Featherstone, 2000).
First, a gravimetric quasigeoid model was computed, following an almost identical computational strategy as AUSGeoid98, using the latest release of the Australian gravity database (~1.4 million measurements), more up-to-date altimetry derived gravity anomalies -without including the shipborne data (Featherstone, 2009), and using EGM2008 to degree 2160 (Pavlis et al., 2012(Pavlis et al., & 2013 as the reference model in the remove-compute-restore process. The gravity data were gridded at 1 by 1 arc minute and a range of geoid models were computed by varying the integration cap radius and, this time, also the Stokes kernel modification degree. The optimal parameters were determined to be an integration cap of 1 degree and modification degree of M=40. Importantly, the primary indirect effect was not included in the computation so as to compute a "quasi"geoid model (Molodensky et al 1960, Moritz, 1968. Featherstone and Kirby, 1998, which is more consistent with the normal-orthometric height system adopted for the AHD (Featherstone and Kuhn 2006, Filmer et al. 2010. A second layer was then added to the gravimetric quasigeoid model which captured the differences between the "gravimetric-only" quasigeoid model and GPS-derived GDA94 (ICSM, 1996) ellipsoidal heights minus AHD levelled heights. This layer was calculated by gridding the difference using least squares collocation, after removing and then subsequently restoring a tilted plane that was predominantly north-south because of the tilt in the AHD (Feather atone and Filmer, 2011).
After the inclusion of this additional layer, the standard deviation of the ellipsoidal -AHD heights and quasigeoid differences dropped from ±222 mm to ±30 mm. AUSGeoid09 was then used as the national standard for converting GPS measurements to AHD heights for the following seven years.
AUSGeoid2020 (Fig 1.1.1 a) is the most recently computed model used for calculating AHD heights from ellipsoidal heights . AUSGeoid2020 was < A c c e p t e d M a n u s c r i p t > North/North/East direction accounting for the ~ 7 cm/yr continental drift. Importantly, GDA2020 ellipsoidal heights also differ from GDA94 heights by up to 9 cm (ICSM, 2000). The geometric layer (Fig 1.1.2 b) was determined by gridding the difference between GDA2020 ellipsoidal heights minus AHD heights and the gravimetric quasigeoid -for data points at 7624 locations -by least squares collocation .

13
(ii) The AUSGeoid2020 model has an accompanying map of uncertainty estimates (Fig 1.1.1 b).
The uncertainty estimates were calculated by propagating uncertainty estimates of the raw data through each stage of the processing. Gravity anomaly uncertainties were calculated from error estimates of the gravity data and positioning information, propagated through Stokes integral and combined with the published EGM2008 uncertainty grids in a removecompute-restore scheme . These were then augmented with the error propagation from the least squares collocation gridding of the geometric component.
The uncertainty values provide a way for users of the model to calculate how precise their GPS/AUSGeoid2020 derived AHD heights may be.

Patching an error in Port Phillip Bay, Victoria, SE Australia
Following the release of the AUSGeoid2020 model, a number of independent tests -using GDA2020 ellipsoidal heights and levelling data that were not included in AUSGeoid2020 -were conducted by State Government departments to ensure the model was precise and "fit for purpose". However, a test in Victoria identified a spike in the AUSGeoid2020 model in the vicinity of Melbourne and Geelong (Figure 2.1.1).

< A c c e p t e d M a n u s c r i p t >
Melbourne and Geelong are predominately built up around Port Phillip Bay -this is a Bay which covers an area of ~2000 2 with an average depth of only 8 m and opens into the Bass Strait through a very narrow channel. Offshore, gravity data in the underlying AGQG2017 gravimetric component of AUSGeoid2020 is determined entirely from satellite altimetry data -including within this bay. Altimetry data are notoriously unreliably over shallow water and in the near-shore zone (e.g., Deng and Featherstone, 2006).
The altimetry derived gravity anomaly used for AGQG2017 contained a > 50 mGal spike, that carried onshore during the evaluation of Stokes's integral and resulted in a > 1 m error in AGQG2017 (Fig 2.1.1). To remedy this error, we created a patch for the model by computing a new quasigeoid over a 5 by 5 degree area, centred over the bay, replacing the Sandwell et al.

Fine-scale parameter sweeps
The computation of AGQG2017 evaluated Stokes's integral using the FEO Stokes kernel modification , with a range of modification degrees of [0, 40, 60, 80, 100, 120 and 140] and integration cap radii of 0 to 5 degrees in increments of 0.5 degrees. In total, we computed 77 individual quasigeoid solutions. The choice of parameters (identified to be a modification degree of M=40 and integration cap radius of 0.5 degrees on comparison to GPS ellipsoidal minus AHD heights) was somewhat towards the edges of the explored parameter ranges.
To investigate whether any improvements could be made to these parameter choices, we conducted a fine-scale parameter sweep, exploring modification degrees range from 0 to 360 in increments of 10 degrees and integration cap radii between 0 and 5 degrees in increments of 0.1

< A c c e p t e d M a n u s c r i p t >
degrees. We further explored the use of additional deterministic kernel modification schemes, namely the unmodified spherical Stokes kernel, Wong and Gore (1969), Meissl (1971), Heck and Gruninger (1987) and Vaníček and Kleusberg (1987) [also see Vaníček and Sjöberg (1991)] modification schemes. In total, this produced more than 10,000 separate quasigeoid solutions.
Each quasigeoid model was then compared to (a) GDA2020 ellipsoidal heights minus AHD heights, and also (b) GDA2020 ellipsoidal heights minus a readjustment of the AHD levelling height data which were minimally constrained to just a single point and had normal corrections applied (as described in . The intention of using the second ground truth dataset was to further mitigate any issues associated with using "official" AHD heights that are contaminated by a slope originating from the ocean's mean dynamic topography (Featherstone and Filmer, 2011).
The results (Fig 2.2.1 to 2.2.6) indicate that (i) the parameter choices made for AGQG2017 was near optimal in terms of the agreement of the model to GDA2020 ellipsoidal minus AHD levelling heights (approx. ±185 mm). (ii) The Featherstone et al. (1988) and Vaníček and Kleusberg (1987) kernels become unstable for high modification degrees (>150) and larger integration cap radii (>2 degrees). It is not entirely clear where the instability arises from. Both these kernel modification schemes require that a system of linear equations are solved . One candidate for the instability is the matrix inversion required to solve these equations (which increases in size with the modification degree), but this needs further investigation. Future computations should exercise care to ensure the instability does not influence results. The minimum value here for these two kernel modifications actually lies within this area of instability and so it should be interpreted with caution.

Planar terrain corrections versus the Molodensky G series
The height anomaly / quasigeoid is the solution to the Molodensky boundary value problem and the solution can be obtained by using Stokes's integral with the Molodensky gravity anomaly on the Earth's surface (Molodensky et al. 1960Heiskanen and Moritz 1967

< A c c e p t e d M a n u s c r i p t >
The classical planar terrain correction, Δ (Eq. 2.3.2) was used in place of the 1 term for the computation of AGQG2017 as an approximation. To investigate the effect of this approximation on the accuracy of the quasigeoid model, we computed 1 terms and recomputed the series of quasigeoid models, exploring the same ranges of Stokes modification parameters and integration cap radii as used for AGQG2017 . However, using 1 value rather than Δ resulted in a worse fit of the quasigeoid to GDA2020 ellipsoidal heights minus AHD heights. The 1 value and gravimetric terrain corrections used in the computation of AGQG2017, were calculated at the same resolution as the digital elevation model (DEM), a high resolution 1 arc second grid (~30 m). For the 1 terms, this required a grid of free-air gravity anomalies to also be available at the same resolution as the DEM. This made it necessary to grid the gravity data at an unrealistically high resolution given the ~4 2 on-average gravity observation spacing. We suspect this to be the cause of the unexpected worsening of the fit of the quasigeoid model to the GDA2020 ellipsoidal heights minus AHD heights. To investigate this issue, we further recomputed both the terrain corrections and 1 terms using a block-averaged DEM and 1 arc minute grid (equivalent to the resolution of the AGQG2017 model) and produced another suite of quasigeoid models. While in this instance the Molodensky gravity anomaly with the 1 terms gave the superior fit to the GDA2020 ellipsoidal heights minus AHD heights, using the 1 arc second planar terrain correction still gave the best result (Fig 2.3.1).

< A c c e p t e d M a n u s c r i p t >
Note that this conclusion is not definitive because of errors in all the data used. For instance, the < A c c e p t e d M a n u s c r i p t > 24 formal error budget for the GPS-levelling data is ~40 mm on average , which is larger than the change in standard deviations in Figure 2.3.1.

Potential improvements to be offered by airborne gravimetry data
In 2011, an airborne gravity survey was conducted over the Gippsland Basin in southern Victoria, SE Australia. The survey covered the littoral zone up to 40 km offshore and 10-20 km onshore and was conducted as a part of the Victorian Government's CarbonNet Project for geological modelling of subsurface structure for potential carbon dioxide sequestration. Over this offshore area, the gravity data used for AGQG2017 is derived from satellite altimetry -and are generally unreliable close to the coast (e.g., Deng and Featherstone 2006).

< A c c e p t e d M a n u s c r i p t >
We have compared these airborne gravity data (Fig 2.4.1 a) to the gravity anomalies used to compute the AGQG2017 model and see differences with a range of around 9 mGal and a near-zero mean. The comparison was made by interpolating the 1 arc minute gridded AGQG2017 gravity anomaly to the location of the airborne data and calculating the difference. The interpolation was carried out using least squares collocation, with a 3D logarithmic covariance function (Forsberg, 1987), in a remove-compute-restore procedure with EGM2008 to degree 2170 acting as the reference field. The differences (Fig 2.4.1 b) are typically largest onshore in areas where the terrestrial gravity data thin out, and offshore where the altimetry derived gravity anomalies are unreliable. Wu et al. (2019) explored the benefit that these data offer to the quasigeoid modelling in this coastal region. However, there is only a limited number (< 30) of GDA2020 ellipsoidal and levelled heights onshore in this region -so these data offer very little statistical power for assessing the precision of any local quasigeoid model. Instead, Wu et al. (2019) compared quasigeoid models computed with and without the airborne gravity data to independent altimetry derived geoid heights. The comparison to the quasigeoid determined without the airborne data was ±28 mm, including the airborne gravity data improved the agreement by 5 mm -down to ±23 mm. This constituted a potential 20% improvement, however, these data have not yet been assimilated into the Australian quasigeoid model.

Future plans and the Australian Vertical Working Surface (AVWS)
The AHD is a normal-orthometric height system and was realised at a series of differential levelling benchmarks from a constrained and staged least squares adjustment (Roelse et al., 1971). In principle, it is a levelling-only datum . The last 20 years of geoid modelling in Australia has revealed the 1971 realisation of the AHD introduced features that are not related to the Earth's gravity field (e.g.,  and alternatives to AHD have been discussed (e.g., . The geometric layer in AUSGeoid09/2020 attempts to ameliorate these features, but this datum deficiency remains problematic. Firstly, the geometric component has no physical interpretation offshore and so the model must be clipped at the coastline (Fig 1.1.1). This causes difficulties in instances where a seamless onshore/offshore datum is needed, e.g. to align bathymetric and topographic elevation models. Secondly, they are difficult to estimate where no AHD levelling data are available and must be interpolated from nearby values (sometimes > 100 km away).
The inclusion of the geometric layer in AUSGeoid2020 distorts the quasigeoid to fit these features of the AHD and, in doing so, degrades the quality of the model to act as a reference surface for truly physically meaningful heights. For this reason, AUSGeoid2020 will most likely be the last national model of the AHD reference surface. The AVWS (Australian vertical working surface) is an alternative vertical reference frame comprising a normal height system with a quasigeoid model (ICSM, 2021). The model currently underpinning the AVWS is AGQG2017, the gravimetric component of the AUSGeoid2020 model -aligned with GRS80 by the application of a zero degree term (ICSM, 2021). Much like AUSGeoid2020, the AGQG2017 has a companion grid of uncertainty estimates -and these are generally half the size of their AUSGeoid2020 counterparts (cf. Fig. 1.1.1 b & Fig. 3.1 b). The AVWS is directly compatible with GDA2020 ellipsoidal heights, < A c c e p t e d M a n u s c r i p t > since both reference GRS80, however it is also still compatible with heights established via differential levelling when connected to a location with a known GDA2020 ellipsoidal height.
The focus for quasigeoid modelling in Australia is now to refine the AGQGyyyy model by improving the underlying data used to compute it, as well as investigating improved processing routines. This will be achieved by collecting airborne gravimetry observations that densify the existing terrestrial gravity data coverage onshore and improve the quality of the gravity data across the littoral zone. Data collection is currently planned over Greater Adelaide (GA, Fig 3.2 a) and Greater Melbourne (GM) and the Eastern Victoria Highlands (EVH, Fig 3.2 b). GA and GM are two major urban areas in Australia that are located within 50 km of the coast where the quasigeoid model is likely to be impacted by erroneous altimetry-derived gravity anomalies. The EVH has been prioritised because the terrestrial gravity data coverage is particularly sparse (mainly due to inaccessibility) and corresponds to an area of large formal uncertainty in the AGQG2017 model. The planned flight lines are spaced 5 km apart over GA, 2 km over GM and 500 m over the EVH and data will be collected at a maximum height of around 160 m above the terrain, where it is safe to do so.
The intention is to ultimately drive the uncertainty in the model down to 2-3 cm -starting with high impact, urban, areas where the model is of utility to most users.

Concluding remarks
We have provided a historical review of regional geoid and quasigeoid modelling over the Australian continent. The earliest models were produced in the late 1960s and were used to reduce geodetic surveying measurements to the Australian National Spheroid associated with the nongeocentric Australian Geodetic Datum. Following the introduction of the AHD in 1971, and the introduction of GPS in the mid 1980s, the main reason for modelling the geoid was altered so as to allow for ellipsoidal heights (determined from GPS) to be transformed into AHD heights.
With the computation of AUSGeoid98, a number of flaws in the 1971 realisation of AHD heights became evident. To enable users to determine AHD heights more precisely, these flaws were assimilated into two subsequent geoid model computations, AUSGeoid09 and AUSGeoid2020.
However, it is now recognised that distorting the models into better agreement with the AHD significantly degrades their accuracy and their ability act as a reference surface for meaningful physical heights. The underlying "gravimetric-only" layer of the AUSGeoid2020 has been released officially to support users who want to determine physical heights more preciseely and reliably. This new model underpins an alternative vertical height system called AVWS (Australian vertical working surface).
Experiments have been conducted to investigate whether there is scope to improve the computational process used to produce the AGQG model. In particular, (i) a fine-scale parameter sweep of Stokes kernel modification types, modification degrees and integration cap radii has been used to produce more than 10,000 quasigeoid models, and (ii) the effect of approximating the Molodensky series using planar terrain corrections has been scrutinised. The results indicate that the process used to produce AGQG2017 is adequate, given the availability of current gravity data and the limitations of the GPS-AHD heights used to assess the experimental quasigeoid models. The < A c c e p t e d M a n u s c r i p t > future direction of Australian quasigeoid modelling is to improve the underlying gravity data used to produce it. This will be accomplished by collecting airborne gravimetry data to regularise the existing data coverage onshore and improve the quality of the existing gravity data over the littoral zone.
Featherstone WE (1994) An explanation of the Geocentric Datum of Australia and its effects upon future mapping. Cartography 23 (2)