On accurate geoid modeling: derivation of dirichlet problems that govern geoidal undulations and geoid modeling by means of the finite difference method and a hybrid method

The geoid is the reference surface used to measure heights (orthometric). These are used to study any mass variability in the Earth system. As the Earth is represented by an oblate spheroid (Ellipsoid), the geoid is determined by geoidal undulations (N) which are the separation between these surfaces. N is determined from gravity data by Stokes's Integral. However, this approach takes a Spherical rather than an Ellipsoidal Earth. Here it is derived a Partial Differential Equation (PDE) that governs N over the Earth by means of a Dirichlet problem and show a method to solve it which precludes the need for a Spherical Earth. Moreover, Stokes's Integral solves a boundary value problem defined over the whole Earth. It was found that the Dirichlet problem derived here is defined only over the region where a geoid model is to be computed, which is advantageous for local geoid modeling. Moreover, the method eliminates several of the sources of uncertainty in


INTRODUCTION
Present day needs for highly accurate geoid models have driven many attempts to modify Stokes's Integral and to compute ellipsoidal corrections for it (MARTINEC and GRAFAREND, 1997;ARDESTANI and MARTINEC, 2003;HIPKIN, 2004;NAJAFI-ALAMDARI et al., 2006) which can be as big as 1m in some places (MARTINEC and GRAFAREND, 1997;NAJAFI-ALAMDARI et al., 2006).Current high resolution geoid models do not consider these corrections and show decimeter big uncertainties, e.g.JGEOID2008 with 17cm for Japan (KUROISHI, 2009;ODERA et al., 2012) and USGG2009 with as much as 32.1cm for US (WANG et al., 2011).Users of geoid models require at least a decimeter quality as for oceanographers and geophysicists (KUROISHI, 2009) or a 1cm quality as for geodesists (SANSO and SIDERIS, 2013;TAPLEY et al., 2004) which is hardly attainable without methods of improved accuracy.However, the ellipsoidal corrections do not eliminate the error from the spherical approximation but attenuate it.Moreover, they still have the following errors: limiting the area of integration to a small spherical cap (truncation error), approximating the normal derivative by a radial or other derivative, Earth's topography that cause to exist masses outside the geoid, discretization of the input gravity and propagation of the uncertainties of the input gravity.
The Boundary Value Problem (BVP) from which Stokes's Integral is derived is defined over the whole space on and outside the geoid and over the whole geoid.Consequently, its solution requires gravity data over the whole Earth to provide N at a single point.The Dirichlet problem derived here is defined only over the region of computation.Thus, besides computing N directly on the ellipsoid, its solution requires gravity data only over the region of interest.In solving a Dirichlet problem in the Ellipsoid these 3 major sources of error will be eliminated: the spherical approximation, the truncation of the integral and the Earth's topography.
A Dirichlet problem whose PDE governs N over the geoid is derived here and it is shown how to numerically solve it by using the Finite Difference Method (FDM).However this method can handle considerably large matrices which increases its computational cost.Therefore, it is also proposed a modification on FDM (FDM with subgrids, FDM2) that allows the computation of geoid models of large regions by keeping the matrices involved at a small and constant size.It will be shown here that FDM has large errors due to discretization so it is also proposed a Hybrid method to alleviate those uncertainties.To derive it it will be needed a Spherical Dirichlet Problem which will also be derived here.The Hybrid method as opposed to FDM yields subcentimetric differences from Stokes's Integral.Moreover, an error analysis will be developed which indicates that the Hybrid method may well outperform Stokes's Integral.

DERIVATIONS OF STOKES'S INTEGRAL AND DIRICHLET PROBLEM
Stokes's Integral is derived from the fundamental Equation of Physical Geodesy (HOFMANN-WELLENHOF et al., 2006): where T is the disturbing potential, ∂T/∂h is its vertical gradient, γ is the normal gravity, ∂γ/∂h is its vertical gradient and ∆g is the gravity anomaly (reduced to remove the effect of masses outside the geoid).Eq. ( 1) is used as a boundary condition for the PDE given by Laplace's equation.This constitutes the third BVP of potential theory.The BVP is solved for T and N is derived from it by Brun's formula (HOFMANN-WELLENHOF et al., 2006): Stokes's Integral makes a spherical approximation on Eq. (1) to then compute an integral solution for the BVP considering the whole space outside the geoid.To derive the Dirichlet problem for FDM it is considered the generalized Poisson's Equation in ellipsoidal-harmonic coordinates (u, υ, λ) for the gravitational and the normal potentials to derive the Laplacian of T, which inside the Earth is as follows (SANSO and SIDERIS, 2013) (υ is the complement of the reduced latitude β, i.e. β = π/2 − υ , and λ is the geocentric longitude): where u is lesser than or equals b, E=(a 2 −b 2 ) 1/2 is the linear eccentricity, with a being the semi-major axis and b the semi-minor axis of the Ellipsoid of revolution, G is the gravitational constant and ρ is the Earth's density.

The Dirichlet Problem PDE
The PDE that describes T over the Earth, and thus, solely as a function of T and its partial derivatives in υ and λ can be derived from Eqs. ( 1) and (3) after expressing Eq. ( 1) in ellipsoidal-harmonic coordinates u, υ and λ, for details refer to Appendix A. This PDE is as follows: where BN(υ), CN(υ), FN(υ), GN(υ, λ) are as follows: with u=b, F(υ)=−(1/γ)•(∂γ/∂h) (F(υ) is given explicitly by Eq. (A.2)).Furthermore, for the Earth 4πGρ is very small and can be neglected.
This PDE together with its boundary condition, which is to know T or N (geoidal undulation) at the boundary of the surface where it's been applied, defines an Elliptic BVP with a Dirichlet boundary condition.This, in turn, defines a Dirichlet problem, which can be solved numerically by e.g.FDM.

NUMERICAL SOLUTION OF THE DIRICHLET PROBLEM BY FDM AND FDM2
To obtain T by FDM the partial derivatives of T must be expressed by their equivalent form in finite differences (DIEGUEZ, 2005).In other words, to make the following substitutions-considering data on a grid (j, i), 1<j<J−1 and 1< i< I−1: where ∆λ and ∆υ are the grid spacing.
T can be determined from Eq. ( 12) if gravity anomalies ∆g are provided on a grid and T or N are provided on the boundary of this grid, i.e. at j=1, j=J, i=1, i=I.The boundary values of T or N can be determined by: spirit levelling (HOFMANN-WELLENHOF et al., 2006;SANSO and SIDERIS, 2013), a global geopotential model (PAVLIS et al., 2012;PAIL et al., 2011PAIL et al., , 2010)), satellite altimetry if the region is over the seas (HOFMANN-WELLENHOF et al., 2006;SANSO and SIDERIS, 2013) or the remove-compute-restore technique (HOFMANN-WELLENHOF et al., 2006;SANSO and SIDERIS, 2013;ODERA et al., 2012) which can be based on ellipsoidal corrections (MARTINEC and GRAFAREND, 1997;ARDESTANI and MARTINEC, 2003;HIPKIN, 2004;NAJAFI-ALAMDARI et al., 2006).This way, T is determined from the linear system derived from Eq. ( 12) by solving it iteratively by e.g.Jacobi's method.However, for small grids one can get very large systems, e.g. a 1min grid of a 4 o by 4 o region delivers a linear system on more than 57 thousand variables which in double precision requires more than 24.3 Gigabytes of Random Acess Memory (RAM).This can be big for a personal computer to handle.Therefore it is proposed a Finite Difference Method with subgrids (FDM2) to compute geoid models of regions as big as desired while keeping the linear system as small as desired.FDM2 is best suited for high resolution geoid modeling.For further details on FDM2 refer to Fig. 1.

Short Remark on Remove-Compute-Restore Technique and FDM
In FDM, terrain reductions are not mandatory because the Dirichlet problem is derived from Poisson's equation as opposed to Stokes's Integral.Moreover for FDM and FDM2 the long wavelengths need not be removed because the Dirichlet problem is defined only over the region of computation as opposed to the BVP from which Stokes's Integral is derived which is defined over the whole Earth.Therefore, FDM does not require implementation of a remove-compute-restore technique.
Figure 1 -How FDM2 works.The grid (big square) is divided in subgrids (small squares, as many as desired, here 16 for the sake of simplicity).FDM is applied individually to each subgrid instead of a single run on the grid.In addition, T or N must be known at the boundary of each subgrid.
Figure 2 -Region used for assessment of geoid modeling methods.The region where the geoid models were computed is represented by the solid black square.
The island in the center corresponds to Japan.

COMPARISON OF STOKES'S INTEGRAL WITH FDM AND FDM2
Using FDM and FDM2 it was computed geoid models on a 1degree and a 1min grid respectively using ship-borne gravity data from Japan Oceanographic The differences between FDM and Stokes's Integral are displayed in Table 1 for the 1deg and the 1min grids.It can be observed that an increase in resolution improves the quality of the grid as given by the standard deviation (sd).In both grids, the differences were very large, much greater than the decimeter.However, it is known that the sd of Stokes's Integral over Japanese lands is subdecimetric as measured by Kuroishi (2009) and Odera et al. (2012).Thus, this large deviation is mainly due to the errors of discretization of input and propagation of input uncertainties in FDM, see Eq. (A.5), Appendix B.
Importantly, for the 1min grid the computation of such grid for N was not possible by FDM due to the large number of variables but only by FDM2 (FDM2 was designed for the case when FDM cannot be applied).FDM2 was applied with squared subgrids sizing 1deg.In total, 9 subgrids were required.The 1deg grid has too few points to yield a reasonable implementation of FDM2.Therefore, the 1deg grid was computed only by FDM and the 1min grid was computed only by FDM2.
The map of the grid concerning the input gravity anomalies is displayed on Fig. 3.The geoid model computed in a 1min grid using FDM2 is displayed on Fig. 4. The map of the differences between FDM2 and Stokes's Integral in the 1min grid corresponds to Fig. 5.In the difference map it can be observed significant deviations between the methods.The deviation is null at the boundary of each subgrid and generally increases towards their center.Moreover, a plot of the histogram of the differences is displayed on Fig. 6.As given by a -1.2 kurtosis, it indicates that the differences conform to a uniform distribution.
To alleviate the large uncertainties in FDM it will be proposed here a Hybrid Method.Prior to its formulation it is necessary to derive a Spherical Dirichlet problem.

DERIVATION OF A SPHERICAL DIRICHLET PROBLEM
Considering Poisson's equation in spherical coordinates and the fundamental equation of physical geodesy in spherical approximation it can be derived a Spherical Dirichlet problem that will be used to alleviate the uncertainties in the Ellipsoidal Dirichlet problem.The equations read as follows: The Spherical Dirichlet problem can be derived in an analogous way to that of the Ellipsoidal Dirichlet problem.Its Partial Differential Equation (PDE) reads: error due to the topographical masses outside the geoid and ε T is the error due to truncation of Stokes's Integral.ε Topo and ε T are alleviated because of the implementation of RCR technique.ε E and ε 1 are considered equal or very similar for both N RCR and N FDMS because the same approximation is done in both cases.

ALLEVIATION OF THE UNCERTAINTIES IN FDM BY AN OPTIMAL COMBINATION OF TECHNIQUES
The uncertainty due to discretization of input and propagation of input uncertainties can be alleviated if N is computed by the following formula: Note that Eq. ( 22) gives N as a combination from N computed by several techniques.This optimal combination shifts the error due to discretization to that of Stokes's Integral, which is much smaller.Moreover, ε T , ε Topo and ε′ 1 are small, with ε T and ε Topo being alleviated by RCR's implementation and the same being possible for ε′ 1 in future.So, the error of this new method is mainly due to discretization of input and propagation of the uncertainties of the input in Stokes's Integral, which are much smaller compared to FDM.Hereafter, I will refer to the N given by this method as N Hybrid , i.e.   2. Kurtosis is -1.20 which indicates that the differences again conform to a uniform distribution.Moreover, the difference may be expressed as follows: so that the difference is a measure of the uncertainty due to the Earth's ellipticity that is usually not accounted for in RCR.From Eqs. ( 9) and ( 10) N Hybrid 's uncertainty is expected to be inferior to that of NRCR as ε 1 and ε′ 1 may be considered as approximately equal.N Hybrid has the advantage over Stokes's Integral of eliminating ε E and possibly eliminating ε′ 1 in future.The maximum observed difference is very small in magnitude and σ is below centimeter which is remarkable.Moreover, the differences were again null at the boundary of each subgrid and tended to increase towards their center.The same trend was observed here for the Hybrid method.The geoid computed by the Hybrid method is displayed on Fig. 7 and represents the geoid over the region much more accurately than the sole FDM.The difference map is displayed on Fig. 8.The histogram of the differences is displayed on Fig. 9.

CONCLUSIONS AND FUTURE WORK
In this article the formulations of a Dirichlet problem to compute geoidal undulations and its numerical solution by FDM were proposed.The main advantages are the elimination of the errors due to the neglection of Earth's flattening, truncation of Stokes's Integral and the Earth's topography.The major drawback of this approach is the error due to the limited density of the input gravity.Moreover, an optimal combination of techniques was proposed here and named Hybrid method.
Boundary points in FDM/FDM2 can always be determined by Stokes's Integral in a standard remove-compute-restore procedure or a global geopotential model.The remaining points are the vast majority and can be determined by FDM or FDM2.Moreover, for large regions the computation of a geoid model by FDM may become difficult thus highlighting the usefulness of FDM2.
However, FDM and FDM2 are not as yet alternatives for highly accurate geoid determinations in this time when a 1cm geoid is desired.This is mainly due to the error due to discretization of input and the propagation of uncertainties in the input which were shown to be much larger in FDM and FDM2 than in Stokes's Integral.Therefore, future work should concern with the treatment of such uncertainties in order to make FDM and FDM2 comparable or even better than Stokes's Integral just in the way that was done here by means of the Hybrid method.Afer all, an increase in resolution decreased the uncertainty of FDM, so it is expected that in future, with more gravity data of higher quality and denser coverage available, the overall uncertainty of FDM will decrease.Moreover, FDM and FDM2 eliminate 3 major sources of uncertainty in Stokes's Integral.
The error analysis was derived here and it was found that the uncertainty in the Hybrid method is mainly due to the discretization of input and propagation of input uncertainties in Stokes's Integral because other uncertainties cancel each other or are comparably small.This led to very small uncertainties and a remarkable agreement with Stokes's Integral as opposed to FDM.
The differences between the Hybrid method and Stokes's Integral are due to the Earth's ellipticity so that it is expected that the use of ellipsoidal corrections in future work may further improve N RCR and consequentely N Hybrid as well.Moreover, these differences show that for a centimetric geoid a Spherical Earth should not be used.
Further suggestions for future work: if the required data is available, it would be good to compute a geoid model in land in order to compare results with geoidal undulations determined by spirit levelling.Also good would be to compare FDM/FDM2 with Stokes's Integral on ellipsoidal correction in order to improve accuracy of boundary points and slightly improve overall accuracy of geoid model.To derive a PDE for the Dirichlet problem based on a more accurate approximation of the vertical gradient of T, e.g. using the normal to the Ellipsoid, is the final suggestion for future work.
Furthermore, once the error due to the Earth's ellipticity is similar for N RCR and N FDMS techniques, from the error analysis developed here it is expected that the Hybrid method will outperform Stokes's Integral which would be a step towards a centimetric geoid.

Figure 3 -
Figure 3 -Input gravity anomalies used.The input data spans a larger area than the one to be modeled because boundary points were computed by Stokes's Integral which requires a spherical cap, here 3 arc-deg.Contour interval 1.0mgal.

Figure 6 -
Figure 6 -Histogram of differences between FDM and Stokes's Integral (cm).The frequency is represented by means of probability, i.e. their sum equals 1.

Figure 7 -
Figure 7 -Geoid model computed by the Hybrid method in a 1min grid.

Figure 9 -
Figure 9 -Histogram of differences between the Hybrid method and Stokes's Integral (cm).The frequency is represented by means of probability, i.e. their sum equals 1.