Use of GNSS and a refined GGM (XGM2019e) for determining normal heights in the Imbituba Brazilian Vertical Datum and in the International Height Reference System

: This paper aimed to evaluate the use of Global Navigation Satellite System (GNSS) ellipsoidal heights in conjunction with height anomalies provided by Global Geopotential Model (GGM) XGM2019e, refined by Residual Terrain Modelling (RTM) technique, to obtain normal heights in Brazil, referred to the Imbituba Brazilian Vertical Datum (IBVD) and the International Height Reference System (IHRS). For this purpose, a local modelling approach has been analyzed in contrast to the national modeling one on the reference geopotential value. For this, a methodology based on geopotential space was adapted. In the local modeling, two study subregions were defined using the spatial clustering analysis of IBVD and GGM/RTM height anomalies differences outliers. The parameters have been estimated using three different configurations. In the parameters validation step, Root Mean Square Errors (RMSE) of the discrepancies between transformed and Brazilian official normal heights were calculated. In both subregions more accurate results have been obtained with the local modeling. In the SP1 subregion the accuracy increased tenfold (0.97m to 0.10m) and SP2 improved from 0.39m to 0.17m. For the linkage to the future realization of IHRS, the accuracy analysis was not possible. However, discrepancies between calculated normal heights and Brazilian official normal heights have been analyzed.


Introduction
In the wake of a more global view of the Earth's intrinsic processes related to Geodesy, Drewes (2006) stated that the classic goals of Geodetic Sciences were added to others that place Geodesy at the center of the discussion on the challenge of living on a dynamic planet. This more recent interpretation of Geodesy's role as an avantgarde science brought about the need not only for measurement, but also for an understanding of physical entities directly linked to the way the Earth's surface is organized and modified over time. This more global view considers the whole as the Earth System, needing continuous and accurate monitoring over time.
In this context, initiatives promoted by the International Association of Geodesy (IAG) such as the creation of the Global Geodetic Observing System (GGOS) and the search for a global connection of geodetic observations, confirm this purpose. In terms of global unification, the state of the art also points to the definition of a Height Reference System (HRS), which was indicated, in terms of definition and implementation, by IAG Resolution No. 01 of 2015 (Ihde et al. 2017). This is the IHRS and it is future realization International Height Reference Frame (IHRF).
Regionally, as a first step towards the realization of a unified vertical frame, according to de Freitas (2015), the Working Group III of Geodetic Reference System for the Americas (SIRGAS), has operated in the area for more than two decades, seeking to standardize rules, procedures and activities among Latin American countries, aiming to meet international requirements.
In the case of Brazil, after the Readjustment of Altimetric Network with Geopotential Numbers, carried out in 2018 (REALT-2018 -Reajustamento da Rede Altimétrica com Números Geopotenciais in Portuguese), normal heights have been adopted for the High Precision Altimetric Network (HPAN) of the Brazilian Geodetic System (BGS). According to IBGE (2019), this type of height was adopted for being the most adequate to modern concepts and scientific methods and it is in line with the resolutions established by the IAG with aim to define the IHRS. Additionally, the adjustment was carried out using geopotential numbers, in order to prepare the HPAN for future actions under the IHRS; and new spirit leveling and gravimetry observations, as well as correction of eventually detected inconsistencies have been included. However, the new normal heights remain referenced to the Brazilian Vertical Datums (BVD) of Imbituba and Santana, current in Brazil. These datums, in each case, were defined from a single value of the Mean Sea Level, which was classified with data collected in only one tide gauge station, therefore, of a local approach. Thus, given the impossibility of connecting the two leveling networks, it appears that there are two vertical datums in the country. Also, according to IBGE (2019), it is targeted that the future realization of the vertical component of the BGS, in a few years, will already be related to the IHRS/IHRF.
To obtain the normal height value referred to any of the BVD, at a given point in the Brazilian territory, it is possible to carry out spirit levelling from a benchmark (RN -"Referência de Nível" in Portuguese) station and interpolate the real gravity value from Brazilian gravimetric network data. However, with the expressive technological development experienced in recent decades, new methods for obtaining altimetric data have gained projection, supported by the improvement of related observations and techniques. In this context innumerable researches have proposed the use of GNSS ellipsoidal heights in conjunction with data provided by Global Geopotential Model (GGM), which can be refined by Residual Terrain Modelling (RTM) technique using Digital Terrain Models (DTM) data, to obtain physical heights. For example, Hirt, Featherstone, and Marti (2010)like EGM2008, is not capable of representing the high-frequency components of Earth's gravity field. This is known as the omission error. In mountainous terrain, omission errors in EGM2008, even when expanded to degree 2,190, may reach amplitudes of 10 cm and more for height anomalies. The present paper proposes the utilisation of high-resolution residual terrain model (RTM combined residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data, and Castro, Jamur, and de Freitas (2013)considering their distribution and variations in density, have a measurable influence on Functionals of the Gravity Field (FGF explored of the residual topography effect on some functionals of the gravity field.
In the context of normal heights obtained from GNSS/GGM/RTM integration, there is the advantage of avoiding the costly and time-consuming spirit levelling procedure from some benchmark, especially in regions with lower station density. However, there are two issues to consider. The first is the normal height accuracy. While GNSS ellipsoidal height accuracy can be obtained at millimeter level or better, the higher order modern GGM height anomalies, even if refined using DTM, are accurate in the order of magnitude of centimeter or decimeter. The second is the question of the datum, associated with the GGM/RTM zero-height surface, which has a global definition approach.
For the first question, the solution is to limit the use of such a methodology to certain types of applications in which the orders of magnitude of the accuracies may be permissible. In the second question, the solution is, somehow, to estimate a transformation parameter or set of transformation parameters between Height Reference Frames (HRFs), and, therefore, establish the linkage of the normal height obtained from GNSS/GGM/RTM to one of the BVD. For this purpose, some methods have been proposed over the last years. For example, the methodology based on geopotential space as indicated in Ihde et al. (2017) and the geodetic boundary value problem (GBVP) method proposed initially by Rummel and Teunissen (1988).
Given the above, this research aimed to evaluate the use of GNSS ellipsoidal height in conjunction with height anomalies provided by a GGM, refined by RTM technique, for the determination of normal heights in Brazil, referred to the IBVD and to the future IHRS. For the linkage, an adaptation of the method based on geopotential space indicated in Ihde et al. (2017) has been used. In the case of linkage at the IBVD, local estimated parameters and the one obtained by Sánchez and Sideris (2017) at national level have been analyzed. In the second case, the IHRS geopotential value at the geoid W 0 has been used for the linkage.
It is noteworthy that, according to IBGE (2019), in Brazil, six stations belonging to the Brazilian Network for Continuous Monitoring of GNSS Systems (PPTE, IMBT, CUIB, MABA, CEFT and BRAZ) have been pre-selected to integrate the IHRF, being the first in the study region of this research and within the set of stations used in this research. In this sense, studies are necessary about the densification of gravimetric and the linkage of vertical datums. This paper aims to contribute to the second topic.

Data and Study Region
Normal heights obtained from GNSS/GGM/RTM are referred to a zero-height surface associated with the zerolevel geopotential W 0 value of the GGM/RTM solution used. In order to link such heights at the BVD, as mentioned, a transformation between HRFs must be performed (Rummel and Teunissen 1988). As will be presented in item 2.2, in the methodology indicated in Ihde et al. (2017), the zero-level geopotential value in a local datum is estimated as parameter. For this purpose, the differences between BGS and GGM/RTM height anomalies are used, considering a set of points. A given BGS height anomaly, at a given RN/Global Positioning System (GPS) connection station, can be obtained subtracting the BGS ellipsoidal height and the BGS normal height, referred to one of the BVDs. On the other hand, a GGM/RTM height anomaly can be obtained from the International Centre for Global Earth Models (ICGEM, http://icgem.gfz-potsdam.de/, last access: 6 August 2021) (Ince et al. 2019) and calculation of the RTM effect using one or more DTMs. It should be noted that currently there are several options for obtaining global DTMs for free on the web.
Considering all available RN/GPS connection stations in the Brazilian territory at the time of the study, Sánchez and Sideris (2017) obtained the mean (W 0 ) values for IBVD and Santana Brazilian Vertical Datum (SBVD), whose values were 62,636,849.61±0.18 m 2 /s 2 and 62,636,852.93 ± 1.51 m 2 /s 2 respectively. From these values, it is possible to carry out the transformation associated with each datum for the entire Brazilian territory. However, in this research, with the purpose of investigating the problem of transformation at the local level, the value of (W 0 ) was calculated considering a smaller study region, using the subset of available RN/GPS data. The defined study region comprises the state of São Paulo, an administrative region with 102 RN/GPS stations distributed throughout the entire territory. In the Figure 1 it is shown the spatial distribution of the RN/GPS stations. All benchmarks in this region have normal height referred to IBVD, materialized through the station code "4X", which is the fundamental benchmark of the IBVD. For the experiments, all 102 stations have been used, being a part separated to estimate the transformation parameter (zero-level geopotential value in local datum), and the complementary part for it is evaluation, as will be seen below. All data were obtained free of charge from the BGS Geodetic Database (BGD).
For linking the normal heights determined by GNSS/GGM/RTM to the future IHRS, a transformation between HRFs is also required. However, by using the same transformation strategy for the previous case, an estimate of W 0 is not necessary. This is because the IHRS zero-level geopotential value can be directly used. In this case, only the GGM/RTM height anomalies have been used to calculate the geopotential value at the point. It is noteworthy that, for both linkage cases, to compare the values of normal heights obtained from GNSS/GGM/RTM, the same RN/GPS stations set in the study region has been used.
During the process of choosing the GGM, the study developed by Nicacio, Dalazoana, and Freitas (2018) was taken as reference, which verified the employability of some of the most recent GGMs for the extraction of normal-orthometric heights in Brazilian territory. It is noteworthy that, this research was carried out in 2017, when normal-orthometric heights were in use. According to IBGE (2019), the differences between previous and current HPAN heights in the study region vary almost entirely between -0.01 m and 0.01 m. Therefore, the results obtained by Nicacio, Dalazoana, and Freitas (2018) remain valid after REALT-2018. As GGM XGM2016 (Pail et al. 2018) had the best performance and has been updated to XGM2019 (Zingerle et al. 2020), this has been selected for the development of this research. It is noteworthy that, this GGM is more accurate than its predecessor, by adding topographically derived gravimetric information over the terrain. The ICGEM calculation tool was used with data from the GGM XGM2019e up to degree and order 2190. The zero-degree term was selected. The reference ellipsoid used was the GRS80 and the mean tide system was selected.
To minimize the omission error in the GGM data and in order to increase the spectral resolution, the topography data can be used. This information is extracted from DTMs. As mentioned, this procedure proposes a refinement in the GGM height anomalies. For this, the RTM technique has been applied, by using the prisms method, with a radius of approximately 200 km, in agreement with Hirt, Featherstone, and Marti (2010). It should be pointed out that harmonic corrections have been not taken in account. For a given point P, after calculating the residual height anomaly (ζ RTM P ) by the RTM technique using a grid of rectangular prisms and planar approximation (MacMillan 1930) (Heck and Seitz 2007)the topographic surface of the Earth is often divided with respect to geographical or map-grid lines, and the topographic heights are averaged over the respective grid elements. The bodies bounded by surfaces of constant (ellipsoidal, the refinement of the GGM height anomaly ( GGM P ζ ) can be done by: (1) is the refined GGM height anomaly.
The selected DTMs were those that had the best cost/benefit ratio based on factors such as accuracy, computational cost and resolution in previously tests performed by the authors. It was chosen the ETOPO1 (Amante and Eakins 2009) model with 1' spatial resolution and the SRTM15_PLUS (Tozer et al. 2019) 15" spatial resolution, obtained, respectively, from the ICGEM and the ERDDAP data server of the National Oceanic and Atmospheric Administration. For the application of the RTM technique, the ETOPO1 model has been used as reference surface, while SRTM15_PLUS model has been used as detailed surface. Both models cover oceanic areas derived from bathymetric data.
In the data integration, the differences between permanent tide effect have been taken in account. While the BGS ellipsoidal heights are referred to tide free system, the BGS normal heights and DTM heights are referred to mean tide system and zero tide system, respectively. For the compatibilization, the mean tide system has been adopted and the Rapp (1989) formulation has been used, as recommended by the IAG in the conventions for the adoption of the IHRS (Ihde et al. 2017).

Strategies for normal heights linkage
After defining the study region, obtaining and performing the compatibilization of all data, the linkage of the GNSS/GGM/RTM normal heights to IBVD and to the future IHRS were carried out. For this, the methodology presented by Ihde et al. (2017) to define gravitational potential values for regional datums and their respective differences from a global reference has been adapted. This methodology is based on obtaining zero-level geopotential value (W 0k ) of a Local Vertical Datum (LVD) k, and it is difference (∆W 0k ) in relation to the zero-level geopotential value (W 0 ) of a global HRS. Using the HRS unification principle using geopotential source values (Ihde et al. 2017), quantities related to a LVD and a LVD can be obtained. The proposed adaptation in this research is to use the value of directly W 0k to obtain the normal height already referenced to LVD.
As previously mentioned, for linkage to IBVD, besides analyzing the value of _ 0 IBVD SS W obtained by Sánchez and Sideris (2017) for the entire Brazilian region, a more adjusted local model for the study region was also searched. Thus, it was necessary to obtain the mean value of 0 IBVD W from the height anomaly data subset.
In this context, in a given P point, BGS RN/GPS station, the anomalous potential can be obtained from the GGM ( GGM P T ), and the h P ellipsoidal height can be obtained from BGS/BDG, referenced to SIRGAS2000 and transformed from tide-free system to the mean tide system. With these quantities available and knowing the zerolevel spheropotential (U 0 ) of the adopted Reference Ellipsoid (GRS80) (Moritz 1980), the geopotential value for point P (W P ) can be calculated by: where U P is the spheropotential value for point P; On the other hand, the same value of W P is also related to the 0 IBVD W value and to the geopotential number value of point P referred to IBVD ( IBVD P C ) by: or yet: From the combination of Equations (2) and (4), a value of 0 IBVD W can be obtained for point P by: or yet: where N IBVD P H is the normal height referred to the IBVD at point P, obtained from BGS/BGD; and ζ BGS P is BGS the height anomaly.
Knowing the values of , 0 IBVD P W for a set of n points, the mean value of the geopotentials can be assumed as the local benchmark: Finally, the normal height obtained from GNSS/GGM/RTM, at a given point Q, referred to IBVD, is: where W Q can be calculated with the form of Equation (2) with the data of point Q. It is noteworthy that, in this strategy, the quantity 0 IBVD W is considered as a transformation parameter between the HRFs.
In the second approach of this research, for the case of linkage to the future IHRS, the same functional mathematical modeling has been used. The difference, however, is that there is no 0 W parameter to be estimated, since IHRS zero-level geopotential value ( 0 IHRS W ) has already been defined by IAG Resolution 01/2015. In this case, the normal height can be obtained by:

Validation of linkage and detection of outliers
In the context of linkage to IBVD, in order to estimate the local parameter 0 IBVD W , from Equation (8), a BGS height anomalies dataset is required. Subsequently, after estimating the parameter, it must be evaluated, as well as the _ 0 IBVD SS W parameter calculated by Sánchez and Sideris (2017). Therefore, from the 102 RN/GPS stations in the study area, 68 stations have been separated to estimate the parameters, being called reference points (r); and 34 stations have been separated to evaluation of the parameter, being called check points (t). The ratio of 2/3 and 1/3 has been defined empirically. The designation of the points in each one-third has been made looking for a homogeneous spatial distribution throughout the study region, but without a quantitative criterion.
According to the three one-thirds procedure, an initial one-third configuration of points was called Alpha (A). Subsequently, two other adopted one-thirds were named Bravo (B) and Charlie (C), using exchanges between the reference and check points of the one-third part A. In this way, all stations were, in some of the three one-thirds, considered reference and check points. Thus, the question to test the use of RN/GPS stations for both purposes, allowing transformation parameters in each one-third to be obtained, validated from all stations and compared. was verified. First, in the detection step, a Jarque-Bera (JB) normality test has been performed (Miot 2017), with a 95% confidence level.
After the test, with the identification of outlier points, in the adaptation process, the hypothesis of spatial correlation was verified, that is, if the outliers points were spatially grouped in a subregion. Therefore, a spatial data clustering analysis technique has been used, based on the standard deviation value. This technique considered that spatially correlated outliers can express a different local bias from that verified in the rest of the data in the sample set, rather than simply the presence of gross error. Thus, for the study region, there was the possibility of carrying out of linkage studies for one or more areas, called study subregions, with different local parameters.
In the case of linkage to the IHRS, there is no _ 0 IBVD SS W parameter to be estimated. To calculate the discrepancies were used only the check points. Furthermore, in the case of linkage to the future IHRS, there is no possibility of calculating the accuracy since there are no values to be taken as a reference.

Definition of study subregions and one-third parts for the case of local modeling
For the case of local modeling, after obtaining the discrepancies between the height anomalies ζ BGS P and GGM P ζ , calculated for the 102 GPS/RN connection stations in the study region, a histogram was generated, which is shown SP in Figure 2.
As can be seen, most of the discrepancies were located between 0 and 0.4 m, with a mean equal to 0.259 m and standard deviation equal to 0.343 m. However, occurrences of values above 0.7 m can also be verified, moving the histogram away from a normal distribution.
When analyzing this question, the JB normality test was applied to the dataset of the study region. The test result showed an index higher than the defined critical value for 95% confidence level.
The presence of outliers was verified using the Data Snooping method. After the process, the points considered outliers were those with the greatest discrepancy between the calculated height anomalies, located in the right part of the histogram. After applying the spatial data clustering analysis, it was verified that the outlier points were all clustered in the southeast area of the study region, indicating a spatial correlation. This corroborates the result seen in the right part of the histogram. Since the subregion of the outlier points is in a coastal and mountainous area, the difference in the systematic component may be related GGM/RTM inaccuracies. The outlier points subregion was named SP2 and the other SP1, without excluding any of the points in this classification.
After ranking the points in each subregion, the normality of the associated datasets was verified. The JB normality test was applied to the dataset of the SP1 study subregion and the Shapiro-Wilk normality test was applied to the dataset of the SP2 study subregion, as it is a sample of less than 50 observations. In both cases, the datasets were within the confidence level and can be considered normally distributed, although in SP2 the sample is small for conclusions and the test results may not be significant.
The histograms of discrepancies between the height anomalies ζ BGS P and / ζ GGM RTM P in the SP1 and SP2 study sub-regions are shown in Figure 2. Through these results, the idea of estimating different local parameters in each subregion has been adopted. Regarding the question of evaluating the estimated parameters, the procedure of one-thirds was applied with the stations available in each subregion, alternately Alpha, Bravo and Charlie, as mentioned in item 2.3. The criterion used was the search for three groups with points distributed evenly in space and in discrepancies. In Figure  3, the check points in each one-third of the subregion are presented.

Analysis of the linkage results
In the case of local modeling, using the reference points defined for each one-third in each subregion, it was possible to determine the 0 IBVD W local transformation parameters, according to Equation (8). The calculated values are shown in Table 1, as well as their respective standard deviations. For comparison purposes, the value of _ 0 IBVD SS W obtained by Sánchez and Sideris (2017) is also presented.  The estimated parameters values indicate satisfactory agreement between the defined one-thirds. While in subregion SP1 the discrepancy values are in the order of 10 -2 m 2 /s 2 or ~10 -3 m, in subregion SP2, in the worst case (Alpha/Charlie) the value is 0.14 m. This indicates that the local parameters can be estimated regardless of the location of the points, especially in the SP1 subregion, which has the greatest number of points. Regarding the difference for the value of _ 0 IBVD SS W , the mean for subregion SP1 was 9.409 m 2 /s 2 or ~1 m and for subregion SP2 it was 3.410 m 2 /s 2 or ~ 0.35 m, indicating differences that can be significant for some applications, mainly in the SP1 subregion. The smaller difference between the values of 0 IBVD W and _ 0 IBVD SS W in SP2 indicates greater adherence and better applicability of the national parameter in this subregion.
After estimating the transformation parameters, using the r reference points in each one-third of the study subregion, the determination of normal heights _ N trf IBVD Q H for the t check points have been carried out, according to Equation (9), for the local case. Then, the sets of differences In Figures S1 and S2 (supplementary files) the differences for all points are presented, by one-thirds, in the two study subregions.
For the linkage to IBVD, as can be seen, for both subregions, in all one-thirds, the differences  Table 1 with equal displacement values. Consequently, it is possible to verify that there was a better adjust of the data in the case of local modeling. In general, it is also possible to verify that the choice of one-thirds did not change the results. The differences in this topic did not show any spatial correlation with the topography.
In the case of linkage to the future IHRS, in subregion SP1 the difference values _ trf IHRS δ were all positive with a mean value around 0.55 m. On the other hand, in subregion SP2 with occurrence balanced of positive and negative values, bringing the mean of the differences to close to zero. In this case, while in subregion SP1 there is a certain systematic bias, in subregion SP2 there is a possible presence of random errors or random errors with systematic errors of low magnitude.
In order to statistically verify whether systematic bias are significant or not, at a confidence level of 95%, statistical tests with Student's t distribution (t ≤ 30) were performed in the two subregions, using the three onethirds, with the following hypotheses: where is the mean of a given dataset of differences for a given one-third part/subregion. Theoretical statistics for subregions SP1 and SP2 are ±2.05 and ±3.18, respectively. The values of the calculated t student statistics are shown in Table 2. δ cannot be considered free of significant systematic biases at a 95% confidence level only in the SP1 subregion.
As presented in item 2.3, to obtain the accuracy of the transformation in each one-third of each subregion, the RMSE of each dataset of differences _ δ trf IBVD and _ _ δ trf IBVD SS has been calculated. Additionally, the mean RMSE have been also calculated. The results are shown in Table 3. As mentioned, in the case of linkage to the future IHRS, it is not possible to calculate the accuracy, since there are no height values to be taken as reference. As can be seen, the agreement in the RMSE values between the one-thirds is predominantly at the level of mm, with a maximum difference of the order of cm.
In subregion SP1, the accuracies were practically ten times better when using the local modeling approach, regardless of the one-third part of reference points used. The mean difference was 0.90 m. The values obtained with the local parameter were approximately 0.10 m.
On the other hand, in subregion SP2, the differences in accuracy between the local and national modeling approach were smaller, in the order of 0.22 m, with the local modeling presenting more accurate values. In general, the values were less accurate in relation to the SP1 subregion. The more accurate values with the national modeling approach in subregion SP2 in relation to subregion SP1 is a direct consequence of the greater proximity between the values of 0 IBVD W and _ 0 IBVD SS W in this sub-region, as indicated above with the data in Table 1.

Conclusions and outlook
The use of GNSS positioning integrated with GGM data, refined by DTMs, to determine normal heights in Brazil, referred to a LVD such as IBVD or IHRS, configures an alternative to the traditional spirit leveling method. However, there are two issues to be considered, the lower accuracy and the difference in terms of vertical datum. While the first question is solved by limiting the use of such methodology to certain types of applications, in the second, the solution is to carry out the transformation between HRFs.
Considering the second question, the objective of this work was to evaluate the use of GNSS ellipsoidal heights, in conjunction with height anomalies provided by GGM, refined by RTM technique, to obtain normal heights in Brazil, referred to the IBVD and future IHRS. In the case of linkage to IBVD, local and national modelling approach have been analyzed using the parameter obtained by Sánchez and Sideris (2017). The study region was the state of São Paulo.
For the linkage, the methodology indicated in Ihde et al. (2017) was adapted. The adaptation was made to directly provide the normal height values already transformed. In the case of linkage to IBVD with local modeling, the presence of outliers in the dataset of differences between the BGS and GGM/RTM height anomalies was firstly investigated. After the spatial clustering analysis, the spatial correlation in the outliers was detected and two subregions were defined for modeling.
Thus, in the linkage to IBVD, it was found that the differences between the local and national parameters approach for the SP1 and SP2 subregions were ~1 m and ~0.35 m, which may be significant for some applications, especially in the SP1. This, together with the histogram analysis of the differences, corroborates the idea of using a parameter better adjusted to the local case.
In the parameters validation step, for the case of linkage to IBVD, in both subregions, regardless of the configuration of reference and check points used, more accurate results have been obtained with local modeling approach than with national modeling approach. Also in this context, a significant systematic bias (95% confidence level) in the differences of normal heights referred to IBVD has been verified using the national modeling approach. For the case of linkage to the future IHRS, only in subregion SP1 was there a significant systematic bias (95% confidence level) in differences from normal heights.
With this research, the perspective is broadened of determining normal heights of points of interest in a study region in Brazil, with a fieldwork of only a GNSS survey of the point. This survey integrated with a post-processing, using the transformation parameter, allows obtaining the normal height of a given point referred to IBVD and the future IHRS with a decimetric precision. Although the calculated accuracies are below those desired by the IAG/ GGOS, the local modeling approach is more accurate than a national modeling approach, with possibilities for adjustments in procedures and part of the improvement process that is sought at this level of research.
The methodology comprises a modeling with open data with wide access, provided by government agencies and institutions and widely accepted in the scientific community. Thus, the same procedure performed here can be extended to other states and regions of Brazil. Additionally, it is recommended analysis of: -other GGM/RTM solutions; -improvements in the calculation of the RTM effect, using tesseroids and harmonic correction; -studies of other linkage strategies, such as the one based on the solution of the GBVP (Rummel and Teunissen 1988) and; -more rigorous analysis of clusters to obtain smaller subregions, to improvements of accuracy.
Finally, it is worth noting that since the geodetic stations have vertical movement, because of temporal variations of physical heights, the parameters calculated in the research are valid for a given time and region (Ferreira et al. 2019)suggesting that South America's hydrological cycle is changing, which impacts its water availability and, consequently, the Earth's surface due to its elastic response to the surface mass loading/unloading. Therefore, we analyzed 3 to 15 years of vertical crustal displacements (VCDs. Due to the lack of an estimate of vertical velocity and acceleration, due to the difficulties of modeling it is causative factors, vertical coordinates may change over time at the centimetric level (Blewitt, Hammond, and Kreemer 2018).Thus, consideration is given to the need to update the parameters, within what is recommended by the methodology, whenever there is a change in epoch or region.

AUTHOR'S CONTRIBUTION
All authors contributed equally.