PHYSICAL CONNECTION BETWEEN BVRF SEGMENTS BASED ON LEVELING ASSOCIATED WITH GRAVIMETRY

Considering the efforts to establish Global Reference Systems linked to the geopotential space, new alternatives are sought to address the problems found in the classic national vertical networks. The Brazilian Vertical Reference Frame (BVRF) was materialized in two different segments with independent datums (Imbituba and Santana tide gauges) due to the terrain difficulties for conventional leveling. The 2018 BVRF realization, in the geopotential space, still remains without interoperability between its segments. We analyze alternatives for physical connection based on the new precepts of the International Association of Geodesy (IAG) involving the geopotential space. Some proposed solutions for physical connection based on GPS leveling associated with gravimetry are presented. These solutions were developed with the aim of evidencing the discrepancy between the two BVRF segments, now carried out in terms of geopotential numbers and normal heights. The results indicate differences ranging from about 45 cm to 140 cm between the two segments depending on the strategy employed. Comparisons with previous determinations based on indirect strategies and involving previous BVRF realizations are made.


Introduction
Earth System change monitoring activities demand accuracy and reliability in realized spatial position referencing.The requirements for quantifying the dynamical process within the planet are strongly dependent on modeling the mass flux inside the Earth's System based on changes in its geometry, rotation and gravity field (Kutterer, Neilan, and Bianco 2012).In this sense, the quality of Geodetic Reference Systems (GRS), involving positions, gravity, and geodynamical aspects, must support accuracy requirements tending to 1 ppb in the forthcoming years (Plag et al. 2009).These aspects refer to the main present challenges in the global Geodesy.
Considering the need for an effective Earth Observing System for monitoring global changes, in 2003 the International Association of Geodesy (IAG) established the Global Geodetic Observing System (GGOS) as a central element for integrating all global geodetic infrastructure composed by systems involved at diverse levels of observation, conventions, databases, and models as well as new geodetic methods and services.GGOS Theme 1 -Global Unified Height System -was established in 2010 as a fundamental element for integrating geometry and gravity field aspects.These aspects were reinforced by United Nations General Assembly Resolution A/RES/69/266 in February 2015.This resolution established "A Global Geodetic Reference Frame for Sustainable Development", within the "United Nations Global Geospatial Information Management (UN-GGIM)" by recognizing the coordinated approach of IAG with regard to GGOS.IAG Resolution 1/2015 defined the International Height Reference System (IHRS) with the meaning of a global equipotential surface as a common global vertical reference with geopotential value W0 = 62 636 853.4 m 2 s-2 and the geopotential number CP = W0 -WP, as the primary vertical coordinate (IAG 2015a).IAG Resolution 2/2015 related to the "Global Absolute Gravity Reference System" (GAGRS) intended to replace the "International Gravity Standardization Net 1971 -IGSN71" (IAG 2015b).As a consequence of this, the Global Geodetic Reference System (GGRS) and its realization (GGRF) were defined, where each point P is described by its geometric coordinate vector  ⃗  (X, Y, Z) referring to the space-oriented Earth involved in the realization of the International Terrestrial Reference System (ITRS) , and also by the geopotential WP, the physical height H(CP), and the gravity vector  ⃗ involved in IHRS.GGRS thus includes the geometry, gravity field, and the Earth's orientation in space (IAG 2016).
Several efforts are currently underway for realizing the global IHRS (IHRF) frame.There are still several limitations for establishing IHRF, such as heterogeneities in standards for realizing National Vertical Reference Systems (NVRS) and the lack of geodetic information with adequate metadata for supporting the solution of the Geodetic Boundary Value Problem (GBVP) independent of local reference frames and with qualification of data and elimination of outliers (Santacruz and De Freitas. 2015).These backgrounds are fundamental for the determination of geopotential numbers with a minimum accuracy of 1x10 -2 m 2 s -2 coherent with IHRS and GGOS purposes (Drewes et al. 2016).
The coordination of activities related to the IHRS realization is done by the "IAG/GGOS WG 0.1.2on the Strategy for the Realization of the IHRS" in which the SIRGAS Working Group III (Vertical Datum) takes part.In this context the following SIRGAS WGIII protocols are distributed to member countries: • Strategies for the establishment of vertical networks by physical heights [HP = f(CP)]; • Link of national vertical networks to the SIRGAS GNSS Continuous Stations; • Integration of national vertical networks of members countries in the geopotential space; • Approaches for referring the SIRGAS Vertical Frame to the W0 value of IHRS; • Association to a specific epoch by considering the realization epoch and temporal variations in coordinates; • Planning of activities for establishing a GGRF/IHRF station profile in the SIRGAS region.
Argentina, Brazil, and Uruguay achieved partial relevant results related to the SIRGAS WG III protocols by realizing their NVRS in the geopotential space although still linked to local Vertical Datums.However, several problems remain when IHRS standards are considered.
The Brazilian Geographic and Statistics Institute (IBGE) is working to surpass some structural problems for modernizing the Brazilian Vertical Reference Frame (BVRF).One of these problems is related to the realization of the network in two different segments.Most benchmarks are linked to the southern Brazilian Vertical Datum of the Imbituba tide gauge (BVD -I) and about five percent of stations are linked to the northern Vertical Datum of the Santana tide gauge (BVD -S).It is almost impossible to link the two segments by conventional spirit leveling because of the difficulties in surpassing the Amazon River mouth and large inaccessible portions of the Amazon Rainforest involving distances over three hundred kilometers of inaccessible wet regions.
The studies conducted described in this manuscript aim to provide solutions to several problems based on data coming from different modern geodetic techniques.A fundamental step is to restrict data to the same Geodetic Reference System and permanent tide System.These conditions are essential in order to observe the interoperability condition.The reference frames mentioned above are followed in this work and it is in accordance with IAG recommendations for realizing IHRS and for modernizing NVRS (Sánchez et al. 2017;Sánchez and Sideris 2017).IBGE's strategies and current activities in cooperation with the Federal University of Paraná (UFPR) for solving the connections between the two BVRF segments in the geopotential space and the results obtained are presented in this paper.

Amazonian Connection Area and description of involved Database
During the course of the historical cooperation between IBGE and UFPR related to BVRF modernization (De Freitas, Ferreira, and Luz 2018), an area 10° x 10° was delimited in the Amazonian region involving the states of Amapá and Pará with the aim of establishing a physical connection between the two BVRF segments (Fig. 1).New strategies for leveling associated with gravimetry were considered.In order to enable a better understanding, the data set was divided into two segments as per Figure 1 and Table 1.Figures 2 and 3 present the reference surfaces of the initial points of the leveled segments by considering the geoid height as well as height anomaly computed from the Global Geopotential Model (MGG) XGM2016 (Pail et al. 2016) in its maximum degree and order (d/o).We justify the use of this model because of its performance in the region, even considering its medium resolution (Nicacio, Dalazoana and De Freitas 2018).
IBGE's Geodesy Coordination carried out the surveys related to the two connection segments mentioned.In segment 1 (A -B in Fig. 1), which relates to BVD -I conventional spirit leveling is impossible.The survey of level differences in this segment was therefore carried out based on GPS and gravimetry along a stretch of about 269 km, considering the maximum distance of 15 km between two neighboring GPS stations because of the uneven characteristics of region involved and requirements for gravity frequency measurements at 10 km to 15 km intervals in undulating topography regions according to Ramsayer (1963) andDe Freitas andBlitzkow (1999).Spirit leveling associated with gravimetry and sparse GPS positioning was realized along about 172 km of segment 2 which is linked to the BVD -S.Besides the densification of geodetic observations in the region, this new set of information was established with the aim of improving some studies on connecting the Santana and Imbituba BVD conducted by LARAS (Laboratory of Geodetic Reference Systems and Satellite Altimetry of UFPR).These studies are described in Montecino and De Freitas (2014), Moreira andDe Freitas (2016) andDe Freitas et al. (2016).These previous works were based on the least square estimation of differences between modeled height anomalies, and Bulletin of Geodetic Sciences, 25(1): e2019006, 2019 GPS/Lev.differences.Spectral decomposition, Global Geopotential Models (GGMs) enhancement and RTM were applied without observing true physical level differences.The results are synthesized in Table (3) which already includes corrections to previous different BVRF realizations with the aim of making them compatible with the 2018 BVRF realization.

Database Analysis
The data obtained for heights and gravity from the two leveled segments were preanalyzed as to their consistency.Several actions were necessary for correcting above all changes of nomination, duplicity, and several inconsistencies related to different numerical forms for expressing the same kind of measurements (e.g.owing to lack of significative numbers).Analyses were also performed based on GGMs with the aim of detecting possible outliers because of the complex characteristics of the area involved and new techniques employed.It was also necessary to qualify previously available data in several different databases relating to different reference frames and epoch.All the data obtained were reduced to the mean permanent tide system (see e.g.Tenzer et al. 2011).
With regard to segment 1, the GPS survey results of 23 points, the associated gravity observations and the normal height of the initial point were analyzed.It must be emphasized that initial point A is not part of the official BVRF because of the large destruction of BMs in the region.The starting point used (BM104, Table 2) is a second order point derived from the 2011 BVRF.
In segment 2 it was necessary to recover raw leveling differences referring to BVD -S and in some cases to interpolated gravity information for computing geopotential number differences with a severe discrimination of wrong information.

Methodology and Results
In accordance with the conventions for IHRS established in IAG Resolution 1/2015 (IAG 2015a), the global reference zero level is that associated with the equipotential surface with the conventional value W0 = 62 636 853,4 m 2 s -2 (Sánchez et al. 2016).In this system, the primary vertical coordinates of a point P with geopotential WP is the geopotential number, already referred in the Introduction, given by (more details in Hofmann-Wellenhof and Moritz 2005, pp.55): Considering a particular NVRS, the geopotential number refers to the particular equipotential surface W0i and is given by: Where the geopotential number can be computed with sufficient approximation by the sum of the products of the mean gravity   value in a leveling section and the observed level differences ∆ if the leveling section is not too large (see practical conditions e.g. in De Freitas and Blitzkow 1999).Because the characteristics of the two segments of leveling used in the connection are quite different (GPS + gravimetry in segment 1 and spirit leveling + gravimetry + sparse GPS in segment 2) some strategies were developed for testing the approach used.In segment 1 the ellipsoidal height differences and gravity measurements were considered as the main subject while the raw spirit leveling differences associated with gravity values in each section were the main subject in segment 2. Different solutions were employed for each segment as set out in the following section.

Strategies for computing normal heights and geopotential numbers.
The leveling line of each segment started from a point where normal height was determined.At starting point A (BM104, Table 2) of leveling segment 1, a discrepancy of about -216 mm was considered between the normal heights in the 2018 BVRF adjustment and the normal-orthometric height in the 2011 adjustment.The reason is that in the region around point A the interval of the discrepancy detected by the 2018 BVRF adjustment, which was concluded by IBGE on July 30 th , 2018, ranges from -215 mm to -217 mm.In leveling segment 2 we obtained the normal height at starting point C by considering (Hofmann-Wellenhof and Moritz 2005, pp.326): Where   and   are respectively the normal-orthometric and the normal heights in the local system, (  −  ) is the difference between the geoid height and height anomaly (see Fig. 2 and Fig. 3),   is the Bouguer anomaly (computed using GGM XGM2016 in this case) and   is the mean normal gravity value between the telluroid and reference ellipsoid (see Fig. 2 and Fig. 3).In either case (local or global system) normal height is related to the referred mean normal gravity value and the geopotential number according to this equation: If the normal height of the initial point    is known, then it is possible to obtain the normal heights (and geopotential numbers) along the leveling line if the normal height differences   between neighboring points are known.Provisional values are obtained as follows: Where  (1−2)  is the observed level difference.It is possible to obtain corrected values for normal height if corrected normal height differences    are known using an iterative approach by applying the following equation, for example: Where the iteration is realized by recomputing the mean normal gravity value and the corrected normal height coming from the corrected normal height difference until there is stabilization in an iterative process, remembering that the mean normal gravity   value between the ellipsoid and the telluroid can be computed by using a simple or a rigorous method.
By using a simple approximation, which is useful for many cases, such as in uneven regions, the simple method is: Or by using a rigorous method such as the one used in this study (Gemael 1999): Considering that , ,  are constant parameters of the level ellipsoid associated with the GRS under consideration (see numerical values for GRS80 in Hofmann-Wellenhof and Moritz 2005, pp.86).
The determination of    can be done by applying the Clairaut theorem or the Somigliana formula.Both are functions of geodetic latitude and are based on geometric and physical parameters of the level ellipsoid (see, e.g.Hofmann-Wellenhof and Moritz 2005, pp.71-72 and 86).

The corrected normal height 𝐻 𝑐
for each point in the segment is obtained after the convergence of the iterative process (usually two or three loops are sufficient).
The determination of normal heights and geopotential numbers can be based on rigorous Normal Correction (NC) to be applied on the observed level difference between two adjacent points in each leveled segment in substitution of Eq. 6 in the following form (Hofmann-Wellenhof and Moritz 2005, pp. 168): Its solution considers  0 =  45°= 9.806199203 ms -2 to be the normal gravity value for latitude 45°,   obtained by the previously mentioned iterative process and    of each point.Then    is obtained by considering e.g. between two points (1 and 2) the addition of  (1−2) to the observed level difference value ∆ (1−2) .This approach results in:

Proposed Solutions
Solution 1: This solution took into consideration gravity correction on the observed level differences in both leveling segments according to Eq. ( 6).Even considering that point A is not directly part of BVRF, it was considered that around point A the mean systematic difference previously referred to is observed: It must be emphasized that the condition expressed in Eq. ( 12) was applied in all subsequent solutions.Solution 2: Heights and geopotential numbers by considering the rigorous NC (Eq.10) and initial conditions expressed in Solution 1.
Observation: It must be noted that GPS leveling involves open questions when applied to determine level differences with physical significance.The physical significance is obtained by considering the correction coming from the deflection of the vertical in a direct form on the GPS + gravimetry level differences or in an indirect form by using height anomaly differences derived from high-resolution GGMs, in a relative form, as explored below.This approach can be useful for improving the connection computation based on GPS leveling differences for obtaining normal height differences or geopotential number differences even in the local BVRF.This must be considered because the critical distances involved in the GPS + gravimetry leveling when comparing the usual sections of spirit leveling and the possibility of existing strong local physical effects coming from local crustal density anomalies.Solution 3: Determination of normal heights in leveling segment 2 only based on GPS leveling associated with gravimetry and rigorous NC.In leveling segment 1 we take the same value for Solution 2. The aim of this approach was to test the GPS leveling procedure in heterogeneous conditions, involving level differences over heterogeneous distances ranging from 0.12km to 33km.This aspect must be explored in a future more rigorous investigation.Sciences, 25(1): e2019006, 2019 Solution 4: As per the observation made in Solution 2 and because of the lack of external control over previous solutions based on GPS + gravimetry leveling, we introduced height resolution enhanced GGM as part of the solution.The basic selected functional was the differences of height anomaly in each GPS leveled section by high resolution enhanced GGM GECO 2015 with harmonic development up until degree and order 2190 (Gilardoni, Reguzzoni and Sampietro 2016).It must be stated that according to these authors this model involves an optimal least-squares combination between the EGM2008 (Pavlis et al. 2008) and GOCE spherical harmonic coefficients.Comparisons with other recent combined models, such as EIGEN-6C4 (Förste et al. 2014), and a local geoid/quasigeoid based on new gravity datasets show that the proposed combination, weighting the different input contributions not only on a global basis but also according to some local error information, can perform even better than other more sophisticated combinations in areas where the input global error description is not reliable enough.Solution 4 is more performant than the solution given by De Freitas et al. ( 2016) because the referred local error was not used in the enhancement of combined model in that solution as described in Ferreira, De Freitas and Heck (2015).

Bulletin of Geodetic
It is possible to observe in Fig. 4 (a, b, c and d) that in the region involved there are severe changes in the height anomaly.Note that all the considered models repeat similar behavior and patterns related to a strong variation of height anomaly in the region of leveling segment 1.It must be emphasized that the model satellite only based GOCO 05s (Mayer-Gürr 2015) is independent of some local reference frames and observations.Thus, the strong variation mentioned is real and reflects even on the satellite orbit of GOCE.This fact is not associated with observation errors.This anomaly evidently affects GPS + gravimetry leveling.We, therefore, introduced such variations as an indirect way for considering the discrepancies between the level ellipsoid and the geoid/quasigeoid.Solution 5: Same as Solution 4 but considering GGM XGM2016, with harmonic expansion up until degree and order 719.This model was adopted following suggestions by Nicacio, Dalazoana and De Freitas (2018) relating to the relative good behavior of this model in the region and the new approach used for its computation by considering its regional error behavior and data availability.
Each solution presented had the purpose of testing techniques for determining physical heights and geopotential numbers based on independent GBVP solution approaches.Table 4 presents the computed discrepancies for    and   for solutions coming from each leveling segment respectively linked to BVD -I (adopting different solutions at point A) and BVD -S (adopting point C as the origin).The five solutions explained are linked to level differences based on GPS ellipsoidal heights associated with gravimetry considering the BVRF 2018 realization.GPS leveling associated with gravimetry is the only current possibility in solving the physical connection in segment 1.We can note that among the solutions developed, Solution 4 has more agreement with previous solutions summarized in section 2 (see Table 3).We emphasize that this solution is based on GPS observed level differences in segment 1 corrected by the difference of computed height anomalies from local observations and those modeled by GGM GECO.
Figure 5 shows the normal height discrepancy (between segment 1 related to BVD -I and segment 2, associated with BVD -S) at point B of all solutions performed in this research, which varies between 0.4512 m and 1.4195 m.With the aim of checking the consistency of solutions based on GPS leveling and gravimetry, the recomputed height anomalies were considered in view of the corrections to the GPS level differences.Fig. 6 shows the height anomalies for Solution 2 (almost constant on the Fig. 6 scale) and Solutions 4 and 5 by considering the physical variations coming from regional existing crustal heterogeneities as shown in Fig. 4 by different GGMs.

Remarks, Conclusions and Outlook
Five solutions for establishing the first leveling connection, with physical significance, point to point, between the two parts of BVRF were tested.Two leveling lines were made in the connection region (point A to point B, and point C to point B, as per Fig. 1).Segment A to B is about 270 km long and GPS associated with gravimetry was performed because of the land characteristics in the region.The other segment is about 170 km long and spirit leveling line associated with gravity was used.The solutions were compared with previous connection attempts based on indirect regional approaches in the geopotential space.
The previous indirect solutions in regional approaches have different characteristics to the new proposed point to point leveling solutions.
The previous solutions involved different realizations of BVRF, distributed GPS/Lev points with irregular distribution and little information on the BVD Santana side.An allowance was made for some error estimates based on a least square adjustment without, however, considering the possible error sources in the spirit leveling, GPS observations and commission and omission errors in the GGMs and DEMs used.
The five solutions proposed here depend strongly on point-based observations, mainly at points A, B, and C. It must be emphasized that GPS leveling aiming level differences with physical meaning is still an open subject of investigation.It was evidenced that solutions 1, 2 and 3, based on conventional gravity and normal corrections to level differences by GPS did not appear to address strong lateral heterogeneities in the Earth´s gravity field along the leveling lines evidenced in Fig. 4. Solutions 1 and 2 must be tested for short baseline among leveling stations.This is a fundamental step to be accomplished for understanding the viability of GNSS leveling associated only with local gravity corrections.
Solutions 4 and 5 integrate in a relative way the regional characteristics of crust based on two GGMs.Solution 4 is based on the high-resolution GECO enhanced model and Solution 5 is based on medium resolution GGM XGM2016, a previous test-model of the future EGM2020.Both enhanced models involve an optimum combination of harmonic coefficients of existing GGMs and the weighting of local errors for reducing local commission errors.This represents the cutting edge of development for generating a new generation of GGMs.Given the characteristics of the GECO model and comparisons with other recent models referred to, we consider Solution 4 to be the most controlled one.
It must be emphasized that Solution 4 shows good agreement with the previous solutions presented in Table 3.We emphasize the previous solution developed by De Freitas et al. (2016) in the determination of the offset between the two parts of BVRF, respectively 1.237 m (Solution 4, Table 4) and 1.200 m (Table 3).However, there is no clear possibility for quality control of these because they are based on open leveling lines.
A real possibility for establishing control on the connecting leveling lines is to determine geopotential numbers (or height anomalies) related to IHRS/IHRF and to check their discrepancies related to the height anomalies computed in the five solutions proposed.This approach assumes the resolution of the GBVP based on Molodensky gravity anomalies (also known as free-air surface gravity anomalies) or based on surface gravity disturbances.A minimum control could be established by resolving GBVP at least at points A, B, and C.

Figure 1 :
Figure 1: Connection Area, existing Data and new Surveys between the segments linked to BVD -I and BVD -S (referred to the BVRF).

Figure 2 :
Figure 2: Point A position and related physical and geometric reference surfaces.

Figure 3 :
Figure 3: Point C position and related physical and geometric reference surfaces.

*
In these solutions, we can consider only error estimation of about 0.11m at the initial point A related to the 2018 BVRF adjustment.Error propagation along segment 1 between points A and B as well as along segment 2 between points C and B cannot be controlled (open leveling lines).Each leveled point in segment 1 has an implicit error coming from the GPS ellipsoidal height determination of about 0.05 m.There is no possibility of control in segment 2 between points B and C. ** In these solutions there is the possibility of R.M.S estimation only based on the GGM used local determination of commission errors exists.Unfortunately, there are only estimates of such errors based on other GGMs.

Figure 5 :
Figure 5: Normal height discrepancies at point B.

Table 1 :
Details of the segments

Table 3 :
Different computed offsets between BVD-I and BVD-S The offset calculated byMontecino and De Freitas (2014)considering a solution involving GOCE+EGM2008+RTM is 1.420 m (0.56m).This offset shows that the Imbituba datum is located below the Santana datum.This offset value took into consideration the discrepancy of 0.15 m at point A (Fig.1) between AAGP (Ajustamento Altimétrico Global Preliminar -Preliminary Global Altimetric Adjustment) and the 2011 adjustment (IBGE 2011 pp.52), and also the difference of -0.216 m between the 2011 and 2018 adjustments(IBGE 2018, pp.35).There is no control over the error estimation of heights related to AAGP in the region.We then propagated only the error estimation of about 0.14 m around point A from the 2011 adjustment(IBGE 2011, pp.50)and that from the 2018 adjustment (0.11 m) according to IBGE (2018 pp.34) even considering that there is a discrepancy involved between AAGP and the 2011 adjustment.
* ** The offset calculated by Moreira and De Freitas (2016) is 1.30 m (0.11 m) and that calculated by De Freitas et al. (2016) is 1.416 m (0.12 m).These values took into consideration the difference of -0.216 m at point A between the 2011 and 2018 BVRF realizations.The error estimation of the adjusted heights in the point A region relating to 2011 and 2018 adjustments was propagated for both solutions.

Table 4 :
Normal Heights and Geopotential Numbers of point B using different solutions