PL-geoid2021: A quasigeoid model for Poland developed using geophysical gravity data inversion technique

This paper presents the results of research and analyses related to the development of a new quasigeoid model fitted to GNSS/levelling data for the area of Poland (PL-geoid2021). The model was determined employing two procedures based on the Geophysical Gravity data Inversion technique (GGI method): procedure A consisted of the determination of the gravimetric quasigeoid model in the first step and its subsequent fitting to GNSS/levelling data in the second step, and procedure B consisted of a one-step determination of the model fitted to GNSS/levelling data. Both models were developed using the global geopotential model SGG-UGM-2 and gravity data covering the area of Poland, and slightly extend beyond Poland's southern and northern borders. The average model was adopted as the final model. It was demonstrated that the accuracy of the gravimetric quasigeoid model had a very low dependence on the reference topographic mass density model used. On the basis of this model, the GNSS/levelling datasets were also assessed and outliers were identified. The estimated accuracy of the gravimetric model, determined based on four GNSS/levelling datasets, was in the range of ± 1.2 to ± 1.7 cm, in terms of the standard deviation of the differences between the measured and model-determined height anomalies. Due to partial lack of gravity data just beyond the Polish border, the edge effect was also analysed. The accuracy of the final quasigeoid model (estimated in the same way as the gravimetric model) ranges from ± 1.0 to ± 1.2 cm. It should be noted, however, that this assessment is not fully independent because three of the four sets of GNSS/levelling points used for it, were also used to build the final model.


Introduction
Currently, in the determination of heights by geodetic measurements, spirit levelling or GNSS measurements (with the use of a geoid or quasigeoid model) are most often used. In order to ensure the compliance of the heights determined by both methods, the geoid or quasigeoid model used should be compatible with the height datum defined by the levelling network and the reference frame used in GNSS measurements. Therefore, the simple dependencies should be met (e.g., Heiskanen and Moritz 1967): N = h − H o if orthometric heights ( H o ) are used and = h − H N if normal heights ( H N ) are used. In the above equations, N and are the geoid height and the height anomaly respectively, determined from the appropriate model, and h is the ellipsoidal height determined by GNSS measurements.
The geoid and quasigeoid models based only on gravity data (gravimetric models) represent a reference surface independent of the levelling network. Unfortunately, usually these kinds of models are biased or biased and tilted in relation to the GNSS/levelling data, that are caused by systematic errors of levelling and GNSS data as well as the long-wavelength errors of geoid or quasigeoid model. This also necessitates the use of GNSS/levelling data (in addition to gravity data) in the modelling process. Thus, in general, the construction of a geoid or quasigeoid model that can be used for GNSS levelling can follow one-or twosteps procedures: • One-step procedure assumes use of a modelling solution that combines both gravity and GNSS/levelling data and allows the determination of a model fitted to GNSS/levelling data. • In a two-steps procedure the gravimetric geoid or quasigeoid model is determined in the first step and then the model is fitted to the GNSS/levelling data in the second step.
This paper presents the results of research and analyses related to the use of both of the above-mentioned procedures to build a new quasigeoid model fitted to GNSS/levelling data for the area of Poland. The calculations were performed using a method based on geophysical gravity data inversion (GGI).
It should be emphasised that regional geoid modelling for the area of Poland has a long history. The first model was developed based on astro-gravimetric deflections of the vertical (Bokun 1961), and then over twenty different geoid and quasigeoid models have been developed using various techniques and methods. These models were determined using gravity data, deflections of the vertical, GNSS/levelling data and various Global Geopotentials Models (GGMs). A review of the models developed up to the middle of the second decade of the twenty-first century can be found in Kryński (2014). This overview can be supplemented by two other quasigeoid models based on EGM2008. The first model, by Kuczynska-Siehien et al. (2016), was developed using the KTH method. The model accuracy was estimated at the level of ±2.2cm . The second model, named GDQM-PL19 (Kryński et al. 2023), was developed applying the remove-compute-restore technique with the LSC method. The accuracy of this model was estimated at the level of ±1.5cm.
It should be noted that among these models, the GUGiK'2001 quasigeoid model (Pażus et al. 2002) was the first model recommended for use in geodetic practice by the Polish Head Office of Geodesy and Cartography (GUGiK). The model was developed by fitting the gravimetric quasigeoid model Quasi97b (Łyszkowicz 1998, 2012) to GNSS/levelling points using a ninth-order polynomial, resulting in an estimated accuracy of ±2cm (Pażus et al. 2002). The second quasigeoid model to be recommended was the PL-geoid-2011 model (Kadaj 2013). The model was developed through the three-dimensional (seven-parameter) conformal transformation of the EGM2008 quasigeoid model to the GNSS/levelling points. The transformation procedure took into account the local Hausbrandt post-transformation corrections (Kadaj 2013). The accuracy of the PL-geoid-2011 model was estimated to be ±1.5cm , and it was related to the current height data and valid reference frame.
The third model recommended by the GUGiK is the quasigeoid model PL-geoid2021, which is described and analysed in this study. A partial description of the analyses concerning the model developed can also be found in the technical report titled "Technical description of the existing quasigeoid model PL-geoid2021 in the PL-EVRF2007-NH system" (in Polish), available at: http:// www. gugik. gov. pl/ bip/ prawo/ modele-danych. In addition to the analyses contained in the report, this study also includes the edge effect analysis and the influence of the density of topographic masses on the accuracy of the determined, gravimetric quasigeoid model.
Concluding this section, let's pay attention to the indicated by Szelachowska and Kryński (2014) difficulties with the evaluation of geoid and quasigeoid models. In Kryński (2014), it is emphasised that the accuracy of gravimetric models significantly exceeds the accuracy of GNSS/levelling determinations. For example, the accuracy of the gravimetric quasigeoid model GDQM-PL13 (Szelachowska and Kryński 2014)) is estimated to be better than ±1cm Kryński (2014), while a comparison of this model with GNSS/levelling data provided an estimated accuracy of ±1.4 − ±2.2cm depending on the adopted GNSS/ levelling dataset. This is also a problem to consider when fitting the gravimetric model to GNSS/levelling data.

The GGI approach
The quasigeoid modelling approach based on the geophysical inversion of gravity data used in this study has been described in various papers (e.g., Trojanowicz 2012a). In this method, calculations can be carried out according to various procedures. Therefore, it is essential to provide a more detailed description of the applied algorithm.
In general, the GGI solution is based on the construction of a local disturbing potential model (T) determined by three types of masses: • Topographic masses covering a limited area of the study, denoted as volume Ω; • Disturbing masses lying between the Moho discontinuity surface and the geoid surface in the study area, denoted as volume ; • Disturbing masses located outside of the volumes Ω and .
Hence, the disturbing potential at point P located on the terrain surface T P can be written as a sum of three components: where T Ω is the component produced by the topographic masses contained in the volume Ω , T is the component produced by the disturbing masses in the volume and T E (the external potential) represents the influence of the topographic and disturbing masses lying outside the volumes Ω and , respectively.
The T Ω and T components are determined using Newton's integral formulas: (1) where G is Newton's gravitational constant, and are functions of the density distributions in volumes Ω and , respectively, dV Ω and dV are elements of volumes Ω and , respectively, and l is the distance between the attracting mass and the point P.
It is assumed that the T E component describes trend-like changes and is represented by low-degree harmonic polynomials: where X P , Y P are the local coordinates of the point P and H P is the normal height.
The model of the disturbing potential defined by Eq. (1) makes it possible to formulate the following task: find the density distribution functions and and the polynomial coefficients a 1 , a 2 , a 3 , a 4 , a 5 to satisfy Eq. (1) for the measured data.
The formulated problem is solved using the linear inversion of gravity data (Blakely 1995) through the discretization of the functions and . The volumes Ω and are divided into blocks of a constant density. The volume Ω is defined by the regular grid of the digital elevation model (DEM). In order to limit the number of determined unknowns, DEM blocks are grouped into zones of constant, determined densities. The volume is defined as a plate with a thickness approximately equal to the depth of the compensation level and consists of one or more layers of constant density blocks.
Consequently, Eqs. (2) and (3) can be written as follows (Trojanowicz et al. 2020): where n is the number of constant density zones of the DEM; m k is the number of prisms of the DEM in zone k; k is the constant density of zone k; s is the number of prisms defining the volume ; j is the density of the rectangular prism j; and K i and K j represent the coefficients of the solutions of Newton's integrals for the DEM block i and prism j, respectively. So far the calculations have been performed in the local Cartesian coordinate system. The Z axis of the system is directed toward the geodetic zenith at the system origin and the X and Y axes are in the plane of the horizon and are directed to the north and east, respectively. The origin of the coordinate system is located near the centre of the study area. The K i and K j solutions in such a system can be found in, e.g. Nagy (1966) or Nagy et al. (2001).
In order to solve the problem of the ambiguity of the gravity inversion, a deep weighting function proposed by Li and Oldenburg (1998) was used, with the constant coefficients defined in Trojanowicz (2012b).
The unknown parameters of the model (polynomial coefficients a 1 , a 2 , a 3 , a 4 , a 5 and discrete density functions and ) are determined using the least squares method based j GK j on the following input data: GNSS/levelling height anomalies (converted into disturbing potential values) and gravity data, i.e. gravity anomalies or gravity disturbances at points on the terrain surface.
The calculations are performed by assuming certain reference density models for the volumes Ω and , denoted as 0 and 0 , respectively. For the 0 model, a constant value is assumed, while the discrete model 0 is defined by the following equation: where H i is the mean height of the constant density zone i of the volume Ω , which is located exactly above the constant density block j of the volume with a height h j .
Under such assumptions, the vector of unknowns, denoted as , contains the coefficients of the polynomial defining the T E potential and corrections to the reference densities 0 and 0 j . The observation equations for the disturbing potential ( T ) and the gravity disturbances ( g ) are: where v T and v g are adjustment errors, and z are known vectors of the parameters determined by Eqs. (4)-(6) and the approximate values T 0 and g 0 are derived from the reference density models 0 and 0 .
Previously conducted analyses (e.g., Trojanowicz 2015Trojanowicz , 2019 showed that the use of GGMs in the calculations significantly increases the accuracy of the quasigeoid models determined by the GGI method. In this version, the calculations are performed using the remove-compute-restore ( RCR ) procedure. When GNSS/levelling and gravity data are used, the calculations are carried out as follows: 1. The global component is removed and the residual values of the data T r and g r are determined: where g is the gravity value, W is the gravity potential and T GM and W GM are calculated from the GGM. 2. Based on the residual values, the residual disturbance potential model is determined using the GGI solution described above: where the components T E , T Ω and T are the residual components of T E , T Ω , and T (after the parts contained in the GGM are removed). 3. For new points, the residual disturbance potential from the GGI model and the value of the gravity potential from the GGM are determined. They are used to determine the gravity potential values from the GGI method W GGI , where U is the normal potential at the determined point and Q is the normal gravity acceleration on the telluroid.
In this case, observation Eqs. (8) and (9) are rewritten for the residual values: The quasigeoid model developed in this way is fitted to the GNSS/levelling height anomalies used and is not a gravimetric quasigeoid model. In Trojanowicz (2021), a solution of the GGI method has been proposed that also enables the determination of the gravimetric quasigeoid. Such a model, as mentioned in the introduction, represents a reference surface for height determination, independent of the levelling network. In the case of the determination of the gravimetric model using the GGI method, the GNSS/ levelling height anomalies in the computation procedure are replaced with height anomalies determined from the GGM. Since the disturbing potential from the GGM can be determined at any point, a regular grid of points extending beyond the study area is assumed for calculations. The points of this grid are assumed to be points with known height anomalies from the global model, so, as a consequence, we obtain a quasigeoid model fitted to the GGM adopted for calculations. Because in this case we assume that W = W GM ( T r = 0 ), observation Eq. (15) for the value of the disturbing potential is written as: Equation (16) Ihde et al. (2006)).
Finally, we can write If a quasigeoid model fitted to GNSS/levelling data is developed, the gravimetric model should be fitted to these data using an appropriate function.

Data used
The analyses used the following data provided by the Head Office of Geodesy and Cartography of Poland ( Fig. 1a): 1. 68,663 gravity points and free-air gravity anomalies: • 55,497 points from the collection of the Polish Geological Institute-National Research Institute (PGI-NRI) (light blue colour), • 2083 free-air gravity anomalies from the area of the Czech Republic (blue colour), • 8779 free-air gravity anomalies from the area of Slovakia (dark blue colour), • 516 points from the collection of the Warsaw University of Technology (red colour), • 1788 points from the collection of the Gdańsk University of Technology (green colour).
All gravity values were defined in the IGSN71 system. Among the aforementioned collections, detailed analyses are available only for the largest one, provided by PGI-NRI. The gravity values of this collection were determined in the years 1957-1989 and come from a gravity database of more than 800,000 gravity points (Królikowski 2006). The accuracy of gravity of this data is estimated at the level of ±0.075mGal (Kryński 2007). Horizontal coordinates of the points were determined based on topographical maps with accuracy at ±50m and their heights were determined by spirit levelling with accuracy estimated at ±4cm (Królikowski 2006). The accuracy of the data provided by the Warsaw University of Technology is estimated at ±0.015mGal (based on the list of points).
2. GNSS/levelling data with known ellipsoidal and normal heights, in the form of point sets that are part of the three Polish geodetic networks: ASG-EUPOS, POLREF and EUVN.
• The Polish Active Geodetic Network (ASG-EUPOS) constitutes the fundamental horizontal geodetic network for the area of Poland. ASG-EUPOS network points form the sets composed of station and two eccentric points located nearby. For calculations and analyses only the stations and one of two eccentric points for each station with known normal heights, were selected. The normal heights of the eccentric points were determined by precise, spirit levelling based on a levelling network. In relation to the eccentric points, the heights of stations were determined by satellite levelling. Due to the different methods of determining the normal heights of the stations and eccentric points, the points were separated into two sets, with slightly different accuracies of the determination of the GNSS/levelling height anomalies expected in these sets. These sets are hereinafter referred to as ASG_Stat (for stations) and ASG_Ecc (for eccentric points). • The POLish REference Frame (POLREF) network was established in the 1990s as a densification of the EUREF-POL network, which is an extension of the European reference frame ETRF89 to Polish territory (Bosy 2014, Zieliński 1993. The set of used points will be referred to as POLREF. • Points that are the Polish part of the European Vertical Reference Network (EUVN), hereinafter referred to as EUVN points.
The ellipsoidal heights of the ASG-EUPOS and POLREF network points were the official coordinates of the Polish basic horizontal geodetic network, and they were determined in ETRF2000 for epoch 2011.0. These heights are referred as PL-ETRF2000-GRS80h. The ellipsoidal heights of the EUVN network points were taken from report titled "Integration of the basic geodetic network in the area of Poland with reference stations of the ASG-EUPOS system stage IV develop and adjustment of GNSS observations" (not available) and are given in ITRF2005 for epoch 2011.0. For this reason, these points cannot be directly used in the process of the determination of a quasigeoid model fitted to GNSS/levelling data, and they are used only for comparison purposes.
The normal heights in the Polish height system (marked as PL-EVRF2007-NH) for the main points of the EUVN network were determined during the international adjustment of the EUVN and EUVN DA network in the European height system-EVRF2007-NH. The heights of the remaining EUVN points were determined in the basic national height network adjustment.
The normal heights of the POLREF network points were determined by readjusting the measurements carried out in 1995-1997 using spirit levelling in reference to the basic levelling network.
It is clear from the presented description that the GNSS/levelling height anomalies determined for individual sets do not constitute uniform data. Hence, in further work, it was assumed that the analyses would be performed independently for each set. Generally, therefore, the following GNSS/levelling datasets have been adopted at this stage ( Fig. 1b): • ASG_Stat-94 points; • ASG_Ecc-95 points; • POLREF-297 points; • EUVN-42 points.
The calculations also used the global geopotential model SGG-UGM-2 up to degree and order of 2190 (Liang et al. 2020). The Moho depth model for the European plate (Grad et al. 2009) was used to define the volume (Fig. 2b), and the digital elevation model DEM SRTM v4.1 (Jarvis et al. 2008) was used to define the volume Ω (Fig. 2a).
Based on the SRTM model, three numerical elevation models were prepared and directly used in the calculations: • 100 × 100-m-resolution model, which covers the mountainous areas of southern Poland up to the latitude 52 • N; • 500 × 500-m-resolution model, which covers the lowlands of northern Poland (latitudes above 52 • N); • 1000 × 1000-m-resolution model, which covers the GGI model construction area. Models with resolutions of 100 × 100 m and 500 × 500 m were used for calculations in the vicinity of measurement points (up to 10 km distance). For calculations at further distances, the 1000 × 1000-m-resolution model was used.

The conducted analyses
As mentioned above, the GGI method enables the development of a quasigeoid model fitted to GNSS/levelling data in two ways. In the first version the gravimetric model is determined using the GGI method (fitted to the GGM used), and then it is transformed into GNSS/levelling data. In the second version, the GGI method directly determines the model fitted to the GNSS/levelling data. It is worth noting that in this version, possible errors in the GNSS/levelling data are more difficult to detect and directly distort the quasigeoid

Gravimetric quasigeoid model development
The gravity data used to develop the gravimetric quasigeoid model are presented in Fig. 1a.
The applied solution assumes the use of a regular grid of points for which observation Eq. (17), related to the disturbance potential, will be formed. The points of this grid are marked with grey dots in Fig. 1b. The horizontal extent of volumes Ω and is marked in Fig. 1a and b with a black, dashed line. Although the final regular grid of the quasigeoid model with a resolution of 0.01 • × 0.01 • was developed for an area between 48 • and 56 • N and 13 • and 26 • E, the values obtained from the GGI model can be significant only in the area covered by the gravity data. Hence, in Fig. 1a, the area inside which the height anomalies were determined using the full model values T r GGI (Eq. (12)) is marked with a black continuous line. At a distance of about 10 km from the indicated line, outside this area, a transition zone was defined in which the T r GGI values are gradually attenuated to zero. For the remaining outer area, height anomalies from the GGM were used ( T r GGI = 0 ). The border of Poland, marked in Fig. 1a and b with a solid grey line, is entirely within the area for which the quasigeoid values were determined using the GGI method. It is worth noting that the gravity data only partially cover the area beyond the border of Poland, which is unfavourable due to the edge effect. However, in the mountainous area (near the southern border), gravity data coverage (although with a lower resolution) is also provided beyond the border. Analyses regarding the edge effect are carried out later in this study. In Trojanowicz et al. (2021), it was shown that the accuracy of the gravimetric quasigeoid model developed using the GGI method depends on the adopted reference model of the density 0 for the mountainous area. In order to determine if such relationships occur for the lowland area, which represents most of the territory of Poland, a number of models of the residual potential ( T r GGI ) and the heights of the gravimetric quasigeoid model at GNSS/levelling points were determined according to Eq. (18). Calculations were performed for values of 0 in the range 0 − 2400 kg m 3 , with a step size of 200 kg m 3 , as well as for the values 2500 kg m 3 and 2670 kg m 3 (the density 0 was calculated according to Eq. (7)). For each group of points, the differences Δ grav and their standard deviations, which represent the accuracy parameter, were calculated: where Grav is the height anomaly calculated from Eq. (18), and GNSS∕lev is the height anomaly calculated using GNSS/levelling data. Figure 3 shows the dependence of the determined values on the assumed reference densities.
The graphs presented in Fig. 3 indicate a very weak relationship between the accuracy of the model and the values of the reference densities 0 . The accuracies change only in the range of a few tenths of a millimetre, and the smallest variability occurs for the sets with the highest accuracy (ASG_Ecc and ASG_Stat). Based on this, the model developed for a reference density 0 = 1000 kg m 3 was adopted for further work. Let us note, however, that with such small differences, any value from the entire density range studied could be utilised without substantially changing the final result.
Due to the occurrence of outliers in the POLREF set, Fig. 3 shows the standard deviations of Δ grav for the POLREF* set, which represents the POLREF set after the outliers have been removed. Analyses related to this procedure are described later in this study.
For the selected reference density ( 0 = 1000 kg m 3 ), the model of the residual potential ( T r GGI ) and the gravimetric model of the quasigeoid were determined for the 0.01 • × 0.01 • grid and in the above-mentioned range. The residual potential allows us to calculate the differences between the gravimetric quasigeoid model developed with the GGI method and the quasigeoid model determined on the basis of the SGG-UGM-2 model. A map of these differences is presented in Fig. 4a. Figure 4b presents the developed gravimetric quasigeoid model.

GNSS/levelling data quality analysis and gravimetric quasigeoid model accuracy assessment
The developed gravimetric quasigeoid model is a reference surface independent of the levelling data; therefore, it can be used to assess the uniformity of GNSS/levelling datasets. On the other hand, a verified and internally consistent set of points with known GNSS∕lev values enables the reliable estimation of the accuracy of the developed gravimetric model. Therefore, the assessment of the developed model accuracy was carried out, along with the analysis of the occurrence of gross errors (outliers) in the GNSS/levelling datasets. The analysis is based on the calculated differences Δ grav (Eq. (19)) for each of the analysed groups of GNSS/levelling points. Table 1 presents the most important statistics concerning these differences as well as shows the maximum deviation of Δ grav from their mean value ( Δ mean grav ): It was assumed that for outliers, the deviation from the mean value in a given set ( v Δ grav ) exceeds three times the standard deviation of the Δ grav values: It was found that in the sets ASG_Stat, ASG_Ecc and EUVN, such observations do not occur (for all points of these sets, relation (21) was satisfied). However, in the POLREF set, at least one such observation was found (the maximum deviation max   Table 1 is almost 10 times greater than the standard deviation). Therefore, an iterative procedure for identifying and eliminating outliers was carried out for this set. In each iteration, deviations were calculated using Eq. (20), and the point with the largest max | | | v Δ grav | | | deviation, which did not satisfy relation (21), was removed. This operation was repeated until relation (21) was satisfied by all points. In this way, 8 GNSS∕lev points were identified as outliers and removed from the POLREF set, leaving a set of 289 points, which is referred to as POLREF * . The most important statistics concerning the Δ grav values for this set are presented in the last row of Table 1.
The results contained in Table 1 also allow a preliminary assessment of the compatibility of the developed gravimetric quasigeoid model with the GNSS/levelling data. The parameters describing this compliance include the following: 1. The standard deviation of the differences Δ grav , which describes differences of a random nature. This quantity is also considered one of the accuracy parameters of the determined model (in addition to the root mean square error (RMSE)). 2. The vertical bias and tilt of the gravimetric quasigeoid relative to the GNSS/levelling quasigeoid, which describes differences of a systematic nature.
Some of the values of the indicated parameters are included in Table 1. In terms of systematic differences, there is a shift of the gravimetric quasigeoid surface in relation to the GNSS/levelling data for individual groups of points with a value of −0.3 to −1.2cm ( Δ mean grav values). In relation to the ASG-EUPOS network sets, the shift is −0.3cm (ASG_Stat) and −0.5cm (ASG_Ecc). The shift relative to the EUVN and POLREF* points is slightly larger ( −1.1 and −1.2cm , respectively). While the shift of the GNSS/levelling quasigeoid in relation to the gravimetric quasigeoid was expected and should not affect the quality of the model fitted to the GNSS/levelling data, the differences in the shifts for individual sets pose a certain problem. Note that the differences in the pairs of sets ASG_Stat and ASG_ Ecc and EUVN and POLREF* are small ( 2 − 3mm ) and can be considered acceptable. However, the differences between these pairs are about 8 mm. Such differences cannot be neglected and are a problem with the selection of a set of points that can be used to build a fitted quasigeoid model.
The standard deviations of the differences Δ grav shown in Table 1 indicate the accuracy of the developed gravimetric quasigeoid model. Since they also depend on the accuracy of the GNSS∕lev values, they allow a certain evaluation of these data. The lowest standard deviation value occurs for the points of the ASG_Ecc set (1.2 cm), which suggests that the points in this set have the highest accuracy. It is worth noting, however, that the standard deviations shown in Table 1 may be overestimated in the case of other than bias systematic errors (e.g. the tilt of the analysed surfaces relative to each other). The horizontal distributions of the differences Δ grav for individual groups of points are shown in Fig. 5.
The distributions of Δ grav values for the points of the ASG_Stat and ASG_Ecc sets shown in Fig. 5 are very similar, which was expected. For all sets, there are visible changes in the differences, from clearly negative in the western part of the map towards positive in the eastern part. These changes are most pronounced for the EUVN and POLREF* sets. Changes in the north-south direction are also clearly visible for all sets. Especially in the eastern part of the country, the analysed differences are clearly positive in the north and south areas of the country, while in the central part of the country, these values are negative. This suggests that the differences between the GNSS/levelling data and gravimetric quasigeoid model are also of a local nature. This should be taken into account when selecting the function that transforms the gravimetric quasigeoid model into GNSS/levelling data.

Selection of the function that transforms the gravimetric model into GNSS/ levelling data and determination of the fitted quasigeoid model
The differences Δ grav discussed in the previous subsection will be used to determine the analytical model of differences Δ mod , which will be used to calculate the quasigeoid model fitted to GNSS/levelling data based on the following equation: Before selection of the form of the Δ mod difference model, it was assumed that the function defining this model should be a function of low degree with a small number of parameters. The functions given in Table 2 were selected for analysis. For individual groups of GNSS/levelling points, the transformation models given in Table 2 were built and the Ai GNSS∕lev values according to Eq. (22) were determined (i is an indicator of the transformation model). The standard deviations of the differences v Ai = GNSS∕lev − Ai GNSS∕lev were calculated, which are presented in Table 2. Among the analysed transformation models, for models e) and g) the lowest standard deviations of the differences v A i for all sets are presented. The model with fewer number of parameters is model e), and this model was chosen to transform the gravimetric quasigeoid model into a fitted model. Thus, finally, we write: It should be noted that the points of the ASG_Stat and ASG_Ecc sets are the points of the Polish fundamental horizontal geodetic control network (the most accurate geodetic network), and the ASG_Ecc set shows the best fit to the developed gravimetric quasigeoid model in all analyses. Taking this into account, the gravimetric quasigeoid model was fitted according to Eq. (23) to points of the ASG_Stat and ASG_Ecc sets.

(23)
A GNSS∕lev = Grav + a o + a 1 x + a 2 y + a 3 x 2 + a 4 xy Table 2 Standard deviations of the differences v Ai = GNSS∕lev − Ai GNSS∕lev for different transformation models and for particular groups of points (i is an indicator of the transformation model) (f) a o + a 1 x + a 2 y + a 3 y 2 + a 4 xy 1.29 1.04 1.14 1.37 (g) a o + a 1 x + a 2 y + a 3 x 2 + a 4 y 2 + a 5 xy 1.24 1.05 1.13 1.27 (h) a o + a 1 cos cos + a 2 cos sin + a 3 sin + a 4 sin 2 1.32 1.11 1.20 1.29 The calculated corrections that are used to transform the gravimetric model into GNSS/ levelling data are presented in Fig. 6.
The developed model was compared with height anomalies for individual sets of points. Table 3 contains the most important statistics concerning the differences

Development of a quasigeoid model fitted according to procedure B ( B GNSS∕lev )
As presented in the description of the GGI method, in accordance with Eq. (14), the quasigeoid model fitted to the GNSS/levelling data is directly determined. Therefore, the selection of a reliable set of GNSS/levelling points is an important part of this approach. In the previous paragraph, the quality of the available sets of these points was analysed based on the gravimetric quasigeoid model, resulting in verified GNSS/levelling datasets. Due to the fact that as many GNSS/levelling points as possible should be used to build a quasigeoid model fitted in accordance with procedure B, it was decided that the points of the POLREF* set should be used to build the model, while the points of the other sets will be used to assess the quality of the model. The same GGM (SGG-UGM-2), gravity data and DEMs that were used in the previous calculations were utilised. Similarly to the previously determined model, this model was compared with height anomalies for individual sets of points. The basic statistics concerning the differences Δ B = GNSS∕lev − B GNSS∕lev are presented in Table 4.

Edge effect analysis
As mentioned above, the gravity data coverage for the land area of Poland is very high and completely sufficient. It is very unfavourable, however, that there is a lack of gravity data outside the eastern, western and northern parts of the land border and very little data coverage  Table 4 The basic statistics of the differences in the area of the Baltic Sea. This may be related to the slightly lower accuracy of the developed model in border areas due to the so-called edge effect. In order to evaluate this effect, another gravimetric quasigeoid model (created according to procedure A) and a model fitted to GNSS/levelling data (created according to procedure B) were determined, limiting the area of the gravity data and using points from the POLREF* set. The data used in these calculations are shown in Fig. 7. The Ω and volumes were left unchanged. To determine the gravimetric model, the previously used grid of points was adopted for the construction of Eq. (17). The heights of the gravimetric quasigeoid ( EEGrav ) and the model fitted to GNSS/levelling data ( B EEGNSS∕lev ), calculated based on limited datasets (at regular grids of points, as in the previous calculations), were compared with the corresponding height anomalies developed previously (using all data) to obtain the following differences: The unfavourable influence of the edge effect should be most noticeable near the border of the area covered by the gravity data. Thus, the EEGrav and B EEGNSS∕lev values for points close to the border of the gravity data used to calculate these values (solid black line in Fig. 7) will be much more distorted by the edge effect than the Grav and B GNSS∕lev values calculated at the same points but using a wider range of data (this data range is shown in Fig. 1). Therefore, it was assumed that in these analyses, the edge effect is represented by the differences calculated using Eqs. (23) and (24). First, the points located close to the border of the area covered by the gravity data were analysed. In this area, nine zones consisting of 5 − km-wide, narrow strips located from 15 km outside the gravity data border (blue dashed line in Fig. 7) to 30 km inside the area covered by the gravity data have been defined. This internal border is marked with a dashed red line in Fig. 7. The tenth zone consists of points located inside this border. The standard deviations and variation ranges of the differences calculated using Eqs. (23) and (24) were calculated for the points located in the individual zones. In order to assess the viability of the GGI model in the zones close to the gravity data border, basic statistics were also calculated for the following differences: where GGM is the height anomaly calculated from the global model used.
The standard deviations of the differences calculated using Eqs. (23)-(26) are presented in Fig. 8a, while their maximum absolute values are presented in Fig. 8b. In these diagrams, the distances from the border of the gravity data for the zones lying outside this border are negative, while the distances for the zones lying inside the gravity data border are positive. The values of statistics for a single zone are constant.
When the variability of the statistics presented in Fig. 8 is assessed, it is clear that the edge effect reduces the accuracy of both analysed models. However, it is equally clear that right at the border of the gravity data, the influence of the GGI model is positive. The statistics of the differences calculated from Eqs. (23) and (24) presented in Fig. 8 are distinctly smaller than the standard deviations of the differences calculated from Eqs. (25) and (26). For the standard deviation values, this trend is visible already in the first zone next to the border and for the maximum magnitudes of the absolute values in the zone 5 − 10km from the border. Outside the limit of the gravity data, the GGI model in a narrow band can be considered neutral in relation to the global model, and at further distances, it has a negative effect. This means that outside the gravity data border, the GGI model should not be used, and it is better to leave the values calculated from the global model uncorrected. This justifies the use of the procedures adopted in this study.
In addition, the horizontal distributions of the absolute values of the differences calculated from Eqs. (23) and (24) are presented in Fig. 9 ( Fig. 9b). As can be seen on the presented maps, the greatest differences occur in the vicinity of the southern border of the elaboration area. This is the area closest to the mountainous areas (cf. Fig. 2a). Therefore, in mountainous areas, it should be expected that the edge effect will have a more unfavourable impact. Since the final quasigeoid model, due to its use as a correction surface in GNSS measurements, should best match the GNSS/levelling data, it was assumed that the final model should be built on the basis of both developed models. The A GNSS∕lev model was adopted as the base model, and the final height anomalies (marked as GGI ) at points in the area of Poland were calculated as the average of the values from models A and B according to the following equation: This value takes into account aforementioned the average bias between models A and B equal 0.0085 m. Beyond the borders of Poland, the final height anomalies were taken from model A: A map of the GGI height anomalies is presented in Fig. 10. The most important statistics concerning the developed model are presented in Table 5. The differences Δ A GGI = A GNSS∕lev − GGI determine the effect of model B on the final model values. A map of these differences is shown in Fig. 11, and the most important statistics concerning these differences are also presented in Table 5.
The final model was compared with height anomalies for particular sets of points. The most important statistics concerning the differences Δ GGI = GNSS∕lev − GGI are presented in Table 6.
From among the sets shown in Table 6, only the EUVN set was not directly utilised in the development of the final model in any stage of the work. Therefore, only for this set, the values presented in Table 6 constitute an independent assessment of the model accuracy. Fig. 10 The final, fitted quasigeoid model developed using the GGI method ( GCI ) . (Color figure online) The GGI model was accepted by Polish Head Office of Geodesy and Cartography as the current quasigeoid model for the area of Poland and is available to users under the name PL-geoid2021. The evaluation of this model provided by the GUGiK, based on 24 and 62 GNSS/levelling points indicated its accuracy at ±1.9cm and 2.6cm respectively (Kalinczuk-Stanałowska and Mielczarczyk 2022; Kryński et al. 2023).

Conclusions
The analyses conducted lead to a number of important conclusions concerning the application of the GGI method.
First, with regard to the gravimetric quasigeoid model, a low correlation (which can even be considered insignificant) between the accuracy of this model and the reference density 0 is indicated. On the basis of this model, GNSS/levelling datasets were also assessed to identify outliers. The accuracy of the gravimetric model, estimated using four GNSS/ levelling datasets, was in the range from ±1.2 − ±1.7cm (in terms of the standard deviation of the differences between the measured and model-determined height anomalies).
Two quasigeoid models fitted to GNSS/levelling data were developed according to two procedures, which provided very similar results. The average model was adopted as the final model. Due to the availability of the gravity data, the values obtained by the GGI method, for both procedures, basically cover only the area within the land borders of Poland, slightly exceeding this area from the north and south. The final regular grid of  the quasigeoid model, with a resolution of 0.01 • × 0.01 • , was developed within the regular limits between 48 • and 56 • N and 13 • and 26 • E , but the values obtained from the GGI model are meaningful only in the area covered by the gravity data. For the remaining area, the height anomalies from the GGM were utilised. The accuracy of the final model (estimated in the same way as the gravimetric model) ranges from ±1.0cm (for the ASG_Ecc set) to ±1.2cm (for the ASG_Stat, EUVN and POLREF* sets). However, it should be emphasised that this assessment is not fully independent, because only the EUVN set was not directly used in the development of the final model in any stage of the work. Provided by the GUGiK the evaluation of this model based on 24 and 62 GNSS/levelling points indicated its accuracy at ±1.9cm and 2.6cm respectively (Kalinczuk-Stanałowska and Mielczarczyk 2022; Kryński et al. 2023). These discrepancies confirm the indicated by Szelachowska and Kryński (2014) difficulties in evaluation of the geoid and quasigeoid models. In order to estimate the impact of the lack of gravity data beyond the western, eastern and part of the northern boundaries on the quality of the quasigeoid model, analyses of the edge effect were performed. In both analysed cases (for the gravimetric quasigeoid model and for the model fitted according to procedure B), the edge effect turned out to be similar and should be considered small. In both cases, next to the gravity data border, there is a clear improvement in the quality of the quasigeoid model determined by the GGI method in relation to the global model used. At a distance of about 5-10 km from the border, the standard deviation of the differences between the models distorted and undistorted by the edge effect decreased to less than ±1.0cm , and at a distance of 25 − 30km , this standard deviation decreased to less than ±0.5cm.