## 2 Theoretical Basis  
### The GGI Model  
The GGI method is based on a local model of disturbing potential, which consists of three components (Figure 1) [11; 12]. Figure 1: Components of the GGI disturbing potential model (reproduced from [16]). The first component is the potential \(T_{\Omega}\) of the topographic masses, which lies above the geoid, included in volume \(\Omega\), with density distribution function \(\rho\). The second component is the potential \(T_{\kappa}\), of the disturbance masses occurring between the geoid's surface and the compensation level, included in volume \(\kappa\). The density distribution function of the masses will be marked as \(\delta\). Volumes \(\Omega\) and \(\kappa\) 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 \(T_{E}\) (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 \(T_{E}\) should be of the trend type and can be modelled using harmonic polynomials of a low degree. Finally, we can write:  
\[T=T_{\Omega}+T_{\kappa}+T_{E} \tag{1}\]  
Components \(T_{\Omega}\) and \(T_{\kappa}\) are given by [PERSON]'s integrals:  
\[T_{\Omega}=G\iiint_{\Omega}\frac{\rho}{r}dV_{\Omega} \tag{2}\]  
\[T_{\kappa}=G\iiint_{\kappa}\frac{\delta}{r}dV_{\kappa} \tag{3}\]  
where \(r\) is the distance between the attracting masses and the attracted point \(P\), \(G\) is [PERSON]'s gravitational constant, and \(dV_{\Omega}\) and \(dV_{\kappa}\) are elements of volumes. In our calculations we use component \(T_{E}\) in the form:  
\[T_{E}=a_{1}+a_{2}X_{P}+a_{3}Y_{P}+a_{4}X_{P}Y_{P}+a_{5}H_{P} \tag{4}\]  
where \(X_{P}\), \(Y_{P}\), \(H_{P}\) are coordinates of the point \(P\) and \(a_{1},\dots,a_{5}\) are determined coefficients. We can now formulate an inversion problem: Find the density distribution functions \(\rho\) and \(\delta\) in defined volumes \(\Omega\) and \(\kappa\) and find the coefficients of the polynomials modelling the potential \(T_{E}\) to satisfy Equation (1) for given data. Using linear inversion [17], the solution of the task requires discretisation of the continuous 3D functions \(\rho\) and \(\delta\). Volumes \(\Omega\) and \(\kappa\) are divided into finite volume blocks and a constant density is assigned to each of the blocks. In previous studies, we used volume \(\Omega\), defined by digital elevation model (DEM) in the form of rectangular blocks that are grouped into zones of constant, searched density. The \(\kappa\) 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 \(\Omega\) and \(\kappa\) 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 \(\Omega\) and \(\kappa\) volumes in the form of rectangular prisms for which the solutions of [PERSON]'s integrals are presented by [18; 19].