Refinement of a gravimetric geoid model for Japan using GOCE and an updated regional gravity field model

Abstract

We developed a refined gravimetric geoid model for Japan on a 1 × 1.5 arc-minute (2 km) grid from a GOCE-based satellite-only global geopotential model and a regional gravity field model updated in this study. First, we have constructed a regional gravity field model for Japan using updated gravity datasets together with a residual terrain model: 323,431 land gravity data, 77,389 shipborne marine gravity data, and Sandwell’s v28.1 altimetry-derived global marine gravity model. Then, the geoid was determined with the gravity field model. The methodology for gravimetric geoid determination was based on the remove–compute–restore technique with Helmert’s second method of condensation of topography (Stokes–Helmert scheme). Here, the hybrid Meissl–Molodensky modified spheroidal Stokes kernel was employed to minimize the truncation error under an appropriate combination of different kinds of gravity data. In addition, a high-resolution GSI-DEM on a 0.4 × 0.4 arc-second (10 m) grid, together with the SRTM-DEM on a 7.5 × 11.25 arc-second (250 m) grid, was utilized for precisely applying terrain correction to the regional gravity field model. Consequently, we created a gravimetric geoid model for Japan, consistent with 971 GNSS/leveling geoid heights distributed over the four main islands of Japan with a standard deviation of 5.7 cm, showing a considerable improvement by 2.3 cm over the previous model (JGEOID2008). However, there remain some areas with large discrepancies between the computed and GNSS/leveling geoid heights in northern Japan (Hokkaido), mountainous areas in central Japan, and some coastal regions. Since terrestrial gravity data are especially sparse in these areas, we speculated that the largeness of the geoid discrepancies there could be partly attributed to the insufficient coverage and accuracy of gravity data. The Geospatial Information Authority of Japan has started airborne gravity surveys to be covered over the Japanese Islands, and in future, we plan to develop a geoid model for Japan further accurately by incorporating airborne gravity data to come.

Introduction

The geoid is the equi-geopotential surface best fitted to the global mean sea level. Though classical geoid determination is restricted to regional computations using the astro-geodetic techniques (Drews and Adam 2019), by grace of the enhancement of space geodesy techniques, the two approaches described below are mainly used to determine the figure of the geoid globally and regionally as well in recent decades.

One approach adopts direct measurements of geoidal undulations at points on the Earth’s surface using the combination of spirit leveling and gravity measurements, and GNSS positioning. The former provides the physical orthometric height between the point and geoid, whereas GNSS positioning provides the geometrical ellipsoidal height of the point. The geoid height can be obtained by subtracting the orthometric height from the ellipsoidal height (e.g., Heiskanen and Moritz 1967). This approach is termed as GNSS/leveling, enabling a high-precision measurement of the geoid height over a short distance. However, it is difficult to cover long distances because spirit leveling is a time-consuming process and is susceptible to cumulative errors with increasing distance.

The other is the so-called gravimetric approach, based on boundary value problems on the Earth surface, which can be applied to computing the geoid model from gravity anomaly as the boundary value. Even in the computation of the geoid for a local area, the global integration of gravity anomalies is essentially required, but it can be done in these days using the mixture of data from satellite, airborne, and terrestrial gravimetry. This approach is more reliable than GNSS/leveling for a high-resolution geoid modeling in a regional or global scale, because gravity anomalies are comparatively easy to measure over a wide area without cumulative errors in distance.

Once a precise geoid model becomes available for the target area, the model is commonly used to transform the GNSS-derived ellipsoidal height to the orthometric height. Recent advances in GNSS techniques make it possible to measure ellipsoidal heights with an accuracy of a few centimeters. Therefore, if a well-established geoid model is available, orthometric heights can be derived precisely by GNSS positioning, and accordingly a reference height system can be established via the geoid model, under which GNSS positioning does realize the height determination instead of spirit leveling (and gravity measurements).

Most of the current national height reference datums, including that of Japan, are established and maintained by spirit leveling, the zero levels of which are conventionally defined as the local mean sea level at single or multiple tide-gauge stations. Because of the considerable work required for spirit leveling to cover all the benchmarks, many countries are making increasing efforts to establish/reconstruct national geoid-based height reference systems. For example, Canada constructed a high-precision gravimetric geoid model and established a geoid-based height reference system in 2013 (Véronneau and Huang 2016). New Zealand conducted nation-wide airborne gravity surveys to refine their gravimetric geoid model and shifted to a geoid-based height reference system in 2016 (Amos 2016). The United States launched a project named GRAV-D in 2012 to realize a geoid-based height reference system by 2022 for 10 countries around North America (Smith and Roman 2010).

The Geospatial Information Authority of Japan (GSI) is also considering establishing its own height reference system based on the geoid. To achieve this, it is essential to develop a highly accurate geoid model. With an emphasis on accuracy and precision over land, a regional gravimetric geoid model for Japan has been developed by GSI (e.g., Kuroishi 1995). The latest model, named JGEOID2008, was constructed by applying the generalized Stokes’ integral in a remove–compute–restore (RCR) technique under the assumption of Helmert’s second method of condensation (Stokes–Helmert scheme) (Kuroishi 2009). This model incorporates a global geopotential model (GGM) from Gravity Recovery And Climate Experiment (GRACE) satellites (Tapley et al. 2005), 267,805 land gravity data, 579,186 marine gravity data collected by shipborne gravimetry, and an altimetry-derived marine gravity model, KMS2002 (Andersen et al. 2005). Consequently, JGEOID2008 is consistent with GNSS/leveling geoid heights at 971 benchmarks with a standard deviation of 8 cm.

Following the successful operation of a dedicated satellite gravity mission, called—Gravity field and steady-state Ocean Circulation Explore (GOCE)—during 2009–2013, a number of GGMs with high accuracy and high spatial resolution have been released by research institutes and universities. Odera and Fukuda (2014) first employed a GOCE-based GGM for the determination of the geoid undulation over the Japanese archipelago and demonstrated that utilization of the latest GOCE-based GGM could improve the accuracy of the gravimetric geoid model for Japan.

This study attempted to refine a gravimetric geoid model for Japan, set the same area of 20°–50° N latitude and 120°–150° E longitude as previous models, by utilizing a GOCE-based GGM and an updated regional gravity field model for Japan. We reconstructed the regional gravity field model using updated gravity datasets together with a residual terrain model (RTM) (e.g., Forsberg 1984). The RTM provides the short-wavelength gravity information implied by topographic masses, which is useful to fill-in gaps of terrestrial gravity data as well as to improve the accuracy of gravity data interpolation. We computed the gravimetric geoid model for Japan based on the RCR technique with Stokes–Helmert scheme, which is the same methodology as that used by Kuroishi (2009) and by Odera and Fukuda (2014). Here, we employed the modified Stokes’ integral kernel proposed by Featherstone et al. (1998) instead of that proposed by Meissl (1971) and used in the previous studies. This is expected to enable a more accurate reduction of the truncation error caused by Stokes’ integration over a limited spatial area, as well as to allow a more optimal combination of terrestrial gravity data with the GGM. In addition, we used a high-resolution digital elevation model (DEM) with 0.4 × 0.4 arc-second (10 m) grid spacing, which enables more accurate computation of the terrain correction. Finally, we evaluated the accuracy of the newly developed gravimetric geoid model by comparison with GNSS/leveling geoid heights at 971 benchmarks installed in the four main islands of Japan (Hokkaido, Honshu, Shikoku, and Kyushu).

Data

GNSS/leveling data

The Japanese GNSS/leveling data were reanalyzed in 2014 in response to the large-scale crustal deformations and gravity changes after the 2011 Tohoku-Oki earthquake, of moment magnitude 9.0 (e.g., Ozawa et al. 2011; Matsuo and Heki 2011; Miyahara et al. 2014). The geocentric coordinates on the benchmarks were recomputed using GPS data observed over 6 h, based on the ITRF2008 frame (epoch: May, 2011) for eastern area of Honshu and the ITRF1994 frame (epoch: January, 1997) for other regions. The Helmert orthometric heights on the benchmark were also recomputed using spirit leveling data and gravity data following Kuroishi et al. (2002) and Hiyama et al. (2011). The height difference cause by the epoch difference between GNSS data and spirit leveling data was corrected using a semi-dynamic parameter model developed by GSI (Tanaka et al. 2007). Consequently, the GNSS/leveling geoid heights used here are expected to have an accuracy of better than 3 cm on average (Miyahara et al. 2014). Thus, they are sufficiently accurate to validate a gravimetric geoid model for Japan. A geographical distribution of GNSS/leveling geoid heights is presented in Fig. 1. The GNSS/leveling benchmarks are installed at 170 sites in the Hokkaido region, 565 sites in the Honshu region, 66 sites in the Shikoku region, and 170 sites in the Kyushu region.

Fig. 1

Geographical map of GNSS/leveling geoid heights at 971 benchmarks. The total tide effect is eliminated (tide-free system) following Kuroishi et al. (2002)

Global geopotential model (GGM)

The International Centre for Global Earth Models collects all existing GGMs and makes them publicly available on their web site (http://icgem.gfz-potsdam.de/ICGEM/) (Barthelmes and Köhler 2016). They are provided in the form of spherical harmonic coefficients (Stokes’ coefficients) and their accuracy and spatial resolution vary from model to model.

To choose the best GGM for our geoid modeling, we tested all the three latest GOCE-based satellite-only GGMs, complete to degree and order 300: GOCONS_GCF2_DIR_R6 (Bruinsma et al. 2014), GOCONS_GCF2_TIM_R6 (Brockmann 2014), and GOCO06s (Kvas et al. 2019). Using each GGM with different values for the modification parameters explained in the next sentence, we computed the gravimetric geoid model and calculated the model-derived geoid heights at 971 GNSS/leveling benchmarks, based on the methodology described in the next section and compared them with the GNSS/leveling geoid heights. In the geoid computation, a set of various values for the two modification parameters of Stokes’ integral kernel, namely the integration radius ψ₀ and the degree of Molodensky modification L, are tried to find out an optimum solution (described later in detail).

Table 1 lists the statistics (standard deviation) of the difference between the resulting model-calculated and GNSS/leveling geoid heights. Though the difference is minor, the GOCONS_GCF2_TIM_R6 model shows the best results among the three GGMs. The same performance is confirmed in terms of gravity by Odera and Fukuda (2017). Comparing GOCE-derived gravity anomalies with first-order gravity data over Japan, they concluded that the GOCE time-wise solution (GOCONS_GCF2_TIM model) performed better than the other solutions. Thus, we adopted the GOCONS_GCF2_TIM_R6 model as the foundational GGM.

Table 1 Standard deviation of difference between the computed and GNSS/leveling geoid heights at 971 benchmarks with the use of the GOCONS_GCF2_DIR_R6 model (Bruinsma et al. 2014), the GOCONS_GCF2_TIM_R6 model (Brockmann 2014), and the GOCO06s model (Kvas et al. 2019) (unit: centimeter)

Land gravity data

We collected land gravity data from the databases constructed by research institutes and universities in Japan (Shichi and Yamamoto 1994; Kuroishi 1995; Yamamoto et al. 2011; Honda et al. 2012; Miyakawa et al. 2015; Yahagi et al. 2018). Over the target area are retrieved 355,657 land gravity data. These gravity data, except for several dozens of data, were obtained with relative gravimeters, and different gravimeters and software were used in the respective organizations. Besides, most of the gravity data were collected before the 1990s and some geographical coordinate data recorded (longitude, latitude, and orthometric height) may be inaccurately determined because GNSS had been unavailable at the time of measurement. Therefore, we should suspect that there exist some systematic errors and outliers among the data. To identify and remove gravity data of suspicious quality due to such systematic errors and outliers, we assessed the data quality based on K-nearest neighbor’s prediction method (Saleh et al. 2013) as described later. Further, we took the average of gravity in duplicate within the cell of 10 × 10 m for reducing aliasing effects. In total, there remain 323,431 land gravity data, which are 58,406 more than those used for JGEOID2008. The geographical distribution of the land gravity data is shown as the red dots in Fig. 2.

Fig. 2

Geographical distribution of terrestrial gravity data used in this study. The red dots are the land gravity data (323,431), the dark blue dots the shipborne gravity data (77,389), the light blue dots the altimetry-derived gravity data (2,412,685), and the light green dots EGM2008/RTM fill-in gravity data (109,315), respectively

Marine gravity data

For the oceanic areas, we relied mainly on an altimetry-derived marine gravity model, because most of the shipborne marine gravity data around Japan were obtained between the 1960s and 1980s, and these gravity data are known to suffer from significant systematic errors in some areas (Kuroishi 2001; Kuroishi and Keller 2005).

This study employed the global marine gravity version 28.1 model (V28 model) developed by Sandwell et al. (2014). V28 model was created from altimetry data acquired by the JASON-1, JASON-2, Cryosat-2, and SARAL/Altika satellites. The accuracy of this model is evaluated at approximately 2 mGal on average, approximately 2.5 times better than that of the KMS2002 model used for JGEOID2008. The altimetry-derived marine gravity model covers the whole oceans; however, the quality along the coastlines is still poor because the altimetry data are contaminated by spurious reflections from land. Therefore, we used the shipborne marine gravity data instead for the coastal regions.

For JGEOID2008, 579,186 shipborne marine gravity data are used, most of which are provided by the Bureau Gravimétrique International (BGI). Because the BGI data are particularly old and they appear to contain significant systematic errors, we excluded the BGI data and used the data provided by Koizumi et al. (1994), the Geological Survey of Japan (Miyakawa et al. 2015), and the Japan Oceanography Data Center (https://www.jodc.go.jp/jodcweb/index.html). The observation for these data was mostly conducted between the 1980s and 2000s, and they have already been processed (i.e., drift correction, Eӧtvӧs correction, and tidal correction) by the data providers. We first removed outliers from the datasets according to the same procedures as applied to the land gravity datasets and then extracted the data in the coastal regions within a 40 km from the coast on the four main islands of Japan. In total, there remain 77,389 shipborne marine gravity data. The geographical distribution of the altimetry-derived and shipborne marine gravity data is shown as the light and dark-blue dots, respectively, in Fig. 2.

Fill-in gravity data

The computation of a gravimetric geoid model by the RCR technique implements the regional integration of terrestrial gravity over some limited spatial area. If there exist any substantial data gaps of gravity within the integration radius, the integration will lead to inaccurate results. Therefore, to compute the gravimetric geoid precisely and accurately, one must compensate for the gaps with relevant estimates via interpolation, for example.

We compensated for the missing pieces of gravity information in the data gaps on land in two ways; Least-squares collocation (LSC) interpolation to the gaps on the four main islands and a fill-in method to the gaps for the other areas. In the latter, we computed the fill-in gravity from the Earth Gravitational Model (EGM) 2008 (Pavlis et al. 2012) in the combination with RTM, following a spectral enhancement approach presented by Vu et al. (2019). The EGM2008 free-air gravity anomaly (FGA) was derived using the Stokes’ coefficients up to degree 2190 and order 2159, and the RTM gravity was calculated using a 7.5 × 11.25 arc-second (250 m) grid Shuttle Radar Topography Mission (SRTM) V4.1 DEM (Jarvis et al. 2008) and the Earth 2014 topography model up to degree 2190 and order 2159 (Hirt and Rexer 2015). The geographical distribution of the fill-in gravity data is shown as the light green dots in Fig. 2.

Digital elevation model (DEM)

A DEM is used for terrain correction of observed gravity data, evaluation of primary indirect effect, and computation of RTM. Since these processes are the key to the precise gravimetric computation of the geoid, it is important to use a high-precision and high-resolution DEM. A Japanese DEM is created by GSI from aerial photography and airborne laser altimetry. A 7.5 × 11.25 arc-second (250 m) grid DEM was used in JGEOID2008, and a 1.5 × 2.25 arc-second (50 m) grid DEM was used by Odera and Fukuda (2014). Here, we used a 0.4 × 0.4 arc-second (10 m) grid DEM (GSI-DEM) on the Japanese islands, which has a vertical accuracy of approximately 5 m (http://www.gsi.go.jp/kiban/faq.html). In areas where GSI-DEM is unavailable, we appended a 250-m grid SRTM-DEM and prepared the combined DEM for the area of 15°–55° N latitude and 115°–155° E longitude, with five degrees wider on each side than the target area for determining the gravity field and geoid.

Regional gravity field modeling for Japan

We constructed a regional gravity field model in terms of FGA on a 1 × 1.5 arc-minute (2 km) grid for the area of 20°–50° N latitude and 120°–150° E longitude. The free-air gravity anomalies ($\Delta g_{{\text{FGA}}}$) over the target area were derived from the datasets described in the previous section using the following equation (e.g. Featherstone and Dentith 1997: Kuroishi 2001):

$$\Delta g_{\text{FGA}} = g - \gamma_{Q} + \delta g_{\text{FC}} + \delta g_{\text{AC}}$$

(1)

where $g$ is the observed gravity data, $\gamma_{Q}$ is the normal gravity on the ellipsoid, $\delta g_{\text{FC}}$ is the free-air gravity correction, and $\delta g_{\text{AC}}$ is the atmospheric gravity correction. $\gamma_{Q}$, $\delta g_{\text{FC}}$, and $\delta g_{\text{AC}}$ can be computed as follows:

$$\begin{aligned} \gamma_{Q} & = \gamma_{a} \frac{{1 + k\sin^{2} \phi }}{{\sqrt {1 - e^{2} \sin^{2} \phi } }}, \\ \delta g_{\text{FC}} & = \frac{{2\gamma_{a} }}{a}\left( {1 + f + m - 2f\sin^{2} \phi } \right)H_{P} - \frac{{3\gamma_{a} }}{{a^{2} }}H_{P}^{2} , \\ m & = {{\omega^{2} a^{2} b} \mathord{\left/ {\vphantom {{\omega^{2} a^{2} b} {GM}}} \right. \kern-0pt} {GM}},\quad f = \left( {{{\left( {a - b} \right)} \mathord{\left/ {\vphantom {{\left( {a - b} \right)} a}} \right. \kern-0pt} a}} \right), \\ \delta g_{\text{AC}} & = 0.8658 - 9.727 \times 10^{ - 5} H_{P} + 3.482 \times 10^{ - 9} H_{P}^{2} \;\left( {\text{mGal}} \right) \\ \end{aligned}$$

(2)

where γ_a are the normal gravity at the equator, k is the normal gravity constant, e² is the square of the first numerical eccentricity, $\phi$ is the geographical latitude. a and b are the semi-major and the semi-minor axes, respectively, f is the flattening, $H_{P}$ is the orthometric height, ω is the angular velocity, GM is the product of the gravitational constant and the mass of the Earth. This study used the Geodetic Reference System 1980 (GRS80) (Moritz 2000) as the reference ellipsoid.

The irregularly distributed FGAs were interpolated onto the 2-km grid over the target area by LSC with a second-order Markov model (Kasper, 1971). In this gridding process, the short-wavelength gravity information implied by topographic mass was augmented with RTM, following the method proposed by Omang and Forsberg (2000). The gridding was performed according to the following steps:

• 1. The RTM gravity was computed at each of the data points and of the nodes of the grid from the GSI/SRTM-DEM with reference to the Earth2014 topography model up to degree and order 60 using Eqs. (3) and (4) to follow.

• 2. The RTM gravity was subtracted from the observed FGAs at the point, yielding a RTM-reduced FGA.

• 3. The RTM-reduced FGAs were interpolated onto the nodes by LSC, yielding a RTM-reduced FGA grid.

• 4. The RTM gravity at each node was added back to the RTM-reduced FGA grid, resulting a FGA grid.

The RTM gravity was computed as follows:

$$\delta g_{\text{RTM}} = 2\pi G\rho \left( {H_{P} - \left. {H_{P} } \right|_{0}^{\text{nref}} } \right)$$

(3)

$$\left. {H_{P} } \right|_{0}^{\text{nref}} = \sum\limits_{n = 0}^{\text{nref}} {} \sum\limits_{m = 0}^{n} {\left( {\bar{C}_{nm}^{\text{Topo}} \cos \left( {m\lambda } \right) + \bar{S}_{nm}^{\text{Topo}} \sin \left( {m\lambda } \right)} \right)} \bar{P}_{nm} \left( {\cos \theta } \right)$$

(4)

where G is the gravitational constant, ρ is the density of topographic mass (set here at 2670 kg/m³), $\left. {H_{P} } \right|_{0}^{\text{nref}}$ is the elevation of the smooth reference surface, $\bar{C}_{nm}^{\text{Topo}}$ and $\bar{S}_{nm}^{\text{Topo}}$ are cosine and sine components of spherical harmonic coefficients of the topographic height model, respectively, and nref is the maximum degree for the smooth reference surface.

The RTM reduction intends to take medium- to short-wavelength variations in the gravity field, which mostly come from the topographic masses, out of the gravity data distributed irregularly, and to yield the interpolated RTM-reduced gravity at dense nodes with high accuracy. In doing so, we set nref to 60, for the reduction keeps only the long-wavelength features of topography in the resulting RTM-reduced FGA.

In LSC, we empirically assumed the observation errors of the land and shipborne gravity data to be 1 mGal and 3 mGal, respectively, and used the gravity error of V28 model given in the original dataset as such. These values were used as a weighting function in the LSC interpolation. The FGA grid is mapped in Fig. 3.

Fig. 3

Geographical map of the free-air gravity anomalies in and around Japan. The total tide effect is eliminated (tide-free system)

Gravimetric geoid modeling for Japan

Basic concept

The RCR technique with the Stokes–Helmert scheme was applied to construct a gravimetric geoid model. A detailed description of this scheme can be found in many published papers (e.g., Heiskanen and Moritz 1967; Hofmann-Wellenhof and Moritz 2006). The basic equation is as follows:

$$N = N_{\text{GGM}} + N_{\text{co}}^{\text{res}} + \delta N_{\text{PIDE}}$$

(5)

where N is the geoid height, $N_{\text{GGM}}$ the approximated geoid height derived from GGM, $N_{\text{co}}^{\text{res}}$ the residual cogeoid height to the reference GGM geoid grid, and $\delta N_{\text{PIDE}}$ the primary indirect effect. The parameter $N_{\text{GGM}}$ reflects the long- to medium-wavelength components of the geoid undulation, and $N_{\text{co}}^{\text{res}}$ represents the residual component that the GGM cannot recover because of deficiencies in its spatial resolution. The parameter $\delta N_{\text{PIDE}}$ is a correction term for the conversion of the cogeoid height to the geoid height.

GGM-derived geoid height

The GGM-derived geoid height $N_{GGM}$ with respect to the GRS80 ellipsoid and a geoidal potential value W₀ can be calculated at each computation point from the Stokes’ coefficients using the equation presented in Smith and Roman (2001):

$$\begin{aligned} N_{\text{GGM}} & = \frac{{GM_{1} - GM_{2} }}{{R\gamma_{Q} }} - \frac{{W_{0} - U_{0} }}{{\gamma_{Q} }} \\ & \quad + \frac{{GM_{1} }}{{r\gamma_{Q} }}\sum\limits_{n = 2}^{{n_{\text{max} } }} {\left( {\frac{{a_{1} }}{r}} \right)^{n} } \sum\limits_{m = 0}^{n} {\left[ {\left( {\bar{C}_{nm} - \left( {\frac{{GM_{2} }}{{GM_{1} }}} \right)\left( {\frac{{a_{2} }}{{a_{1} }}} \right)^{n} \bar{C}_{2n}^{*} } \right)\cos \left( {m\lambda } \right) + \bar{S}_{nm} \sin \left( {m\lambda } \right)} \right]} \bar{P}_{nm} \left( {\cos \theta } \right) \\ \end{aligned}$$

(6)

where GM₁ and GM₂ are the product of the gravitational constant and the mass of the Earth given by GGM and GRS80, respectively, R is the mean radius of the Earth (6371 km), W₀ the geopotential value of the geoid, U₀ the normal potential of GRS80 on the ellipsoid, r the geocentric distance of the approximate geoidal surface under the computation point, a₁ the scale factor of GGM, and a₂ the equatorial radius of GRS80. n and m are the degree and order of Stokes’ coefficients, n_max is the maximum degree of GGM, $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are cosine and sine components of Stokes’ coefficients, respectively, $\bar{C}_{2n}^{*}$ is the even zonal Stokes’ coefficient of GRS80, λ is the geographical longitude, $\bar{P}_{nm}$ is nth degree and mth order fully normalized Legendre function, and $\theta$ is the geocentric co-latitude. Here, we used the conventional W₀ value of 62,636,853.4 m²/s² from Sánchez et al. (2016) as recommended in the International Association of Geodesy (IAG) resolution 2015. The r value was recursively estimated following Barthelmes (2013). The maximum degree n_max of the GGM, i.e., the GOCONS_GCF2_TIM_R6 model, is 300 as previously described.

Residual cogeoid height

The residual cogeoid height $N_{\text{co}}^{\text{res}}$ can be computed from the Stokes’ integral of the residual Faye gravity anomaly given in the following equation:

$$N_{\text{co}}^{\text{res}} = \frac{R}{{4\pi \gamma_{Q} }}\int {\int\limits_{0}^{{\sigma_{0} }} {\Delta g_{\text{Faye}}^{\text{res}} \tilde{S}\left( \psi \right)} } {\text{d}}\sigma$$

(7)

where σ₀ is the integration area, $\Delta g_{\text{Faye}}^{\text{res}}$ the residual Faye gravity anomaly, $\tilde{S}\left( \psi \right)$ a modified form of Stokes’ integral kernel, ψ the angular distance between the computation point and the running point, and dσ the surface element.

The calculation of $N_{\text{co}}^{\text{res}}$ was performed by direct numerical integration with the closed Newton–Cotes formula of degree 1. We should note that the integrant $\tilde{S}\left( \psi \right)$ may diverge as the angular distance ψ approaches zero. To avoid this problem, we used, for the area within a radius 5 km from the computation point, the approximation formula of Stokes’ integral kernel with a circular ring system given in Hofmann-Wellenhof and Moritz (2006).

The residual Faye gravity anomaly $\Delta g_{\text{Faye}}^{\text{res}}$ can be computed from the free-air gravity anomaly $\Delta g_{{\text{FGA}}}$ as follows:

$$\Delta g_{\text{Faye}}^{\text{res}} = \Delta g_{{\text{FGA}}} + \delta g_{\text{TC}} + \delta g_{\text{SIDE}} - \Delta g_{\text{GGM}} - \delta g_{\text{aliasing}}$$

(8)

where $\delta g_{\text{TC}}$ is the terrain correction (TC), $\delta g_{\text{SIDE}}$ is the secondary indirect effect (SIDE), $\Delta g_{\text{GGM}}$ is the GGM-derived gravity anomaly, and $\delta g_{\text{aliasing}}$ is the variation of gravity due to the topographic masses at half-wavelengths shorter than the grid resolution, which induce the aliasing error in geoid computation. These terms can be computed as follows (e.g., Wichiencharoen 1982; Forsberg 1984; Smith and Roman 2001):

$$\delta g_{\text{TC}} = \frac{{G\rho R^{2} }}{2}\iint {\frac{{ (H - H_{P})^{2} }}{{l_{0}^{3} }}}{\text{d}}\sigma$$

(9)

$$\delta g_{\text{SIDE}} = 2\gamma_{Q} {{\delta N_{\text{PIDE}} } \mathord{\left/ {\vphantom {{\delta N_{\text{PIDE}} } R}} \right. \kern-0pt} R} \cong - {{2\pi G\rho H_{P}^{2} } \mathord{\left/ {\vphantom {{2\pi G\rho H_{P}^{2} } R}} \right. \kern-0pt} R}$$

(10)

$$\begin{aligned} \Delta g_{\text{GGM}} & = - \frac{{GM_{1} - GM_{2} }}{R^{2}} + \frac{2}{R}\left( {W_{0} - U_{0} } \right) \\ & \quad + \frac{{GM_{1} }}{{r^{2} }}\sum\limits_{n = 2}^{{n_{\text{max} } }} {\left( {\frac{{a_{1} }}{r}} \right)^{n} \left( {n - 1} \right)} \sum\limits_{m = 0}^{n} {\left[ {\left( {\bar{C}_{nm} - \left( {\frac{{GM_{2} }}{{GM_{1} }}} \right)\left( {\frac{{a_{2} }}{{a_{1} }}} \right)^{n} \bar{C}_{2n}^{*} } \right)\cos \left( {m\lambda } \right) + \bar{S}_{nm} \sin \left( {m\lambda } \right)} \right]} \bar{P}_{nm} \left( {\cos \theta } \right) \\ \end{aligned}$$

(11)

$$\begin{array}{*{20}c} {\delta g_{\text{aliasing}} = 2\pi G\rho \left( {H_{P} - H_{P} |_{{2\;{\text{km}}}} } \right)} \\ \end{array}$$

(12)

where H is the orthometric height at the running point, l₀ (= 2Rsin(ψ/2)) the horizontal distance between the computation point and the running point, and $H_{P} |_{{2\;{\text{km}}}}$ the orthometric height of the topographic surface smoothed by averaging the GSI/SRTM-DEM over the cell of 2 × 2 km.

The TC term plays a role in condensing the topographic masses outside the geoid onto an infinite layer on the geoid as well as in applying downward continuation under a planar Earth approximation with a constant topographic density (e.g., Wang et al. 2016). The SIDE term is to upward (or downward) continue the reduced terrestrial gravity on the geoid to the cogeoid. The subtraction of $\delta g_{\text{aliasing}}$ from the FGA is hypothetically intended to smooth the computed gravimetric geoid model with a 2 km spatial resolution and accordingly to reduce the aliasing errors arising from high-frequency undulations (Kuroishi 2001).

The calculation of TC was performed by direct numerical integration, which was an extremely time-consuming task. We calculated the TCs by splitting the integration area into two zones: an inner zone and an outer zone. The TCs in the inner zone and the outer zone were calculated using 10 and 250 m-grid DEMs, respectively. Concerning the integration radii, Hwang et al. (2003) investigated the best radii for the calculation of TCs in Taiwan. They confirmed that the integration radii of 20 km for the inner zone and 200 km for the outer zone were the best to achieve the sufficient accuracy (0.1 mGal) as well as to reduce the calculation time. Because the topographic features of Japan are similar to those of Taiwan, we used the same integration radii as Hwang et al. (2003) for our calculations.

It is known that the TC derived using Eq. (9) may diverge at a point around which the topographic gradient exceeds 45 arc-degrees. In this regard, we replaced it by the one interpolated from nearby (McCubbine et al. 2017). The largest TC of approximately 124 mGal was confirmed in our case around Mt. Fuji (elevation 3776 m), which is the highest mountain in the country.

In theory, the global integration of residual Faye gravity anomaly using the Stokes’ integral kernel yields the residual cogeoid height rather than the residual geoid height. However, because of the incomplete coverage of terrestrial gravity and the limitation of computational capability, it is difficult to execute numerical integrations for the entire Earth. In practice, numerical integrations are performed within a spherical cap of limited spatial extent around each computation point. Such omission of the integration area leads to a so-called “truncation error,” which hampers the precise gravimetric determination of the geoid using Stokes’ formula (e.g. Hagiwara 1972). The effect of the truncation error on the geoid determination can be largely reduced by employing the RCR technique with a modified form of Stokes’ integral kernel.

In this study, we introduced the hybrid Meissl–Molodensky modified spheroidal kernel proposed by Featherstone et al. (1998), which is the so-called FEO kernel. Though both Kuroishi (2009) and Odera and Fukuda (2014) used a modified Stokes’ integral kernel proposed by Meissl (1971), we applied the FEO kernel because of its superiority over Meissl’s modification as described below.The FEO kernel is defined as:

$$\begin{aligned} \widetilde{S}\left( \psi \right) & = S_{p}^{L} \left( \psi \right) - S_{p}^{L} \left( {\psi_{0} } \right), \\ S_{p}^{L} \left( \psi \right) & = S_{p} \left( \psi \right) - \mathop \sum \limits_{k = 2}^{L} \frac{2k + 1}{2}t_{k} P_{k} \left( {\cos \psi } \right), \\ S_{p} \left( \psi \right) & = S\left( \psi \right) - \mathop \sum \limits_{n = 2}^{p} \frac{2n + 1}{n - 1}P_{n} \left( {\cos \psi } \right), \\ S\left( \psi \right) & = \csc \left( {\frac{\psi }{2}} \right) - 6\sin \left( {\frac{\psi }{2}} \right) + 1 - \cos \psi \left( {5 + 3\log \left[ {\sin \left( {\frac{\psi }{2}} \right) + \sin^{2} \left( {\frac{\psi }{2}} \right)} \right]} \right) \\ \end{aligned}$$

(13)

where ψ₀ is the spherical cap distance of the truncation radius, L is the degree of Molodensky modification (Molodensky et al. 1962), which should be set smaller than n_max, p is the degree of Wong and Gore modification (Wong and Gore, 1969), which is set equal to L in practice, t_k is the kth Vanicek and Kleusberg modification coefficient (Vanićek and Kleusberg 1987), P_k is the kth degree Legendre function. The application of Wong and Gore modification replaces the low-degree contributions from terrestrial gravity data with those from the GGM, whereas those of Molodensky modification and Vanićek and Kleusberg modification minimize the L₂ norm of the error kernel for the selected ψ₀ and L (Li and Wang, 2011). Therefore, unlike the Meissl kernel (namely L = 0), the FEO kernel not only reduces the truncation error but also optimally combines terrestrial gravity data with the GGM when the appropriate modification parameters (ψ₀ and L) are selected.

Ideally, these modification parameters should be selected based on the performance of the GGM and the systematic errors of terrestrial gravity data. The performance of the GGM can be examined from the spectral ratio between the signal degree variance and the error degree variance. However, the systematic errors of terrestrial gravity data are generally unknown because they are caused by various factors, such as the usage of different types of gravimeter, different specifications of data processing, and the influence of time-variable gravity signals (e.g., tides, fluid circulations, and crustal deformations).

To determine the appropriate modification parameters, we utilized GNSS/leveling geoid heights as reference. The modification parameters were selected by conducting a grid search for the optimum values that minimize the standard deviation of the difference between the computed and GNSS/leveling geoid heights at 971 benchmarks. Table 1 summarizes the statistics (standard deviation) of the difference between the computed and GNSS/leveling geoid heights with different sets of the two modification parameters. The results show that the modification parameters of ψ₀= 90 km and L = 240 are consistent best with the GNSS/leveling geoid heights. Thus, we chose these modification parameters.

Primary indirect effect (PIDE)

The Stokes–Helmert scheme does not directly provide the geoid, but rather the cogeoid, because of the mathematical condensation of the topographic masses outside the geoid onto an infinite layer on the geoid. This discrepancy between the geoid and cogeoid is called the primary indirect effect (PIDE) and can be computed using DEM. We computed the PIDE based on a planar approximation of the Earth with a constant density of topographic mass presented by Wichiencharoen (1982), using a DEM of 2 km resolution. The DEM of 2 km resolution was obtained by smoothing the GSI/SRTM-DEM over the cell of 2 × 2 km. The equation of PIDE is given as follows:

$$\delta N_{{{\text{PIDE}}}} = -\frac{{\pi G\rho }}{{\gamma _{Q} }}H_{P} ^{2} - \frac{{G\rho R^{2} }}{{6\gamma _{Q} }}\iint {\frac{{H^{3} - H_{P} ^{3} }}{{l_{0} ^{3} }}{\text{d}}\sigma }$$

(14)

The integration radius was set as 200 km, which is the same as in the computation of TC. The largest negative PIDE, of about − 0.7 m, was confirmed around Mt. Kiso (elevation 2956 m) in central Japan.

Assessment of the land and shipborne marine gravity data

The assessment of the land and shipborne marine gravity data was carried out according to the following steps:

• 1. The data were converted into the FGA using Eqs. (1) and (2).

• 2. The RTM gravity at the point was computed from the orthometric height data with the Earth2014 topography model up to degree and order 300 using Eqs. (3) and (4).

• 3. The GGM gravity at the point was computed from the foundational GGM model up to degree and order 300 using Eq. (11).

• 4. The residual Bouguer gravity anomaly (RBA) at the point was computed by subtracting the RTM and GGM gravity from the FGA.

• 5. RBA at the point was predicted by LSC with a second-order Markov model using the 10 nearest neighbor RBAs computed from the data around the point. Note that the RBA computed at the point should not be included in the 10 nearest RBAs in prediction. In LSC, the errors of all RBA data were assumed to be 1 mGal.

• 6. In the case that the difference between the computed and predicted RBA exceeds 3 mGal and 6-σ LSC estimation error, the data were regarded as having suspicious quality and then removed.

Eventually, 1309 land gravity data and 335 shipborne marine gravity data were removed as outlier. Then, we took the average of the gravity data in duplication within the cell of 10 × 10 m for reducing aliasing effects.

Results and discussion

In the Stokes–Helmert scheme, a gravimetric geoid model is composed of three components ($N_{\text{GGM}}$, $N_{\text{co}}^{\text{res}}$, and $\delta N_{\text{PIDE}}$) expressed in Eq. (5). Summing up these three components readily provides us with a gravimetric geoid model. Figure 4 displays the geographical map of the gravimetric geoid model developed. Hereinafter, we refer to the model as JGEOID2019.

Fig. 4

Geographical map of the gravimetric geoid model developed in this study

The statistics of the difference between the JGEOID2019 and GNSS/leveling geoid heights over the four main islands of Japan are compiled in Table 2. For comparison, similar results with EGM2008 and JGEOID2008 are also shown. The geoid heights from EGM2008 were computed using the Stokes’ coefficients up to degree 2190 and order 2159 following Eq. (6) of this paper and Eq. (69) of Barthelmes (2013), in which the GRS80 parameters and the geoidal potential W₀ of 62,636,853.4 m²/s² were used.

Table 2 Statistics of difference between gravimetric geoid models and GNSS/leveling geoid heights at 971 benchmarks for EGM2008, JGEOID2008, and this study

As shown in Table 2, EGM2008 has standard deviations (SDs) of 5.2 to 7.4 cm for the four areas. JGEOID2008 has SDs of 5.5, 5.8, and 5.4 cm for Honshu, Shikoku, and Kyushu, but has a large SD of 10.4 cm in Hokkaido. JGEOID2019 exhibits a significant improvement over these models, having SDs of 3.5 to 4.2 cm for the four areas. As for the whole area (Hokkaido, Honshu, Shukoku, and Kyushu), EGM2008, JGEOID2008, and JGEOID2019 have SDs of 8.0 cm, 8.0 cm, and 5.7 cm, correspondingly. Namely, the new model improves the consistency with the GNSS/leveling geoid heights by 2.3 cm over JGEOID2008. Figure 5 represents the geographical distribution of the geoid difference between the models and the GNSS/leveling at all 971 benchmarks. In each plot, the mean value of geoid differences was subtracted. The histograms of the geoid differences are also shown into a bin width of 5 cm in the lower right inset of the figure. The relative frequencies of the geoid difference within ± 5 cm are 49.99% for EGM2008, 56.33% for JGEOID2008, and 81.15% for JGEOID2019. The ranges of the differences are from − 22.21 to + 32.84 cm, from − 25.95 to + 55.69 cm, and from − 15.03 to + 20.59 cm for EGM2008, JGEOID2008, and JGEOID2019, correspondingly.

Fig. 5

a–c Geographical map of difference between gravimetric geoid models and GNSS/leveling geoid heights. The mean value of each geoid difference was subtracted from each model. The graphs inside each figure are histograms of the geoid difference from the mean value, whose horizontal axis represents the geoid difference and vertical axis shows the frequency of the geoid difference for each bin width of 5 cm. d Geographical map of the topography of Japan, in which the locations with large geoid discrepancies were indicated

According to Fig. 5b for JGEOID2008, the largest geoid difference of + 55.69 cm is observed within the area spread of notoriously large discrepancies in eastern coastal areas of Hokkaido. We speculated that the area with these largest discrepancies was contaminated by the systematic errors at medium wavelengths in the gravity field model used, because of the areal scarcity of relevant marine and land gravity data in the fore-arc zones above the subducting Pacific Plate.

This study, on the other hand, used a GOCE-derived GGM and up-to-dated altimetry-derive marine gravity model, resulting in improving the gravity field modeling at medium wavelengths greatly in the above-mentioned area; however, there still remains comparatively large geoid discrepancies in eastern coastal areas of Hokkaido, to which the quality degradation of altimetry along the coastlines is considered to be attributed. Similar characteristics of large geoid discrepancies are also found around other coastal areas, e.g., Tsugaru Strait located between Hokkaido and Honshu, Tokyo Bay, Seto Inland Sea between western Honshu and Shikoku, Ariake Sea located in western Kyushu, Kagoshima Bay in southern Kyushu. The areas with the large geoid discrepancies are indicated in Fig. 5d. To improve the accuracy of the gravity field and geoid, it is essential to add new gravity data of high quality in such coastal seas obtained by state-of-the-art shipborne or airborne gravimetry.

We also noted that there exist large geoid discrepancies distributed in central Honshu (35°–37° N latitude and 137.5°–139° E longitude). High mountain ranges are outspread there, with the mean elevation of approximately 1000 m including the highest peak of 3776 m at Mt. Fuji in Japan. We considered that there may be two main candidates for the reasons of large geoid discrepancies. One is the data gaps of land gravity. As shown in Fig. 2, the land gravity data distribute quite sparsely there. This is because it is considerably hard to conduct gravity and positioning surveys densely in mountainous regions. To overcome such difficulties, it would be one efficient way to implement airborne gravimetry.

The other candidate is the inaccuracy of downward continuation applied to the land gravity data. This study estimated the downward continuation through the computation of terrain corrections expressed in Eq. (9), in which the planar approximation and the linear relationship between gravity anomaly and topographic elevation are assumed. Such a procedure degrades the computation accuracy of downward continuation especially in high topographic areas, resulting in large geoid errors. One way to tackle with the problem is applying analytical downward continuation (e.g., Moritz 1980). Matsuo and Forsberg (2019) confirmed that such an application to the case of the Colorado Plateau with the mean elevation of approximately 2000 m reduces the errors in the gravimetric geoid computation significantly when compared to the DEM-based downward continuation in planar and linear approximations. Therefore, it is expected to alleviate the geoid discrepancy in central area of Honshu by introducing the improved method of downward continuation and terrain correction. We will work on it in future.

Ideally, the mean difference represents the deviation of the reference (zero-height) surface realized by the vertical datum (for example, the mean sea level of Tokyo Bay for Japan) from the geoid with a geoidal potential W₀. In reality, it is also influenced by the systematic errors inherent to GNSS/leveling data, the long-wavelength errors and the omission error of gravimetric geoid models, and the uncertainty of geodetic parameters such as the geocentric mass constant of the Earth. Therefore, we should recognize that the mean geoid differences obtained here are inclusive of the mixed effects of them.

Looking at the mean difference in geoid heights between gravimetric models and GNSS/leveling data, EGM2008, JGEOID2008, and JGEOID2019 show the offsets of + 2.3 cm, − 22.5 cm, and + 4.8 cm, respectively. Note that the value for JGEOID2008 is quite different from those for EGM2008 and JGEOID2019. This is simply because of the difference in the geoidal potential used in each model; JGEOID2008 employs the value of 62,636,856 m²/s² from the IAG recommendation in 1999 (Bursa 1995), whereas EGM2008 and JGEOID2019 employ 62,636,853.4 m²/s² from the IAG recommendation in 2015 (Sánchez et al. 2016). According to the Bruns formula, the difference in these two W₀ values yields the relative bias of approximately 26.53 cm, almost identical to the mean difference for the entire area between JGEOID2008 and JGEOID2019 (27.28 cm).

It is also noticeable in Table 2 that the mean geoid difference in Hokkaido deviates commonly by approximately 10 cm from that in Honshu. We speculated that the deficit of the land gravity data and the systematic biases of GNSS/leveling data were attributable to such a bias.

As vivid in Fig. 2, the land gravity data distribute especially sparsely in Hokkaido. Hokkaido is one of the challenging areas to conduct land gravity surveys because of its high coverage of mountain forest areas and snowy climate in winter. The estimates of propagation error of EGM2008 indicate the lacks of gravity data. The average propagation error of EGM2008 is approximately 5 mGal in Hokkaido, approximately 1.5 times larger than those in the other areas of Japan. For clarifying the impact of the data deficit on the geoid bias, we need to have more gravity data.

The GNSS/leveling data used would be contaminated by significant biases in Hokkaido. Such a bias is mainly coming from the cumulative errors of orthometric heights determined by spirit leveling. The origin of the vertical datum is located in Tokyo, whose benchmark is connected to a single tide–gauge at Aburatsubo tidal station near Tokyo Bay. Spirit leveling is consecutively conducted in the national networks along the First-order leveling routes starting at the origin and the determination of orthometric heights at benchmarks over the Japanese islands uniquely refers to the datum origin. Imakiire and Hakoiwa (2004) estimated that the cumulative error associated with spirit leveling is expected to reach up to 4.5 cm in Hokkaido and 3.5 cm in Kyushu.

In addition, large measurement errors are also suspected because the connection measurements are quite weak between Hokkaido and Honshu, and between Honshu and Kyushu islands. The connection surveys of spirit leveling between Hokkaido and Honshu were carried out in the Seikan undersea tunnel with a route length of approximately 23.2 km. The spirit leveling in the tunnel was quite harsh because of its dark, humid, and narrow environment. Otaki (2005) reported that a measurement error of 2.5 cm was estimated to occur in the Seikan undersea tunnel. Therefore, the majority of the geoid bias between Hokkaido and Honshu could be attributed to the error in the spirit leveling data.

To evaluate the gravimetric geoid models without the influence of the systematic errors of the GNSS/leveling data, we estimated the planar fitting trend of the geoid difference by three-parameter regulation using a least-squares adjustment (e.g., Chen and Luo 2004). The tilt rate and azimuth of the maximum inclination direction of the fitted plane are listed in Table 3. The planar fitting of the geoid difference shows the azimuth of a nearly north–south direction in EGM2008 and JGEOID2019, and a northeast–southwest direction in JGEOID2008. The tilt rates are 0.14 ppm, 0.08 ppm, and 0.11 ppm for EGM2008, JGEOID2008, and JGEOID2019. Since it is based on a GOCE-based GGM, JGEOID2019 could be expected to demonstrate the best performance at long wavelengths among the geoid models; however, the tilt rate of JGEOID2019 is not smaller than that of JGEOID2008. This may imply that the GNSS/leveling data are not valid in evaluating the long-wavelength accuracy of geoid models because of the existence of long-wavelength features in the systematic errors. We need to check the quality of the GNSS/leveling data in the future. Table 3 also presents standard deviations (SDs) of post-fit residuals of the geoid difference: EGM2008, JGEOID2008, and JGEOID2019 show SDs of 6.56 cm, 6.85 cm, and 4.06 cm, correspondingly. Our model shows the considerably better result than the other two in terms of short wavelengths.

Table 3 Planar fitting of gravimetric geoid models to GNSS/leveling geoid heights at 971 benchmarks for EGM2008, JGEOID2008, and this study

Though we put emphasis on the performance over Japanese islands as mentioned in introduction, that of JGEOID2019 over the oceanic areas is also of interest to us. Sasahara et al. (2008) developed a marine geoid model around Japan and evaluated it in terms of sea surface dynamic topography (SSDT). Since SSDT can be derived from altimetry-derived sea surface heights and the geoid model and be also calculated with the Conductivity–Temperature–Depth (CTD) data, comparison of SSDT shows the consistency between the two results. We will conduct such an evaluation in future.

Finally, we discuss the performance of analysis strategy on the improvement of gravimetric geoid modeling. First, we discuss the modification of Stokes’ integral kernel based on Table 1. We employed the FEO kernel with an integration radius ψ₀ of 90 km and the Molodensky modification degree L of 240, and obtained the SD of geoid difference of 5.71 cm. In the case of the Meissl kernel (namely L = 0) with the integration radius ψ₀ of 110 km, the resulting geoid model matches to the GNSS/leveling geoid heights with an SD of 6.20 cm. It indicates that the introduction of the FEO kernel contributed to an improvement by about 0.5 cm.

Next, we investigate the impact of a high-resolution DEM with 10 m-grid spacing on the gravimetric geoid modeling. If we replace the TC terms by those computed using a 250 m-grid DEM, the resulting geoid model using the FEO kernel gave an SD of the geoid difference of 5.78 cm against the GNSS/leveling geoid heights, and the model using the Meissl kernel an SD of 6.32 cm. In either case, an improvement of approximately 0.1 cm is observable. This may indicate that the inclusion of topographic gravity signals at spatial half-wavelengths shorter than 250 m has little impact on the recovery of the geoid undulations on average. The remaining contribution to the model improvement is considered to mostly come from the incorporation of a GOCE-based GGM and the update of a regional gravity field model for Japan.

Conclusions

The Japanese gravimetric geoid model has been greatly improved by incorporating a GOCE-based GGM, updating the regional gravity field model, and adopting the FEO kernel. The new gravimetric geoid model (JGEOID2019) was found to be consistent with GNSS/leveling geoid heights with an SD of 5.71 cm, which demonstrates an improvement by 2.29 cm over the previous version (JGEOID2008). The existence of areally distributed large discrepancies between the new model and GNSS/leveling geoid heights are noticeable in Hokkaido, mountainous areas in central Japan, and some coastal regions. It may indicate the deficiency of gravity data and the room for improvement in geoid computation for mountainous areas. In the future, we will further improve our modeling by incorporating airborne gravity data, introducing a laterally variable topographic density model, and refining our methodology.

An improved gravimetric geoid model could facilitate the implementation of a geoid-based height reference system in Japan that could be maintained more efficiently than the present system supported by spirit leveling. Furthermore, a geoid-based height reference system can help realize the global unification of the height reference systems within the framework of the Global Geodetic Observing System, as stated in the IAG resolution 2015 and 2019.

Future prospects

For further improving a gravimetric geoid modeling for Japan, the implementation of the following analysis strategies is expected to be effective.

• 1. Focus should be put on the enhancement of both the coverage and the accuracy of gravity data. The distribution of land gravity data remains sparse in some mountainous regions, e.g., Hokkaido and central areas of Honshu. In addition, the accuracy of marine gravity data near the coasts is known to be poorer than in offshore regions. An effective measure to cope with these problems would be to conduct airborne gravimetry surveys. The GSI have started to implement airborne gravimetry over the Japanese islands in 2019 and plan to complete it over the whole country by 2022.

• 2. The introduction of a laterally variable topographic density model will also be of importance. Recently, Sheng et al. (2019) developed a two-dimensional global topographic density model at a 30 × 30 arc-second resolution. Considering such a heterogeneous density distribution would enhance the accuracy of the terrain gravity correction, the indirect effect, and the RTM.

• 3. The application of an alternative methodology for gravimetric geoid modeling, i.e., the UNB approach (e.g., Vanićek et al. 2013; Huang and Véronneau 2013), should be considered and compared in performance. In this approach, terrestrial gravity data are analytically continued onto the cogeoid under the spherical Earth approximation (e.g. Moritz 1980; Heck 2003). This would improve the accuracy of gravimetric geoid models, especially in high topographic areas (Matsuo and Forsberg 2019).

• 4. Finally, instead of a deterministically modified Stokes’ kernel (Featherstone et al. 1998), the adoption of a stochastically modified Stokes’ kernel (Sjöberg 2003) should be evaluated. This approach is advantageous because it minimizes the truncation error and reduces the errors of the GGM and/or terrestrial gravity data as well in a least-squares sense. A GGM and terrestrial gravity data are optimally combined, if their error degree variances are reliably known or can be well estimated, to develop a more robust gravimetric geoid model (e.g., Ellmann 2005).

In the future, we will investigate these strategies and develop a more accurate gravimetric geoid models for Japan.