Gravity field modelling in mountainous areas by solving the nonlinear satellite-fixed geodetic boundary value problem with the finite element method

The paper concerns a novel concept based on the iterative approach applied for solving the nonlinear satellite-fixed geodetic boundary value problem (NSFGBVP) using the finite element method (FEM). We formulate the NSFGBVP that consists of the Laplace equation defined in the 3D bounded domain outside the Earth, the nonlinear boundary condition (BC) prescribed on the disretized Earth’s surface, and the Dirichlet BC given on a spherical boundary placed approximately at the altitude of GOCE satellite mission and additional four side boundaries. The iterative approach is based on determining the direction of actual gravity vector together with the value of the disturbing potential. Such a concept leads to the first iteration where the oblique derivative boundary value problem is solved, and the last iteration represents the approximation of the actual disturbing potential and the direction of gravity vector. As a numerical method for our approach, we implement the FEM with triangular prisms. High-resolution numerical experiments deal with the local gravity field modelling in the Andean, Himalayan and Alpine mountain range, where we focus on the contribution to the disturbing potential solution by solving the nonlinear geodetic boundary value problem in comparison with the solution to the oblique derivative geodetic boundary value problem. The obtained results showed the maximal contribution of the presented approach 0.0571 [m2s- 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}$$\end{document}] (Andes), 0.0702 [m2s- 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}$$\end{document}] (Himalayas), 0.0066 [m2s- 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text{m}}^{{\text{2}}} {\text{s}}^{{{\text{ - 2}}}}$$\end{document}] (Alps) in the disturbing potential, that is located in the areas of the highest values of the deflection of vertical.


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 Extended author information available on the last page of the article 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. (2006Díaz et al. ( , 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 largescale 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.

Formulation of the nonlinear satellite-fixed geodetic boundary value problem and the iterative procedure
Let us consider the bounded domain depicted in Fig. 1, similarly as it was defined in Fašková et al. (2007Fašková et al. ( , 2010. The domain is placed in the external space above the Earth and is bounded by the bottom surface 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. In the domain, we suppose the nonlinear satellite-fixed geodetic BVP (NSFGBVP) for the disturbing potential T(x) that was defined in Macák et al. (2016) and is formulated in the following form where U(x) denotes the normal gravity potential (Hofmann-Wellenhof and Moritz 2006), g(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 and insert (4) into (2) to obtain If we denote by v(x) the unit vector i.e., the vector defining the direction of the gravity vector, the Eq. 5) becomes and after some rearrangement we get Since the vector v(x) is unknown and depends on T(x) , BC (8) is still nonlinear, but its form allows to use an iterative approach for determining v(x) and T(x) such that BVP defined as (1) -(3) is fulfilled. Then the iterative procedure for solving NSFGBVP is defined as follows , 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 for n = 0, 1, 2, … , where We start the iterations by choosing T 0 (x) = 0 and consequently for v 0 (x) we get where s(x) represents the direction of the normal gravity vector. In this way, in every iteration we solve the geodetic BVP for T n+1 (x) with the prescribed oblique derivative vector v n (x) , while in the first one it is given by where (x) = |∇(U(x))| denotes a magnitude of the normal gravity vector and g(x) stands for the gravity disturbance.
In the following iterations, we improve the direction of the unit vector v(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 where is a user-specified small real number. Then the last iteration represents the approximation of the disturbing potential T(x) and direction of gravity vector v(x) in (1) -(3), and the sum T n+1 (x) + U(x) represents the approximation of actual gravity potential W n+1 (x) in every point of the computational domain .

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 To solve (16) -(18), we apply the FEM, see (Brenner and Scott 2002) or (Reddy 2006). Since the boundary , i.e., the discretized Earth's surface is complicated and irregular, it is natural to choose triangular prisms, namely the finite pentahedral element e with six nodes and five faces. We discretize the whole computational domain into n 1 , n 2 , n 3 elements in latitudinal, longitudinal and height direction, respectively. Then to specify the position of an element 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 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 e where n denotes the unit normal to e . Due to the oblique derivative BC (17) prescribed on the bottom boundary , 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 v into one normal, n , and two tangential, t 1 , t 2 , components, and we insert this decomposition into (17) From (20) we express the normal derivative where ∂T/∂t 1 and ∂T/∂t 2 denote derivatives of the disturbing potential in the direction of t 1 and t 2 , respectively, moreover we assume that c 1 ≠ 0 , and we insert (21) to (19) to get After some rearrangement, we obtain where n is the normal vector and t 1 , t 2 are tangent vectors to e ⊂ e ⊂ R 3 , where e denotes the bottom boundary of an element e .
Then for a pentahedral element e with six nodes we can write i.e. we approximate the unknown value T as T e using a linear combination of basis functions j with coefficients T e j , j = 1, ..., 6 . We substitute (24) into the weak formulation (19) and consider i for weight function w. We obtain the i th equation in the form where ∂ψ/∂x, ∂ψ/∂y, ∂ψ/∂z are directional derivatives, q n = ∇T ⋅ n denotes the projection of the vector ∇T along the unit normal n.
For the row of elements 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 = i , we obtain the i th equation in the form where index j = 1, ..., 3 refers to nodes of the element e that lie on the bottom boundary of the computational domain . Now we can write (25) and (26) in a compact matrix form where K e = [K ij ] denotes the element stiffness matrix, T e = (T 1 , ..., T 6 ) is a column vector of unknowns and Q e denotes the right-hand side vector.
To evaluate K e and Q e we proceed in the following way. We choose one basis function 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 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 − . We have Fig. 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 obtained the global linear system of equations with a column vector of unknown global nodal values T where the matrix K is sparse and positive definite, and entries of the column vector Q are zeros except for the nodes with the prescribed BC (17).

Testing numerical experiment
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 = 6378 [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 • S and 80 • N , longitude 0 • and 160 • 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 × 41 × 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 (28) KT = Q, the values obtained with the exactly calculated vector v(x) (see eq. (17)), so this value includes only the discretization error.

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  (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 , 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).

The Andean mountain range
The computational domain for the Andean mountain range is bounded by longitude 278 • and 300 • , and by latitude 56 • S and 12 • N . The meshing was set to the number of division: 1100 × 3400 × 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 • × 0.02 • . 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 [ m 2 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 [ m 2 s −2 ] in the disturbing potential which corresponds approximately to 5.7 [mm] in the quasigeoidal height, see Table 2.

The Alps mountain range
The computational domain for the last experiment located in Alps region is bounded by longitude 4.799 • E and 17.503  Table 3, we stopped the computation after the third iteration. The highest values of differences in the disturbing potential have been approximately − 0.0066 [m 2 s -2 ] , i.e. approximately to 0.6 [mm] in the quasigeoidal height, see Table 3 and Fig. 5.

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 Table 2 Andean mountain range: Statistics of residuals [ m 2 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 . STD denotes the standard deviation of residuals