Geodetic Boundary Value Problem Research Articles

Modeling Gravity Gradients from Surface Gravity Anomaly Data

Gravity gradients are useful to characterize near mass anomalies since they are much more sensitive to short wavelength anomalies than gravitational accelerations. Estimating gravity gradients from surface gravity data is based on numerical implementations of solutions to geodetic boundary value problem for determination of disturbing potential. One of methods to solve this problem is least-squares collocation which is basically based on data and a defined covariance function. This study deals with estimating gravity gradient tensor components from along track surface gravity anomaly data. The Least-Squares Collocation solution is based on a stationary local covariance function defined for the disturbing potential which allows upward continuation of the observations to a desired altitude. The modeling method is evaluated in using Earth Gravitational Model 2008(EGM2008) and real airborne gravity gradiometry data collected over Southern Texas, Oklahoma region. The results show that modeled gravity gradients estimated in both on the ground and at a certain altitude have basically good agreement with EGM08 gradients. Modeled gradients including horizontal components in the east-west direction exhibit some discrepancies in comparison to the airborne gradiometry data, which may be attributed to some measurement errors in the gradient data.

Connecting the Brazilian Vertical System to the International Height Reference Frame by estimating the vertical datum parameters

One of the challenges in the unification of the global vertical datum is to determine the vertical offset between the local vertical datum and the global vertical datum. Since 2015, the scientific community has been working on the unification of the global vertical datum, following the publication of the resolution for the definition and realization of an International Height Reference System. This paper aims to estimate the vertical offset in Brazil (Imbituba and Santana Data) by combining a regional gravity field model, using the Geodetic Boundary Value Problem approach applying the Remove-Compute-Restore procedure, and GNSS/leveling data. In the Imbituba vertical datum 1271 stations were used and 66 stations were connected to the Santana tide gauge. The results showed that the geopotential value of Imbituba was 62,636,849.87 ± 0.224 m2s-2, and the vertical offset concerning the reference surface (W0) was 0.358 ± 0.023 m. In Santana, the geopotential value was 62,636,836.88 ± 4.108 m2s-2, and the vertical offset was 1.689 ± 0.420 m. The offset between Imbituba and Santana was computed once the physical connection between Brazil's two local vertical data is a challenge. The estimated value was 1.331 ± 0.421 m. This means the Santana vertical datum is 1.331 m higher than the Imbituba vertical datum.

Uncertainties associated with integral-based solutions to geodetic boundary-value problems

Physical geodesy applies potential theory to study the Earth’s gravitational field in space outside and up to a few km inside the Earth’s mass. Among various tools offered by this theory, boundary-value problems are particularly popular for the transformation or continuation of gravitational field parameters across space. Traditional problems, formulated and solved as early as in the nineteenth century, have been gradually supplemented with new problems, as new observational methods and data are available. In most cases, the emphasis is on formulating a functional relationship involving two functions in 3-D space; the values of one function are searched but unobservable; the values of the other function are observable but with errors. Such mathematical models (observation equations) are referred to as deterministic. Since observed data burdened with observational errors are used for their solutions, the relevant stochastic models must be formulated to provide uncertainties of the estimated parameters against which their quality can be evaluated. This article discusses the boundary-value problems of potential theory formulated for gravitational data currently or in the foreseeable future used by physical geodesy. Their solutions in the form of integral formulas and integral equations are reviewed, practical estimators applicable to numerical solutions of the deterministic models are formulated, and their related stochastic models are introduced. Deterministic and stochastic models represent a complete solution to problems in physical geodesy providing estimates of unknown parameters and their error variances (mean squared errors). On the other hand, analyses of error covariances can reveal problems related to the observed data and/or the design of the mathematical models. Numerical experiments demonstrate the applicability of stochastic models in practice.

Global Geopotential Models assessment in Ecuador based on geoid heights and geopotential values

Abstract Since the 1960s, the analysis of disturbed satellite orbits to infer Earth’s gravity field functionals has been an important element in determining the Earth’s gravitational field. The long wavelengths of the gravitational field are recovered through the analysis of non-Keplerian variations in the orbital path of artificial satellites, from their tracking from ground stations (Satellite Laser Ranging, Doppler Orbitography and Radiopositioning Integrated by Satellite, and Precise Range And Range-Rate Equipment), from satellite-to-satellite tracking, or by microwave interferometry. In addition, differences in gravitational acceleration in three mutually orthogonal dimensions can be determined by employing a differential accelerometer carried on artificial satellites (satellite gravity gradiometry, SGG). Satellite gravimetry provides global information (long wavelengths) of the Earth’s gravitational field, which is the fundamental basis for the implementation of Global Geopotential Models (GGMs). The GGMs are one of the key tools for the representation of the Earth’s gravity field and, therefore, for the establishment of a Global Height System (i.e., International Height Reference System), whose fundamental reference surface is defined in terms of a geopotential value. In this study, the validation of high-resolution GGMs (coefficients up to degree 2190) was performed based on their performance in Ecuador by comparing geoid heights estimated by the GGMs with the corresponding values derived from Global Navigation Satellite System/leveling records. Furthermore, geopotential values from the GGMs are compared with the corresponding value obtained for the Ecuadorian Vertical Datum by solving the fixed geodetic boundary value problem. The obtained results indicated that the precision of the high-resolution GGMs does not reach the established requirements for the geopotential computation in the International Height Reference Frame fundamental stations.

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 [{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) in the disturbing potential, that is located in the areas of the highest values of the deflection of vertical.

Analysis of a Relative Offset between the North American and the Global Vertical Datum in Gravity Potential Space

The accurate estimation of the zero-height geopotential level in a local vertical datum (LVD) is critical for linking traditional height reference systems to a global height system. In this paper, we investigate the theoretical and practical challenges involved in determining the offset between the North American vertical datum (NAVD) and the global vertical datum (GVD). Drawing on the classical theory of the vertical system in physical geodesy, we define the vertical datum offset and derive rigorous formulas for its calculation. We examine various factors that affect the determination of the offset, including the global gravitational models (GGMs), geodetic reference system, tide system, tilt error, and omission error. Using terrestrial gravity data and gravity anomalies from multiple GGMs in conjunction with Global Navigation Satellite System (GNSS) and orthometric heights, we estimate the vertical offset between the NAVD and GVD. Our results indicate that the geopotential difference approach and the geodetic boundary value problem (GBVP) approach yield consistent results. When the normal gravity geopotential of the geodetic reference system is selected as the gravity geopotential of the global height datum, the NAVD is approximately 0.04 m higher than the GVD relative to the GRS80 ellipsoid, and 0.97 cm higher than the GVD relative to the WGS84 ellipsoid. When the Gauss–Listing geopotential value is chosen as the gravity geopotential of the global height datum, the NAVD is roughly 1.45 m higher than the GVD relative to the GRS80 ellipsoid, and approximately 0.52 m higher than the GVD relative to the WGS84 ellipsoid.

Linkage of Normal Heights Obtained from GNSS and Refined XGM2019 GGM to Brazilian Vertical Data Using Different Approaches

This work aimed to analyze the use of different approaches to link normal heights obtained via Global Navigation Satellite System (GNSS)/Global Geopotential Model (GGM) refined by the RTM technique to the Brazilian Vertical Data (Imbituba Brazilian Vertical Datum – IBVD and Santana Brazilian Vertical Datum – SBVD). Specifically, it analyzed approaches based on the weighted mean of discrepancies between height anomalies, the zero- level geopotential value, the Geodetic Boundary Value Problem (GBVP) solution, and the use of parametric modeling of a plane with a scale factor. For the numerical tests, two different study regions have been used, the first with heights referenced to IBVD and the second to SBVD. Using the first three approaches, the local modeling idea has been investigated in both regions. In this context, spatial cluster analysis of the outliers of differences between local and global height anomalies defined the sub-regions. In the fourth approach, the treatment of local modeling was initially considered. In the accuracy analysis of linkages, it has been verified that approaches based on the mean of the discrepancies between height anomalies and using zero- level geopotential value propose practically the same results. On the other hand, there were improvements at the centimeter level with the use of the GBPV solution-based approach compared to the first two, except for two worsening cases. With the approach based on parametric modeling, the accuracy results were mainly worse considering the approaches with local modeling. The most significant differences reached the decimeter level.

A detailed quasigeoid model of the Hong Kong territories computed by applying a finite-element method of solving the oblique derivative boundary-value problem

Abstract New gravity and precise levelling measurements have been performed throughout the Hong Kong territories to modernize a vertical geodetic datum that is currently realized by heights of levelling benchmarks defined in the Hong Kong Principal Datum (HKPD). Modernization of the HKPD involved delivering various products, including new detailed geoid and quasigeoid models and newly determined orthometric and normal heights of levelling benchmarks. In this study, we present the result of gravimetric quasigeoid modelling. The method used to compute a detailed gravimetric quasigeoid model is based on the finite-element method to solve the geodetic boundary-value problem with oblique derivative boundary conditions considered directly at computational nodes on the discretized Earth’s topography. The result of a gravimetric quasigeoid modelling shows a good agreement with a geometric quasigeoid model at the Global Navigation Satellite System (GNSS)-levelling benchmarks. The standard deviation of differences between the gravimetric and geometric quasigeoid heights of ±3.3 cm is compatible with the expected accuracy of gravity, levelling, and GNSS measurements.

On the connection of the Ecuadorian Vertical Datum to the IHRS

Abstract In this work, the determination of the discrepancy between the Ecuadorian Vertical Datum (EVD) and the International Height Reference System (IHRS) is presented. The vertical offset was estimated at the EVD based on the fixed geodetic boundary value problem approach. The focus of the experiment was the determination of the anomalous potential in the EVD, which in turn enable determination of the respective geopotential value. Taking a geopotential space-based approach, two estimates of the EVD offset with respect to the IHRS were obtained that amount to −1.51 and −1.61 m2/s2.

Estimation of Geopotential Value W0 for the Geoid and Local Vertical Datum Parameters

Unification of the global vertical datum has been a key problem to be solved for geodesy over a long period, and the main challenge for unifying the global vertical datum system is to determine the geopotential value W0 of the geoid and to calculate the vertical offset between the local vertical datum and the global vertical datum W0. The geopotential value W0 can be calculated using the grid mean sea surface (GMSS) data and the global geopotential model (GGM). In this study, this GMSS data was measured with adjustment methods and 24 years of merged multi-satellite altimetry data. The data of HaiYang-2A (HY-2A) and Jason-3 were first used to calculate W0. The geopotential value W0 was determined to be 62,636,856.82 m2s−2 by combining the EIGEN-6C4 (European Improved Gravity Model of the Earth by New Techniques) and the GMSS data. Then, the geopotential difference approach and geodetic boundary value problem (GBVP) approach were used to determine the vertical datum parameters in this study. To compensate for the omission error of the GGM, this study utilized the remove–compute–restore (RCR) technique and the residual terrain model (RTM)-recovered high-frequency gravity signals. Finally, as a result of the GBVP solution, the geopotential value of the Australian Height Datum (AHD) was 62,636,851.935 m2s−2, and the vertical offset of the AHD relative to the global vertical datum W0 was 0.4885 m. As a result of the geopotential difference approach, the geopotential value of the Chinese Height datum was 62636861.412 m2s−2, and the vertical offset of the Chinese Height datum was −0.4592 m.

Решение задачи Молоденского по чистым аномалиям силы тяжести с относительной ошибкой порядка квадрата сжатия

The solution of the geodetic boundary value problem for determining the anomalous potential from gravity measurements in the spherical approximation, taking into account the relief and compression, was developed in sufficient detail in 1960 and is also based on the results of G. G. Stokes. Flattening of the reference surface was taken into account by D. V. Zagrebin in several works dated 1940–1970; in 1956 M. S. Molodensky proposed a simpler method for an oblate ellipsoid, based on the already known one for a sphere. It turned out that at the very beginning of developing the mentioned solution, an error was made when obtaining the derivative of the reciprocal distance between two points on the ellipsoid in the direction of the outward normal; later it was repeated in issue 131 of the Proceedings of the TsNIIGAiK. That inaccuracy did not influence the conclusions of V. V. Brovar in 1996 and 1999, where the relief was added to the consideration of compression. This solution of Molodensky’s problem, taking the compression into account, was seldom used, and that is probably why the mistake went unnoticed. V. V. Popadyev saw it and found the correct expression for the mentioned derivative, I. I. Mezhenova adapted subsequent conclusions for the gravity disturbances.

Computation and Analysis of Geopotential Number in São Paulo, Brazil

In recent decades, important steps have been taken to implement the physical concepts of Geodesy in practice, con- cerning height systems. Despite the difficulties involving gravity field modeling, with the establishment of conventions, standards, and computation strategies, the realization of the International Height Reference System (IHRS) is well underway. For a global system, there are constraints for some countries, especially for those with sparse gravity data, mountain regions, and vast areas. In terms of methodology, the computation can be performed directly using the Global Geopotential Models (GGM), recovering existing geoid models, or determining pointwise the gravity potential using integral formulas. In general, the regional gravity modeling is given by numerical integration or least-squares collocation and more recently adopting the spherical radial basis functions. The first approach allows determining the earth’s gravity component at a specific point and adjusting the integral formula according to the gravity coverage. Since so far there is no common sense about the best methodology, computation strategies are been analyzed. In this con- text, the paper aims to contribute to IHRF, computing the geopotential number in the scope of IHRF, using numerical integration to solve the Geodetic Boundary Value Problem and an existing recent quasi-geoid model in four stations in São Paulo state, Brazil. The first approach was performed considering two cases: a radius of 210 km and 110 km of gravimetric data coverage and the Global Geopotential Model GOCO05S truncated at 100 and 200, respectively. The results between solutions have shown a maximum difference of 94 cm, and a minimum difference of 10 cm.

Comparison of kriging and least-squares collocation – Revisited

Abstract The well-known to physical geodesists method of least-squares collocation and the geostatistical method of kriging probably known to the broader audience are compared. Both methods are rooted in Wiener–Kolmogorov’s (W–K) prediction theory; but, since necessity is the mother of invention, the W–K foundations have been extended to satisfy the needs of particular applications. The paper presents a link or rather an equivalence of the two methods as far as their basic forms are considered (specialization to geodetic boundary-value problems, covariance propagation between functionals and nonlinear geostatistical methods are excluded from this comparison). Only scalar random fields (univariate case) and the assumption of a second-order structure of a random function are considered. Due to the equivalence of their basic formulas, both techniques share the same advantages and disadvantages. The paper also shows the difference as to the predicted values and prediction variances in case of exact and filtered (noise reduction) prediction models. This theoretical comparison of the methods has practical implications because of readily available geostatistical software that, in local as well as global applications, can be used for predictive problems occurring in geodesy and surveying.

Determination of the geopotential value on the permanent GNSS stations in Vietnam based on the Geodetic Boundary Value Problem approach

SUMMARY In the realisation of the International Height Reference System, the determination of the geopotential value and its variations in time plays an important role. In this study, the geodetic boundary value problem approach is applied for direct determination of the gravity potential value using a GOCE global gravity field model enhanced with terrestrial gravity data. This determination is carried out on the Global Navigation Satellite System-Continuously Operating Reference Stations (GNSS-CORS) stations towards the realisation of the International Height Reference System in Vietnam. First, the effects of the GOCE global gravity field model omission error, the indirect bias term on the disturbing potential and the systematic cumulative errors in levelling data are estimated. These errors affect the estimated geopotential value. The results calculated on the GNSS/levelling points show that the effect of the GOCE DIR-R5 (up to degree/order 260) omission error on the offset potential value is quite significant. This effect was eliminated using high-resolution terrestrial gravity data using the remove-compute-restore technique. The indirect bias term on the disturbing potential can be safely neglected by using a GOCE global gravity field model for degrees higher than 60 for this study region. The systematic cumulative errors in levelling data can be modelled and removed using a third-order polynomial model. Then, the mean zero-height gravity potential of the Vietnam local vertical datum is estimated equal to ${\rm{W}}_0^{{\rm{LVD}}}$ = 62 636 846.69 m2 s–2 with standard deviation of 0.70 m2 s–2 based on the proposed approach. Finally, the geodetic boundary value problem approach was used to determine the geopotential on the surface of three GNSS-CORS stations in Vietnam. Based on time-series of the vertical component derived from the GNSS observations as well as InSAR data, temporal variations in geopotential are also estimated on these permanent GNSS stations. The purpose is to monitor deformation of the vertical datum. The results indicate that the geopotential value needs to be monitored and determined with the time-dependent component on the three Vietnamese permanent GNSS stations for a vertical datum. These stations may contribute to increase the density of reference points in the International Terrestrial Reference Frame, which is being researched and implemented by the International Association of Geodesy.

Estimation of vertical datum offset for the South African vertical datum, in relation to the international height reference system

The vertical offset and the geopotential value over South Africa is estimated on the four fundamental benchmarks in relation to the international height reference system (IHRS). It is estimated to obtain discrepancies between the South African local vertical datum (W_P) and the global vertical datum (W_0). A single-point-based geodetic boundary value problem (GBVP) approach was used following Molodensky theory for estimating the height anomalies from the disturbing potential (T_P) using Bruns’s formula. The gravity potential at each tide gauge benchmark (TGBM) in South Africa deviates from the potential of the global reference surface by 0.589,-1.993,-2.593 and 2.154 m2s-2 for Cape Town, Port Elizabeth, East London and Durban, respectively. The corresponding vertical datum offsets between the international height reference system and the four fundamental benchmarks over South Africa are 6.013, -20.347, -26.478, and 21.996 cm for Cape Town, Port Elizabeth, East London and Durban, respectively. These offsets can be used for the unification of the South African vertical datum at the four tide gauge benchmarks in a manner that is consistent to the international height reference system.

Estimation of Vertical Datum Parameters Using the GBVP Approach Based on the Combined Global Geopotential Models

Unification of the global vertical datum has been a key problem to be solved for geodesy over a long period, and the main challenge for a unified vertical datum system is to determine the vertical offset between the local vertical datum and the global vertical datum. For this purpose, the geodetic boundary value problem (GBVP) approach based on the remove-compute-restore (RCR) technique is used to determine the vertical datum parameters in this paper. In the RCR technique, a global geopotential model (GGM) is required to remove and restore the long wavelengths of the gravity field. The satellite missions of the GRACE (Gravity Recovery and Climate Experiment) and GOCE (Gravity field and steady-state Ocean Circulation Exploration) offer high accuracy medium–long gravity filed information, but GRACE/GOCE-based GGMs are restricted to medium–long wavelengths because the maximum degree of their spherical harmonic representation is limited, which is known as an omission error. To compensate for the omission error of GRACE/GOCE-based GGM, a weighting method is used to determine the combined GGM by combining the high-resolution EGM2008 model (Earth Gravitational Model 2008) and GRACE/GOCE-based GGM to effectively bridge the spectral gap between satellite and terrestrial data. An additional consideration for the high-frequency gravity signals is induced by the topography, and the residual terrain model (RTM) is used to recover the omission errors effect of the combined GGM. In addition, to facilitate practical implementation of the GBVP approach, the effects of the indirect bias term, the spectral accuracy of the GGM, and the systematic levelling errors and distortions in estimations of the vertical datum parameters are investigated in this study. Finally, as a result of the GBVP solution based on the combined DIR_R6/EGM2008 model, RTM, and residual gravity, the geopotential values of the North American Vertical Datum of 1988 (NAVD88), the Australian Height Datum (AHD), and the Hong Kong Principal Datum (HKPD) are estimated to be equal to 62636861.31 ± 0.96, 62653852.60 ± 0.95 and 62636860.55 ± 0.29 m2s−2, respectively. The vertical offsets of NAVD88, AHD, and HKPD with respect to the global geoid are estimated as −0.809 ± 0.090, 0.082 ± 0.093, and −0.731 ± 0.030 m, respectively.

Effect of great lakes on gravity reduction and geoid determination caused by unclassified DTMs: case study for Lake Victoria, Africa

The determination of the gravimetric geoid is based on the magnitude of gravity observed at the topographic surface of the Earth. In order to satisfy Laplace’s equation, the masses between the surface of the Earth and the geoid must be removed or shifted inside the geoid. Then the gravity values have to be reduced to the geoid, forming the boundary values on the boundary surface. Gravity reduction techniques with unclassified height models usually presume that positive elevations are reserved for positive rock topography. Great lakes, however, are filled with water and may be situated partially or fully above sea level. In case of Lake Victoria, the whole lake including its bed floor is situated above sea level, i. e., having positive elevation (orthometric height). This leads to an obvious error in the topographic-isostatic reduction using, for example, the TC program (Forsberg in A study of terrain reductions, density anomalies and geophysical inversion methods in gravity, 1984; Forsberg and Tscherning in Sanso F, Rummel R (eds) Geodetic boundary value problems in view of the one centimeter geoid, Lecture notes in earth sciences. vol 65, pp 239–272, 1997. https://doi.org/10.1007/BFb0011707 ) by assuming rock topography filling the lake instead of water. The aim of this paper is to determine the effect of Lake Victoria on gravity reduction and geoid computation in Africa, as a prototype of the effect of great lakes on gravity reduction and geoid determination. The results prove that the masses of Lake Victoria have significant effect both on the reduced anomalies and on the computed geoid, which then have to be considered for precise geoid determination with correct density values.

Colorado geoid modeling at the US National Geodetic Survey

The geoid computation method currently used at the US National Geodetic Survey is based on the analytical downward continuation solution of the Molodensky geodetic boundary value problem. The quasigeoid is computed first in 1′ × 1′ grids for the whole area, and then, the geoid is obtained by adding the geoid–quasigeoid separation term based on the simple Bouguer anomalies. A few variations in data combination and topographic/gravity reductions are tested, and three models are selected for the Colorado experiment. Historical GPS/leveling data are used as a guide in the computations. The standard deviation of geoid differences between the historical GPS/leveling data and geoid models is 5.2, 5.0 and 4.9 cm for the models of the initial run, the first and second iteration, respectively. After the models were submitted to the Colorado experiment, they were validated using the GSVS17 (Geoid Slope Validation Survey 2017) GPS/leveling data. The standard deviation of the geoid differences on 223 GSVS17 marks is 3.2, 2.6 and 2.3 cm for the three runs, respectively. In addition to the (quasi)geoid models, the geopotential numbers at 223 GSVS17 marks are also computed and submitted to the Colorado geoid computation experiment and it is evaluated in Sanchez et al. (2020 this issue).

Computational optimization in solving the geodetic boundary value problems

<p style='text-indent:20px;'>The finite volume method (FVM) as a numerical method can be straightforwardly applied for global as well as local gravity field modelling. However, to obtain precise numerical solutions it requires very refined discretization which leads to large-scale parallel computations. To optimize such computations, we present a special class of numerical techniques that are based on a physical decomposition of the computational domain. The domain decomposition (DD) methods like the Additive Schwarz Method are very efficient methods for solving partial differential equations. We briefly present their mathematical formulations, and we test their efficiency in numerical experiments dealing with gravity field modelling. Since there is no need to solve special interface problems between neighbouring subdomains, in our applications we use the overlapping DD methods. Finally, we present the numerical experiment using the FVM approach with 93 312 000 000 unknowns that would not be possible to perform using available computing facilities without aforementioned methods that can efficiently reduce a numerical complexity of the problem.

The Finite Element Method as a Tool to Solve the Oblique Derivative Boundary Value Problem in Geodesy

Abstract In this paper, we propose a novel approach to approximate the solution of the Laplace equation with an oblique derivative boundary condition by the finite element method. We present and analyse diverse testing experiments to study its behaviour and convergence. Finally, the usefulness of this approach is demonstrated by using it to gravity field modelling, namely, to approximate the solution of a geodetic boundary value problem in Himalayas.