A contribution for the study of RTM effect in height anomalies at two future IHRS stations in Brazil using different approaches, harmonic correction, and global density model

1 Introduction

The definition and realization of the International height reference system (IHRS) have been proposed in the Resolution No. 01/2015 of the International Association of Geodesy (IAG – International Association of Geodesy, Resolution No. 1, 2015). The goal is to provide a global unified high-precision gravity field-related height system and an international standard for the precise determination of physical heights worldwide (Drewes et al. 2016, Ihde et al. 2017, Sanchéz et al. 2021). This is important not only to provide most adequate, unified, and precise altimetric coordinates for the society activities context, but also for scientifical purposes as determination and investigation of global changes in Earth, and for the detection of sea level change. The definition and realizations of the IHRS is an integral part of the Global Geodetic Reference System proposed by the United Nations for Sustainable Development in 2015 (A/RES/69/266).

In this context, since 2011, the Global geodetic observing system (GGOS) Joint working groups (JWG) 0.0.1, 0.0.2, and 0.0.3 have developed works and investigations for an accurate, consistent, and well-defined realization of the IHRS. For example, definition of the zero-level geopotential value (W ₀), adoption of criteria for the network design and establishment of the future stations, evaluation of different strategies for the determination of physical coordinates at the reference stations, proposals for standards, procedures and conventions needed for the realization of the reference system, and the strategy for the implementation of the IHRF (e.g., Sánchez 2012, Dayoub et al. 2012, Sánchez et al. 2014, Cunderlík 2015, IAG – International Association of Geodesy, Resolution No. 1 2015, Sánchez et al. 2015, Sánchez et al. 2016, Ihde et al. 2017, Mäkinen 2017, Sánchez and Sideris 2017, Sánchez et al. 2018, Wang et al. 2018, IAG – International Association of Geodesy, Resolution No. 3 2019, Sánchez et al. 2019, Barzaghi et al. 2020, Sánchez and Barzaghi 2020, Mäkinen 2021, Sanchéz et al. 2021, Wang et al. 2021).

In the context of evaluation of different strategies for the determination of physical coordinates, the Colorado Experiment (Sánchez et al. 2018) has been conducted and performed by GGOS JWG 0.0.2, IAG JWG 2.2.2, IAG Sub-commissions 2.2 (SC 2.2), and Inter-Commission Committee on Theory Joint Study Groups 0.15 (ICCT JSG 0.15). In this experiment, using exactly the same input data provided by NOAA’s National Geodetic Survey and a set of basic standards, 15 research teams proposed their individual solutions obtained from their own procedures and using their own software (Wang et al. 2018, Wang et al. 2021). From the analysis of the solutions, it was observed that the discrepancies between the different solutions are highly correlated with the topography.

According to Sanchéz et al. (2021), for the determination and accuracy assessment of IHRF coordinates, there are three existing resources: high-resolution Global Geopotential Model (GGM), Regional Gravity Field Model (RGFM), and existing physical height systems. In the first two cases, in order to reduce the omission errors and propose an augmentation or spectral extension on the GGM or RGFM solutions (remove-restore procedure), a short-scale gravity forward modeling with the widely used and well-known residual terrain modeling (RTM) technique can be used (Forsberg and Tscherning 1981, Forsberg 1984). This is to obtain the high-frequency effects of the gravity field.

This technique is based on the gravitational attraction of residual topographical masses between a detailed topographical surface and a smoothed reference topographical surface, in which the spectral information is assumed to be already contained in the GGM or RGFM. Both surfaces can be obtained from regional, national, or global digital elevation models (DEM). For the computation of the gravitational potential, the well-known Newton’s volume integral (NI) is used. As examples of works with RTM, investigations over the last two decades can be cited: Omang and Forsberg (2000), Nagy et al. (2000), Seitz and Heck (2001), Heck and Seitz (2007), Hirt and Flury (2008), Hirt et al. (2010), Omang et al. (2012), Tsoulis (2012), Bucha et al. (2016), Uieda et al. (2016), Hirt et al. (2019), Lin and Denker (2019), Yang et al. (2020), and Yang et al. (2022).

In the spatial domain methods context (Hirt and Kuhn 2014), for the practical evaluation of NI, the residual continuous topographic masses are discretized using elementary mass bodies associated with density information. Subsequently, in the forward modeling, the gravitational effects at the computation point are evaluated by superimposing the contributions from all elementary bodies. Currently, there are four body types associated with different approaches for the solution of NI: rectangular prisms with flatted or inclined top and bottom faces (MacMillan 1930, Anderson 1976, Forsberg 1984, Nagy et al. 2000, Heck and Seitz 2007, Tsoulis et al. 2009, Wu and Chen 2016, Jiang et al. 2017), tesseroids (Seitz and Heck 2001, Heck and Seitz 2007, Grombein et al. 2013, Deng et al. 2016, Shen and Deng 2016, Lin et al. 2020), point-masses (Heck and Seitz 2007), and Polyhedron (Okabe 1979, Götze and Lahmeyer 1988, Tsoulis 2012, D’Urso 2014, Benedek et al. 2018, Yang et al. 2020).

Other issues investigated over the last decades in the RTM gravity and gravitational potential computation are the Harmonic correction (HC) and the use of global or regional density models. The HC must be applied when computation points reside inside the smoothed reference topographical surface (Forsberg 1984). In this case, the RTM gravitational potential cannot be considered harmonic and NI cannot be used. This is considered today as a major theoretical problem with the RTM technique (Hirt et al. 2019, Sanchéz et al. 2021). Recent works have investigated the non-harmonicity problem and have proposed numerical methods to HC evaluation in terms of gravity anomaly and geoid height, for example, via regularized analytical downward continuation using Taylor series expansion (Omang et al. 2012, Bucha et al. 2016) or classical condensation method (Forsberg and Tscherning 1981, Forsberg 1984, Omang et al. 2012, Yang et al. 2022). For the specifical case of height anomalies, Omang et al. (2012) proposed a HC from continuation theory, using classical condensation method, assuming the horizontal gradient of potential as equal to zero and using Poisson’s equation.

For the issue of density values, generally the standard Harkness’s crystalline rock constant density value of 2,670 kg m⁻³ is used. However, as indicated by Yang et al. (2018), although this standard value is a reasonable estimate in terrain composed of crystalline rock, 75% of continental surface is underlain by sedimentary rocks. Nowadays, more realistic density values can be obtained from global or regional density models, obtained by means of seismic velocities or mass density. As example of more recent results, one has the CRUST series published by US Geological Survey and the Institute for Geophysics and Planetary Physics at the University of California. For example, the CRUST 1.0 Global Crustal Model (Laske et al. 2013) with 1° × 1° spatial resolution and eight layers: three sediments, three crusts, water, and ice.

According to Sanchéz et al. (2021), considering the learnings from the Colorado experiment, in which different elevation models (SRTM V4.1 and Earth2014) and different parameters (prism size, reference topography resolution, HC, integration radii, and computation approach) have been applied, for the moment, there is no conclusion about how the RTM effect must be computed in IHRS context. Due to the large variety of options to be considered in the RTM computation, Sánchez et al. (2021) indicate that is not possible to establish a standardization at this moment and recommend “to identify the best RTM configuration by contrasting the results with supplementary, independent data, being precise leveling co-located with GNSS the best option”.

Considering that in terms of strategies for the IHRS the RTM computation is an open problem, in this research, contributing for the study of RTM effect in height anomalies at PPTE and BRAZ IHRS future stations in Brazil has been proposed. This is based on the use of different approaches, with HC, and use of the CRUST 1.0 global crustal model to provide mass density values. In order to identify the most accurate configuration, RTM height anomaly values have been compared with height anomalies values derived from Brazilian geodetic system (BGS).

2 Methods

2.1 Computation of the RTM gravitational potential using different approaches

In this research, for the case of a flat-topped rectangular prism with geodetic latitude (φ) and longitude (λ)-dependent density (lateral variation only) ρ and dimensions Δx = x ₂ – x ₁, Δy = y ₂ – y ₁, and Δz = z ₂ – z ₁, with coordinate axes parallel to the edges of the prism (Figure 1(a)), the NI has been evaluated by the analytical solution proposed by Nagy et al. (2000). In this approach, the RTM gravitational potential for a prism i (V _{RTM i} ) is given by:

(1) V RTM i = G ρ ( λ , φ ) ∣ ∣ ∣ ( x P − x ) ( y P − y ) ln ( z P − z + l ) + ( y P − y ) ( z P − z ) ln ( x P − x + l ) + ( x P − x ) ( z P − z ) ln ( y P − y + l ) − ( x P − x ) 2 2 arctan ( y P − y ) ( z P − z ) ( x P − x ) l − ( y P − y ) 2 2 arctan ( x P − x ) ( z P − z ) ( y P − y ) l − ( z P − z ) 2 2 arctan ( x P − x ) ( y P − y ) ( H P − z ) l x 1 i x 2 i y 1 i y 2 i z 2 i z 2 i ,

Figure 1

(a) Geometry of flat-topped rectangular prisms and (b) geometry of tesseroid. Source: Heck and Seitz (2007).

with

(2) l = ( x P − x ) 2 + ( y P − y ) 2 + ( z P − z ) 2 ,

where G is the Newton’s constant of gravitation and l is the Euclidean distance between the computation point P(x, y, z) and the geometrical center Q of the prism (Figure 1a).

An important remark to be taken in account using flat-topped rectangular prisms approach is the Earth curvature effect for the application in large areas. The topocentric cartesian coordinate system attached to the direction of the vertical at the point P is not parallel to edge system of the prism. In this case, a coordinate transformation was necessary as indicated in Heck and Seitz (2007) and a Local Geodetic System has been used.

In the case of tesseroid approach, for the NI evaluation, the approximated solution based on Taylor series expansion (TSE) has been used (Heck and Seitz 2007, Grombein et al. 2013, Deng et al. 2016). This solution is simpler and faster than the solution based on numerical quadrature rules (e.g., Gauss–Legendre quadrature – GLQ – Wild-Pfeiffer 2008, Li et al. 2011, Shen and Deng 2016, Uieda et al. 2016), and as indicated in Lin et al. (2020), and can produce accurate results at low latitude, as for the case of the BRAZ and PPTE stations (between ∼14°S and ∼24°S).

As shown in Figure 1(b), the tesseroid is an elementary body bounded by two spherical surfaces with constant ellipsoidal heights h ₁ and h ₂, two meridional planes with λ ₁ and λ ₂, and two coaxial circular cones with φ₁ and φ₂, (Anderson 1976, Heck and Seitz 2007). Using the solution based on TSE with the Taylor expansion point at the geometrical center Q of the tesseroid (r ₀ , ϕ₀, λ ₀), neglecting terms of order four and higher in Δr, Δφ, and Δλ, and using the band-wise approach proposed by Smith (2002) to take into account the Earth’s ellipticity, the RTM gravitational potential for a tesseroid i can be obtained by (Heck and Seitz 2007)

(3) V RTM i ( r , λ , φ ) = G ρ ΔλΔ φ Δ r K 000 + 1 24 ( K 002 Δ λ 2 + K 020 Δ φ 2 + K 200 Δ r 2 ) + O ( Δ 4 ) ,

where

(4) Δ λ = λ 2 − λ 1 ,

(5) Δ φ = φ 2 − φ 1 ,

(6) Δ r = r 2 − r 1 = ( N 2 + h 2 ) − ( N 1 + h 1 ) ,

and

(7) λ 0 = ( λ 1 + λ 2 ) 2 , φ 0 = ( φ 1 + φ 2 ) 2 , r 0 = ( r 1 + r 2 ) 2 ,

where N ₁ and N ₂ are the prime vertical radius for h ₁ and h ₂, and K ₀₀₀, K ₀₀₂, K ₀₂₀, and K ₂₀₀ are second-order coefficients which have been found in Heck and Seitz (2007).

For the point-mass potential approach, one has the idea to concentrate the total mass of the tesseroid at P ₀. For this, only the zero-degree term in TSE is considered, resulting in equation (8), which has been used in this research. Generally, this approach is used for large integration radii (Heck and Seitz 2007), due to the Newton’s Law and for the simplicity of the equation.

(8) V RTM ( r , λ , φ ) = G ρ ΔλΔ φ Δ r [ K 000 + O ( Δ 2 ) ] .

In the case of polyhedron approach, for the NI evaluation, the analytical solution with rectangular coordinates proposed by Tsoulis (2012) has been used as done by Yang et al. (2020) in the TGF software. For this, the Fortran code polyhedron written by the mentioned author has been accessed from MATLAB environment through MEX-function, as reported by Yang et al. (2020). As in the prism approach, a polyhedron coordinate system is used. Consequently, a coordinate transformation also was necessary to compute the gravitational potential in a local topocentric coordinate system attached to the computation point and avoid the Earth curvature effect (Figure 2(a)).

Figure 2

(a) Geometry of polyhedron and (b) special cases for polyhedron. Source: (b) Yang et al. (2020).

A very important issue which has been considered is the special cases in which there is no analytical solution for NI. As indicated by Yang et al. (2020), this occurs when the tops’ corners are located on different sides of the bottom surface (Figure 2(b)). In order to detect and solve this problem, the authors developed a strategy based on the subdivision of polyhedron in two parts and involving a surface of constant height −1,000 m, as shown in Figure 2(b). This strategy has been implemented in TGF software (Yang et al. 2020) and also used in this research.

An important remark is the approach performance in terms of mass integration radius zone. In previous works, poor results have been verified when the prisms and polyhedrons are located far away from the computation point (e.g., Heck and Seitz 2007, Tsoulis et al. 2009, Wu and Chen 2016, Jiang et al. 2017). This is due to the rounding errors in numerical calculations. On the other hand, using tesseroids approach, better results have been verified when the integration radius is larger (e.g., Heck and Seitz 2007, Grombein et al. 2013).

In view of this, in this research, all approaches have been tested individually in order to identify the integration zone with better performance and the use of the approaches in a more optimized way. For this, for every 1 km interval in integration radius, the comparison with the reference BGS height anomaly value has been calculated. The point-mass approach has also been tested to replace tesseroid approach in order to verify simplifications in the calculations.

Regarding the HC, it must be applied when computation points reside inside the smoothed reference topographical surface. In this research, the traditional mass condensation technique (Forsberg 1984) together with the Poisson’s equation as proposed in Omang et al. (2012) has been used. In this approach, the HC is obtained from continuation theory, considering that the horizontal gradients of gravitational potential equal zero and a constant vertical potential gradient. In terms of gravitational potential, it can be obtained by equation (9):

(9) V HC = − 4 π G ρ ( H P − H ref P ) 2 ,

where H _{ref P} and H _P are the heights at the computation point extracted from smoothed reference surface DEM and detailed surface DEM, respectively.

As previously mentioned, the gravitational potential at the computation point is evaluated by superimposing the contributions from all elementary bodies. Furthermore, to obtain RTM height anomaly and HC in meters from RTM gravitational potential, the Bruns equation has been used. Thus, for both the stations, considering the HC, one has

(10) ζ RTM P = ∑ i = 1 nlin ∑ j = 1 ncol V RTM i j ( r i j , λ i j , φ i j ) γ P + ζ HC , H ref P > H P 0 , H ref P ≤ H P ,

where γ _P is the normal gravity in m/s²; and nlin and ncol are, respectively, the number of lines and columns of the residual DEMs with RTM heights obtained from the subtraction of the detailed DEM by the smoothed reference DEM (H − H _ref). The maximum integration radius for RTM height anomaly calculations in each station is extended up to 210 km according to Sánchez, Barzaghi and Vergos (2019).

Except for the polyhedron approach, the prism, tesseroid, and point-mass approaches have been coded in MATLAB environment. All codes were validated by comparing the results of RTM height anomaly contribution obtained for the Zugspitze area in Germany, with the results obtained by Hirt et al. (2010) for the same area. In Figure 3(b), the result obtained with the prism approach used in this research is presented as an example. In Figure 3(a), the result obtained by Hirt et al. (2010) in the Zugspitze area is presented.

Figure 3

(a) RTM height anomaly contribution obtained by prism approach for the Zugspitze and (b) RTM height anomaly contribution obtained by Hirt et al. (2010) for the same area.

2.2 Study area and datasets

The BRAZ and PPTE IHRS future stations in Brazil are located as shown in Figure 4. Both the stations are part of the International Terrestrial Reference Frame and have ellipsoidal height referenced to BGS and normal heights referenced to Imbituba Brazilian Vertical Datum (IBVD) defined from tide gauge in Imbituba-SC. So, after compatibilization of tide system, reference height anomalies can be obtained and used to evaluate the RTM height anomalies provided by the tests. The topographic reliefs of both the regions around the stations are slightly rough, with variations of just over 1,000 m.

Figure 4

Study areas at (a) BRAZ and (b) PPTE IHRS future stations in Brazil. The red circles delimit the mass integration zones. The red dots are the stations.

For the detailed topographical surface, Multi-Error-Removed Improved Terrain DEMs (MERIT DEMs – Yamazaki et al. 2017) have been used. This is a 90 m spatial resolution quasi-global DEM (90°N–60°S), referenced to EGM96 geoid and WGS84 Terrestrial Reference System, which was developed by removing multiple error components (absolute bias, stripe noise, speckle noise, and tree height bias) from the existing spaceborne DEMs (SRTM3 v2.1 and AW3D-30m v1). The tree canopy signals have been corrected from U-Maryland Landsat forest cover data and NASA Global Forest Height Data and filtering techniques. After the error removal, land areas mapped with 2 m or better vertical accuracy were increased from 39 to 58%. For more detail, refer Yamazaki et al. (2017).

The smoothed reference topographical surfaces have been obtained from down-sampling of MERIT DEMs to 4,000 m, using simple block-means as in Hirt et al. (2019). This was done with the purpose to guarantee spectral consistency with the future scenario of 4 km in spatial resolution for gravity data around both the stations, as recommended by Sánchez, Barzaghi and Vergos (2019). Subsequently, in order to maintain the detailed DEM grid size, the smoothed DEMs have been up-sampled to 90 m spatial resolution using nearest neighbor interpolation. As a result, after subtraction of the detailed DEMs by the smoothed DEMs (H − H _ref), the residual DEMs with RTM heights have been maintained with 90 m spatial resolution in all masses integrations zones around the computation points.

The obtention of density values from CRUST 1.0 Global Crustal Model (Laske et al. 2013) has been done using the MATLAB codes developed by Dr Michael Bevis at Ohio State University. As a result, the code provides density and layer thickness values. An example for a given position around the PPTE station is shown in Figure 5. In this research, only lateral density variations in the layer sediments 1 have been taken into account. In other words, the change in density values in terms of terrestrial crust layer thicknesses (vertical density variations), considering the RTM heights, have not been taking into account.

Figure 5

Example of CRUST 1.0 values at a given point around the PPTE station.

For the calculations over the entire integration zones around the computation points at BRAZ and PPTE stations, grids of density have been generated. Since CRUST 1.0 global grid have spatial resolution of 1°, they have been up-resampled to 90 m using nearest-neighbor interpolation as reported by Yang et al. (2018). For comparison purposes, tests using the Harkness constant density value have also been performed.

In the evaluation step, the RTM height anomalies values have been compared with the reference height anomalies values. For this, BGS ellipsoidal heights (h _BGS) in tide-free system were transformed to the mean tide system. The normal heights referenced to the IBVD are already in the mean-tide system. Thus, the discrepancies Δζ have been calculated as

(11) Δ ζ = ζ RTM − ( h BGS − H IBVD N ) = ζ RTM − ζ BGS ,

where H IBVD N is the normal height referenced to IBVD.

3 Results

After the subtraction of the detailed DEMs by the smoothed reference DEMs, the residual DEMs with 90 m spatial resolution have been obtained for integration zones around the BRAZ and PPTE stations. They are shown in Figure 6(a) and (b), respectively. In Figure 7(a) and (b) are shown the CRUST 1.0 density grids for the layer sediments 1 at BRAZ and PPTE stations, respectively. At BRAZ station, two values of density have been verified, while at PPTE station only one. Such results were already expected due to the low spatial resolution of the CRUST 1.0 model.

Figure 6

Residual DEM around (a) BRAZ and (b) PPTE stations.

Figure 7

CRUST 1.0 density grid for sedimentary 1 layer at (a) BRAZ and (b) PPTE stations.

For the individual evaluation of each approach using CRUST 1.0 densities, as mentioned, at each interval of 1 km in integration radius, the RTM height anomalies have been compared with the reference BGS height anomaly values. Graphics with the discrepancies (Δζ) at BRAZ and PPTE stations are shown in Figures 8 and 9, respectively. It should be noted that with polyhedron approach, the analysis of the discrepancies was performed up to 50 km of integration radius, since results are very close when compared to three other approaches and the computation time is increased dramatically, as previously indicated by other authors.

Figure 8

Differences between RTM and BGS height anomalies (Δζ) on the integration radius intervals of each approach at BRAZ station using densities provided by CRUST 1.0.

Figure 9

Differences between RTM and BGS height anomalies (Δζ) on the integration radius intervals of each approach at PPTE station using densities provided by CRUST 1.0.

The HC for height anomaly has been calculated only at BRAZ station, since at this station the smoothed reference topography surface has been found above detailed surface. The calculated correction value is ζ _HC = –4.53 × 10⁻⁵ m. In view of its order of magnitude, this correction can be disregarded without prejudice to the calculation of the RTM height anomaly at BRAZ station.

As can be seen, at both stations, all approaches have proposed very close results. The prism, tesseroid, and point-mass approaches are in agreement at the 10⁻⁴ m level in most of the results. The differences of these results with respect to the results obtained using the polyhedron approach reach a maximum level of 10⁻⁴ m. It can also be noted that with the increase in the integration radius, the amplitudes of discrepancy variations decrease. At the BRAZ station, after 50 km of integration radius, the amplitude of the variation is around 1 mm. At the PPTE station, the amplitude is smaller, approximately 0.5 mm. This is because the RTM gravitational potentials and height anomalies are inversely proportional to the distance from the computation point (Newton’s law of gravitation).

As mentioned, for comparison purposes, tests using the Harkness constant density value have also been performed. The differences in relation to the results obtained with the CRUST 1.0 densities along the integration radii at both stations, using tesseroid approach as example, are shown in Figure 10. As can be seen, the use of densities provided by CRUST 1.0 proposed improvements in the results at the 10⁻⁴ m level. The average improvements were 0.0005 and 0.0002 m at BRAZ and PPTE stations, respectively.

Figure 10

Differences between RTM height anomalies using CRUST 1.0 densities and Harkness constant density value, using tesseroid approach as example.

In the sequence, the spatial contributions of RTM height anomaly on the maximum integration radius (210 km) at both stations are presented in Figure 11(a) and (b). The results are presented using only the tesseroid approach as example. Table 1 shows the BGS height anomaly values and RTM height anomaly values considering all the contribution of the entire integration area (integration radius of 210 km) using prism, tesseroid, and point-mass approaches, with densities provided by CRUST 1.0 and Harkness constant density, at both the stations. In the same Table 1, the discrepancies in relation to the reference height anomalies are also presented. As can be seen, the RTM height anomaly values, using CRUST 1.0 densities, reached 1.5 mm at the BRAZ station and 1.4 mm at the PPTE station.

Figure 11

Example of RTM height anomaly contributions using tesseroid approach and CRUST 1.0 at (a) BRAZ and (b) PPTE stations.

4 Discussion

In the individual analysis of each approach can be observed a proximity in the results, with discordance at the submillimeter level, not in agreement with previous research works such as Heck and Seitz (2007), Tsoulis et al. (2009), Tenzer et al. (2010), Wu and Chen (2016), and Jiang et al. (2017). This may be related to the characteristics of the topographic relief of the regions. Both are slightly rough, with a variation of just over 1000 m. Furthermore, as can be seen in Figure 6(a) and (b), most RTM heights have values less than 50 m in both stations. This can also be seen in the histogram shown in Figure 12(a) and (b). It should be noted that the RTM height anomaly values at the millimeter level are also related to this.

Figure 12

Histograms of RTM heights at (a) BRAZ and (b) PPTE stations.

Another issue that may have contributed to the proximity of the results is the residual DEM spatial resolution of 90 m for the entire mass integration area, rather than decreasing the spatial resolution with the increase in the integration radius as practiced for example by Yang et al. (2020). In this case, the differences in geometry between prism, tesseroid, and polyhedron may not have been significant. It should be noted that the work by Ferraz and Souza (2021), using a Digital Terrain Model with 2 m spatial resolution for up to 20 km of integration radius, also did not find significant differences between the prism, tesseroid, and point-mass approaches.

Specifically, the proximity of the results obtained from the use of tesseroid and point-mass approaches indicates that, for the study areas, with the dataset here used here, the deviations of the tesseroid mass bodies from point-masses are insignificant, i.e., the approximated solution by TSE, considering only the zero-degree term, is already sufficient (Heck and Seitz 2007). Since point-mass approach is simplest and fastest and has also proposed results with discrepancy at the 10⁻⁵ and 10⁻⁴ m level in relation to prism and polyhedron approaches, it can be adopted as an optimal solution. It should be noted that in this research, no strategy has been used to improve the results with tesseroid approach in the vicinity of the computation point, as rotation or subdivision of them (Li et al. 2011, Grombein et al. 2013, Uieda et al. 2016, Marotta and Barzaghi 2017).

An important remark is that the polyhedron is a prismatic mass body with two inclined triangle tops, as shown in Figure 3a. Furthermore, the corners are coincident with the pixel-centers of the residual DEM. Consequently, the top corners have different height, while the bottom corners are averaged from the four height values extracted from smoothed reference DEM. However, since the RTM heights have little variation in the zones of integration near the computation points (Figure 6(a) and (b)), the geometry changes from prism to polyhedron are small. This may also have contributed to the proximity of the results.

In the context of density values, as can be verified in Figure 8, the use of the CRUST 1.0 model has proposed most realistic values different from the Harkness constant density value. The use of these values has proposed slight improvements at both stations, in agreement with Tziavos and Featherstone (2001), Kuhn (2003), Sjöberg (2004), Tenzer et al. (2015), Chen et al. (2018), Yang et al. (2018), and Lin et al. (2020), which obtained improvements in the modeling of the gravity field using density models. It is obvious that with the increase in the integration radius, the differences in relation to the use of constant density value also increases, since there are accumulations of errors in the summations. It is also important to point out that the slight improvements may be related to the use of the sediment layer 1 density values only. That is, the vertical variation in density was not considered. Therefore, more studies must be realized in order to analyze the use of linear and nonlinear vertical density variation along the elementary mass bodies. Previous works have investigated these strategies, e.g., Wu and Chen (2016), D’Urso and Trotta (2017), Jiang et al. (2017), Chen et al. (2018), Ren et al. (2018), Lin and Denker (2019), and Lin et al. (2020).

Regarding the HC, the value here obtained at BRAZ station in terms of height anomaly, using the approach proposed by Omang et al. (2012), can be considered as insignificant, since its order of magnitude is 10⁻⁵ m. As can be seen in equation (10), the HC value is directly proportional to the square of the RTM height and inversely proportional to γ_P. Thus, the smaller the RTM height is, smaller is the HC value. In the case of BRAZ station, the RTM height at computation point is −9.0836 m, which has contributed for the smaller HC value. This result is in agreement with those obtained by Yang et al. (2022) in terms of HC for RTM geoid height. From the obtained results, these authors have indicated that HC is more significant for regional geoid determination over mountain areas, achieving cm-level. Consequently, HC could be ignored over ∼99% of continental areas. Also, according to Klees et al. (2022), in general, HC for height anomaly is recommended only in areas of strong topographical variations.

5 Conclusion and future works

For the dataset used here, at both stations, the point-mass approach has been indicated as an optimal solution. This is because it proposed results with maximum discrepancy at 10⁻⁴ m level in relation to the other three approaches, it is less complex, and has a shorter processing time.

For the HC, the traditional mass condensation technique and Poisson’s equation with horizontal gradients equal to zero idea have been used. This correction was only applicable at BRAZ station and proposed a slight improvement in the results, at 10⁻⁵m level. Therefore, it can be disregarded without prejudice.

Regarding the use of densities provided by the CRUST 1.0 model, improvements at the mm level were verified in relation to the use of the Harkness constant density. The low improvement may be related to the disregard of the vertical variation in the density values.

As future works, one has to expand the studies to other future IHRS stations in Brazil and the contribution to the study of RTM effects in gravity disturbance using different approaches and HC. In addition, tests will be performed using Digital Terrain Models obtained from topographic maps at 1:50,000 scales provided by the Brazilian Institute of Geography and Statistics. For density, values will be tested using two options: (i) the geological map of Brazil available by the Brazilian Institute of Geography and Statistics; and (ii) the laterally varying topographical density model with 900 m spatial resolution developed by Sheng et al. (2018), at University of New Brunswick (UNB_TopoDensT_2v01).