Bathymetric effect on geoid modeling over the Great Lakes area

Bathymetry data over lake areas are not included in the current and previous NGS (National Geodetic Survey) geoid models. Lake surfaces are simply treated as land surfaces during the modeling regardless of the apparent density difference between water and rock, resulting in artificial masses that distort the model from the actual gravity field and the corresponding geoid surface. In this study, compiled high-resolution bathymetry data provided by National Centers for Environmental Information are used to identify the real volume of water bodies. Under the mass conservation principle, two strategies are deployed to properly account the water body bounded by the mean lake surface and the bathymetry indicated lake floor into the current NGS geoid modeling scheme, where the residual terrain modeling method is used to account for topographic effects. The first strategy condenses water bodies into equivalent rock masses, with the cost of changing the geometrical shape of the water body. The second one keeps the shape of the water body unchanged but replaces the water and rock densities inside each topographical column bounded by the geoid surface and the mean lake surface by an averaged density. Both strategies show up to 1-cm geoid changes when compared with the previous geoid model that does not consider bathymetric information. All three geoid models are evaluated by local GNSS/Leveling benchmarks and multi-year-multi-mission altimetry indicated mean lake surface heights. The results show that both strategies can improve the geoid model precision. And the second strategy yields more realistic results.


Introduction
The NGS (National Geodetic Survey) geoid models are computed based on a modified Molodensky approach (Li et al. 2019).First a high-degree global reference model, such as the EGM2008 model (Pavlis et al. 2012), is used as an approximation of the external long-wavelength gravity field of the Earth.The residual high-frequency gravity field signal is modeled via the residual terrain modeling (RTM) (e.g., Forsberg 1984;Yang et al. 2018;Hirt et al. 2010Hirt et al. , 2019aHirt et al. , 2019b)).Through removing the gravity field contributions of both reference model and RTM from gravity observations, the residual observations are transformed into the residual height anomalies via the Stokes' integral using a modified kernel (Wong and Gore 1969).These three parts are aggregated together for modeling the entire external gravity field of the Earth and hence generating the height anomalies on the Earth's surface.Finally, the classical geoid and quasigeoid transformation is employed to transform height anomalies into geoid heights (Heiskanen and Moritz 1967;Rapp 1997).
Within the NGS geoid modeling scheme, the topographic information is disassembled into two parts: a global reference section defined by the reference field and a residual part that is fluctuating around the reference surface.This decomposition reduces the amplitude of the disturbing potential a magnitude smaller.The resulting residual signals are largely free of systematic signals and hence easier to model by the least-squares collocation (e.g., Heiskanen and Moritz 1967;Moritz 1980).However, the artificial reference topography can break the free space assumption of the exterior space of the real topography.Harmonic corrections (e.g., Forsberg 1984;Denker 2013;Hirt et al. 2019a;Klees et al. 2022Klees et al. , 2023;;Yang et al. 2023) have to be applied to overcome this problem.In addition, the gravity forward modeling of the topography as well as the condensation approaches, which are normally used in the classical Stokes-Helmert geoid modeling scheme (e.g., Vaníček et al. 1999;Huang andVéronneau 2005, 2013;Matsuo and Kuroishi 2020), have to be modified to account for mass anomalies, especially in the areas that have large water bodies totally or partly above the mean sea level, such as the Great Lakes area in North America.
Previous NGS geoid models approximated lake surfaces as land surfaces.Due to the density difference between water and rock, this approximation apparently will cause distortions of the true gravity field by artificially including extra masses into the model.Bathymetry data are needed to correct this approximation.NOAA's NCEI (National Centers for Environmental Information) is responsible for preserving, monitoring, assessing, and providing public access to geophysical data and information.Within NCEI, the Digital Elevation Model (DEM) Team develops integrated bathymetry and topography data products depicting the Earth's surface.It has partnered with other federal agencies and organizations to create an integrated framework for DEMs of U.S. coasts that supports a sustainable, seamless development of DEMs ranging from the near-shore to the global scale.These DEMs are a foundational dataset for tsunami modeling, hazard mitigation, spatial planning, habitat research, coastal change studies, and Earth visualization.Albeit some new higher resolution (less than 10 m) bathymetry data are available in the Great Lakes area (Great Lakes Observing System 2021), the current bathymetry data from NCEI are used in this study to match up the NGS DEM.
In comparison with the single use of the DEM, the combined use of the DEM and bathymetry data can provide more precise mass information for the Earth's topography, especially in areas having large lakes due to the ability of identifying water bodies and rocks.By applying gravity forward modeling either in the spatial domain (e.g., Kuhn and Seitz 2005;Heck and Seitz 2007;Wild-Pfeiffer 2008;Tsoulis et al. 2009;Kuhn and Hirt 2016;Grombein et al. 2013Grombein et al. , 2016) ) or in the spectral domain (e.g., Rummel et al. 1988;Ramillen 2002;Kuhn and Seitz 2005;Tenzer 2005;Hirt andKuhn 2012, 2014;Ince et al. 2020), these mass anomalies can determine more realistic Earth's topographic potential and its derivatives in the scale from local to global.This study is focused on how to put bathymetric information in the current NGS geoid modeling scheme.It is useful to mention that the treatment of these mass anomalies differs with the topographic reductions (e.g., Tziavos et al. 2010) that are adapted to different geoid modeling schemes.Current relevant studies are mainly concentrated on the classical Helmert's method of condensation used in the Stokes-Helmert geoid modeling scheme and the topographic-isostatic reduction used in the geoid modeling scheme based on the window remove-restore technique (Abd-Elmotaal andKühtreiber 1999, 2003).For example, Martinec et al. (1995) investigated the impact of the water masses of Lake Superior on the computation of topographic effects in the Stokes-Helmert scheme as well as its impact on geoid modeling.The results revealed that the correction to the geoid height due to the direct topographic effect on gravity lies between − 1.1 cm and 1.3 cm, the correction due to the primary indirect effect on potential ranges from − 0.24 cm to 0.0 cm, and the correction due to the second indirect effect on gravity is orders of magnitudes smaller that can be neglected.Abd-Elmotaal et al. (2018) estimated the impact of considering the water bodies of Lake Nasser in Egypt on topographic reductions and geoid heights.The results illustrated that the effect of Lake Nasser on both topographic-isostatic gravity anomalies and geoid undulations is limited to the lake area.Over the whole research area, the impact of Lake Nasser on gravity anomalies ranges from 0.0 mGal to 1.4 mGal, while the impact on geoid heights is between 0.3 cm and 4.3 cm.Most recently, Abd-Elmotaal et al. (2020) carried out a deep investigation on the effect of Lake Victoria on computing the direct topographic effects, the isostatic effects based on the Airy-Heiskanen model, and the topographic-isostatic effects by using two different gravity forward modeling approaches, as well as its impact on the geoid determination.The results showed that the two approaches provide similar topographic effects and corrections to geoid heights.The effect of Lake Victoria on topographic-isostatic gravity anomalies behaves locally, reaching up to about 4 mGal.Its impact on geoid heights is regional, with the maximum value of about 28 cm in the lake.The above studies clearly demonstrate that the inclusion of water bodies can cause centimeter-level geoid changes that cannot be ignored in precise geoid modeling over the area with large lakes.However, these water body treatments cannot be directly adopted into the current NGS geoid modeling scheme in which the RTM method is used for the topographic reduction.
There are not too many studies about the RTM considering the bathymetry and other information.The traditional RTM method is based on the assumption that the topography has a homogeneous density distribution, so that the reference topography can be directly obtained by filtering the original detailed topography, yielding the residual terrain model that consists of the geometry and density information for residual topographic masses.However, the homogeneous density assumption breaks if there are other kinds of masses in the topography, such as water and ice masses.In this case, it is inconvenient to apply the traditional RTM method for obtaining a reasonable residual terrain model, requiring us to modify the current method to meet the request of its application in those areas with more complex topography, such as in the Great Lakes area.For this purpose, Hirt (2013) proposed to modify the traditional RTM method by employing the concept of the rock-equivalent topography, allowing to use a single uniform constant density in the RTM forward modeling.By using this approach together with the high-resolution DEM/DBM (Digital Bathymetric Model), he successfully augmented the gravity field based on the EGM2008 model in the coastal zones of Greece and Canada.This approach is easy to implement, but changes the geometrical shape of the original topography with considerable mass displacements.In order to preserve the original topography shape, Schwabe et al. (2014) developed a new RTM method to take ice and water masses into account, in which three reference topographies were individually computed based on the ice-surface, ice-bottom, and lake-bottom topographies.The total RTM effects were obtained by means of summing the contributions from the water, ice, and bedrock.Using this RTM approach, they successfully computed the geoid model in the area of the subglacial Lake Vostok, Antarctica.In principle, this RTM approach can rigorously treat ice, water, and rock masses in the geoid modeling.Thus, the residual terrain model only considering water and rock masses can also be modeled by using the same concept.However, the construction of such a residual terrain model is much more complicated than that of using the rock-equivalent topography, and its implementation needs a significant modification of the currently available software such as the "tc.for" program (Forsberg 1984).Therefore, this concept is not considered in this investigation.Alternatively, we employ another simple concept of averaging the densities of different kinds of masses in the same topographical column, resulting in a new topography with a laterally varying density distribution but the unchanged geometrical shape.The residual terrain model having the same density distribution is then generated by filtering the new topography.To this end, the main objective of this study is to design ad hoc algorithms adapted to the widely used RTM-based geoid/ quasigeoid computation scheme to incorporate water bodies more properly.
The remainder of this paper is organized as follows."Methods" section explains the methods for RTM effect computation, RTM modeling, and geoid computation used in this study."Bathymetry and gravity data" section describes the sources and precision of the bathymetry and gravity data."Results and analysis" section firstly shows RTM effects and geoid models, analyzes the differences caused by considering bathymetry data, and then evaluates geoid models by independent datasets.Finally, some conclusions derived from numerical results are given in "Conclusions" section.

RTM effect computation via tesseroid-based gravity forward modeling
In the RTM method, the gravitational effects induced by the residual topographic masses, which are modeled as the mass layer between the reference and real topographic surfaces, need to be evaluated accurately.The real topographic surface is commonly represented by a high-resolution DEM such as the SRTM (Shuttle Radar Topography Mission) data (Farr et al. 2007).The reference topographic surface represents the long-wavelength component of the real topographic surface.It is normally determined either by a moving average of the real topography or by a high-degree spherical harmonic topography such as the DTM2006.0 model (Pavlis et al. 2007).And the former approach is used in this study.The resulting residual terrain model actually contains real masses with positive densities (i.e., mass surpluses) and filled masses with negative densities (i.e., mass deficits).
The principle of discretized Newtonian integration (e.g., Grombein et al. 2013;Kuhn and Hirt 2016) is employed to calculate the gravitational effects implied by the residual topographic mass distribution (i.e., RTM effects), where the Newton's integral is evaluated by means of dividing the mass distribution into a set of regularly shaped mass elements whose gravitational field can be calculated analytically or numerically.The gravitational effects of the complete mass distribution are then obtained by superposition of individual effects.For the gravitational potential V at the point P , generated by the residual topographic mass distribution ⊂ R 3 , it can be obtained through the approximation of the Newton's integral by where G is the gravitational constant, ρ(Q) is the den- sity at the running integration point Q ∈ , ℓ(P, Q) is the Euclidean distance between the computation point P and the running integration point Q , ρ i Q ′ denotes the den- sity at the running integration point Q ′ inside the ith dis- cretized mass element � i (i = 1, 2, . . ., K ) and is usually set as a constant, ℓ i P, Q ′ denotes the Euclidean distance between the computation point P and the running inte- gration point Q ′ , and V i (P) represents the individual grav- itational potential of i .Using the same principle, the first-and second-order derivatives of the gravitational potential can be derived by means of differentiating Eq. ( 1), and their expressions can be seen, e.g., in Kuhn and Hirt (2016), Lin and Li (2022).
For practical evaluation of the gravitational potential shown in Eq. ( 1) as well as its first-and second-order derivatives, we use the tesseroids (Anderson 1976) as the mass elements in this study because they can be directly generated from the DEMs without any approximations in geocentric spherical coordinates and allow for geodetic and geophysical applications in any scale.In the past decades, various approaches have been developed for precise evaluation of the gravitational effect implied by a tesseroid (e.g., Asgharzadeh et al. 2007;Heck and Seitz 2007;Wild-Pfeiffer 2008;Li et al. 2011;Grombein et al. 2013;Shen and Deng 2016;Uieda et al. 2016;Fukushima 2017Fukushima , 2018;;Marotta and Barzaghi 2017;Lin and Denker 2019;Zhong et al. 2019;Lin et al. 2020;Qiu and Chen 2020).Here we choose the two-dimensional (2D) Double Exponential quadrature approach (Fukushima 2017(Fukushima , 2018)), the three-dimensional (3D) Gauss-Legendre quadrature approach along with the 2D adaptive subdivision technique (Lin and Denker 2019;Lin et al. 2020), and the Taylor series approach up to the zeroth and second order (Heck and Seitz 2007) to precisely compute gravitational effects of residual topographic masses.To take the advantage of each tesseroidal approach and to achieve the balance between accuracy and efficiency, a strategy of combining these three approaches is proposed by Lin and Li (2022) and applied here.It is worth noting that the used combination approach also allows to include a lateral topographical density model.We also want to emphasis that the investigation of the effect of using different tesseroidal approaches or their combinations on RTM effect computation and geoid modeling is out of the (1) scope of this study.Relavent studies can be seen, e.g., in Lin and Li (2022).
The mass redistribution associated with the RTM reduction leads to the situation that some obervation points are inside residual topographic masses.In this case, the reduced gravity field quantities (e.g., reduced gravity anomalies) at these points are inside the masses, where the associated potential function is not harmonic.To solve this problem in approximation, the RTM effects on gravity at the points with h sur < h ref are corrected by the harmonic correction (Forsberg 1984): where h sur and h ref denote the height of the observation point on the Earth's surface and the height of the corresponding point on the reference topographic surface, respectively.After removing the corrected RTM effects of Eq. ( 2) from the observed gravity anomalies, the reduced values are consistent with a harmonically downward continued gravitational potential of the RTM reduced Earth (Klees et al. 2022(Klees et al. , 2023)).

Residual terrain models over the Great Lakes area
The RTM method is applied in the NGS geoid determination over the Great Lakes area.In traditional RTM effect computation by using the "tc.for" program and unclassified DEMs, the positive heights are reserved for rock masses.However, the Great Lakes are filled with water and might be situated totally or partly above the geoid.Thus, the filling of rock masses instead of water masses in the Great Lakes will introduce an obvious error in the computed RTM effects, and hence in the geoid modeling.How to properly treat water masses in the RTM scheme is therefore of high value in pursuing a precise geoid model over the Great Lakes area.
In the RTM scheme, the data are reduced for the gravitational effect of the residual topographic masses between the real Earth's topography and a smoother reference topography.The resolution of the reference topopgraphy usually corresponds to that of the global geopotential model (GGM) used in the remove-compute-restore (RCR) technique in order to avoid subtracting those parts of the signal twice that are already included in the GGM.If the real topography only consists of rock masses with the same density, it is straightfoward to generate the reference topographic surface by a moving average of the real topographic surface and then compute the gravitational effects caused by residual topographic masses.When the real topography contains not only rock masses but also other kinds of masses (e.g., water and ice masses), the reference topographic surface obtained by the moving average approach is not reasonable.This (2) is because the density of the real topography varies both laterally and vertically, making the density distribution of the reference topography quite ambiguous and hard to determine.In this case, if we simply treat the density of the reference topography to be homogeneous, the resulting residual terrain model should be inaccurate.Obviously, such a problem exists in the geoid determination over the Great Lakes area.To reduce the errors caused by the incorrect treatment of water masses in traditional RTM effect computation, two dedicated residual terrain models based on a certain mass transformation are developed and investigated in this study.For brevity, they are denoted as the RTM-Condensed model and the RTM-Denstity model, respectively.In addition, the traditional residual terrain model that directly replaces water masses by rock masses is denoted as the RTM-Base model.
To identify water masses in the Great Lakes, highresolution DEM and DBM are necessary.Both models deviate in the lakes but coincide in the land.Figure 1 illustrates the sketches of real topographic models in the Great Lakes area while considering the cases of water masses being totally (Case 1) or partly (Case 2) above the geoid.Let Land denote the land surface, Lake the lake surface, and ′ the horizontal position described by the latitude ϕ ′ and the longitude ′ .In both cases, the land and water surface heights H s ′ are provided by the DEM and the water floor heights H f ′ are given by the DBM.The lake depths d ′ are then equal to H s ′ − H f ′ .In Fig. 1, if the water masses with the density ρ w are substituted by the rock masses with the density ρ 0 , the real topographic models then become the topographic models corresponding to the RTM-Base scheme.It is clear that extra masses are included in this scheme, which is not a true case.In the following, the RTM-Condensed and RTM-Density schemes will be introduced.
The basic idea of the RTM-Condensed scheme is to trasnform the real topographic model including both rock and water masses into a new topographic model consisting of only rock masses.To achieve this goal, the rock-equivalent principle is applied here, in which water masses are compressed into equivalent rock masses while keeping total masses unchanged.Figure 2 shows the sketches of the topographic models corresponding to the RTM-Condensed scheme.In comparison with the real topographic model, the main change is that the water masses (the blue part in Fig. 1) are replaced by the equivalent rock masses (the orange part in Fig. 2).Since the total masses are presumed to be unchanged and the rock density ρ 0 is larger than the water density ρ w , it is easy to infer that the condensed equivalent surface is lower than the real water surface.Let H ′ s ′ be the height of the condensed equivalent surface at a horizontal position ′ and d ′ ′ be the thickness of the corresponding equivalent rock mass column, in planar approximation, they are calculated by and (3) In the case of the water column being totally above the geoid, all water masses between the lake surface and lake floor are condensed.The formulas of calculating H ′ s ′ and d ′ ′ are those satisfying the condi- Eq. ( 4).In the case of the water column being partly above the geoid, through considering the fact that only the masses above the geoid are needed to be removed or shifted in the Stokes' integral, those water masses below the geoid remain as they are and the masses between the lake surface and the geoid are compressed.The relevant formulas are the ones holding the condition of Eq. ( 4).It should be noted that the replacement of water masses by equivalent rock masses results in a changed geometry with considerable mass displacements.Due to the distance dependency, this has also an impact on the computation of topographic effects.The outcome of this approach is a new topographic model with a homogenous density distribution.The heights of the land surface of this model are the same as those of the initial DEM, while the heights of the lake surface are replaced by the heights of the condensed equivalent surface.Since the masses bounded by the new topographic surface and the geoid have the same rock density ρ 0 , the corresponding reference topographic surface can be directly obtained by a moving average of the new topographic surface, and then the RTM effects are computed accordingly.
In contrast to the RTM-Condensed scheme in which a topographic model with the changed geometrical shape and the homogeneous density distribution is generated, the basic idea of developing the RTM-Density model is to obtain a topographic model having the unchanged geometrical shape and the laterally varying density distribution.There exist land and lakes in the Great Lakes area, implying that the density varies laterally and vertically.
Here we aim at transforming the real 3D density model into a 2D density model in which the density at the horizontal position ′ does not change from the surface to the geoid.To do this, we express the density ρ � ′ , assumed to be varying only laterally, as an average of the actual 3D density ρ r, � ′ along the topographical column of the height H s ′ .In the land area, it is easy to obtain ρ � ′ = ρ 0 because the topographical column only contains rock masses.In the lake area, since the topographical column consists of both water and rock masses, the computation of ρ � ′ in planar approximation follows: In the above equation, the formula corresponding to the first condition is suited for the case that the water column is totally above the geoid, while the one for the second condition is dedicated to the case that the water column is partly above the geoid.Figure 3 depicts the topographic models corresponding to the RTM-Density scheme.It is clear that each topographical column with the vertically varying density is replaced by the same topographical column but with the homogeneous density in the lake area (the colored part in Fig. 3).To summarize, the averaged density is ρ � ′ = ρ 0 for the topographical column satisfying ′ ∈ Land .If ′ ∈ Lake , the averaged (5) ′ is calculated by Eq. ( 5).In comparison with the topographic model for the RTM-Condensed scheme, the topographical geometry does not change in the RTM-Density scheme while remaining total masses unchanged.The outcome of this approach is a new topographic model bounded by the geoid and the real topographic surface described by the initial DEM, however, the density of this model varies laterally.In this case, it is also straightforward to generate the reference topographic surface by means of applying the moving average approach to the real topographic surface, and the computation of RTM effects accounts for the lateral density variations of residual topographic masses (Lin and Li 2022).

Geoid modeling over the Great Lakes area
The geoid model over the Great Lakes area is computed by using exactly the same approach as used in determining the xGEOID19A model (Li et al. 2019).By applying the RCR technique, the geoid model is computed as follows: with and In the above equations, g OBS is the observed grav- ity anomalies over the Great Lakes area; g ATM and ( 6) Bouguer gravity anomaly, γ the mean normal gravity, and H the orthometric height (e.g., Heiskanen and Moritz 1967;Rapp 1997).

Bathymetry and gravity data
The NCEI Great Lakes Bathymetry datasets (NOAA National Geophysical Data Center 1996, 1999a, 1999b, 1999c, 1999d) are used in this study.These datasets were created from 120 years of soundings and more recent multibeam bathymetry measurements.From 1996 to 1999, scientists from NCEI led an effort to compile existing bathymetry data collected by U.S. agencies, the Canadian Hydrographic service, and other partners to develop geophysical maps.Source data for these maps range in age from the early 1900s to the late 1990 and consisted of digital data sounding, digitized sounding points from National Ocean Service (NOS) smooth sheets, nautical chart data, and hand drawn contours.Bathymetry was compiled using the entire array of historical hydrographic soundings collected in support of nautical charting over a 120-year period by the NOAA NOS and its predecessor agency for Great Lakes surveying, the Army Corps of Engineers; see Holcombe et al. (2005).The datasets were converted into a consistent vertical reference frame with the xGEOID20 DEM (Krcmaric 2022), which, itself, is a combined model from Tan-DEM-X (Wessel et al. 2018), MERIT (Yamazaki et al. 2017), and the USGS 3D Elevation Program (3DEP) dataset (United States Geological Survey 2021).Spatial resolution of the bathymetry grid is 3 ′′ .Accuracy for the xGEOID20 DEM component of the bathymetry grid is estimated to be 1 m in flat terrain and 2-3 m in mountainous terrain, while the accuracy of the NCEI bathymetry grids is hard to access because the NOS hydrographic survey database is a historical set of surveys.Recently (after 1965) rigid standards were maintained for survey control.These accuracy standards and requirements are stated in the NOS Hydrographic Surveys Specifications and Deliverables publication, which focus on the minimal criteria needed for both in-house and contract NOS surveys.
The terrestrial gravity dataset used in this study is the same as that used for the computation of the xGEOID19 model (Li et al. 2019).The dataset is made up of terrestrial gravity observations from NGS, NGA (National Geospatial-Intelligence Agency), NRCan (Natural Resources Canada), and DTU (Technical University of Denmark), consisting of 1650169 gravity points.For more details about this dataset, please refer to Li et al. (2019).Since the geoid modeling is limited to the area covering the Great Lakes in this study, only those gravity points located inside the research area are selected for the computation.

Results and analysis
The contribution of considering bathymetric information in the geoid determination over the Great Lakes area will be displayed by means of comparing the RTM effects and geoid models computed by using three different residual terrain models (i.e., RTM-Base, RTM-Condensed, and RTM-Density) and of validating geoid models against independent data such as the ground GNSS/Leveling benchmarks and the multi-year-multi-mission averaged mean lake surface heights measured by satellite altimetry.The inter-comparison of RTM effects and geoid models can give an area-based perspective of the bathymetric effect (i.e., the effect of properly treating the water masses of the Great Lakes instead of incorrectly treating them as rock masses) on geoid undulations, while the external validation can provide a solidate answer of whether the inclusion of bathymetry data improves the geoid model precision or not.For brevity, the geoid models computed by using the RTM-Base, RTM-Condensed, and RTM-Density models are denoted as the Geoid-Base, Geoid-Condensed, and Geoid-Density models, respectively.

RTM effect comparison
The research area is bounded by 40.5 • N ≤ ϕ ≤ 50.5 • N and 266.5 • E ≤ ≤ 285.5 • E , covering all of the five Great Lakes and containing 227764 terrestrial gravity observations that are used for the geoid determination.Figure 4 displays the high-resolution SRTM v4.1 topography (Jarvis et al. 2008) of the area with a slightly larger extent bounded by 40 • N ≤ ϕ ≤ 51 • N and 266 • E ≤ ≤ 286 • E .It shows us that the lake sur- faces are almost flat.Among the five lakes, the surfaces of Lake Superior, Michigan, Huron, and Erie are nearly at the same level with the height of about 180 m.The surface of Lake Ontario is much lower, at the height of about 75 m.The NCEI bathymetry data provide the lake floor heights that are depicted in Fig. 5.It is easy to see that, besides Lake Erie that is totally above the geoid, the other four lakes are partly above the geoid.Among all lakes, Lake Superior is the deepest one.
When computing the RTM effects associated with the RTM-Base scheme, the SRTM v4.1 DEM shown in Fig. 4 is directly used to calculate the reference topography.The water masses in the Great Lakes are obviously treated as the rock masses in this scheme.When computing the RTM effects associated with the RTM-Condensed scheme, the condensed equivalent surface over the lake area is obtained by using the lake surface heights provided by the SRTM v4.1 DEM data and the lake floor heights from the NCEI bathymetry data (see Eqs. 3 and 4).The condensed equivalent surface heights over the lake area and the topographic surface heights over the land area are then merged, yielding a new 3 ′′ resolution DEM that is shown in Fig. 6.Comparing Figs. 4, 5 and 6, the condensed equaivalent surface is lower than the real lake surface and exhibits height variations.When computing the RTM effects associated with the RTM-Density scheme, the lateral topographical density model is calculated by using both DEM and bathymetry data as well as Eq. ( 5). Figure 7 depicits the estimated densities over the whole research area.As expected, the density varies laterally in the lake area and remains as 2.67g cm 3 in the land area.The averaged density for Lake Erie is bigger than the others.This is because Lake Erie is much shallower than the other four lakes, resulting in a higher percentage of rock masses in each topographical column.
Figure 8a shows the RTM gravity anomaly differences between the RTM-Condensed scheme and the RTM-Base scheme, which are computed at 227764 gravity observation points.The differences range from about -6 mGal to 6 mGal with a standard deviation (SD) of about 0.3 mGal.Because two residual terrain models have the same geometry in the land area but differ from each other in the lake area, gravity anomaly differences with the values around zero are observed in the land and larger differences are seen in the lakes and their shores.
Figure 8b presents the RTM gravity anomaly differences between the RTM-Density scheme and the RTM-Base scheme.The differences are in the range of about -0.4 mGal and 3.1 mGal and have a SD of about 0.04 mGal.In comparison with Fig. 8a, the magnitudes of differences shown in Fig. 8b are smaller.This is because the RTM-Density and RTM-Condensed models are developed by different concepts, and they have different geometries and density distributions.When looking at the distribution of differences, it is interesting to see that the differences inside the lakes are almost zero while their values in the shores and land are slightly larger.This phenomenon looks weird because the discrepancies of the used densities for the RTM-Density and RTM-Base schemes are significant in the lakes (see Fig. 7).Considering the fact that the RTM-Density and RTM-Base models have the same real topographic surface, the lake surfaces are almost flat, and the reference topographic surface is a moving average of the real topographic surface, it is easy to infer that the reference topographic surface inside the lakes should be very close to the real topographic surface, resulting in rather small resdiual topographic masses.As a consequence, the different density distributions for the RTM-Density and RTM-Base models in lake areas almost do not influence the computation of RTM effects.In contrast, the reference and real topographic surfaces deviate from each other evidently in the lake shores and land, and hence the different density distributions do affect the RTM effect computation, leading to slightly larger differences with the values mostly within ± 0.06 mGal.
The differences of the RTM gravity anomalies between the RTM-Density scheme and the RTM-Condensed scheme are illustrated in Fig. 8c.Large differences are mainly distributed inside the lakes and along the shores.The minimum, maximum, and SD of the differences are about -5.9 mGal, 7.7 mGal, and 0.30 mGal, respectively.This implies that the two RTM schemes that consider water masses would have different contributions to the geoid modeling.
The restoring RTM height anomalies are computed at the nodes of a 1 ′ × 1 ′ grid that covers the whole research area.Figure 9a shows the differences between the RTM-Condensed scheme and the RTM-Base scheme.The differences mainly appear in the lakes and their surroundings, with the minimum, maximum, and SD of -16 mm, 19 mm, and 1 mm, respectively.More specifically, both negative and positive differences are observed inside the lakes and only positive differences are seen along the lake shores.Figure 9b displays the differences between the RTM-Density scheme and the RTM-Base scheme, with the minimum, maximum, and SD of -1 mm, 18 mm, and 1 mm, respectively.The behavior of these differences is rather different from that shown in Fig. 9a.Significant positive differences are only located along the lake shores, while the differences are close to zero in the other areas.Figure 9a, b suggests that the effect of water masses on the height anomaly are quite different when using the RTM-Density and RTM-Condensed models.For a clear comparison, the height anomaly differences between both RTM schemes are depicted in Fig. 9c.We find that the differences are concentrated on the lake area, with the values ranging from − 14 mm to 27 mm and a SD of 2 mm.

Geoid model comparison
Figure 10a shows the geoid differences between the Geoid-Condensed model and the Geoid-Base model, reaching to about 25 mm over the whole area.Because of the typical behavior of the gravitational potential, the differences decrease slowly from the lake towards the land.Specifically, positive and negative differences close to zero are observed in most of the land area.This is due to the fact that the residual topographic masses corresponding to the RTM-Condensed and RTM-Base models are almost the same in the land area.Significant differences with positive and negative values mainly exist in the Fig. 6 New DEM for the Great Lakes area corresponding to the RTM-Condensed scheme lakes, while the differences on the lake shores are almost positive.This implies that the bathymetric effect on geoid undulations are not limited to the lakes but also propagated into the lake shores.The difference behavior varies among the five lakes.The variation of differences in Lake Superior is the strongest, while the smallest difference variation is found in Lake Erie due to its shallow depth and small lake surface area.The distribution of the geoid differences between the Geoid-Condensed model and the Geoid-Base model acts approximately as a Gaussian distribution.Since the topographic model used in the RTM-Base scheme directly replaces water masses by rock masses, the RTM-Base scheme actually has extra masses than the RTM-Condensed and RTM-Density schemes.
If water masses are properly treated in the RTM-Condensed or RTM-Density scheme, the resulting geoid differences with repsect to the Geoid-Base model should not be distributed in a Gaussian form.One possibility is that the geometrical shape change caused by condensing water masses into equivalent rock masses introduces additional errors into the Geoid-Condensed model.
Figure 10b shows the geoid differences between the Geoid-Density model and the Geoid-Base model.In comparison with Fig. 10a, the behavior of differences looks quite different in this panel.The lake area is dominated by positive differences with the magnitudes less than 1 mm, for example, over the whole Lake Erie and part of the other four lakes.Larger differences are found Fig. 7 Estimated lateral topographical density model for the Great Lakes area corresponding to the RTM-Density scheme along the shore of Lake Superior, in the north-east part of Lake Huron and Lake Michigan, and in the west part of Lake Ontario.Special attention should be paid on the north shore of Lake Superior where significant differences ranging from about 4 mm to 18 mm are observed.As expected, most of the land area is featured by small differences fluctuating around zero.The geoid differences between the Geoid-Density model and the Geoid-Base model over the lake area act as a trend surface, and their distribution is not in the Gaussian form.According to previous explanations, this may imply that water masses are better treated in the RTM-Density scheme.
Finally, the geoid differences between the Geoid-Density model and the Geoid-Condensed model are pictured in Fig. 10c.The differences range from − 19 mm to 34 mm over the whole area, and very small differences are seen in the land area, as expected.Significant positive differences are observed inside the lakes, and negative differences are seen along the lake shores.These evident differences clearly demonstrate that the use of the

Some remarks
Previous numerical analysis showed that the bathymetric effect on computing the RTM effect and geoid model is small compared to the surface area of the Great Lakes.This might be due to the property of the RTM reduction itself, and in the following, we try to explain why this happens.
In the RTM reduction, the magnitudes of RTM effects mainly rely on the total volume of the residual terrain model.The RTM-Condensed scheme condenses water masses into equivalent rock masses, yielding the uneven condensed equivalent surface in the lakes (see Fig. 6).The residual terrain model for this scheme has a different geometrical shape in the lakes and their shores when comparing to the one for the RTM-Base scheme, while the density distrbution for both residual terrain models is the same over the whole area.The differences between the RTM effects for both schemes reflect the bathymetric effect on RTM effect computation and attribute to the geometrical shape difference.In this study, the maximum degree of the reference topographic surface is chosen as 2160 (i.e., the 5 ′ resolution).This makes the residual ter- rain model for each scheme has a small volume, and thus resulting in small magnitudes of RTM effects and bathymetric effect.The RTM-Density scheme takes the water masses of the Great Lakes into account by means of averaging the densities of water and rock masses in each topographical column, resulting in a laterally varying density distribution in the lakes (see Fig. 7).The residual terrain model for this scheme has a different density distribution in the lakes when comparing to the one for the RTM-Base scheme, while the geometrical shape for both residual terrain models is the same over the whole area.Therefore, the bathymetric effect on RTM effect computation attributes to the density difference.Since the reference topographic surface is a moving average of the real topographic surface and the lake surface is flat, it is easy to infer that both topographic surfaces would be close to each other in the lakes.This will result in a quite small volume of the residual terrain model in the lake area for each scheme, yielding small magnitudes of RTM effects and bathymetric effect.The use of a high-resolution reference topographic surface is another reason for generating small RTM effects and bathymetric effect.Upon the above analysis, we may summarize that, when using the RTM reduction for the geoid modeling over the Great Lakes area, the magnitude of the bathymetric effect may depend on (1) the lake depth, (2) the degree of the lake floor surface variation, (3) the lake surface area, and (4) the resolution of the used reference topographic surface.
Since the RTM reduction is not the only topographic reduction used for the geoid modeling, the magnitude of the bathymetric effect might be different if using other topographic reductions, such as the topographicisostatic reduction (Abd-Elmotaal et al. 2020) and the Helmert's method of condensation (Martinec et al. 1995).In comparison with the RTM reduction, both reduction schemes use larger volumes of topographic masses for computing topographic effects.Lin and Li (2022) showed that the magnitudes of the RTM effects obtained by using a 5 ′ resolution reference topograohic surface are much smaller than those of topographic-isostatic effects.Therefore, when using either the Helmert's method of condensation or the topographic-isostatic reduction, the bathymetric effect might also be larger than that for the RTM reduction, especically when the RTM reduction uses a high-resolution reference topographic surface.Which topographic reduction is the best for the geoid modeling over the Great Lakes area would be of interest to be investigated in future studies.

Geoid model validation
Previously, we gave an area-based view of the bathymetric effect on geoid undulations through inter-comparison between the Geoid-Base, Geoid-Condensed, and Geoid-Density models.The resulting geoid differences only show the geometrical distribution of the bathymetric effect as well as its magnitude but cannot prove whether the effect improves the geoid model precision or not.Therefore, an external validation of the computed geoid models by independent data is needed.In the Great Lakes area, two kinds of validation data are available.They are GNSS/Leveling benchmarks located along the lake shores and multi-year-multi-mission averaged altimetry data inside the lakes.Figure 11 displays the distribution of all available validation data, and only part of them are used here.It should be noted that, the precision of used GNSS/Leveling data is difficult to access because they are old and hard to check the year of measurement, especially for the leveling data that were used to define the NAVD88 datum (Zilkoski et al. 1992).Eventhough, the GNSS/Leveling dataset is the only dataset we have to evaluate the geoid model over the lake shore area and should provide meaningful information about the geoid model precision.
Before applying the Stokes' integral to transform residual gravity anomalies into residual height anomalies, the original scattered residual gravity anomalies should be interpolated into grid values.In this study, we employ the LSC method for the gridding, which requires the error information about gravity observations.Because all gravity observations are old and measured at different historical time periods, it is very hard to know their true error information.To investigate whether the geoid model improvement independent of the knowledge of the true error information about gravity observations or not, three gravity observation noise levels are artificially set up, namely the low, middle, and high levels.The Stokes' truncation degree is also not fixed to allow a direct view of the responses in different frequency bands.Because of the existence of NAVD88 errors in the conterminous United States (Li 2018), the geoid models cannot be fairly evaluated by using all GNSS/Leveling benchmarks available in the Great Lakes area.Alternatively, GNSS/Leveling data are separated into state by state, and then are used for geoid validation separately.
Figure 12 shows the results by comparing geoid models with 155 GNSS/Leveling benchmarks located in the State of Minnesota (MN) on the west side of the Great Lakes.Multiple parameters are used to make the comparisons more independent and solid.It is clear that the SDs of geoid differences as a function of Stokes' truncation degree perform similarly.The SDs converge toward the high truncation degree and large observation noise as expected.With the increasing of the truncation degree, the SDs firstly decrease, then reach their smallest values, and finally increase slightly.The optimal truncation degree is located between 700 and 800 for all geoid models and observation noise levels.When the truncation degree is smaller than 600, the Geoid-Base model has a slightly higher precision than the Geoid-Condensed model.This turns to be opposite when the truncation degree is larger than 600.For all observation noise levels and truncation degrees, the precision of the Geoid-Density model is higher than the other two geoid models.A geoid model precision improvement of about 2 mm can be obtained when using the RTM-Density scheme instead of the RTM-Base scheme.This improvement corresponds well to the significant geoid differences on the west shore of Lake Superior shown in Fig. 10b.In general, the results verify that the inclusion of bathymetric information has a positive impact on geoid modeling over Lake Superior.More specifically, the RTM-Density scheme outperforms the RTM-Condensed scheme.
In addition to the geoid model evaluation by using ground GNSS/Leveling benchmarks, the geoid models, in particularly over the lake and marine areas, can be validated by altimetry data.The method for evaluating the geoid model by altimetry data follows Li et al. (2016).Figures 13 and 14 show the evaluation results along Track01 and Track02 passing over Lake Michigan.Track01 is located in the south part of the lake, and Track02 is in its middle part (see Fig. 11).Both figures reveal that the SDs of geoid differences decrease as the truncation degree increases, and finally converge towards the high truncation degree.A geoid model precision improvement of about several centimeters can be obtained by means of selecting an optimal truncation degree.The truncation degree yielding the smallest SDs is between 1200 and 1300 for Track01 and between 650 and 700 for Track02.For all truncation degrees, the SDs for the Geoid-Density and Geoid-Base models coincide quite well, but are about half milimeter smaller than that of the Geoid-Condensed model.The results demonstrate that, on the one hand, the inclusion of bathymetric information has a minor impact on the geoid model over Lake Michigan, on the other hand, the use of the RTM-Condensed scheme has a slightly negative impact on geoid modeling.On the basis of the comparisons among Figs.12, 13 and 14, we see that the bathymetry data help to improve the geoid model precision and the RTM-Density scheme is recommended for the geoid determination over the Great Lakes area.

Conclusions
In the current and previous NGS geoid determination over the Great Lakes area where the RTM method is used to account for topographic effects, the Great Lakes filled with water are incorrectly treated to be filled with rocks (i.e., RTM-Base).This results in extras masses that distort In this investigation, through considering the high-resolution bathymetry data of the Great Lakes provided by NCEI, the water bodies of the Great Lakes bounded by the mean lake surface and the lake floor can be accurately identified.According to bathymetric information and the hypothesis of mass conservation, two different RTM schemes have been developed to properly treat the water bodies in the RTM effect computation.In the first scheme (i.e., RTM-Condensed), water masses are compressed into equivalent rock masses, resulting in a residual terrain model having a different geometrical shape from that for the RTM-Base scheme in the lakes and their shores but the same density distribution over the whole area.In the second scheme (i.e., RTM-Density), the density of each topographical column from the lake or land surface to the geoid is taken as an average of the densities of different kinds of masses included in this column, resulting in a residual terrain model having a different density distribution from that for the RTM-Base scheme in the lakes but the same geometrical shape over the whole area.
Obviously, the difference in the RTM effects computed by the RTM-Condensed scheme and the RTM-Density scheme attributes to the geometrical shape difference and density difference for the residual terrain models.Both RTM schemes are used to compute the geoid models (i.e., Geoid-Condensed and Geoid-Density) over the Great Lakes area, which are then compared with the one (i.e., Geoid-Base) obtained without considering bathymetric information.The area-based inter-comparison between geoid models shows up 1-cm geoid changes caused by the inclusion of bathymetry data.The geoid differences between the Geoid-Condensed and Geoid-Base models look quite different from those between the Geoid-Density and Geoid-Base models.The formers are approximately distributed in a Gaussian form but the latters are not.In order to further prove the significance of considering bathymetric information in the geoid modeling, all three geoid models are validated by local GNSS/Leveling benchmarks and mean lake surface heights averaged from multi-year-multi-mission altimetry data.The validation results show that, when the Stokes' truncation degree is properly selected, both Geoid-Condensed and Geoid-Density models can provide better precisions than the Geoid-Base model regardless of the used gravity observation noise level.Furthermore, the Geoid-Density model looks more realistic than the Geoid-Condensed model.Therefore, the RTM-Density scheme is recommended for future precise geoid modeling over the areas having large lakes if high-resolution bathymetry data are available.We also would like to emphasis that, the density in the land area is assumed to be constant in the current study, which is not the true case.With the availability of a global or regional lateral topographical density model (e.g., Sheng et al. 2019), the lateral density variation in the land area can also be taken into consideration for the geoid modeling over the Great Lakes area.This can be achieved by employing our currently developed software that is able to compute various kinds of topographic effects by using either a constant topographical density model or a lateral topographical density model (Lin and Li 2022).Further investigations will be carried out in our future studies.Current studies are limited to the Great Lakes area, in which there are only water and rock masses.The proposed RTM-Condensed and RTM-Density schemes should be also suited for applications in those areas having more kinds of masses.This is also of interest to be investigated in future works.

Fig. 1
Fig. 1 Sketches of the real topographic model in the case of water masses being a totally or b partly above the geoid.If the water masses with the density ρ w (the blue part) is substituded by the rock masses with the density ρ 0 , the real topographic model then becomes the topographic model corresponding to the RTM-Base scheme

Fig. 2
Fig. 2 Sketches of the topographic model corresponding to the RTM-Condensed scheme in the case of water masses being a totally or b partly above the geoid corrections and the reference gravity anomalies synthesized from the reference model xGeoid16RefA (Li et al. 2019) up to d/o 2160, respectively; g 2161→216000 RTM is the RTM effects on the gravity anomaly associated with different residual terrain models that were introduced before, where the resolution of the real topography and the reference topography is chosen as 3 ′′ and 5 ′ , respectively; S M (ψ) is the modi- fied Stokes' integration kernel; ζ 2→2160 REF and ζ 2161→216000 RTM are the restoring terms from the reference model and the residual terrain model, respectively; N is the classical geoid-quasigeoid separation term with g SB the simple

Fig. 3
Fig. 3 Sketches of the topographic model corresponding to the RTM-Density scheme in the case of water masses being a totally or b partly above the geoid

Fig. 4
Fig. 4 High-resolution SRTM v4.1 DEM of the research area without bathymetry.The lake surface is treated as flat land

Fig. 5
Fig. 5 Bathymetry of the Great Lakes provided by NCEI

Fig. 8
Fig. 8 Differences between the RTM gravity anomalies computed by the RTM-Condensed scheme and the RTM-Base scheme (a), by the RTM-Density scheme and the RTM-Base scheme (b), by the RTM-Density scheme and the RTM-Condensed scheme (c)

Fig. 9
Fig. 9 Differences between the RTM height anomalies computed by the RTM-Condensed scheme and the RTM-Base scheme (a), by the RTM-Density scheme and the RTM-Base scheme (b), by the RTM-Density scheme and the RTM-Condensed scheme (c)

Fig. 10
Fig. 10 Differences between the Geoid-Condensed model and the Geoid-Base model (a), between the Geoid-Density model and the Geoid-Base model (b), between the Geoid-Density model and the Geoid-Condensed model (c)

Fig. 11
Fig. 11 Distribution of all available GNSS/Leveling benchmarks (circles are from U.S. and squares are from Canada) and tracks of multi-year-multi-mission averaged altimetry data in the research area.Here we use part of them for the valiation

Fig. 12
Fig. 12 Standard devations of geoid model differences with the increasing of the Stokes' truncation degree.155 GNSS/Leveling benchmarks located along the shore of Lake Superior close to the state of MN are used for the validation.The geoid models are computed by using three different RTM schemes (RTM-Base: dashed lines; RTM-Condensed: solid lines; RTM-Density: plus symbols) and three gravity observation noise levels