Determination of the Selected Gravity Field Functionals by the GGI Method: A Case Study of the Western Carpathians Area

: In this paper, some features of the local disturbing potential model developed by the GGI method (based on Geophysical Gravity Inversion) were analyzed. The model was developed for the area of the Western Carpathians covering the Polish–Slovak border. A detailed assessment of the model's property was made regarding the accuracy of the disturbing potential values (height anomalies), gravity values, complete Bouguer anomalies (CBA), and differences between geoid undulations and height anomalies (  −  ). Obtained accuracies of the GGI quasigeoid model (in terms of standard deviation of the residuals to the reference quasigeoid models) were at the level of ± 2.2 cm for Poland and ± 0.9 cm for the Slovak area. In terms of gravity, there was shown dependence of the accuracy of the GGI model on the digital elevation model (DEM) resolution, the point height, the density of gravity data used, and used reference density of topography model. The best obtained results of gravity prediction were characterized by an error of approximately 1  . The GGI approach were compared with classical gravity prediction methods (using CBA and topographic-isostatic anomalies supported by Kriging prediction), getting very similar results. On the basis of the GGI model, CBA and differences (  −  ) were also determined. The strong dependence of resolution of the CBA model obtained by GGI approach, on the size of the constant density zones, has been demonstrated. This significantly reduces the quality of such a model. The crucial importance of the topographic masses density model for both determined values (CBA and (  −  )) was also indicated. Therefore, for determining these quantities, all available information on topographic


Introduction
Gravity measurements in geodesy are used in many ways. One of the basics is their use in defining and determining the height of geoid and quasigeoid as a reference surface, as well as determining the height of points in a given height system in classical spirit levelling. For the purpose of developing geodetic levelling networks, gravity values are measured, interpolated, or determined from the appropriate gravity potential model [1,2]. Note that the classic approach to the problem of geoid and quasigeoid modelling using gravity data, basically involves construction of the disturbing potential model in the external space [3][4][5]. Such a model could also be used to determine the gravity values. However, the modelling procedures themselves, which runs using global geopotential models, are usually complicated and require the introduction of various corrections and reductions to the data [6,7], causing that gravity models that could be used, for example, to develop levelling networks, are not directly (during this procedure) determined. Obviously, techniques for building such gravity models e.g. [8][9][10], most often in the form of regular grids, used in classical geoid and quasigeoid modelling procedures, could be also used to determine the gravity values. Nevertheless, we can state that classical geoid and quasigeoid modelling techniques do not allow for the construction of models that allow determination of height anomalies, geoid undulations, or gravity directly from the same model. Individual mentioned characteristics of the gravity field are determined in separate, though often related, procedures.
Models that have the ability to determine from the same model the values of different characteristics of the gravity field are global geopotential models. Their accuracy is currently very high, but particularly in mountainous areas, they are still much lower than the aforementioned classical solutions and requires introduction of the appropriate corrections. Thus, they are primarily used as reference models on the basis of which more precise, local, and regional models of the gravity field are built.
A model that allows to determine various parameters of the gravity field is the GGI model (based on Geophysical Gravity Inversion) [11,12]. The model is built on basis of a rare network of GNSS/levelling height anomalies and a dense network of gravity points. Based on this solution, height anomalies characterized by very high accuracy can be derived, at the level of accuracy of the classical solutions [13,14]. Directly from the GGI model, can be also easily determined other characteristics of the gravity field, such as gravity values, differences between geoid undulations and quasigeoid heights, and complete Bouguer anomalies (CBA). Out of the three characteristics mentioned above, only preliminary analyses of the accuracy of the gravity values determined from the model were carried out [15]. Hence, this study analyzes in detail the accuracy of gravity values and is an attempt to assess the quality of the differences of the geoid undulations and quasigeoid heights as well as the CBA values determined on the basis of the GGI model.

The GGI Model
The GGI method is based on a local model of disturbing potential, which consists of three components ( Figure 1) [11,12].
The first component is the potential of the topographic masses, which lies above the geoid, included in volume , with density distribution function . The second component is the potential , of the disturbance masses occurring between the geoid's surface and the compensation level, included in volume . The density distribution function of the masses will be marked as . Volumes and are horizontally limited to the area of elaboration, while the data used for calculations also contain information about the density distribution outside these volumes. Consequently, the potential (i.e., external disturbing potential) is introduced, whose role is to model the part of the data resulting from masses outside of the study area and to cover long-wavelength errors of gravity data and systematic errors of both levelling and GNSS data. This indicates that potential should be of the trend type and can be modelled using harmonic polynomials of a low degree. Finally, we can write: Components and are given by Newton's integrals: where r is the distance between the attracting masses and the attracted point P, G is Newton's gravitational constant, and and are elements of volumes. In our calculations we use component in the form: where , , are coordinates of the point P and , … , are determined coefficients. We can now formulate an inversion problem: Find the density distribution functions and in defined volumes and and find the coefficients of the polynomials modelling the potential to satisfy equation (1) for given data. Using linear inversion [17], the solution of the task requires discretisation of the continuous 3D functions and . Volumes and are divided into finite volume blocks and a constant density is assigned to each of the blocks. In previous studies, we used volume , defined by digital elevation model (DEM) in the form of rectangular blocks that are grouped into zones of constant, searched density. The volume is defined as a slab with a thickness that is nearly the same as the depth of the compensation level and consists of one or many layers of constant density blocks. In previous and current studies, a one-layer version was used. In the horizontal plane, volumes and exceed the border of data occurrence. Our calculations are performed in the local Cartesian coordinate system. The Z-axis of the coordinate system is directed towards the geodetic zenith at the origin point. The X and Y axes lie on the plane of the horizon and are directed towards the north and east, respectively. The definition of the coordinate system enables the determination of the and volumes in the form of rectangular prisms for which the solutions of Newton's integrals are presented by [18,19].
Equations (2) and (3) can now be written as follows: = − . The model parameters (vector and the coefficients , … , of the polynomial that defines potential ), are determined by the least-squares method. It should also be mentioned that calculations can be performed with or without (version use in the article) a global geopotential model.

Calculation of Selected Functionals from the GGI Model
From this model, we can determine various functionals of the gravity potential. In the following section, we will focus on determining and analyzing the following:  The height anomalies ( ), which are determined based on disturbing potential;  The gravity values;  The complete Bouguer anomalies; and  The difference between the geoid undulations (N) and the height anomalies.

The Height Anomalies
Because model (1) is disturbing potential model, the value of the height anomaly, at a particular point P, can be determined by Brun's formula [3]: where is the normal gravity at point on telluride ( Figure 2).

The Gravity Values
If we neglect the slight differences between the real and normal plumb lines, then by differentiating the disturbing potential (1) along the normal to the equipotential surface at the observation point (hence differentiating with respect to height H) we get the expression defining the gravity disturbance: Now, if we have the position of the point ( , , ℎ), e.g., from the GNSS measurement, we can determine the value of the normal gravity ( ) at the point, and then the gravity value ( ): These values can be used, for example, to determine system corrections for precise spirit levelling instead of performing gravity measurements.

The Complete Bouguer Anomalies
The complete Bouguer anomaly is defined as [20,21]: where Δ is a free-air anomaly on a geoid, is gravity on the geoid, is normal gravity on the ellipsoid, and is topographic reduction. Used in (10), the value is calculated based on gravity at point P ( ) by introducing free-air reduction [4] (Figure 2): where = − is free-air reduction, is orthometric height and is vertical gravity gradient. For many practical purposes, vertical gravity gradient is approximated by normal gravity gradient . For areas with small differences between the geoid and the quasigeoid, the orthometric height can be replace with the normal height ( ). Free-air reduction can thus be defined as: Topographic reduction ( ) is expressed as: where l is the distance between the attracting masses and the attracted point P, is density, and is integration area. The topographic reduction is classically determined based on DEM grids with different resolutions. This resolution is very high near the gravity point and it decreases with distance. Integration of topographic masses is also carried out to a certain, defined distance around a gravity point, often 167 km [22], introducing corrections for further zones only for large-scale and precise geophysical studies [23].
When calculating the CBA based on the GGI model, we can use the classical procedure that was presented earlier. By putting equation (11) into equation (10) and replacing value with (9) and (8), we get: wherein where ℎ is ellipsoidal height.
In the GGI model (1), topographic masses are represented by masses Ω , and topographic reduction ( ) is approximated by the vertical derivative of the potential: By putting equations (16) and (15) into equation (14) we get the simple expression that defines the complete Bouguer anomaly: Note that complete Bouguer anomalies calculated by equation (17), assume approximation (16), in which the integration takes place in relation to the entire volume Ω (not to a certain distance around the gravity point) with densities determined from the GGI model. These parameters are different from those described earlier in the classic approach of determining the topographic reduction.

Figure 2.
Schematic relationship between orthometric ( ), normal ( ), and ellipsoidal (h) height (slight differences between the real and normal plumb lines as well as normal to ellipsoid are neglected). Detailed explanations are given in the paper.

The Difference Between the Geoid Undulations and the Height Anomalies
An exact formula can be used to determine the difference between geoid undulations and height anomalies, resulting from a comparison of normal ( ) and orthometric ( ) heights (Figure 2) e.g [19]: where ̅ is the mean gravity between geoid and point P defined by (19), and ̅ is the mean normal gravity between ellipsoid and telluroid defined by (20).
Although the value ̅ , which is used to determine the orthometric heights, can be calculated by various methods, its determination is difficult. In practice, this value is often calculated on the basis of Poincaré and Prey reduction [4]. In this approach, the mean gravity ( ̅ ) can be calculated at point P', lying in the middle height of point P in three steps: 1. Remove the effect of topography from the gravity value at point P; 2. Introduce the free air reduction moving point P down to point P'; and 3. Restore the effect of topography at point P'.
The calculations can either use the full topographic reduction (13) or only take into account the so-called the Bouguer plate (steps 1 and 3) and use the real or normal gravity gradient in the free-air reduction (step 2).
Using the GGI model, the value of ( ̅ ) can be calculated directly at point P' based on (9). Value ̅ is in practice calculated at the middle height of telluroid (point Q) above the ellipsoid.

Characteristics of the Data and the Area of Elaboration
The area of study included part of the Western Carpathians, which is located on the border of Poland and Slovakia. In our calculations, for definition of the Ω volume, four versions of DEM resolution were used: 0.1 × 0.1 , 0.25 × 0.25 , 0.5 × 0.5 , and 1 × 1 . All DEM grids were determined based on the Shuttle Radar Topography Mission (SRTM) model, which was shared by the CGIAR-CSI Consortium for Spatial Information (http://srtm.csi.cgiar.org/contact-us/) with a resolution of 3′′ × 3′′. Horizontal range of digital terrain model used in the calculations is shown in Figure 3. In all calculations where the DEM model with resolution 0.1 × 0.1 was used, two DEM models were used in order to reduce the computational effort. The 0.1 × 0.1 model was used for a square with side of 10.500 , in the center of which there was a determined point and 0.5 × 0.5 was used for the rest of the area. Volume was defined based on the Moho depth model for the European plate [24]. In the calculations, two constant density values were assumed as the reference density model for the volume Ω: = 2670 ⁄ and = 2200 ⁄ (this value is close to the mean density of topographic masses for lowland areas of Poland [25]) as well as the global topographic mass density model UNB_TopoDens [26]. The reference density model for the volume ( ) was adopted as "negative density", which balanced the topographical masses of the volume Ω. The density of a separate block of volume was calculated based on the equation

= −
, where , are mean height and reference density of zone i of the Ω volume, which are situated directly above block j of volume ; and ℎ is the height of block j. Based on the quasigeoid models, PL-geoid-2011 [27] for the area of Poland and the model GMSQ03C [28] for the area of Slovakia, a regular grid of points with known height anomalies has been designated. These points were divided into two sets: The first was used to build the GGI model and consisted of 95 points (about 1 point for 200 ); while the second set consisted of 81 test points and was used for assessment of the accuracy of the quasigeoid model. In the calculations, we also used a set of gravity points that were divided into two sets: A set of 10,532 known points and a set of 10,513 test points (the approximate density of the data and test points was about 1 point for 2.4 ). For the Polish part of the elaboration area, the gravity points were made available by the Polish Geological Institute National Research Institute (PGI-NRI). For the Slovak part, they were made available by Geofyzika, a.s. Brno, ČGS -Geofond. Both groups of data points (height anomalies and gravity points) evenly cover the entire study area.
Taking into account the various datasets from the Slovak and Polish areas (various quasigeoid models and various gravity datasets), it was decided to use the GGI model that separated the Polish ( ) and Slovak ( ) parts: Both parts ( and ) share parts and , and have independent part . This separation is caused by the assumption that the differences in the used data are the trend type.

Assessment of the Accuracy of the Gravity Values and height Anomalies Determined from the GGI Model
The accuracy analyses that we carried out examined four factors, as follows: DEM resolution, constant density zones dimensions, density of gravity data, and the initial density of topographic masses. The results of the analyses are presented in two parts, which will be described in following sections. Section 4.1.1 covers the first two factors and Section 4.1.2. covers the other two. This Section also includes a comparison with classical methods of gravity prediction.

The Impact of DEM Resolution and Constant Density Zones Dimensions on Determined Gravity Values and Height Anomalies
In order to assess the significance of DEM resolution for the accuracy of the GGI model, the calculations were done using the four resolutions listed in Section 3. Each of the volumes Ω and were divided into 3600 constant density zones with dimensions of about 4 × 4 (7200 estimated constant densities). Calculations were performed assuming = 2670 ⁄ as the reference density model for the volume Ω. The disturbing potential model was built based on the known points and according to the procedure described in Section 2.1. Due to the adoption of various datasets for the Polish and Slovak parts, the accuracy parameters were determined separately for the areas of both countries. For evaluating the gravity prediction, the test points were also separated in terms of their heights, giving the dependence of the accuracy of the model on this parameter. In respect to this, three groups of points were separated. The first group ( ) represents points located in lowland terrain up to 300 m high, the second group ( ) represents points located in the area of low mountains with heights between 300 m and 700 m, and the third group ( ) has points located in the area of higher mountains with a height above 700 m. Differences Δ = − and Δg = g − g were determined for test points, where is the height anomaly determined from the adopted quasigeoid model, g is the measured gravity, and g are height anomaly and gravity calculated from the GGI model according to equations (7) and (9) respectively, wherein the disturbing potential is calculated by equations (21,4,5,6), and the gravity disturbance value is determined by equation (8).
Since the average values of the analyzed differences were close to zero, the standard deviations of these differences turned out to be almost the same as their RMS errors. Thus, Table 1 and 2 show the standard deviations of these differences as indicators of the GGI model fit to the appropriate quasigeoid models and to the measured gravity. The tables also show the mean heights of the test points and the corresponding accuracy parameters of the EGM2008 model [29], accepting Δ = − and Δg = g − g .  Commenting on the above results, we should first point out the lack of dependence between the DEM resolution and the accuracy of the GGI model in terms of the disturbing potential (quasigeoid accuracy). Standard deviations of differences Δ are at a similar level for all analyzed DEM resolutions. Looking at gravity disturbances, the results are slightly different. There is a clear relation between the used DEM resolution and the height of the test points and the accuracy of the modelwhich is as expected. Ultimately, it can be estimated that the average accuracy of the gravity disturbances (and thus gravity) determined from the GGI model for the best versions is about 1 . There are also noticeable differences in standard deviations of Δ for the Polish and Slovak parts. The designated GGI model is better suited to the Slovak quasigeoid model. This is probably related to the methods used for construction of both quasigeoid models. In the process of developing the Polish quasigeoid model, gravity data were not used directly (this is the EGM2008 model fitted to the GNSS/levelling points [27]). Thus, the presented differences may indicate greater accuracy for the Slovak quasigeoid model. However, this is only a suggestion, of which the results are not settled.
In order to investigate the significance of the constant density zones dimensions to the accuracy of the model, calculations were carried out by assuming zones of constant density in used earlier 4 × 4 and three additional versions: 8 × 8 (900 constant density zones), 6 × 6 (1600 constant density zones), and 3 × 3 (6400 constant density zones), with all other parameters left unchanged. Constant density zones dimensions are presented in Figure 4. Table 3 presents the results of these calculations, giving statistics for all test points.  Taking into account the height anomalies, it should be stated that there is no dependency of the accuracy of their determination on the analyzed dimensions of the constant density zones. Again, there are differences in the accuracy of determining gravity. Only for the two versions, the used earlier 4 × 4 and smaller 3 × 3 , these accuracies are the highest (approximately 1 ). Enlarging the zones of constant densities visibly reduces the accuracy.

The Impact of the Density of the Gravity Data and the Initial Density of Topographic Masses Model on the Determined Gravity Values and Comparison with the Classical Methods of Gravity Prediction
The high accuracy of gravity obtained from GGI models shown earlier (Tables 1-3) is also the result of the very high density of gravity data used in the calculations. Along with the decreasing number of data points, the accuracy of the model should be expected to decrease. The estimation of the influence of gravity data density on the accuracy of the model in terms of height anomalies was initially estimated by [13]. We will now focus on assessing the accuracy of the gravity. The analyses were carried out for DEM with resolution of 0.1 × 0.1 and constant density zones dimension 4 × 4 . The research was supplemented with the classical methods of gravity data interpolation, using CBA and Pratt-Hayford topographic-isostatic anomalies (hereinafter referred as Pratt anomalies and marked as Δ ). In the classical approach, the calculations were made according to the scheme further marked as CBA/Kriging and Pratt/Kriging:

Based on CBA
and Δ , we build a dense grid of these anomalies using the Kriging method (we used for this step the Surfer 16, Golden Software). 3. Using the grid specified in the previous point, we interpolate the CBA and Pratt anomalies at test points (CBA ; Δ ). 4. We compare the interpolated and computed anomalies at test points (Δg = CBA − CBA ; Δg = Δ − Δ ), and then determine the accuracy indicators.
In the analyses, CBA has been calculated on the basis of Equation (10), whereas Pratt anomalies were defined as [4]: where is compensation attraction according to Pratt-Hayford system: where D is depth of compensation level and Δ is the density difference between the constant density and the actual density of the compensation model. In the calculations, we adopted a compensation model using the constant density of topographic masses used for Bouguer reduction, for which the mass equality has the form [4]: Hence: It is worth noting that in the GGI model, the initial density model ( ) of the volume is determined from an analogous relationship (see Section 3). This means that in the modelling process, after taking into account the initial density model ( and ), we use a kind of topographic-isostatic gravity disturbance and the component -represents the topographic-isostatic gravity disturbances after optimization of the densities ( and ) with respect to the measured data.
Both, CBA and Pratt anomalies were calculated using normal gravity gradient in free-air reduction, and an integration radius of 167 was used for topographic and isostatic reductions. For isostatic reductions, the standard depth of the compensation level = 100 was also assumed. The calculations used a DEM with a resolution of 0.1 × 0.1 for a square with side of 10.500 , in the centre of which there was a gravity point and 0.5 × 0.5 for the rest of the area. Calculations were also performed for four versions of the density of the known gravity points: 2.4; 4.9; 9.7; and 19.4 . Given that the accuracy of the classical interpolation methods considered here also depends on the density of the topographic masses accepted for reduction [e.g., 9], the calculations were performed for two constant density models: = 2200 ⁄ , = 2670 ⁄ and for the global topographic mass density model UNB_TopoDens. For comparison, a prediction based on free-air anomalies and Kriging interpolation ( Δ /Kriging) was also performed. A description of the analyzed prediction versions is given in Table 4.  Table 4 were adopted as prediction accuracy parameters for each analyzed method. The standard deviations are shown in Tables 5 and 6.  When analyzing the presented results, it should be noted that there is a clear dependence of the prediction accuracy of both: The density of gravity data and the density of topographic masses used in the calculations. The influence of the density of topographic masses is significantly visible for the versions using a lower density of gravity data. However, differences in prediction accuracy are visible for areas of Slovakia (an area with a slightly higher location), even for the highest density of gravity data. The best results for each of the methods were obtained for the version using = 2670 ⁄ . It is surprising to obtain significantly worse prediction results using the UNB_TopoDens model. This suggests that when predicting gravity in a given area, it is advisable to carry out preliminary tests using various density models and to then choose the density model that gives the best results for the final calculations.
By comparing the prediction accuracy of the GGI method and the classical method (using CBA or Pratt anomalies), the results for the highest density of gravity points were obtained at a similar level of accuracy. However, alongside the reduction of the density of gravity data, there is a slight advantage of classical methods at the level of 5 to 10%.
It is also worth emphasizing that there is a huge decrease in the accuracy of the prediction if the topography is not included in the prediction process (FA/ version).

Designation of the CBA and the Distance Between the Geoid and Quasigeoid
With the developed GGI model using density of the gravity points 2.4 , the UNB_TopoDens model as a reference density of topography model, DEM with resolution of 0.1 × 0.1 , and assuming zones of constant density equal 4 × 4 , the complete Bouguer anomalies were determined according to equation (17) while the differences between the geoid undulations and the height anomalies were determined according to equation (18). These quantities will hereinafter be designated as CBA and ( − ) . The values calculated from the GGI model were compared with the classical approach using the UNB_TopoDens model and DEM with resolution of 0.1 × 0.1 for a square with side of 10.500 , in the centre of which there was a gravity point and 0.5 × 0.5 within a radius up to 167 km. The CBA values were calculated according to equation (10). We also use Poincaré and Prey reduction (according to the procedure presented in Section 2.2.4, using full topographic reduction and normal gravity gradient) to calculate the values of ̅ that are required in equation (18). These classically determined quantities are hereinafter marked as CBA and ( − ) . In addition, the following differences were also calculated: CBA = CBA − CBA and ( − ) = ( − ) − ( − ) . The basic statistics of the aforementioned values are presented in Table 7.

A B B A A B
Analyzing the maps shown in Figure 5, it should be pointed out that, in their general course, the determined CBA correspond to CBA values and also correspond to CBA determined in other studies for this area [30][31][32]. Unfortunately, the imperfections of the GGI approach in terms of CBA determination are also very clearly visible. Apart from the differences in the values of the determined CBA, which are shown in detail in Figure 7A, the basic problem of the CBA values concerns its resolution. The CBA values were calculated according to the equation (17), in which this value depends primarily on the potential, defined by the volume , the resolution of which is defined by constant density zones (in this case 4 × 4 ). Therefore, increasing the resolution of the CBA requires a significant increase in the resolution of the volume . In its current form, the CBA calculated by the GGI model may only be use in little detailed considerations.
The map of differences ( − ) is significantly better. Even small structures are clearly visible. The course of the determined values ( − ) corresponds to the values ( − ) , however the differences are also visible.
Referring to the values included in Table  , respectively, and are relatively small and acceptable. However, the extreme values of these differences are already very large. In Figure 6, representing the spatial decomposition of these differences, it can be seen that the highest values occur in the area of the eastern part of the Tatra Mountain and in the higher mountain parts of Slovakia.

Estimation of the Formal Errors of the Complete Bouguer Anomalies and the ( − ) Differences Resulting from Errors of Density Estimation.
We will now consider the differences ( CBA and ( − )) that were presented in the previous Section while taking into account the density distribution of topographic masses from the UNB_TopoDens model ( ) and the standard deviations of these densities ( ) (which are part of the UNB_TopoDens model). These data make it possible to determine, on the basis of the low of covariance propagation, the formal errors of both the CBA ( ) and the ( − ) differences ( ), which result from errors of density estimation; however, these are not the only factor affecting the accuracy of the Bouguer anomalies and the ( − ) differences. The maps of this values are presented in Figure 8 ( and ) and Figure 9 ( and ).
(A) (B)  The errors and may be determined on the basis of topographic reduction (13). First, let us express this reduction in the form of the executable equation. If DEM is expressed in rectangular coordinates, the topographic reduction will take the form: The solutions of Newton's integral (26) are presented by [18,19]. If we have the mean square errors (standard deviations) of densities ( ), we will write on the basis of the law of covariance propagation: Further, we assume that the impact of the density errors on the topographic reduction is equal to the impact of these errors on the CBA, so: Value can be found based on equation (18), where ̅ is determined, as mentioned in Section 2.2.4, by Poincaré and Prey reduction according to the procedure presented: where is free-air reduction from point P down to point ′ (Fig. 2), is given by equation (26), and is topographic redaction calculated for point ′: When topographic reductions and are calculated based on the same DEM, equation (29) can be written in the form: Based on this, with given values , we can calculate the mean square error of the value ̅ ( ) resulting from errors of density estimation: and on the basis of equation (18): their basic statistics are presented in Table 8.  (Fig. 9)). Therefore, it should be considered that in the regions where the differences CBA and ( − ) are greater (Fig. 7), the CBA and ( − ) values are inappropriate (also remember that the CBA values we previously considered low quality).
On the other hand, as was shown previously, the gravity prediction results using the UNB_TopoDens model are clearly not as good as the other analyzed versions. These results do not follow to the principle that with well-chosen densities of topography, the CBA and topographicisostatic anomalies should be smoother, which should result in better prediction. This indicates the difficulty of unambiguously assessing of the demonstrated, largest values of the CBA and ( − ) and thus, the values of CBA and ( − ) determined from the GGI model.

Conclusions
In this paper, the features such as height anomalies, gravity values, CBA, and ( − ) differences determined by the GGI method were analyzed.
Since the influence of various parameters of the GGI model on the determined height anomalies has been studied many times before, the analyses of this functional were now limited and confirmed a small impact on determined height anomalies accuracy of the DEM resolution and the size of the constant density zones. As a result of the analyses, a slightly better fit of the GGI model to the quasigeoid model of Slovakia was also found, which is probably caused by the slightly higher accuracy of this quasigeoid model. The Polish quasigeoid model used was developed without direct use of gravity data (the EGM08 model fitted to GNSS/levelling points).
In terms of gravity values, a dependence of the accuracy of the GGI model on all analyzed parameters was shown: The DEM resolution, the point height, the density of gravity data used to build the model, and also on used reference density of topography model. The best results of gravity prediction were characterized by an error of approximately 1 . These results were compared with classical prediction methods (using CBA and topographic-isostatic anomalies supported by Kriging prediction), getting very similar results to the GGI approach. The sensitivity of the classical prediction methods to the model of topographic mass density distribution used for reduction was also confirmed. However, for the best versions ( = 2670 ⁄ ), the classical methods showed slightly better results than GGI method (5 to 10%), in the versions with poorly selected topographic mass density, slightly more favorable prediction results were obtained by the GGI method.
On the basis of the GGI model, CBA and differences ( − ) were also determined. The strong dependence of resolution of the CBA anomaly model on the size of the constant density zones has been demonstrated, which significantly reduces the quality of such a model. Hence, determining the CBA directly from the GGI model requires a change in the approach of the solution itself (a significant reduction in the size of constant density zones is needed). The crucial importance of the topographic masses density model for both determined values was also indicated. Therefore, for determining these quantities, all available information on topographic mass densities should be used in modelling.