Quasigeoid computation in the Nordic and Baltic region

1 Introduction

Heights on land are historically estimated by geometric levelling to one or more local tide gauges on the coast of the relevant region, hence the subsequent derivation of “height above sea level.” The levelled heights are commonly processed into either the orthometric heights, H , or normal heights, H * . The orthometric height is the distance between a point on the surface and the geoid, measured along the plumb line. However, knowledge of the density distribution of the masses, between the surface and the geoid, is required to determine the direction of the plumb line. To avoid this challenge, the normal height was introduced as the distance between a point on the surface and the quasigeoid, measured along the normal plumb line. Both H and H * refer to local mean sea level as the (quasi)geoid, determined at a tide gauge. It can be affected by the mean dynamic ocean topography and geodynamic effects.

The development of global navigation satellite systems (GNSS) allows for direct height measurements from satellite positioning, without the resource-demanding work of levelling. The GNSS ellipsoidal heights, h , refer to a mathematical reference ellipsoid of the Earth, e.g. the World Geodetic System 1984 (WGS84) or the Geodetic Reference System 1980 (GRS80) ellipsoids, respectively.

To smoothly transform between these references a height transformation surface is a prerequisite, either a (local) geoid for orthometric heights, or a quasigeoid for normal heights. For the geoid, the equipotential level surface in the gravity field of the Earth coinciding with the mean sea level, the transformation is obtained by

where N is the height of the geoid above the reference ellipsoid. Similarly, for the quasigeoid,

where ζ , called the height anomaly, is the height of the quasigeoid above the reference ellipsoid. The quasigeoid is a geoid-like surface, but is not an equipotential level surface in the gravity field, except over the oceans. On land, its deviation from the geoid is dependent on the local topography (Flury and Rummel 2009, Sjöberg 2010). Local and regional models of N or ζ are therefore the tie between the physical height systems in a country and the application of GNSS height measurement. They are also an essential part in the implementation of the global physical height system, the International Height Reference System, and its realisation, the International Height Reference Frame (Sánchez et al. 2021).

The recently concluded Finalising Surveys for the Baltic Motorways of the Sea (FAMOS) project (Schwabe et al. 2020; Liebsch et al. 2023) has resulted in several new marine and airborne gravity datasets in and around the Baltic sea, as well as the inclusion of more gravity data from neighbouring countries into the FAMOS gravity database. The project also resulted in quality control of older gravity surveys in the database and more homogenised GNSS-levelling data on land. The FAMOS region is a sub-region in the region made up of the countries in the Nordic Geodetic Commission (NKG). The geodetic units in the NKG have a tradition of making official NKG gravimetric geoids/quasigeoids at regular intervals of the whole region, starting in 1986 with Tscherning and Forsberg and last in 2015 (Ågren et al. 2016). The NKG members consist of Denmark, Estonia, Finland, Iceland, Latvia, Lithuania, and Sweden. On the background of contributing in the computation of gravimetric quasigeoid for the new Baltic Sea Chart Datum 2000 under the FAMOS project, and in light of the updates to the gravity database, an interesting question is to what degree the new data would improve the quasigeoid in the whole NKG region (not including Iceland). The NKG region is a topographically varied region: from predominantly flat and hilly terrain in the southern and eastern countries to rugged, mountainous terrain in Fennoscandia. Especially, the inclusion of Norway, with its rugged topography, makes the full NKG region different from the FAMOS region. This sort of terrain poses an interesting and challenging task in (quasi)geoid modelling, particularly with the major goal of reaching the 1 cm accuracy, to match the increasing accuracy of ellipsoidal heights from GNSS measurements. There are many different aspects to consider.

The main aim of this study is therefore to investigate improvements to the accuracy of the gravimetric quasigeoid in the NKG region with the inclusion of the new gravity data and the quality-controlled older gravity surveys from the recently FAMOS project. There is a big variety of different parts in the quasigeoid modelling, which influence the resulting gravimetric quasigeoid. This study is limited to the consideration of the following parts, by varying:

1. the resolution of the Gaussian filter in the creation of the mean elevation surface (m.e.s.) for the Residual Terrain Model (RTM) method (Forsberg 1984, 1985), which is used to represent the short-medium wavelength gravity signal from the local terrain.

2. the range of the bounding degrees in the Wong–Gore modification of the Stokes’ kernel, to locate an optimal transition from satellite gravity data to terrestrial gravity.

Regarding the latest point, for previous geoid models computed at DTU Space, the custom has been to use no modification at all (NKG86) or shorter intervals from degree 30 to degree 60 in later geoids. For the last NKG geoid computation of 2014/2015, the chosen Wong–Gore kernel modification range was from degree 90–100 for the DTU Space solution (Ågren et al. 2016). However, with the advancement in satellite gravimetry, the kernel modification ranges need to be revised, to adapt to the more accurate long wavelength gravity signal in the global geopotential models (GGMs). The results from the Colorado geoid computation experiment (Wang et al. 2021) showed that this is applied by the many geoid computation centres, which use kernel modifications in their geoid computation strategies. In the NKG region, a similar study, testing multiple combinations of Wong–Gore modification ranges, has been conducted in regional quasigeoid modelling in Finland (Saari and Bilker-Koivula 2018).

The influence on and the accuracy of the resulting gravimetric quasigeoids are assessed by comparison to two sets of GNSS-levelling height anomalies. The first covers the FAMOS region and the countries involved herein: the NKG member states, Germany, and Poland. The second is an older dataset covering all the NKG countries, which was used for the assessment of the NKG geoid computations for the NKG2015 gravimetric quasigeoid. The fit, in form of a normalised root mean square error (RMSE) of the residuals, is analysed in a twofold approach: first, by direct comparison, and second, by applying a simple tilt correction in the manner proposed by Gruber et al. (2011) and Gruber and Willberg (2019).

A variety of methods can be and are used for the determination of the geoid or the quasigeoid, as was seen in the Colorado geoid computation experiment. In this study, we use the spherical fast Fourier transform (FFT), with a Wong–Gore modified kernel, to compute the gravimetric quasigeoid in a remove–compute–restore (R-C-R) scheme.

The rest of the article is structured as follows. Section 2 describes the methodology of the geoid determination, the computation, and assessment. Section 3 describes the data used to compute the gravimetric quasigeoids, and Section 4 shows and analyses the results of the assessment. The article is concluded with a discussion of the results in Section 5 and a conclusion in Section 6.

2 Methods

The following section will shortly explain the main theory behind the method applied in the quasigeoid computation, as well as the steps in the computation of the quasigeoid for the region.

The anomalous gravity potential T , which is the difference between the actual gravity potential W and the normal gravity potential U , is reduced by removing the known trends, in such a fashion that

where T GGM is the component from the GGM representing the long wavelength part of the anomalous gravity potential, T RTM is the component due to the effect of the terrain, computed with the RTM method, representing the medium-short wavelength part of the gravity potential, and T res is the residual anomalous gravity potential. Analogously, the gravity anomalies Δ g and the geoid N or height anomaly ζ are similarly divided into a long wavelength part, a short-medium wavelength part, and a residual part.

The classical method of determining the Geoidal undulation N involves determining Stokes’ integral

where R is the radius of the Earth, γ 0 is the normal gravity, Δ g is the gravity anomaly, and S ( ψ ) is Stokes’ function (Hofmann-Wellenhof and Moritz 2006). Solving Stokes’ integral requires that the topography is shifted inside the geoid or otherwise removed to meet the requirements of the spherical harmonic field outside the geoid. The formula yields N , by integrating over the whole sphere. In principle, this would require gravity data over the whole sphere, which is near impossible. However, local gravity and height data in combination with a global gravity model can be used to estimate a local geoid, where the GGM represents the long wavelengths and the gravity observations and digital elevation models (DEM) represent the medium-short wavelengths, as mentioned earlier.

The Molodensky geodetic boundary value problem can be solved to compute the quasigeoid instead of the classical geoid, by essentially the same integral, except for a correction to take into account the variations of the terrain. The Molodensky integral (Molodensky et al. 1962) to first order is derived as

where Δ g is the free-air gravity anomaly, but now referring to the terrain level (the “Telluroid”) instead of the geoid, assumed in the Stokes integral (cf. equation 8–23, Hofmann-Wellenhof and Moritz 2006). The term g 1 is the first-order Molodensky correction (cf. equation 8–62, Hofmann-Wellenhof and Moritz 2006). The height anomaly ζ denotes the distance from the topography to the telluroid and can, in modern interpretation, be thought of as the harmonic upward continuation of the classical geoid. The relation between ζ and the anomalous potential T is described by an equivalence of Bruns’ theorem:

The difference between N and ζ can be determined by

where g ¯ is the mean gravity along the plumb line between the geoid and the ground, γ ¯ is the mean normal gravity along the normal plumb line between the ellipsoid and the telluroid, H is the orthometric height, the distance from the geoid to the topography, and H * is the normal height, which is the distance from the ellipsoid to the telluroid. The term g ¯ − γ ¯ is simply approximated by the Bouguer gravity anomaly (cf. equations 8–113, Hofmann-Wellenhof and Moritz 2006). In practice, it is therefore easy to convert between geoid and quasigeoid, since the linear approximation in equation (7) is corresponding to the similar linear approximation in the conversion of levelled geopotential numbers in (Helmert) orthometric heights.

To determine the Molodensky integral, it is rewritten as a spherical convolution integral, in the geodetic coordinates ( φ , λ ) , to use the FFT techniques (Strang van Hees 1990):

where F is the Fourier transform and F − 1 is the inverse Fourier transform, and S ( ψ ) is defined as

with

Furthermore,

where Δ φ = φ − φ P and Δ λ = λ − λ P , for the computation point φ P , λ P . The last term cos ( φ P ) cos φ of equation (11) is approximated as

where φ m is the mean latitude for the area.

The components of equation (3) are described in the following part of Section 2, as they are part of the R-C-R scheme. T GGM is the component from a GGM, representing the long wavelength part of the anomalous gravity potential, also called the reference field. It is used to compute gravity anomalies Δ g GGM and the height anomalies ζ GGM , for the remove and restore steps, respectively. By spherical harmonic expansions, we can express Δ g GGM at point ( φ P , λ P ) as

and, similarly, for ζ GGM , we obtain

where G is the constant of gravity, M is the mass of the Earth, R is the radius of the Earth, γ is the normal gravity, C n m and S n m are the normalised spherical harmonics, and P n m are the fully normalised associated Legendre functions, with degree n and order m to the maximum degree N in the spherical harmonic expansion (Heiskanen and Moritz 1967)

T RTM is the component from the terrain effects representing the medium-short wavelength components of the anomalous gravity potential T . The remove and restore components, Δ g RTM and ζ RTM , respectively, are computed with the RTM method (Forsberg 1984, 1985), where the reduction of the terrain effect is made relative to a low-pass filtered m.e.s. of the topography, with the revisited harmonic correction for points situated below the m.e.s. (Klees et al. 2022). The RTM effect on the gravity anomalies is approximated by the following expression:

where ρ is the density (set to 2.67 g/cm³, h is the height of the topography, h mes is the height of the m.e.s., and t c is the classical terrain correction) (Forsberg 1984). In this study, we did, however, use the full prism integration of the terrain effects. The similar terrain effect for the height anomaly ζ RTM may similarly be approximated in a planar approximation by the following expression:

yielding the terrain effect from the height anomaly, which can be determined with FFT methods (Forsberg 1985, Omang and Forsberg 2000) and is used here.

A Wong–Gore kernel modification (WG-mod.) on the Stokes kernel is applied (Wong and Gore 1969), to keep the accurate gravity field information, in the long wavelengths of the spherical harmonic expansion of the gravity field, from being influenced by the terrestrial gravity data. The WG-mod. is of the form

where the transition coefficient α ( n ) increases linearly from 0 to 1, resulting in a linear tapering from satellite data to terrestrial gravity data, between the degrees N 1 and N 2 , with

for n = 2 , … , N . The optimal choice of transitional degrees, N 1 and N 2 , is both loosely based on the error characteristics of satellite gravity missions and is empirically estimated from tests of fit with the available GNSS-levelling data in the region.

2.1 R-C-R

The R-C-R technique is used to remove the GGM and terrain effects from the gravity anomalies, making them as small as possible in terms of RMS variability, for smoother gridding of the datapoints. Furthermore, by evaluating Molodensky’s integral on the gridded residual gravity anomalies, the g 1 -term is of an insignificant influence once RTM terrain corrections are applied (Forsberg and Sideris 1989), which simplifies the computation of the residual height anomalies. After determining the residual height anomalies, the GGM and terrain effects are restored. All computations are made with respect to the GRS80 reference ellipsoid (Moritz 2000) and in the zero-tide permanent tide system. An overview of the data is given in Section 3.

The remove step consists of representing the free air gravity anomalies by long and medium-short wavelength part and subtracting these from the free air gravity anomaly data, to obtain the residual gravity anomalies. More explicitly this is expressed as

(19) Δ g FA = Δ g GGM + Δ g RTM + Δ g res .

The remove steps for the gravimetric quasigeoids are as follows:

1. The part of the GGM, Δ g GGM , is reduced. This is done in the form of a “sandwich grid” interpolation, i.e. between two Δ g GGM grids, computed at 0 and 2,500 m, enveloping the topography of the region. The Δ g GGM value at the observation height is then subtracted from the gravity observations by three-dimensional grid interpolation.

2. The RTM terrain effect is computed and subtracted by use of a m.e.s. of resolution 30’ which is constructed by a Gaussian filter with a full width at half maximum of corresponding resolution, thus, yielding the residual gravity anomaly, Δ g res , at every observation point. Different filter lengths were tested for the region and the best results were found for a filter length of around 30’, with marginal differences in using a longer filter length, and slightly worse for a shorter filter length.

Table 1 shows the statistics of the remove step, when subtracting the GGM, computed to d/o 719 and a 30’ m.e.s in the RTM method from Δ g . The mean is close to 0, as it should be, when removing what is possible to model from the gravity observations. It seems that the XGM2019e (Zingerle et al. 2020) at d/o 719 minimises to a satisfactory extent, although there remains some variation in the Δ g res .

Table 1

Statistics of the gravity reduction in the remove step

	Mean	Std
Original data	− 0.82	24.71
Free-air - GGM	− 0.91	11.88
Free-air - GGM - RTM	0.42	6.40

No of gravity datapoints: 485301. XGM2019e d/o 719 and 30 ′ m.e.s. in the RTM terrain effect. Units in mGal.

The residual gravity anomalies are gridded over the entire region, with a resolution of δ φ = 0.0 1 ∘ , δ λ = 0.0166 7 ∘ , corresponding to approximately 1 km², using least-squares collocation, with individual weights from the associated standard deviations of the gravity data, and a 15 km correlation length. The specific choice of correlation length was chosen based on empirical testing of several different correlation lengths.

The gridded residual gravity anomalies are evaluated with Molodensky’s integral using the spherical FFT method (equation (8)) with five reference parallels, dividing the region up into five latitudinal bands, overlapping the adjacent bands, during the computations (Forsberg and Sideris 1993). A WG-mod. is used on Stokes’ function, to remove the lower harmonics up to a certain degree followed by a linear tapering in transition to a higher degree. The optimal linear modification range is investigated (Section 4.2) by comparison to GNSS-levelling data in the region.

The restore step consists, similar to the remove step, of restoring the full height anomaly from the computed residual height anomaly,

(20) ζ = ζ GGM + ζ TC + ζ res .

After computing the residual height anomalies, ζ res , the medium-short and long wavelength components are restored in the following steps:

1. RTM terrain effects of the height anomaly, ζ RTM , are computed in a grid for the whole region using the m.e.s. and spherical FFT, with a WG-mod. range corresponding to that used in the evaluation of the Molodensky integral, and then added to ζ res .

2. the GGM component, ζ GGM , is, as in the remove step, also computed in a sandwich grid interpolation to model the 3D variation of the GGM. The height anomalies from the spherical harmonic expansion are computed at elevations of 0 and 2,500 m, and then interpolated onto the fine resolution DEM, to yield ζ GGM on the topography. This grid is then restored by addition to finally yield the gravimetric quasigeoid ζ .

The gravimetric quasigeoid is shown in Figure 3. All computations are made with the GRAVSOFT suite of Fortran programs (Tscherning 1992).

Figure 1

Free air gravity anomalies extracted from the NKG gravity database. Units in mGal.

3 Data

An overview of the gravity data, DEMs, and GGM grids is shown in Table 2. The DEMs and the GGM grids are deliberately covering a larger area than the gravity data to account for computational effects at the edges of the computation region.

3.1 Gravity anomalies

The gravity data are extracted from the NKG gravity database (maintained by DTU Space), which includes the latest measurements from the FAMOS project, and covers the area shown in Figure 1. The excerpt of data from the database consists of ground, marine, and airborne data, yielding 493,447 gravity observations after thinning to approximately 1 observation pr. 1 km², primarily to decimate very densely sampled marine data sources.

Figure 2

DEM of the NKG region.

3.2 DEMs

A DEM, based on the NKG2014 DEM, is used in the quasigeoid computation, in three versions, covering the region ( ϕ ; λ ) = ( 5 1 ∘ to 7 5 ∘ ; − 4 ∘ to 3 8 ∘ ), as shown in Figure 2. The detailed version has a resolution of 0.3 ′ × 0.45 ′ , approx. 500 × 500 m resolution, while a coarser version has a resolution of 0.6 ′ × 1 ′ , approx. 1 × 1  km. The third grid is a m.e.s., which is a low-pass filtered DEM. All three versions are used for computing the RTM effect in the remove and restore steps. Different filter lengths are used for the m.e.s., as mentioned in Section 2.1.

Figure 3

Gravimetric quasigeoid: heights are in metres relative to the GRS80 reference ellipsoid, with contour lines 2.5 m apart.

3.4 GNSS-levelling height anomalies

To assess the computed quasigeoids, two available sets of GNSS-levelling data are used.

The first consists of GNSS-levelling height anomalies, which were used for evaluation in the FAMOS project (FAMOS GNSS-levelling data), covering the FAMOS quasigeoid area, from latitude 53 ∘ to 66.5 ∘ and longitude 8.5 ∘ to 31 ∘ . The location of the FAMOS GNSS-levelling data is shown in Figure 4.

Figure 4

Plots of the residuals from the comparison of the FAMOS GNSS-levelling height anomalies and the gravimetric quasigeoid: (a) normalised residuals, (b) tilt-corrected residuals, and (c) tilt correction.

The GNSS data are in ITRF2008 (Altamimi et al. 2011), with the zero tide system, and a land uplift referring to epoch 2000.0 in the NKG2016LU (Vestol et al. 2019) for the NKG countries.

The second set consists of GNSS-levelling data (NKG2008 GNSS-levelling data) for all the NKG member countries in the region (Hakli et al. 2016). The GNSS data are in ETRF2000, with the zero tide system, and a land upplift referring to epoch 2000.0 in EVRF2007 (Sacher et al. 2009). It was used in the assessment of the latest official NKG geoid computation prior to the release of the official NKG2015geoid. The location of the NKG2008 GNSS-levelling data is shown in Figure 7.

Figure 5

Normalised (circles) and corrected (diamonds) RMSE of the residuals between the FAMOS GNSS-levelling height anomalies and seven gravimetric quasigeoids, with Wong–Gore kernel modification ranges of 20 degrees (start and end degree indicated on bottom axis). The colours of the markers correspond to the countries in the legends.

The FAMOS GNSS-levelling data is a subset of the NKG2008 GNSS-levelling data. There are improvements to this dataset, compared to the NKG2008 dataset, especially in the Baltic countries.

4 Results

The result of the assessment is analysed in the current section. It is divided into four subsections. First the general results of the tilt correction, followed by an evaluation of the optimal WG-mod. range, where two approaches are pursued: best fit with a 20 degree WG-mod. range, and best fit with a modification range from 50 and up to 120 degrees, in the full range from N 1 = 100 being the lowest and N 2 = 240 as the highest. Finally, a comparison to previous NKG gravimetric quasigeoids is made.

4.2 Wong–Gore kernel modification range

Gravimetric quasigeoids, computed with XGM2019e to d/o 719 and an RTM effect from a 30 ′ resolution Gaussian filtered m.e.s., were computed with a variety of different WG-mod. ranges, and compared to the FAMOS GNSS-levelling height anomalies.

The WG-mod. ranges in intervals of 20 degrees are tested from degree 100–120 to 220–240, for the countries individually and for all countries combined. The normalised RMSE of the residuals for the gravimetric quasigeoid with the best-fitting modification range is shown in Table 4, which is for the interval of 160–180 degrees, with an RMSE of 2.06 cm to all the FAMOS GNSS-levelling height anomalies. There is less than a 0.5 mm difference in RMSE between the three best-fitting gravimetric quasigeoids with WG-mod. range intervals of 20 degrees, namely 180–200, and 200–220 degrees that each have best fits to different countries in the region, interchangeably. Figure 5 shows the RMSE for all the different WG-mod. ranges. It illustrates how the RMSE of the residuals varies depending on the bounding degree intervals and that it is at a minimum for most of the countries in the region for the WG-mod. range from degree 160–180.

Table 4

Normalised and corrected RMSE of the residuals between quasigeoids with a different Wong–Gore modification ranges and the FAMOS GNSS-levelling data height anomalies

	Denmark	Estonia	Finland	Germany	Latvia	Lithuania	Norway	Poland	Sweden	All
No. of points	617	131	36	100	250	84	473	29	182	1902
	XGM2019e d/o 719, Wong–Gore modification range degree 160–180
Normalised RMSE	1.67	1.29	1.63	1.39	2.53	1.89	2.66	1.62	1.96	2.06
Corrected RMSE	1.21	1.29	1.48	1.18	1.88	1.83	2.60	1.59	1.92	1.87
	XGM2019e d/o 719, Wong–Gore modification range degree 150–220
Normalised RMSE	1.74	0.98	1.74	1.33	2.20	1.86	2.51	1.57	1.89	1.95
Corrected RMSE	1.34	0.96	1.58	1.15	1.71	1.79	2.46	1.55	1.83	1.78

The best fits for intervals of 20 degrees, and longer intervals are shown, respectively. See A for the RMSE of other setups.

Figure 6

Normalised (circles) and corrected (diamonds) RMSE of the residuals between the FAMOS GNSS-levelling height anomalies and seven gravimetric quasigeoids, with Wong–Gore kernel modification ranges of long intervals (start and end degree indicated on bottom axis). The colours of the markers correspond to the countries in the legends.

The tilt correction improves the fit of all the gravimetric quasigeoids in general, but with difference in magnitude. There is an improvement in the overall fit of around 2 mm for the case of the best-fitting gravimetric quasigeoid model, likely due to the lower RMSE of Denmark and Latvia, where the improvement is most significant, as illustrated by the blue lines in Figure 5.

For the longer intervals of WG-mod. ranges, the modification ranges are varying in ranges from 50 degrees to 120 degrees, from degree 100 to 240, with 11 different intervals tested. The normalised and tilt-corrected RMSE of the best-fitting gravimetric quasigeoids is shown in Table 4. Figure 6 shows that the RMSE of the residuals to all the FAMOS GNSS-levelling data is very similar for all 11 cases. The largest difference between the solutions is smaller than 0.1 cm in comparison with all the GNSS-levelling data. Similar observations are seen for the individual countries, where the differences in general are within 0.2 cm. The improvement from the tilt correction is 0.17 cm to all, for the best-fitting gravimetric quasigeoid. The biggest improvement is for the modification range from degree 100 to 220, of 0.26 cm for all countries and 0.61 cm to the Danish GNSS-levelling network.

Figure 7

Normalised residuals of the gravimetric quasigeoid ζ W G 80 and the NKG2008 GNSS-levelling height anomalies in the full NKG region.

But in general, the RMSE of the residuals from the different WG-mod. ranges gives very similar results, with no significant differences between them.

The results for all the different gravimetric quasigeoids are also given in Tables A1, A2, A4, and A5, with the differences from the tilt correction in Tables A3, and A6 in Appendix I.

5 Discussion

The study has looked at a few aspects in quasigeoid modelling in the NKG region, for optimising quasigeoid determination with the spherical FFT method and using the R-C-R scheme, and the findings will be discussed in the current section.

Firstly, a few different Gaussian filter lengths for the m.e.s. were tested in the RTM method, as mentioned in Section 2.2. Overall, the best fits to the GNSS-levelling data were obtained using either a 60 ′ or 30 ′ resolution Gaussian filter, where the latter was used for the rest of the study. For the large region of the NKG, with its diverse topographic characteristics, there does not seem to be much to gain from a shorter filter length. However, this could be different, if another region was being investigated, e.g. a topographically rugged region with mountains and fjords resulting in a big variation in the high-frequency component of the gravity field, or if a smaller region was investigated. There is also a possibility that a m.e.s. with locally varying resolution of the filter lengths could represent the variety of topography in a region such as the NKG region better and result in an improved regional gravimetric quasigeoid.

In comparing the short and long Wong–Gore kernel modification ranges of 20 degrees to the longer of 50–100 degrees, in Section 4.2, it is observed that the RMSE of the residuals is more sensitive to the choice of bounding degrees for the transition interval, when using a short interval, as can be seen in Figure 5. There is a large variation in the RMSE across the two figures. For WG-mod range intervals of 20 degrees, we observe that the best fits to different national GNSS-levelling data in the region are not obtained by one specific interval, e.g. the best fit for Denmark is for a from degree 160 to 180, for Estonia from degree 180 to 200, and for Norway from degree 200 to 220. This is in contrast to the results in Figure 6. The longer WG-mod. range intervals are, more or less, achieving the same RMSE of the residuals for all the tested bounding degrees. We interpret these results as showing that long intervals, in the kernel modification, are more robust, i.e. it is overall more likely to obtain a good fit of the gravimetric quasigeoid, with long interval bounding degrees. This appears to hold for a region of this type, but it might be different in a more homogeneous type of terrain.

Based on the direct comparison, it seems that there are improvements of some millimetres to be gained by increasing the interval of the WG-mod. range, from 20 degrees and up to around 60–70 degrees. But there also seems to be an upper limit to this range, as seen in the 100–220 degree modification range, which is not overall having a better fit to GNSS-levelling data.

The simple planar tilt correction is a convenient tool for analysing GNSS-levelling networks in the region, and the application gives a strong indication that some countries in the region have inherent systematic tilts, either in their national GNSS-levelling networks or in the gravity data. The latter seems more unlikely, since the quasigeoid model is a combination of several components in the R-C-R computation. The tilt correction pattern is not consistent over the whole region, but individual trends for the countries are observed. The tilt can also be observed from visual inspection of Figure 4(c). The results of the tilt-corrected comparison show an improvement in the fit of 0.15 cm and up 0.45 cm overall, as shown in Section 4.2. This improvement is especially dominated by the reduction of the tilts in the GNSS-levelling networks of Denmark and Latvia, which show big improvements from the tilt correction. Ideally, there would not be a significant improvement in the fit. We observe this for most of the countries in the region. The gravimetric quasigeoid with the longest Wong–Gore kernel modification range (degree 100–220) gives best fit to the whole region, when the planar tilt correction is applied (1.76 cm). However, all the gravimetric quasigeoids with a long modification, shown in Figure 6, give very similar results. This makes the task of determining the most optimal WG-mod. range in the region ambiguous, since there is a variety with equal results.

The second GNSS-levelling dataset, NKG2008, was used for comparing with the last official geoid computation within the NKG. The normalised RMSE of the residuals of the two the best-fitting gravimetric quasigeoids and the NKG2008 height anomalies, called ζ W G 180 − 200 and ζ W G 150 − 230 , to previous NKG quasigeoid models shows a comparable result to the official NKG2015 gravimetric quasigeoid for the whole region. The best fits to the individual NKG countries are varying between the NKG2015 quasigeoid, ζ W G 180 − 200 , and ζ W G 150 − 230 . The NKG2015 quasigeoid is very much superior in Norway, slightly better in Denmark, but worse in Lithuania and Latvia. With these differences in mind, there certainly seems to be a potential of improving the fit in the whole region. It is largely the fit to the NKG2008 height anomalies in Norway that stand out in the fits of ζ W G 180 − 200 and ζ W G 150 − 230 on one side, and the NKG2015 quasigeoid on the other. And it shows that it is possible to obtain a significantly better fit to the NKG2008 GNSS-levelling height anomalies in Norway than has been achieved in this study. The question is then, if this could be done without compromising the fits to the rest of the countries in the region.

The execution and completion of the FAMOS project has resulted in good amount of new marine gravity data in and around the Baltic Sea. In this connection, it should be noted that most of the added gravity data are marine or airborne over the sea, so the improvements in the fit from the added gravity data will, to some part, be the influence of improved marine gravity data and its effect on the geoid/quasigeoid on land. The addition of new gravity data and quality control of old gravity data in the eastern part of the NKG region, and improvements of GGMs, is most likely the reason for the improvement in the fit of the computed gravimetric quasigeoids in this study to the NKG2008 GNSS-levelling data in Latvia and Lithuania, when comparing to the offical NKG quasigeoid of 2015. And with the ongoing project BalMarGrav (https://interreg-baltic.eu/project/balmargrav), the marine gravity data in the Southern and Eastern Baltic should further improve, which could lead to further improvement of the geoid/quasigeoid modelling in the region. Furthermore, with the results of the tilt correction in mind, it should be possible to obtain improved fits between the gravimetric quasigeoid and the GNSS-levelling height anomalies in Denmark and Latvia, the countries showing the largest tilts in the tilt correction, after a quality control of their respective national GNSS-levelling networks.

With the practice of merging national gravimetric geoids with the national GNSS-levelling to make hybrid GNSS-geoids, systematic tilts of this size can easily overshadow the details of local gravity field variations in the gravimetric quasigeoids and geoids, in general. Thus, the possibility of significant systematic tilts in the national GNSS-levelling networks should be investigated and addressed accordingly wherever necessary. As long as local geoid and quasigeoid models are fitted to local tide gauges and national GNSS-levelling networks, to form the fundamental basis of local, national, or regional geodetic infrastructures, the idea of a tilt correction disappears in practice, and a trade-off between the best gravimetric quasigeoid and local GNSS-levelling data will be the result.

6 Conclusion

Gravimetric quasigeoids for the NKG region have been computed with new gravity data in the Baltic Sea from the FAMOS project. The method of quasigeoid determination was the spherical FFT, with a Wong–Gore modification of the Stokes’ kernel, in a R-C-R scheme. The lowest normalised RMSE of the residuals to all the FAMOS GNSS-levelling height anomalies was 1.95 cm. And the lowest normalised RMSE of the residuals to all the NKG2008 GNSS-levelling height anomalies was 2.79 cm. This is comparable to the RMSE of the latest official NKG quasigeoid, and thus not a significant improvement to the overall fit. There is an alternation between the gravimetric quasigeoid solutions of the study and the NKG2015 gravimetric quasigeoid of fitting the individual countries in the region best. Improvements in the fit to Latvia and Lithuania could likely be from the inclusion of the new gravity data from the FAMOS project.

The results, from comparison of the gravimetric quasigeoids to the FAMOS GNSS-levelling height anomalies indicate that a longer transition interval, in the Wong–Gore kernel modification, more robustly results in a low RMSE of the residuals, than a short interval. The optimal compromise between and transition from the long wavelength and the medium-short wavelength gravity field is investigated in the variety of Wong–Gore kernel modification ranges used. With an improving knowledge of the static long wavelength gravity field from satellite gravimetry and in improved GGMs, the range between the bounding degrees, in the WG-mod. range, could be expected to increase in future. Although the best results for the short interval WG-mod. range were comparable to the results of the long intervals, the former show less robustness than the latter for a region with this type of topography. At least, this is observed for the method of quasigeoid determination in this study. The transition degrees could be expected to be similar for computation methods of a similar nature, but this remains to be seen.

Furthermore, the results also indicate that some countries in the NKG region might have significant systematic tilts in their national GNSS-levelling networks. This was shown by performing a planar tilt correction of the normalised residuals from the comparison of the gravimetric quasigeoid solutions with the FAMOS GNSS-levelling data. In Denmark and Latvia, the RMSE of the residuals improved by around 0.4 and 0.5 cm post-tilt correction, respectively.

Finally, it would certainly be interesting to see what accuracy other computational centres in the region could achieve for geoid/quasigeoid models in the NKG region. The gravimetric quasigeoid NKG2015, the last official quasigeoid for the NKG region, is nearing 10 years since its release. We therefore propose the idea to start working on a new official NKG geoid or quasigeoid in the coming years.