1 Introduction

The nonlinear boundary value problems (BVPs) arise in a variety of areas of engineering tasks, and in our study we will focus on those that concern the gravity field modelling. The first application of a nonlinear BVP for the Laplace equation to gravity was given by Backus in Backus (1968). Then Grafarend and Niemeier discussed the free nonlinear BVP exactly solved by metric continuation in Grafarend and Niemeier (1971), and developed and extended this approach in the report published by Grafarend in Grafarend (1989). Koch and Pope (1972) presented a uniqueness proof for the nonlinear geodetic BVP. Few years later Bjerhammar and Svensson (1983) used the general implicit function theorem and gave a solution of the existence and uniqueness problem in the nonlinear case. An interesting study discussing an expansion of the nonlinear boundary condition into a Taylor series, based upon some reference potential field approximating the geopotential, was shown by Heck in Heck (1989). In the same year, Sacerdote and Sansó (1989) further developed the idea used by Bjerhammar and Svensson for an iterative solution and they found explicit convergence conditions. They calculated the respective constant governing the convergence in the ideal case of a spherical boundary. Díaz et al. (2006, 2011) showed the existence and uniqueness of a viscosity solution for the Backus problem. Recently, Macák et al. (2016) published a numerical approach for solving the nonlinear satellite-fixed geodetic boundary value problem (NSFGBVP) by the finite volume method (FVM).

In this paper, we solve the same NSFGBVP as defined in Macák et al. (2016), however, for its solution, we implement the finite element method (FEM) and apply it to a very detailed local gravity field modelling in Andean, Himalayan and Alpine mountain ranges using large-scale parallel computations. Moreover, in this approach we use the triangular discretization of the bottom boundary, which is natural for approximating the complicated Earth’s surface, while in Macák et al. (2016) only ellipsoidal approximation has been applied. In this way, we utilize the main advantages of numerical approaches applied to gravity field modelling, such as a straightforward refinement of the discretization, opportunity to consider real complicated topography as well as feasibility for high-resolution modelling. We expect that taking into account the proposed iterative procedure will improve the solution in the areas of the highest values of the deflection of vertical.

The paper is organized as follows. In the Sect. 2, we formulate the NSFGBVP and derive the iterative procedure for our approach. In the Sect. 3, we present its solution by the FEM using the pentahedral elements (i.e. triangular prisms). The Sect. 5 deals with the high-resolution local gravity field modelling and presents a discussion on the obtained results. The paper ends with conclusion and summary.

2 Formulation of the nonlinear satellite-fixed geodetic boundary value problem and the iterative procedure

Let us consider the bounded domain \(\varOmega\) depicted in Fig. 1, similarly as it was defined in Fašková et al. (2007, 2010). The domain \(\varOmega\) is placed in the external space above the Earth and is bounded by the bottom surface \(\varGamma\) representing a chosen part of the Earth’s surface, the corresponding upper part of the boundary at altitude of the GOCE satellite mission and additional four side boundaries.

Fig. 1
figure 1

Brief sketch of the computational domain \(\varOmega\) (highlighted in blue) for the local gravity field modelling. The boundary \(\varGamma\), visualised by green hatching, represents the approximation of the chosen part of the Earth’s surface, B and L denote ellipsoidal latitude and longitude, respectively

Full size image

In the \(\varOmega\) domain, we suppose the nonlinear satellite-fixed geodetic BVP (NSFGBVP) for the disturbing potential \(T({\textbf {x}})\) that was defined in Macák et al. (2016) and is formulated in the following form

$$\begin{aligned} \varDelta T(\textbf{x})= & {} 0 \quad \textbf{x}\in \varOmega , \end{aligned}$$
(1)
$$\begin{aligned} |\nabla (T(\textbf{x})+U(\textbf{x}))|= & {} g(\textbf{x}) \quad \textbf{x}\in \varGamma , \end{aligned}$$
(2)
$$\begin{aligned} T(\textbf{x})= & {} T_{SAT}(\textbf{x}) \quad \textbf{x}\in \partial \varOmega -\varGamma . \end{aligned}$$
(3)

where \(U(\textbf{x})\) denotes the normal gravity potential (Hofmann-Wellenhof and Moritz 2006), \(g(\textbf{x})\) the magnitude of total gravity vector and \(T_{SAT}\) is, in general, the disturbing potential generated from a satellite-only global geopotential model (GGM) based on the spherical harmonics. The Eq. 2) represents the nonlinear boundary condition (BC).

The derivation of the iterative procedure has been published in Macák et al. (2016), but since it is very crucial for better understanding as well as readability of this paper, we present the main ideas also here.

We rewrite the norm of the gradient on the left-hand side of BC (2) in the form

$$\begin{aligned} |\nabla (T(\textbf{x})+U(\textbf{x}))| = \frac{\nabla (T(\textbf{x})+U(\textbf{x}))}{|\nabla (T(\textbf{x})+U(\textbf{x}))|}\cdot \nabla (T(\textbf{x})+U(\textbf{x})) \end{aligned}$$
(4)

and insert (4) into (2) to obtain

$$\begin{aligned} \frac{\nabla (T(\textbf{x})+U(\textbf{x}))}{|\nabla (T(\textbf{x})+U(\textbf{x}))|}\cdot \nabla (T(\textbf{x})+U(\textbf{x})) = g(\textbf{x}). \end{aligned}$$
(5)

If we denote by \(\textbf{v}(\textbf{x})\) the unit vector

$$\begin{aligned} \textbf{v}(\textbf{x})=\frac{\nabla (T(\textbf{x})+U(\textbf{x}))}{|\nabla (T(\textbf{x})+U(\textbf{x}))|}, \end{aligned}$$
(6)

i.e., the vector defining the direction of the gravity vector, the Eq. 5) becomes

$$\begin{aligned} \textbf{v}(\textbf{x})\cdot \nabla (T(\textbf{x})+U(\textbf{x})) = g(\textbf{x}), \end{aligned}$$
(7)

and after some rearrangement we get

$$\begin{aligned} \textbf{v}(\textbf{x})\cdot \nabla (T(\textbf{x})) = g(\textbf{x})-\textbf{v}(\textbf{x})\cdot \nabla (U(\textbf{x})) \quad \textbf{x}\in \varGamma . \end{aligned}$$
(8)

Since the vector \(\textbf{v}(\textbf{x})\) is unknown and depends on \(T(\textbf{x})\), BC (8) is still nonlinear, but its form allows to use an iterative approach for determining \(\textbf{v}(\textbf{x})\) and \(T(\textbf{x})\) such that BVP defined as (1) - (3) is fulfilled.

Then the iterative procedure for solving NSFGBVP is defined as follows

$$\begin{aligned} \varDelta T^{n+1}(\textbf{x})= & {} 0\quad \textbf{x}\in \varOmega , \end{aligned}$$
(9)
$$\begin{aligned} \textbf{v}^n(\textbf{x})\cdot \nabla (T^{n+1}(\textbf{x}))= & {} g(\textbf{x})-\textbf{v}^n(\textbf{x})\cdot \nabla (U(\textbf{x})) \quad \textbf{x}\in \varGamma , \end{aligned}$$
(10)
$$\begin{aligned} T^{n+1}(\textbf{x})= & {} T_{SAT}(\textbf{x}) \quad \textbf{x}\in \partial \varOmega -\varGamma , \end{aligned}$$
(11)

for \(n=0,1,2,\ldots\), where

$$\begin{aligned} \textbf{v}^n(\textbf{x})=\frac{\nabla (T^{n}(\textbf{x})+U(\textbf{x}))}{|\nabla (T^{n}(\textbf{x})+U(\textbf{x}))|}. \end{aligned}$$
(12)

We start the iterations by choosing \(T^{0}(\textbf{x})=0\) and consequently for \(\textbf{v}^0(\textbf{x})\) we get

$$\begin{aligned} \textbf{v}^0(\textbf{x})=\frac{\nabla (U(\textbf{x}))}{|\nabla (U(\textbf{x}))|}=\textbf{s}(\textbf{x}), \end{aligned}$$
(13)

where \(\textbf{s}(\textbf{x})\) represents the direction of the normal gravity vector. In this way, in every iteration we solve the geodetic BVP for \(T^{n+1}(\textbf{x})\) with the prescribed oblique derivative vector \(\textbf{v}^n(\textbf{x})\), while in the first one it is given by

$$\begin{aligned} \textbf{s}(\textbf{x})\cdot \nabla (T^{1}(\textbf{x})) = g(\textbf{x})-\gamma (\textbf{x}) = \delta g(\textbf{x}), \end{aligned}$$
(14)

where \(\gamma (\textbf{x})=|\nabla (U(\textbf{x}))|\) denotes a magnitude of the normal gravity vector and \(\delta g(\textbf{x})\) stands for the gravity disturbance.

In the following iterations, we improve the direction of the unit vector \(\textbf{v}(\textbf{x})\) and so reduce the linearization error. To stop the iteration process, we use as a criterion the difference between two successive iterations and stop the procedure, if in each point we get

$$\begin{aligned} |T^{n}(\textbf{x})- T^{n+1}(\textbf{x})| < \varepsilon , \end{aligned}$$
(15)

where \(\varepsilon\) is a user-specified small real number. Then the last iteration represents the approximation of the disturbing potential \(T(\textbf{x})\) and direction of gravity vector \(\textbf{v}(\textbf{x})\) in (1) - (3), and the sum \(T^{n+1}(\textbf{x})+U(\textbf{x})\) represents the approximation of actual gravity potential \(W^{n+1}(\textbf{x})\) in every point of the computational domain \(\varOmega\).

3 The FEM solution to the nonlinear satellite-fixed geodetic BVP

We can see that in each step of our iterative process (9) - (11) we deal with the oblique derivative BVP defined as

$$\begin{aligned} \varDelta T(\textbf{x})= & {} 0\quad \textbf{x}\in \varOmega , \end{aligned}$$
(16)
$$\begin{aligned} \textbf{v}(\textbf{x})\cdot \nabla (T(\textbf{x}))= & {} g(\textbf{x})-\textbf{v}(\textbf{x})\cdot \nabla (U(\textbf{x})) = \alpha (\textbf{x}), \quad \textbf{x}\in \varGamma , \end{aligned}$$
(17)
$$\begin{aligned} T(\textbf{x})= & {} T_{SAT}(\textbf{x}) \quad \textbf{x}\in \partial \varOmega -\varGamma . \end{aligned}$$
(18)

To solve (16) - (18), we apply the FEM, see (Brenner and Scott 2002) or (Reddy 2006). Since the boundary \(\varGamma\), i.e., the discretized Earth’s surface is complicated and irregular, it is natural to choose triangular prisms, namely the finite pentahedral element \(\varOmega ^e\) with six nodes and five faces. We discretize the whole computational domain \(\varOmega\) into \(n_1\), \(n_2\), \(n_3\) elements in latitudinal, longitudinal and height direction, respectively. Then to specify the position of an element \(\varOmega ^e\) we use indices k, l, m, where \(k=1,...,\, n_1\) , \(l=1,...,\, n_2\) and \(m=1,...,\, n_3\).

Let us consider an arbitrary element \(\varOmega ^e\) from our finite element distretization with indices \(k=1,...,\, n_1\), \(l=1,...,\, n_2\) and \(m=2,...,\, n_3\). We multiply (16) by a weight function w and using Green’s identity (we omit (x) to simplify the notation in the following equations) we obtain the weak formulation of (16) over an arbitrary above-defined element \(\varOmega ^e\)

$$\begin{aligned} \int \limits _{\varOmega ^e} \nabla T \cdot \nabla w\, \textrm{d}x \textrm{d}y \textrm{d}z = \int \limits _{\partial \varOmega ^e}\nabla T \cdot {\textbf {n}}\, w\,\textrm{d}\sigma , \end{aligned}$$
(19)

where \({\textbf {n}}\) denotes the unit normal to \(\partial \varOmega ^e\).

Due to the oblique derivative BC (17) prescribed on the bottom boundary \(\varGamma\), we derive the weak formulation for elements with indices \(k=1,...,\, n_1\), \(l=1,...,\, n_2\) and \(m=1\) separately. We follow the idea of Macák et al. (2020) or Minarechová et al. (2021). We split the oblique vector \({\textbf {v}}\) into one normal, \({\textbf {n}}\), and two tangential, \({\textbf {t}}_1,\,{\textbf {t}}_2\), components, and we insert this decomposition into (17)

$$\begin{aligned} \nabla T\cdot \textbf{v}=c_1\nabla T\cdot {\textbf {n}}+c_2\nabla T\cdot {\textbf {t}}_1+c_3\nabla T\cdot {\textbf {t}}_2=\alpha . \end{aligned}$$
(20)

From (20) we express the normal derivative

$$\begin{aligned} \nabla T \cdot {\textbf {n}}=\frac{\alpha }{c_1}-\frac{c_2}{c_1}\frac{\partial T}{\partial {\textbf {t}}_1}-\frac{c_3}{c_1}\frac{\partial T}{\partial {\textbf {t}}_2}, \end{aligned}$$
(21)

where ∂T/∂t1 and ∂T/∂t2 denote derivatives of the disturbing potential in the direction of t1 and t2, respectively, moreover we assume that \(c_1\ne 0\), and we insert (21) to (19) to get

$$\begin{aligned}{} & {} \int \limits _{\varOmega ^e} \nabla T \cdot \nabla w\, \textrm{d}x \textrm{d}y \textrm{d}z =\nonumber \\{} & {} \quad = \int \limits _{\varGamma ^e} \left( \frac{\alpha }{c_1}-\frac{c_2}{c_1}\frac{\partial T}{\partial {\textbf {t}}_1}-\frac{c_3}{c_1}\frac{\partial T}{\partial {\textbf {t}}_2}\right) w\,\textrm{d}\sigma +\int \limits _{\partial \varOmega ^e\setminus \varGamma ^e}\nabla T \cdot {\textbf {n}}\, w\,\textrm{d}\sigma . \end{aligned}$$
(22)

After some rearrangement, we obtain

$$\begin{aligned}{} & {} \int \limits _{\varOmega ^e} \nabla T \cdot \nabla w\, \textrm{d}x \textrm{d}y \textrm{d}z + \frac{c_2}{c_1}\int \limits _{\varGamma ^e}\frac{\partial T}{\partial {\textbf {t}}_1}\,w\,\textrm{d}\sigma +\frac{c_3}{c_1}\int \limits _{\varGamma ^e} \frac{\partial T}{\partial {\textbf {t}}_2}\,w\,\textrm{d}\sigma =\nonumber \\{} & {} \quad = \int \limits _{\varGamma ^e} \frac{\alpha }{c_1}\,w\,\textrm{d}\sigma +\int \limits _{\partial \varOmega ^e\setminus \varGamma ^e}\nabla T \cdot {\textbf {n}}\, w\,\textrm{d}\sigma , \end{aligned}$$
(23)

where \(\textbf{n}\) is the normal vector and \({\textbf {t}}_1\), \({\textbf {t}}_2\) are tangent vectors to \(\varGamma ^e\subset \partial \varOmega ^e\subset R^3\), where \(\varGamma ^e\) denotes the bottom boundary of an element \(\varOmega ^e\).

Then for a pentahedral element \(\varOmega ^e\) with six nodes we can write

$$\begin{aligned} T\approx T^{e}=\sum \limits _{j=1}^6 T_j^e \psi _j, \end{aligned}$$
(24)

i.e. we approximate the unknown value T as \(T^{e}\) using a linear combination of basis functions \(\psi _j\) with coefficients \(T_j^e,\,j=1,...,\,6\). We substitute (24) into the weak formulation (19) and consider \(\psi _i\) for weight function w. We obtain the \(i^{th}\) equation in the form

$$\begin{aligned}{} & {} \sum \limits _{j=1}^6 T_j^e\,\,\int \limits _{\varOmega ^e}\left( \frac{\partial \psi _j}{\partial x}\frac{\partial \psi _i}{\partial x} +\frac{\partial \psi _j}{\partial y}\frac{\partial \psi _i}{\partial y} +\frac{\partial \psi _j}{\partial z}\frac{\partial \psi _i}{\partial z}\right) \, \textrm{d}x \textrm{d}y \textrm{d}z =\nonumber \\{} & {} \quad =\sum \limits _{j=1}^6 \int \limits _{\partial \varOmega ^e}q_{n}\, \psi _i \,\textrm{d}x\textrm{d}y, \end{aligned}$$
(25)

where ∂ψ/∂x, ∂ψ/∂y, ∂ψ/∂z are directional derivatives,  \(q_n=\nabla T \cdot {\textbf {n}}\) denotes the projection of the vector \(\nabla T\) along the unit normal \({\textbf {n}}\).

For the row of elements \(\varOmega ^e\) given by indices \(k=1,\,...,\, n_1\), \(l=1,...,\, n_2\) and \(m=1\), we follow the same way and after inserting (24) into (23) and considering \(w=\psi _i\), we obtain the \(i^{th}\) equation in the form

$$\begin{aligned}{} & {} \sum \limits _{j=1}^6 T_j^e\left[ \,\,\int \limits _{\varOmega ^e}\left( \frac{\partial \psi _j}{\partial x}\frac{\partial \psi _i}{\partial x} +\frac{\partial \psi _j}{\partial y}\frac{\partial \psi _i}{\partial y} +\frac{\partial \psi _j}{\partial z}\frac{\partial \psi _i}{\partial z}\right) \, \textrm{d}x \textrm{d}y \textrm{d}z \right] +\nonumber \\{} & {} \qquad +\sum \limits _{j=1}^3 T_j^e\left( \,\,\frac{c_{2}}{c_1}\int \limits _{\varGamma ^e} \frac{\partial \psi _j}{\partial {\textbf {t}}_{1}}\psi _i \,\textrm{d}x\textrm{d}y+\frac{c_{3}}{c_1}\int \limits _{\varGamma ^e} \frac{\partial \psi _j}{\partial {\textbf {t}}_{2}} \psi _i \,\textrm{d}x\textrm{d}y \right) =\nonumber \\{} & {} \quad =\sum \limits _{j=1}^3 \int \limits _{\varGamma ^e}\frac{\alpha _j}{c_1} \psi _i \,\textrm{d}x\textrm{d}y +\sum \limits _{j=1}^6 \int \limits _{\partial \varOmega ^e\setminus \varGamma ^e}q_n\, \psi _i \,\textrm{d}x\textrm{d}y, \end{aligned}$$
(26)

where index \(j=1,..., 3\) refers to nodes of the element \(\varOmega ^e\) that lie on the bottom boundary \(\varGamma\) of the computational domain \(\varOmega\).

Now we can write (25) and (26) in a compact matrix form

$$\begin{aligned} \textbf{K}^e\,\textbf{T}^e=\textbf{Q}^e, \end{aligned}$$
(27)

where \({\textbf {K}}^e=[K_{ij}]\) denotes the element stiffness matrix, \(\textbf{T}^e=(T_1,..., T_6)\) is a column vector of unknowns and \(\textbf{Q}^e\) denotes the right-hand side vector.

To evaluate \({\textbf {K}}^e\) and \(\textbf{Q}^e\) we proceed in the following way. We choose one basis function \(\psi _i\) per vertex \(N^e_i\), see (Reddy 2006) or (Brenner and Scott 2002), and we differentiate the basis functions with respect to a position of each node in cartesian coordinates. To evaluate boundary integrals over a boundary \(\varGamma ^e\) in (26), we approximate tangential derivatives like in the finite difference method. This idea has been applied also in Macák et al. (2020); Minarechová et al. (2021).

Afterwards, we assemble all element equations using two principles (Reddy 2006), namely "continuity of primary variables at the interelement nodes" and "balance of secondary variables in a weighted-integral sense at the interface between two elements", and we take into account the Dirichlet BC (18) for nodes that lie on the \(\partial \varOmega -\varGamma\). We have obtained the global linear system of equations with a column vector of unknown global nodal values \({\textbf {T}}\)

$$\begin{aligned} {\textbf {K}}{} {\textbf {T}}={\textbf {Q}}, \end{aligned}$$
(28)

where the matrix \({\textbf {K}}\) is sparse and positive definite, and entries of the column vector \({\textbf {Q}}\) are zeros except for the nodes with the prescribed BC (17).

4 Testing numerical experiment

Fig. 2
figure 2

The testing numerical experiment - the disturbing potential is obtained as a difference between the gravitational potential generated by the sphere S, representing the simplified Earth and depicted by black colour, and the gravitational potential generated by sphere \(S^*\) representing the normal body and depicted by red colour

Full size image
Table 1 Statistics of residuals between the FEM solution and the exact solution on the bottom boundary in the testing experiment (units: \(m^2s^{-2}\))
Full size table

To test and study the behaviour of the proposed numerical scheme, we performed one theoretical experiment. In this simplified testing case, the disturbing field was generated between two identical homogeneous spheres with the radius \(R=6 378\,[km]\) but with centers mutually shifted by \(100\,[km]\) in the direction of the z-axis, see Fig. 2, where the black sphere, denoted by S, represents the simplified "Earth" and the red sphere, denoted by \(S^*\), the "normal body". Then the computational domain was bounded by latitude \(80^{\circ }\hbox {S}\) and \(80^{\circ }\hbox {N}\), longitude \(0^{\circ }\) and \(160^{\circ }\) and the upper boundary has been placed at the altitude of \(240\,[km]\) above S. The generated gravitational potential, used to obtain the Dirichlet BC and the exact solution of (1), was calculated as GM/r, where GM denotes the geocentric gravitational constant and r denotes the distance from the origin O or \(O^*\). To obtain BC on the bottom boundary, we used the derivative of this exact solution in the form \(-GM/r^2\). We started (see Table  1) with the resolution \(81 \times 41\times 11\), i.e. with the number of nodes in longitudinal, latitudinal and radial direction, respectively. Afterwards, we performed four successive refinements and in each refinement we have calculated four iterations. We can see that with the refinement, we obtained the convergence of the method, and moreover it is obvious that the second iteration of our solution is very close to the "reference" value (see Table  1), so the following iterations bring only small improvement of the results. In this case, the "reference" value is used to designate the values obtained with the exactly calculated vector v(x) (see eq. (17)), so this value includes only the discretization error.

5 High-resolution local gravity field modelling

We applied the developed concept to local gravity field modelling in three areas, namely the Andean, Himalayan and Alpine mountain range, where we expected a significant contribution of our approach due to rugged mountain terrain. In all three experiments, the input surface gravity disturbances were interpolated from Global Gravity Model plus database (GGMPlus) (Hirt et al. 2013) on land areas and complemented by DTU21_GRAV data (Andersen and Knudsen 2019) on marine areas. The bottom surface was the discretized Earth’s surface created from a digital elevation model provided within the GGMPlus database. The upper boundary was placed at the altitude of 230 km above the reference ellipsoid that corresponds to the mean altitude of the GOCE satellite orbits. The disturbing potential on the side boundaries and top boundary was generated from EIGEN-6C4 model (Förste et al. 2014), that is a static combined GGM up to d/o 2190. To solve the nonsymmetric linear system, BiCGSTAB(l) linear solver with \(l=8\), cf. (Sleijpen and Fokkema 1993), was implemented and the stopping criterion was set to \(10^{-8}\). To save memory, we stored only nonzero coefficients. All these large-scale parallel computations were performed on 4 nodes of our cluster with 1.0 TB of distributed memory (each node consists of four 8-core CPUs with 256 GB RAM). Thanks to the NUMA (Non-Uniform Memory Access) architecture of each node, we implemented a hybrid parallelization and the NUMA optimization using 32 MPI processes, each with 4 OpenMP threads (together 128 cores).

5.1 The Andean mountain range

The computational domain \(\varOmega\) for the Andean mountain range is bounded by longitude \(278^{\circ }\) and \(300^{\circ }\), and by latitude \(56^{\circ }\hbox {S}\) and \(12^{\circ }\hbox {N}\). The meshing was set to the number of division: \(1100\times 3400\times 300\) (longitude x latitude x height) creating 1 127 094 801 nodes, including 3 744 501 nodes on the Earth’s surface. These nodes formed 2 244 000 000 pentahedral elements. The corresponding horizontal resolution on the Earth’s surface was \(0.02^{\circ }\times 0.02^{\circ }\). A numerical solution of NSFBVP using the developed FEM approach on such a 3D unstructured mesh took about 15 h of the CPU time. Since the residuals between the third and the second iteration were smaller than 0.00001 [\({\text{m}}^{{\text{2}}} {\text{s}}^{{ - 2}}\)], see Fig. 3d) and Table 2, we stopped our computation after the third iteration. The obtained numerical solution with the differences between two successive solutions is depicted in Fig. 3. One can observe that the highest values of differences are in the areas with the highest values of deflection of vertical and they reach up to \(-\)0.0571 [\({\text{m}}^{{\text{2}}} {\text{s}}^{{ - 2}}\)] in the disturbing potential which corresponds approximately to 5.7 [mm] in the quasigeoidal height, see Table 2.

Table 2 Andean mountain range: Statistics of residuals [\({\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}\)] between the second iteration \(T^{2nd}\) and the first iteration denoted by \(T^{1st}\), and the third \(T^{3rd}\) and the second \(T^{2nd}\) iteration obtained on the bottom boundary \(\varGamma\). STD denotes the standard deviation of residuals
Full size table
Fig. 3
figure 3

The Andean mountain range: a The input surface gravity disturbances [mGal]; b The disturbing potential solution \(T^{1st}\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\); c The contribution of the second iteration, i.e. \((T^{2nd}-T^{1st})\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), where \(T^{2nd}\) denotes the disturbing potential solution after the second iteration; d The contribution of the third iteration, i.e. \((T^{3rd}-T^{2nd})\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), where \(T^{3rd}\) denotes the disturbing potential solution after the third iteration

Full size image

5.2 The Himalayan mountain range

The computational domain \(\varOmega\) for the experiment in Himalayan region is between \(60^{\circ }\hbox {E}\) and \(110^{\circ }\hbox {E}\) longitude, and bounded by latitude \(20^{\circ }\hbox {N}\) and \(50^{\circ }\hbox {N}\). The meshing was set to: \(2000\times 1200\times 320\) (longitude x latitude x height) creating 771 427 521 nodes, including 2 403 201 nodes on the Earth’s surface. In this way, we obtained 1 536 000 000 pentahedral elements. The horizontal resolution on the Earth’s surface has been \(0.025^{\circ }\times 0.025^{\circ }\). The computation (i.e. all three solutions) took about 8,5 h of the CPU time. The obtained solution with the differences between two successive solutions is depicted in Fig. 4. The corresponding statistics is presented in Table 3. One can observe that the highest values of differences are again in the areas of the steepest changes of the disturbing potential, with values up to \(-\)0.0702 [\({\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}\)], see Table 3, which corresponds approximately to 7.0 [mm] in the quasigeoidal height. (Table 4)

Fig. 4
figure 4

The Himalayan mountain range: a The input surface gravity disturbances [mGal]; b The disturbing potential solution \(T^{1st}\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\); c The contribution of the second iteration, i.e. \((T^{2nd}-T^{1st})\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), where \(T^{2nd}\) denotes the disturbing potential solution after the second iteration; d The contribution of the third iteration, i.e. \((T^{3rd}-T^{2nd})\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), where \(T^{3rd}\) denotes the disturbing potential solution after the third iteration

Full size image
Table 3 Himalayan mountain range: Statistics of residuals [\({\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}\)] between the second iteration \(T^{2nd}\) and the first iteration denoted by \(T^{1st}\), and the third \(T^{3rd}\) and the second \(T^{2nd}\) iteration obtained on the bottom boundary \(\varGamma\)
Full size table

5.3 The Alps mountain range

The computational domain \(\varOmega\) for the last experiment located in Alps region is bounded by longitude \(4.799^{\circ }\hbox {E}\) and \(17.503^{\circ }\hbox {E}\), and by latitude \(42.999^{\circ }\hbox {N}\) and \(49.203^{\circ }\hbox {N}\), see Fig. 5. The horizontal resolution on the Earth’s surface was \(0.004^{\circ }\times 0.004^{\circ }\), while the number of divisions in radial direction was set to 600. The total number of nodes in the computational domain was 2 963 353 104, including 4 930 704 nodes on the earth’s surface, and forming 5 911 171 200 pentahedral elements. Again, based on the differences obtained between the third and the second iteration, see Table 3, we stopped the computation after the third iteration. The highest values of differences in the disturbing potential have been approximately \(-\)0.0066 \([ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), i.e. approximately to 0.6 [mm] in the quasigeoidal height, see Table 3 and Fig. 5.

Fig. 5
figure 5

The Alps mountain range: a The input surface gravity disturbances [mGal]; b The disturbing potential solution \(T^{1st}\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\); c The contribution of the second iteration, i.e. \((T^{2nd}-T^{1st})\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), where \(T^{2nd}\) denotes the disturbing potential solution after the second iteration; d The contribution of the third iteration, i.e. \((T^{3rd}-T^{2nd})\,[ {\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}} ]\), where \(T^{3rd}\) denotes the disturbing potential solution after the third iteration

Full size image
Table 4 The Alps: Statistics of residuals [\({\text{m}}^{{\text{2}}} {\text{s}}^{{ - 2}}\)] between the second iteration \(T^{2nd}\) and the first iteration denoted by \(T^{1st}\), and the third \(T^{3rd}\) and the second \(T^{2nd}\) iteration obtained on the bottom boundary \(\varGamma\)
Full size table

6 Conclusion and summary

We developed an iterative approach based on the finite element method applied for solving the nonlinear satellite-fixed geodetic boundary value problem. In this concept, we determined the direction of actual gravity vector together with the value of the disturbing potential. To solve it, we implemented the finite element method with triangular prisms. For the numerical experiments we chose three mountainous locations, namely in Andean, Himalayan and Alpine mountain ranges. As an input data we used the surface gravity disturbances generated from Global Gravity Model plus database (lands) combined with DTU21_GRAV model (seas), and the disturbing potential generated from EIGEN-6C4 model. The obtained maximal contribution to disturbing potential when taking into account the nonlinear boundary condition with comparison to the oblique derivative boundary condition was 0.0571 [\({\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}\)] (Andes), 0.0702 [\({\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}\)] (Himalayas), 0.0066 [\({\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}\)] (Alps). These values correspond to 5.7 [mm] (Andes), 7.0 [mm] (Himalayas), 0.6 [mm] (Alps) in the quasigeoidal heights, and they are located in the areas of the highest values of the deflection of vertical. The obtained results confirm the validity of the presented approach.