Computational accuracy required in the evaluation of global geopotential models

Global geopotential models (GGMs) are mainly used in the remove-compute-restore (RCR) technique applied to gravity field modeling such as geoid determination and height datum unification. The increase in the number and quality of gravity data has led the developers of GGMs to produce models with higher resolution and accuracy. Basically, the long-wavelength coefficients of the gravity field are computed based on satellite data, while the medium- and short-wavelength coefficients are calculated based on terrestrial (land and sea) data. One of the main challenges regarding the evaluation of high-degree GGMs is to compute the associated Legendre functions of the first kind based on the usual recursive formulas. Since most computational softwares use the double-precision format by default, an important question is whether this level of precision is sufficient to numerically evaluate the associated Legendre functions of the first kind? To answer this question, the computation of the associated Legendre functions of the first kind in different degrees and latitudes is studied based on MATLAB software, which uses the double-precision format by default. From the numerical results, we find that the calculation of associated Legendre functions of the first kind up to degree of 2190 (the highest degree of existing GGMs), does not have sufficient accuracy at latitudes between 56^°20^׳ and 78^°33^׳, where the most critical state occurs at the latitude 60^°. We also find that the accuracy of the calculation of associated Legendre functions at the latitude 60^° (the most critical state) significantly decreases for the degrees higher than 2029. These results imply that the usual computational softwares based on the double-precision format are not suitable for calculating the associated Legendre functions in all degrees and latitudes. This is due to the fact that if we consider the associated Legendre functions of the first kind in the form of a matrix with the dimensions corresponding to the degree and order of the functions, as the degree increases, the numbers on the main diagonal approach to the number 10^-308 and thus they are considered zero. In the recursive method, the entries below the main diagonal are calculated from the entries on the main diagonal. Since the entries below the main diagonal become very large as they move away from the main diameter, any error in computing the main diagonal entries leads to a large error in computing the entries below the main diagonal. In this paper, we also study the challenges of using the associated Legendre functions of the first kind in the production of gravity field functionals based on a GGM, utilizing MATLAB software. The results show that the gravity potential computation up to degree of 2190 suffers from very large computational errors at latitudes between 57^°32^׳ and 60^°13^׳. We observe that the safe degrees for the gravity potential computation in all latitudes are degrees less than 2065. The critical latitudes and degrees for the gravity calculation are somewhat different. The results indicate that the gravity computation up to degree of 2190 leads to very large errors at latitudes between 57^°41^׳ and 60^°13^׳. In addition, the maximum degree of expansion that grants sufficient accuracy for the calculation of gravity for all latitudes is estimated to be 2071. Therefore, since the usual computational software based on the double-precision format is not suitable for evaluating the current high-degree GGMs, in this research, a new proposal based on the use of the “long double-precision” format is presented and evaluated. Based on our evaluations, the use of the long double-precision format throughout the computational procedure provides sufficient accuracy to compute the gravity field functionals based on the current high-degree GGMs.