Gravity Field Modeling Using Tesseroids with Variable Density in the Vertical Direction

We present an accurate method for the calculation of gravitational potential (GP), vector (GV), and gradient tensor (GGT) of a tesseroid, considering a density model in the form of a polynomial up to cubic order along the vertical direction. The method solves volume integral equations for the gravitational effects due to a tesseroid by the Gauss–Legendre quadrature rule. A two-dimensional adaptive subdivision technique, which automatically divides the tesseroids near the computation point into smaller elements, is applied to improve the computational accuracy. For those tesseroids having small vertical dimensions, an extension technique is additionally utilized to ensure acceptable accuracy, in particular for the evaluation of GV and GGT. Numerical experiments based on spherical shell models, for which analytical solutions exist, are implemented to test the accuracy of the method. The results demonstrate that the new method is capable of computing the gravitational effects of the tesseroids with various horizontal and vertical dimensions as well as density models, while the evaluation point can be on the surface of, near the surface of, outside the tesseroid, or even inside it (only suited for GP and GV). Thus, the method is attractive for many geodetic and geophysical applications on regional and global scales, including the computation of atmospheric effects for terrestrial and satellite usage. Finally, we apply this method for computing the topographic effects in the Himalaya region based on a given digital terrain model and the global atmospheric effects on the Earth’s surface by using three polynomial density models which are derived from the US Standard Atmosphere 1976.


Introduction
An accurate forward modeling method for the computation of gravitational potential (GP), vector (GV), and gradient tensor (GGT) based on a given mass distribution is highly required in many geodetic and geophysical applications. In geodesy, especially in physical geodesy (Heiskanen and Moritz 1967), it is widely applied in computing mass reductions (e.g., terrain correction, residual terrain correction, isostatic correction, and atmospheric correction) for geoid determination, as the gravitational effects due to the topographic or atmospheric masses are usually removed before the geoid modeling and then added back to the geoid heights after the modeling (e.g. , Forsberg 1984;Denker 2013). In terrestrial gravimetry, the gravity records are disturbed by temporal signals caused by oceanic and atmospheric mass variations. To reveal gravity signals reflecting geophysical processes, the disturbing signals must be precisely evaluated and removed (e.g., Neumeyer et al. 2004;Gitlein 2009). Forward gravity modeling is also commonly used for geophysical interpretations by removing the complete Bouguer effect from terrestrial or satellite measurements to isolate the gravity signals related to the geological structure inside the Earth (e.g., Álvarez et al. 2012;Braitenberg 2015). In the inversion of gravity data, an accurate forward solver is crucial for imaging subsurface density anomalies for resource exploration (e.g., Li and Oldenburg 1998) and investigating crust-mantle structures inside the Earth or other planets (e.g., Liang et al. 2014;Uieda and Barbosa 2017).
Newton's integral serves as the basis for the forward gravity modeling, which requires the geometry and density of the mass distribution. In practice, the mass distribution is usually discretized into a set of elementary bodies associated with density information. The forward modeling aims at evaluating the gravitational effects at the computation point by superimposing the contributions from all elementary bodies.
For local applications, in which the Earth can be approximated as flat and thus the computations are carried out in a Cartesian coordinate system, the elementary body is popularly selected as a flat-topped rectangular prism because its analytical solutions for gravitational effects are available (e.g., Nagy et al. 2000;Tsoulis 2000;Nagy et al. 2002;Heck and Seitz 2007;D'Urso 2012). To improve the computational efficiency of the analytical solutions, an alternative is provided by the approximate solutions of the prism integrals which are computed by the formulas based on a Taylor series expansion of the integral kernel (e.g., MacMillan 1930;Fukushima 2020). Because the use of flat-topped rectangular prisms simply approximates the mass distribution, the rectangular prisms with inclined top and bottom faces are alternatively employed (e.g., Smith 2000;Tsoulis et al. 2003;D'Urso 2015a). Rather than rectangular prisms, polyhedral bodies have better capabilities to describe mass sources with complex geometries. Various analytical solutions have been successfully derived for the gravity field of a polyhedral body (e.g., Okabe 1979; Götze and Lahmeyer 1988;Pohánka 1988;Holstein and Ketteridge 1996;Petrović 1996;Holstein et al. 1999;Holstein 2003;Hamayun et al. 2009;Tsoulis 2012;D'Urso 2013D'Urso , 2014aConway 2015;D'Urso 2015b;D'Urso and Trotta 2017;Ren et al. 2017Ren et al. , 2018a. In addition to classical space-domain analytical approaches, the polyhedral gravity forward problem can also be solved in the Fourier domain (e.g., Wu 2018aWu , b, 2019 or by spherical harmonic expansions (e.g., Chen et al. 2019a, b). We note that the literature about gravity field calculations based on polyhedral bodies is huge. Therefore, we just provided a selection of the most relevant papers. More references can be found in the given references. For applications on either a regional or a global scale, which are usually based on spherical coordinates, the curvature of the Earth cannot be directly taken into account for the rectangular prisms. To solve this, additional coordinate transformations are often applied (e.g., Heck and Seitz 2007;Wild-Pfeiffer 2008). Although the polyhedral representation is valid in the spherical coordinate system (e.g., Benedek and Papp 2009), the approximation of the Earth with curvatures definitely requires a large number of polyhedral bodies. Additionally, since the analytical expressions for the gravity field of the polyhedral body are in Cartesian coordinates, the computation of its effect in a local topocentric coordinate system attached to the computation point also needs coordinate transformations. Alternatively, the elementary bodies are typically chosen as tesseroids rather than rectangular prisms or polyhedral bodies. The tesseroid is a natural representation in the spherical coordinate system (e.g., Heck and Seitz 2007;Wild-Pfeiffer 2008), which is bounded by two spherical surfaces, two meridional planes, and two coaxial circular cones (e.g., Anderson 1976;Grombein et al. 2013). However, analytical solutions for the gravitational effects of the tesseroids expressed in Newton's volume integrals are not available, and one has to employ numerical approaches.
There are two kinds of numerical approaches, which have been investigated extensively for the evaluation of tesseroids, namely the Taylor series expansion (TSE) method and the quadrature method. For the TSE method, Newton's integral is solved by expanding the integral kernel in a Taylor series expansion up to a certain degree, followed by an integration of the subsequent volume integrals (e.g., Heck and Seitz 2007;Wild-Pfeiffer 2008;Grombein et al. 2013;Shen and Deng 2016). The TSE method is fast and can produce accurate results at low latitudes. However, its accuracy decreases rapidly toward the polar regions due to the significant change of the tesseroid surface from the equator (with an almost rectangular horizontal shape) to poles (with nearly a triangular horizontal shape). In the quadrature method, the approximate solutions of Newton's volume integrals are directly solved by using numerical quadrature rules, in which the Gauss-Legendre quadrature (GLQ) is popularly applied, leading to the so-called GLQ method (e.g., Asgharzadeh et al. 2007;Wild-Pfeiffer 2008;Li et al. 2011;Uieda et al. 2016;Deng and Shen 2018;Soler et al. 2019). In some approaches, the volume integral is analytically integrated along the radial direction first. The resulting surface integral is then numerically solved by using the GLQ rule (e.g., Wild-Pfeiffer 2008), or by the split quadrature method using the double exponential quadrature (DEQ) rule, resulting in the so-called DEQ method (e.g., Fukushima 2017Fukushima , 2018. Most recently, Zhong et al. (2019) proposed a new method for computing the gravity field of a tesseroid, in which the original volume integral is first transformed into two surface and four edge integrals, and then the GLQ rule is adapted to evaluate these integrals. In comparison with the TSE method, the quadrature method (e.g., GLQ method) is able to provide more accurate approximations on the cost of more computational time (e.g., Lin and Denker 2019). When the computation point approaches the tesseroids, the integrands become singular. As a consequence, both methods give inaccurate approximations, especially for computing the GV and GGT.
In recent years, many efforts have been made to improve the computational accuracy for computation points near and on the tesseroid. For example, Heck and Seitz (2007), Tsoulis et al. (2009), andTsoulis (1999) suggested to replace the tesseroids near the computation point by the prismatic or polyhedral bodies, so that their analytical solutions can be used. Marotta and Barzaghi (2017) proposed to first rotate the tesseroids into a new coordinate system whose polar axis is the line connecting the computation point and the Earth's center and then analytically calculate the gravity attractions due to the rotated tesseroids. Fukushima (2017, 2018) developed a novel method, which first computes the GP of the tesseroids at an arbitrary point with very high precision by the powerful DEQ rule and then approximates the GV and GGT by numerical partial differentiation of the computed GP. In Li et al. (2011), Grombein et al. (2013), Uieda et al. (2016), Lin and Denker (2019), Soler et al. (2019), and Zhong et al. (2019), the improvement of the approximation is achieved by regularly or adaptively subdividing the tesseroids close to the computation point into smaller tesseroid elements along both horizontal and vertical dimensions or only in the horizontal dimension first, and then summing all effects of the subdivided tesseroid elements computed by the GLQ rule. Among the above-mentioned approaches, the last one is easy to implement and further requires no elementary body conversion, flat Earth approximation, coordinate transformation, tesseroid rotation, and numerical differentiation. Therefore, we follow this approach for tesseroid modeling in this study. Since the numerical comparisons in Lin and Denker (2019) demonstrated the superiority of the GLQ method along with the adaptive subdivision of the tesseroids in the horizontal dimension, it is taken as the basis for the new method presented here. Compared to the method of Lin and Denker (2019), the new method cannot only compute GP and GV for the tesseroid, but also GGT. In addition, a comprehensive analysis of the impact of the vertical and horizontal dimensions of the tesseroid as well as the number of GLQ nodes on the solutions is given. The resulting guidelines for tesseroid modeling are incorporated in the new method, with the aim of ensuring high accuracy in terrestrial, airborne, and satellite applications on a regional or global scale.
So far, in contrast to the significant developments regarding the density modeling for prismatic and polyhedral bodies, the tesseroid-based forward gravity modeling methods mostly deals with homogeneous tesseroids. To fulfill the increasing requirements for applications on a regional or global scale with complex density environments, it is of great interest to extend the constant density model for the tesseroid to a linear or nonlinear one. For this purpose, several efforts have already been made. In the DEQ method proposed by Fukushima (2017Fukushima ( , 2018, the density model for the tesseroids has the form of a polynomial function of the radius up to an arbitrary order. Lin and Denker (2019) presented new formulas for computing the gravity field of the tesseroids with constant density and the density varying linearly along the vertical direction by the GLQ method and the TSE method. However, these formulas are limited to the computation of the GP and GV. Most recently, Soler et al. (2019) proposed a new method, which is able to compute the gravitational effects of the tesseroids whose density varies with depth according to an arbitrary continuous function by the GLQ method. In the computation, the continuous density function for the tesseroid is in fact discretized into a set of smaller tesseroid elements along the radial direction by an adaptive density-based discretization algorithm, while each tesseroid element has a linear density function of the radius. Thus, the tesseroid modeling actually uses the formulas corresponding to the linear density model. Furthermore, the gravitational effects to be computed are also restricted to the GP and GV in their method. Although the DEQ method allows for a direct evaluation of the tesseroids with a polynomial density function of the radius, to the best of our knowledge, there is not such a method using the GLQ rule, which can directly and precisely compute the GP, GV, and GGT due to the tesseroids with density varying nonlinearly along the vertical direction. In this study, an attempt to consider the density model in the form of a polynomial function of the radius up to cubic order is given. Compared to the complex formulas appearing in the DEQ method, the formulas for the new method should be much simpler.
The remainder of the paper is organized as follows. In Sect. 2, we present the method used to evaluate the gravitational effects for the tesseroids with density varying as a polynomial function up to cubic order along the vertical direction. Section 3 first conducts various numerical tests to investigate the influence of the tesseroid dimension, the number of GLQ nodes, and the computation point's position on the solutions, and then gives a comparison of several tesseroid approaches including the DEQ method. In Sect. 4, two applications of the new method are presented to examine its practical applicability. In the first application, the topographic effects in the Himalaya region are calculated by the new method and then compared with those computed by the program TC (Forsberg, 1984). In the second application, global atmospheric effects on the Earth's surface are calculated and compared using three different density models. Finally, some conclusions based on the summary of the numerical results are given in Sect. 5.

Gravitational Effects for Tesseroids
Assuming that the density within a tesseroid is varying along its vertical direction according to an N-th order polynomial function, the density model can be described as where h ′ , h 1 and h 2 denote the altitude of an arbitrary point inside the tesseroid, the bottom and top face of the tesseroid, respectively, and a n are the given density coefficients with respect to order n. In this study, we deal with the case of N = 0, 1, 2, 3 , corresponding to the constant, linear, and polynomial density models up to quadratic and cubic orders. For the sake of conciseness, we call them constant, linear, quadratic, and cubic density models hereafter. For computations in spherical coordinates, the somewhat artificial concept of density as a function of radial distance rather than altitude is required. After introducing the relations r � = h � + R , r 1 = h 1 + R , and r 2 = h 2 + R into Eq. (1), we obtain the density model as a function of radial distance r ′ : The expressions for the polynomial coefficients n based on a n up to cubic order are given in Table 1. Replacing the constant density in Eq. (21) of Grombein et al. (2013) by r ′ , optimized formulas for computing the GP, GV, and GGT due to the tesseroid having a dimension of r 1 , r 2 × 1 , 2 × 1 , 2 and an N-th order polynomial density model in a local north-oriented Cartesian coordinate system with its x-axis pointing toward north, y-axis toward east and z-axis upward can be expressed as

3
where i, j ∈ {x, y, z} , (r, , ) and r ′ , ′ , ′ denote the spherical coordinate of the computation point and the running integration point, G is the gravitational constant that is set as 6.672 × 10 −11 m 3 kg −1 s −2 in this study, K V n , K V i n , and K V ij n are the integral kernels which are specified in Table 2, and (2) based on the density coefficients a n in Eq. (1) up to cubic order Order n Expression  Since the volume integrals in Eq.
(3) comprise elliptic integrals, no analytical solutions are available. Alternatively, they are commonly solved by the GLQ rule (e.g., Stroud and Secrest 1966;Asgharzadeh et al. 2007;Li et al. 2011;Uieda et al. 2016), resulting in the least-squares numerical approximations as given by where The functions f r i ,̂j,̂k are generalizations of the appropriate integral kernels in observation and source point coordinates, w r i , w j , w k are the GLQ weights corresponding to the GLQ nodes -Pfeiffer, 2008), and N r , N , N are the number of GLQ nodes along the radial direction, latitude and longitude, respectively. In practical computation, the GLQ nodes must be scaled to the integration domain r 1 , r 2 × 1 , 2 × 1 , 2 using Equation (5) shows that the least-squares gravitational effects of a tesseroid at each computation point are actually computed by a weighted sum of N r × N × N equivalent point mass effects, where each point mass is located at the source coordinate (r i ,̂j,̂k) inside the tesseroid. If N r = N = N = 1 , the GLQ method is equivalent to the point mass method.
In practical implementation, for those tesseroids near the computation point, the volume integrals in Eq. (3) are only used when the radial distance r of the computation point satisfies r ≥ r 2 or r ≤ r 1 . If r 1 < r < r 2 , the integration domain r 1 , r 2 is divided into two parts, namely r 1 , r and r, r 2 . As an example, Eq. (3a) is rewritten as (4c) cos = sin sin � + cos cos � cos , (4d) C = cos sin � − sin cos � cos , Equations (3b) and (3c) can be expressed in a similar way. It is clear that the original tesseroid is split into two individual tesseroids by the plane passing through the computation point. The gravitational effects for the original tesseroid are now the sum of the effects due to the two divided tesseroids. The reason for this division is to ensure the effectiveness of applying the two-dimensional adaptive subdivision (2DAD) technique (see Sect. 2.2) for improving the computational accuracy. For the tesseroids far away from the computation point, such a treatment is not necessary. This treatment also works if the computation point is inside the tesseroid when computing GP and GV, but not for GGT because it is not defined for points on the boundary surface.

Two-Dimensional Adaptive Subdivision Technique
The adaptive subdivision originates from Li et al. (2011) and Uieda et al. (2016), in which a tesseroid near the computation point is irregularly subdivided into smaller tesseroids along both the horizontal and vertical dimensions, i.e., three-dimensional adaptive subdivision (3DAD) technique. Based on their procedure, Lin and Denker (2019) simplified this strategy by only performing the subdivision along the horizontal dimension and numerically proved its superiority to the 3D case. Here, we follow the 2DAD technique. To determine whether a tesseroid has to be subdivided, the relation is checked for each i ∈ { , } . If the inequality holds for both the latitude and longitude directions, no subdivision is made. Otherwise, the tesseroid is discretized along the direction that failed the condition. In the above equation, 0 is defined as the distance between the geometrical center r 0 , 0 , 0 of the top or bottom face of the tesseroid and the computation point (r, , ) : cos 0 = sin sin 0 + cos cos 0 cos 0 , L i denotes the dimension of the tesseroid along the latitude and longitude: D is a positive value referred to as the "distance-size ratio," which needs to be controlled by the user to achieve a specific requirement. A large D means more tesseroids around the computation point are to be discretized. Some numerical examples on the selection of D can be found in Uieda et al. (2016) and Lin and Denker (2019). Relevant numerical investigations on the optimal choice of D for the new method will be shown in Sect. 3.2. Notice that 0 defined in Eq. (10) is slightly different from that in Lin and Denker (2019). The new definition is more suited for the 2DAD technique based on numerous numerical tests.

Remarks on Numerical Singularity
When observing the integral kernels given in Table 2, we find that the amplitudes of the kernels are inversely proportional to the distance between the computation point and the mass element for GP. As for GV and GGT, the amplitudes are inversely proportional to 3 and 5 , respectively. Since the volume integrals as given in Eq.
(3) are evaluated by the GLQ rule which is in fact a weighted sum of the effects of point masses located inside the tesseroids (see Eq. 5), the numerical singularity, i.e., = 0 , does not exist for exterior points. When the computation point is near or on the tesseroid surface, a very small between the computation point and the GLQ node may exist. Due to the limitations of the computational arithmetic, a rounding off error may occur in the computation of . According to the role of playing in the integral kernels, the error in the computed will definitely be amplified, thus affecting the computational accuracy. For a better explanation, we let =̂ + with , ̂ , and being the computed distance, true distance, and existing error, respectively. For the sake of simplicity, we assume here that the error is positive. We then define the following three ratios to measure the degree of error effect on the integration: with = ∕̂ . The ideal case is that equals zero (i.e., = 0 ), and thus, the ratios reach 1. In practice, a nonzero always exists (i.e., ≠ 0 ). A large deviation between the ratio and 1 indicates a large error effect on the integration, and vice versa. If is a very small value and approaches zero, its impact on the solution is negligible. This often happens when ̂ is much larger than . Therefore, it is easy to see that the numerical result is usually accurate when the computation point is far away from the source mass. However, if ̂ is small, the resulting might not be sufficiently small. As a consequence, its effect on the solution is significant. Figure 1 gives an example of the above three ratios as a function of . Obviously, the degree of the error effect on the integration for three kinds of gravity field quantities follows: GGT > GV > GP. This means that, in the same computational environment, the computation of GGT with high accuracy is the hardest, and then GV and GP.

An Extension Technique
As discussed in Sect. 2.3, too small distances between the GLQ nodes and the computation point may lead to an inaccurate integration due to the significant amplification of the errors in . Accordingly, the evaluation of a tesseroid with a small vertical dimension r = r 2 − r 1 on its surface by the GLQ rule may give an inaccurate solution because there exist one or more GLQ nodes being close to the computation point. In contrast, the tesseroid with a large r might be better approximated. To improve the evaluation of a tesseroid having a small r , we here propose a so-called extension technique. The basic idea is that we replace the evaluation of the tesseroid with a small r by the evaluation of two tesseroids which are related to the original tesseroid and have much larger r , and then make a difference between the gravitational effects due to the two tesseroids to yield the effect of the original tesseroid. Based on the simple illustration in Fig. 2, the extension technique is described as follows: 1. We extend the original tesseroid ( ) along its vertical dimension upward or downward to a predefined value D e , depending on the vertical location of the computation point with respect to the tesseroid ( ). If r ≥ r 2 , then the extension is made downward (see left panel of Fig. 2). If r ≤ r 1 , the upward extension is made (see right panel of Fig. 2); Here the star, tesseroids marked by , , and denote the computation point, original tesseroid, extending tesseroid, and extended tesseroid 2. We compute the gravitational effects caused by the extended tesseroid ( ) and extending tesseroid ( ) by the GLQ rule along with the 2DAD technique; 3. We subtract the gravitational effect of the extending tesseroid ( ) from that of the extended tesseroid ( ), yielding the effect due to the original tesseroid ( ).
Furthermore, some issues need to be clarified when using the extension technique. They are: • The preliminary condition of using the extension technique. It is only used when r ≥ r 2 or r ≤ r 1 , but not for r 1 < r < r 2 . In the case of r 1 < r < r 2 , we first split the tesseroid into two individual tesseroids by Eq. (8), and then consider each of them as a new tesseroid. Because Eq. (8) is only applied to the tesseroids near the computation point in practical computation, the use of the extension technique is also limited to these tesseroids. For those tesseroids far from the computation point, the evaluation of the tesseroid with a small r is sufficiently accurate without using the extension technique. • The definition of et r and D e used in the extension technique. Supposing that a tesseroid with r satisfies the preliminary condition of using the extension technique, we often define a threshold of the vertical dimension et r to judge whether it needs the extension technique or not. If r ≤ et r , we apply the extension technique using a defined D e that is in fact the vertical dimension of the extending tesseroid ( ). Otherwise, no extension technique is used. Since no rules are available for the proper choice of et r and D e , they are empirically determined based on numerical tests. • The density models for tesseroids ( ), ( ), and ( ). In fact, the density model r ′ is only valid for ( ) as can be seen from Eq. (2). But it is also assigned to ( ) in the extension technique, which is actually virtual. As a consequence, the density model for ( ) consists of two parts, namely the real part in ( ) and the virtual part in ( ). Since both density models for ( ) and ( ) are continuous functions, it is straightforward to compute the gravitational effects due to the two tesseroids and then make a difference between them to eliminate the effects caused by the virtual densities in ( ) and ( ) and finally yield the effect of ( ) with the real density model.
As an example of computing the GP of the tesseroid on the point with r ≥ r 2 by applying the extension technique, the volume integral as given in Eq. (3a) is rewritten as When the computation point is located on or below the bottom face of the tesseroid (i.e., r ≤ r 1 ), we then have For the computation of GV and GGT by using the extension technique, the volume integrals are expressed similarly to Eqs. (13) and (14).

Validation Tests
This section aims at examining the computational accuracy of the new method in the double-precision environment. We select a spherical shell as a research object as its gravitational effects can be analytically computed. The spherical shell is located on a spherical Earth with a default radius R = 6378.137 km and has a thickness H s and a polynomial density model r ′ (see Eq. 2). The closed formulas for computing the GP (V), the radial component of GV ( V z ), and the radial-radial component of GGT ( V zz ) of the spherical shell at an arbitrary point are given in Appendix 1. If not differently specified, the computation point is on the shell surface, namely its height h above the spherical Earth equals H s . Notice that V zz becomes indefinite on the shell surface. In this case, we shift the computation point slightly above the shell surface such as h = (1 + e)H s , where e is chosen as 10 −15 here. The computations in the following tests are executed using Fortran codes parallelized by OpenMP with total 40 threads on a server hosting 2 Intel(R) Xeon(R) E5-2660V3 @ 2.60 GHz CPUs with 10 cores per CPU and 2 threads per core. The GP, GV, and GGT have the units of m 2 s −2 , mGal, and Eötvös ( E , 1 E = 10 −9 s −2 ), respectively. Four polynomial density models up to different orders ( N = 0, 1, 2, 3 ) are computed based on the density coefficients a 0 = 10 3 kg m −3 , a 1 = 2 × 10 −2 kg m −4 , a 2 = 2.5 × 10 −5 kg m −5 , and a 3 = 5 × 10 −10 kg m −6 for the following tests. Figure 3 illustrates the density changes with respect to the altitude ranging from 0 m to 10 km above the spherical Earth.

Influence of the Tesseroid Dimension and Number of GLQ Nodes
A number of numerical experiments are carried out in this section, with the aim of investigating the impact of the tesseroid vertical dimension r and horizontal dimension h , as well as the number of GLQ nodes ( N r ∕N ∕N ) on the solutions. Analytical solutions of a homogeneous spherical shell/cap at the computation point on the north polar axis are easily found in the literature (e.g., Heck and Seitz 2007;Grombein et al. 2013). Therefore, the density for the test spherical shell/cap is assumed to be constant if not differently specified in this section. Let us first examine the computational accuracy of the new method without using the extension technique for the modeling of spherical shells with different thickness H s , which varies from 1 m to 10 km. For each spherical shell, the evaluation point is located at the north pole and resides on the shell surface. In the modeling, each spherical shell is regularly divided into 9331200 tesseroids with h = 5 � and r = H s . The gravitational effects of each tesseroid are evaluated by the GLQ rule along with the 2DAD technique. Different combinations of N r ∕N ∕N (i.e., 3/1/1, 3/2/2, 3/3/3, 3/4/4, 3/5/5, 3/10/10, and 5/1/1, 5/2/2, 5/3/3, 5/4/4, 5/5/5, 5/10/10) are compared. Because a large variation between the magnitudes of the gravitational effects of the spherical shells with different H s exists, we prefer to use the relative errors to measure the computational accuracy rather than the absolute errors. Through comparing the solutions obtained by a sum of the gravitational effects of all tesseroids to the analytical solutions of the spherical shells, the resulting relative errors are shown in Fig. 4. It is clear that the relative errors generally decrease as H s increases. When N r is fixed, the increase in N and N improves the computational accuracy. In general, N ∕N = 3∕3 is sufficient for tesseroid modeling in regard to the computational accuracy and efficiency. We also find that V zz is larger than 1 when H s ≤ 100 m , which means the computed V zz is totally unacceptable, representing 100% actual errors.
To investigate the cause of different performances of the method for various shell thicknesses, the computational accuracy of a single tesseroid is analyzed. The single tesseroid is defined as a sector of a spherical zonal band, and thus, its analytical solutions at the computation point, which is located at the north pole, can be computed based on the procedures described in Appendix D.3 of Lin and Denker (2019) and the formulas of a spherical cap given in Heck and Seitz (2007). Additionally, we derived the analytical formula for computing V zz of a spherical cap in Appendix 2. In order to see the dependency of the computational accuracy on the distance between the tesseroid and the computation point, the tesseroid moves from the north pole to the south pole. r is selected as 1 m, 10 m, 100 m, 1 km, and 10 km for comparison purposes. Two horizontal dimensions h = 5 � and 3 ′′ are also considered. The latter h is consistent with a currently high-resolution global digital terrain model (DTM) like SRTM3. The number of GLQ nodes used is set as 3 along each dimension. Let be the typical linear size of the tesseroid, which is computed as the longest distance between two vertices on the top and bottom faces of the tesseroid, be the distance between the tesseroid geometrical center and the computation point, and −1 = ∕ , the relative errors V , V z , and V zz as a function of −1 are illustrated in Fig. 5. With the increase in −1 (i.e., the tesseroid moves away from the computation point), the relative errors decrease rapidly first, then remain stable, and finally increase again. It is clear that the errors mainly stem from the poor evaluation of the tesseroids in the vicinity of the computation point, especially for the nearest one. We observe that a larger r usually results in smaller relative errors. This explains why the relative error decreases as H s increases in Fig. 4. Comparing the results for h = 3 �� to those for h = 5 � , it reveals that a smaller h usually yields larger relative errors. It is now evident that the tesseroid dimension affects the solutions. We also notice that V zz is larger than 1 when evaluating the nearest tesseroid with a small r . Therefore, the evaluation of the nearest tesseroid must be improved for practical usage. To achieve this aim, one possible way is to use an optimal number of GLQ nodes along each dimension for the nearest tesseroid instead of using a fixed number 3 used here. Based on the results shown in Fig. 4, it is sufficient to fix the values of N and N to be 3 for precise tesseroid modeling. We then emphasize the optimal choice of N r and design the following experiment to illustrate the effect of using different N r on evaluating the nearest tesseroid with different r and h .
The single tesseroid is now chosen to be the one that directly connects the computation point located at the north pole. The vertical dimension of the tesseroid varies from 1 m to 10 km. In addition to using h = 5 � and 3 ′′ , another two cases of using 1 ′ and 30 ′′ are considered. For each r and h , six different N r with the values of 1, 3, 5, 10, 30, 50 are selected for the investigation. Corresponding results are shown in Fig. 6. It clearly shows that the performance of using different N r varies with r , h , and the gravity field quantities to be modeled. With the increase in r , the relative errors decrease at first, then increase rapidly for a small N r , or keep on decreasing for a large N r . In general, a small N r is suited for a small r and vice versa. According to Fig. 6, we divide r into several intervals for each h , and an optimal value of N r is selected in each interval. Table 3 summarizes the results. Notice that only four h are analyzed here, and the choice of optimal N r may be different for other h . Alternatively, one may horizontally split the nearest tesseroid having a different h into several tesseroids with one of the above-analyzed h , so that the rule for optimal N r selection given in Table 3 can be directly applied. The cost is that the computational time will be somewhat increased.
We then repeat the experiment shown in Fig. 5. In practical computation, the optimal N r is applied to the tesseroids satisfying −1 ≤ 10 , and N r = 3 is used for the rest. The relative errors for −1 ≤ 100 are shown in Fig. 7. Comparing Fig. 7 with 5, we find that using optimal N r can reduce the relative errors, in particular for V . However, V zz is still larger than 1 for r = 1 m . We then apply the extension technique (see Sect. 2.4) to evaluate V z and V zz of the tesseroid again, which has a r ≤ et r and satisfies −1 ≤ 3 . The values of et r and D e are numerically selected as 1 km and 4 km. When computing V, no extension technique is needed because V is not sensitive to r (Fig. 7). Figure 8 illustrates the relative errors V z and V zz after applying the extension technique. Obviously, they are significantly reduced when comparing them to those in Fig. 7. V zz is now smaller than 1 for all selected r .
Finally, we repeat the experiment shown in Fig. 4. Besides using the optimal N r , we also consider two cases of applying and not applying the extension technique. In addition to the constant density model, the other three density models are also used. The results are shown in Fig. 9, from which it is observed that the use of the extension technique largely improves the computational accuracy. All relative errors are below 1, implying that the computed gravity field is acceptable. To this end, we summarize two numerical findings need to be concerned when using tesseroids for precise gravity forward modeling: 1. The tesseroid vertical dimension r has a major impact on the solution. The smaller r is, the lower is the computational accuracy. Therefore, we recommend to apply the extension technique to evaluate the GV and GGT of the tesseroid, which has a vertical dimension r ≤ 1 km and satisfies −1 ≤ 3 . For the computation of GP, no extension technique is needed. Fig. 7 The same as in Fig. 5, but only displaying the case of −1 ≤ 100 . The optimal N r is used to evaluate the tesseroids satisfying −1 ≤ 10 , while N r = 3 is applied for the rest Fig. 8 The same as in Fig. 7, but only displaying the relative errors V z (a, c) and V zz (b, d). The extension technique is applied to the tesseroids both satisfying r ≤ 1 km and −1 ≤ 3  ( N r ) dimensions is not an optimal choice for evaluating the tesseroids near the computation point. Here we recommend to fix N and N as 3 and use the optimal N r which varies with r and h (Table 3). In practical computation, the optimal N r is applied to the tesseroids satisfying −1 ≤ 10 , while N r = 3 is used for the rest.

Optimal Selection of Parameter D
In the following, we aim to choose the optimal value of the parameter D used in the 2DAD technique by investigating six spherical shells with H s = 1 m , 10 m, 100 m, 1 km, 5 km, and 10 km. Four density models are used for the selection. The relative errors for the computation point at the north pole in terms of D, which varies from 0 to 20 with a step size of 1, are shown in Fig. 10. Here we only show the cases of using the constant and cubic density models because the error curves for the other two density models look similar to the presented curves. We can see that the relative errors generally decrease as D increases and then remain almost unchanged. The optimal value of D is determined as the one after which the relative errors become stable. Table 4 summarizes the chosen optimal values of D for different H s and density models, from which we find that the optimal values vary with H s , density models, and gravity field quantities. Here, we choose the largest value among all cases as the final optimal value of D for each kind of gravity field quantity, resulting in 14, 8, and 7 for GP, GV, and GGT. Finally, the experiment shown in Fig. 9 is  The four values of D correspond to the case of using the constant, linear, quadratic, and cubic density model, respectively 1 m 10 m 100 m 1 km 5 km 10 km V 3/3/3/1 1/1/1/1 3/3/2/1 5/6/2/1 11/11/3/1 14/9/5/1 V z 3/3/3/6 3/3/3/3 2/3/3/3 5/5/5/4 7/7/7/5 8/8/8/5 V zz 3/3/3/3 3/3/3/3 3/3/4/5 4/4/4/5 7/7/7/7 6/6/6/6 revisited by using the optimal D, and the relative errors are shown in Fig. 11. In comparison with Fig. 9, the error reduction is not obvious for V z and V zz . However, V is evidently reduced for H s > 1 km.

Influence of the Computation Point's Height
The computation point is assumed to be on the shell in the previous tests. In this section, we intend to take a deeper insight into the performance of the new method for calculating the gravitational effects of the spherical shell at the point on the shell, near the shell, far away from the shell, and even inside the shell mass. In general, the computation point can be distributed in any direction (e.g., horizontal or vertical). For the sake of conciseness, they are located along a straight radial line passing through an arbitrarily selected point with = 86.5 • E and = 27.5 • N in this test. The test spherical shell has a thickness H s = 10 km . For computing V and V z , the computation points start at the height h = 200 m above the spherical Earth and end at the height h = 300 km . The computation of V zz is only performed at the points exterior to the shell as it cannot be directly calculated on the surface boundary. Again, four different density models are tested. The relative errors as well as the Laplace equation in relative sense computed by the formula |V xx + V yy + V zz |∕ |V xx | + |V yy | + |V zz | (Ren et al., 2018a) are shown in Fig. 12.
In the case of using the constant density model, the relative errors V z and V zz decrease as h increases, with the values in the range between 10 −14 and 10 −8 for V z , 10 −12 and 10 −6 for V zz . V seems to be less dependent on the height of the computation point with the value varying between 10 −15 and 10 −13 . The relative errors V and V z at the points inside the mass are larger than those at exterior points, but still of high accuracy. The relative Laplace equation decreases with the increase in h. Its variation is similar to that of V zz . Similar information can be obtained from the case when using the linear density model. Fig. 11 The same as in Fig. 9, but using the optimal D In the case of using the quadratic and cubic density models, the relative errors are larger than those computed by using the other two density models. This is because the use of the quadratic and cubic density models raises the numerical complexity and hence reduces the computational accuracy. The dependency on the computation point's height is not significant, in particular for V and V z . For the quadratic density model, the relative errors are in the range between 10 −11 and 10 −10 for V , 10 −11 and 10 −7 for V z , 10 −11 and 10 −5 for V zz . The cubic density model produces larger relative errors, with V and V z ranging from 10 −9 to 10 −7 and V zz varying between 10 −8 and 10 −5 . Similar to the constant and linear density cases, the behavior of the relative Laplace equation follows V zz .
Although the performance of using different density models is different, the computational accuracy is sufficient for practical applications, where the computation point can be an arbitrary point exterior to the tesseroid. When computing GP and GV, the precise evaluation inside the tesseroid is also ensured.

Influence of the Computation Point's Latitude
In the previous section, we examined the computational accuracy due to the height variation for the computation point. Here we depict the accuracy due to the latitude variation, in order to confirm whether the change of the geometrical shape of the tesseroid caused by meridional convergence affects the solution. Due to the symmetry of the spherical shell, the gravity field has been computed every 1 • latitude along the meridian 0.05 • E from the equator to the north pole on the shell surface. Two spherical shells with H s = 1 m and 10 km are tested. For each shell, four density models are considered. Figure 13 shows the corresponding relative errors. It is clear that the relative errors remain at similar levels in terms Fig. 12 Relative errors V (a), V z (b), and V zz (c) as well as the relative Laplace equation (d) in log scale as a function of the height of the computation point above the spherical Earth. The cases of using four density models are shown. Notice that the computation points at the heights between 0 and 10 km are inside the shell mass of the latitude without strong jumps. Obviously, the latitude variation for the computation point has a minor impact on the computational accuracy of the new method. Therefore, the new method can be used for the computation at any place on the Earth and it is well suited for regional and global applications.

A Comparison with the Commonly Used Methods
In this section, we carried out a benchmark test to compare the solutions computed by the new method (denoted as GLQ2DADE method) with those obtained from the following commonly used tesseroid methods: the method based on Taylor series expansion of the integral kernel up to the second order (TSE method; Heck and Seitz 2007) and its extension using the regular subdivision technique (TSERD method; Grombein et al. 2013), the method using the GLQ rule along with the 3DAD technique (GLQ3DAD method; Uieda et al. 2016), the method using the GLQ rule along with the 2DAD technique (GLQ2DAD method; Lin and Denker 2019) and its combination with the TSE method (GLQ2DAD_ TSE method; Lin and Denker 2019), and the method approximating tesseroids by rectangular prisms (PR method; Wild-Pfeiffer 2008) and point masses (PM method; Lin and Denker 2019). Considering the fact that most of the above methods deal with a homogeneous tesseroid, only the constant density model is employed.
The test is analogous to the one shown in Fig. 4. The computation point is located at = 0 • E and = 0 • N , where three computation height levels are used: (1) on the shell surface, (2) 3 km, and (3) 250 km above the shell surface. Figure 14 shows the corresponding relative errors. For each method, the relative errors significantly reduce as the computation height increases. For computation points on the shell surface, the GLQ2DADE method outperforms the others, especially for a small shell thickness. Besides the GLQ2DADE method, the GLQ2DAD and GLQ2DAD_TSE methods provide the smallest V z and V zz , where V for the former is smaller than that for the latter. This is because the GLQ2DAD_ TSE method employs the TSE method to evaluate the tesseroids far away from the computation point, while the far-zone effect has a significant impact on the computation of V. For the computation point at the height 3 km above the shell surface, the relative errors are evidently reduced in comparison with the surface case. The relative errors V and V z are smaller than 1 for all methods. For V zz , only the TSE and PM methods fail to produce the values less than 1. All GLQ methods provide similar relative errors which are smaller than those obtained from the other methods. When we move the computation point to the height 250 km above the shell surface, all methods are able to compute acceptable solutions. In this case, the TSE method provides solutions which are comparable to those computed by the GLQ methods. Thus, the TSE method may replace the GLQ method in satellite applications. In general, the GLQ2DADE method is the best in terms of computational accuracy among all mentioned methods, no matter whether the computation point is on the ground or above it.

A Further Comparison with the TSE-Based Method
From the tests in Sect. 3.5, we notice that the TSE method performs well in the case of computation points at the height of 250 km above the shell surface. This is mainly due to the large distances between the evaluation point and source masses, making the approximations accurate and stable. Regarding the higher computational efficiency for the TSE method, it is straightforward to replace the GLQ2DADE method by the TSE method if the distance between the computation point and source mass is larger than a certain value. In  V (a, d, g), V z (b, e, h), and V zz (c, f, i) for different tesseroid methods as a function of the spherical shell thickness H s (from 1 m to 10 km) in log-log scale. Three levels of computation heights are shown: on the shell surface (a-c), 3 km (d-f) and 250 km (g-i) above the shell surface order to briefly answer the questions about the choice of the distance for replacement and the improvement of computational efficiency, we design the following simple numerical tests. The test is analogous to the one shown in Fig. 12, where the homogeneous spherical shell with the thickness H s = 10 km is selected as test model. The horizontal position of the computation point is at = 0 • E and = 45 • N . We then let the computation point move from the shell surface (i.e., h = 10 km ) up to the height of 260 km along the radial direction with a step size of 1 km, resulting in 251 samples. Four computational methods, namely the GLQ2DADE method, the GLQNSE method, the TSE2DAD method, and the TSE method, are used for the approximation. The GLQNSE method is in fact the version for the GLQ2DADE method where the 2DAD technique is not applied, while the TSE-2DAD method is the one where the 2DAD technique is applied to the TSE method. The purpose of using these four methods is on the one hand to see the impact of using the 2DAD technique on the solution and on the other hand to see the performance of the methods themselves. Figures 15 and 16 present the relative error and computational time at each evaluation point.
We start with the relative errors. The use of the 2DAD technique definitely improves the approximations at the sites near the shell surface. However, the improvements become less significant with the increase in the evaluation height. For instance, the GLQ2DADE method is almost equivalent to the GLQNSE method when the evaluation point is 30 km above the shell surface and higher. Comparing the TSE2DAD method to the TSE method, we notice that the application of the 2DAD technique does not always improve the approximations. The sudden jumps for the TSE2DAD method are relevant to the optimally selected values for D (see Sect. 3.2). The GLQ-based method generally outperforms the TSE-based method whether the 2DAD technique is applied or not. We also find that the TSE method provides the approximations with nearly the same accuracy as those computed by the GLQ2DADE method when the computation point is 40 km above the shell surface and higher. From another perspective, we may roughly say that the TSE method can replace the GLQ2DADE method in evaluating the tesseroid whose distance to the computation point is larger than 40 km. Now let us move on to the computational time. Figure 16 clearly shows that the evaluation at sites on or near the shell surface takes slightly more time than those far from the surface. The application of the 2DAD technique has a minor impact on the computational time for the GLQ-based method, but not for the TSE-based method. For a single-point evaluation, the TSE method is the fastest among the four methods and it is about 7-8 times faster than the GLQ-based method. A summary of the computational times is given in Table 5.
Regarding the fact that the GLQ2DADE method can precisely evaluate the tesseroid near the computation point which can be located at any place on the Earth and the TSE method can quickly and accurately evaluate the tesseroid far from the computation point, a combined use of them is preferred in practical applications. On the one hand, the combination can keep the approximation as precise as possible and, on the other hand, it can largely reduce the computational time. However, the issue of the combined approach is beyond the scope of this work. Thus, it will not be discussed in the following, but will be kept for future research.

A Comparison with the DEQ Method
In Sects. 3.5 and 3.6 , we have compared the GLQ2DADE method with several commonly used methods for approximating a homogeneous spherical shell. Here, we examine the computational accuracy of the new method more extensively by comparing it to the DEQ method which could serve as a reference to complement and elaborate existing tesseroid approaches. Two test models will be used here, namely a single tesseroid and a spherical shell. The first test model is used to investigate the performance of the two methods at the evaluation points where the gravity field exhibits nonanalyticity, namely on the surfaces of, along the edges of, at the vertices of the tesseroid, and further inside it. The use of the second test model is to see how precise the two methods can approximate a spherical shell if the computation points are located at different latitudes and altitudes. The implementation of the DEQ method is based on two Fortran programs xqtess.f90 and xtess.f90, which are provided by Fukushima (2018) and conducted in the quadruple-and doubleprecision environments, respectively. Since the two programs deal with the homogeneous tesseroids, the densities for the two test models are assumed to be constant for the sake of convenience. Again, the gravity field quantities to be compared are V, V z , and V zz .

Tests Based on a Single Tesseroid
The single tesseroid is selected as the same as the one used in Figs 5−9 of Fukushima (2018). It is a three-dimensional volume defined as 27 • N ≤ ≤ 28 • N , 86 • E ≤ ≤ 87 • E , and 6340 km ≤ r ≤ 6390 km , where R = 6380 km is the radius of the spherical Earth, i.e., the heights of bottom and top faces of the tesseroid are − 40 km and 10 km. The evaluation points are located on the lines, which are on the surfaces of, along the edges of the tesseroid, and passing through it. Furthermore, the evaluation at vertices of the tesseroid is performed. The relative error computed by |a double ∕a quadruple − 1| is used to measure the computational accuracy, where a double means the value computed by either the GLQ2DADE or DEQ method conducted in the double-precision environment, and a quadruple is obtained from the DEQ method conducted in the quadruple-precision environment as it proves to be sufficiently accurate (Fukushima, 2018).
We begin with V , V z , and V zz , along the tesseroid edges. Figure 17a-c presents the relative errors at three test lines which coincide with the edges, namely Line_A, Line_B, and Line_C. Their definitions are given in Table 6. We omit the results on the other nine edges since there are no significant differences from these three representative cases. When analyzing the relative errors along each test line, they are classified into two groups by their locations. The first group contains the relative errors at the points exactly located on the edge, while the relative errors at remaining points belong to the second group. It is obvious that the relative errors V computed by the DEQ method are smaller than those obtained from the GLQ2DADE method in both groups. In the first group, the DEQ method provides smaller V zz , while the relative errors V z for the GLQ2DADE method are smaller in Line_A and Line_B but larger in Line_C. In the second group, both methods provide similar V z , while the relative errors V zz for the GLQ2DADE method are smaller.
Next come the errors on tesseroid surfaces. Figure 17d-f shows the relative errors at three test lines on the surfaces, namely Line_D, Line_E, and Line_F (Table 6). Again, we omit the results on the other three surfaces since no meaningful differences from these 1 3 three representative cases are observed. The analysis procedure follows that used in the case where the evaluation points are along the edges, and the findings are rather similar to that case.
We then move on to the errors inside the tesseroid. Because there is an infinite number of lines passing through the tesseroid, we only select three as representatives for comparison, namely Line_G, Line_H, and Line_I (Table 6). Since V zz cannot be computed by the GLQ2DADE method if the computation points are inside the tesseroid, the relative errors Table 6 Definition of the test lines used in Fig. 17 When computing V zz , h is shifted to be 10.001 km for Line_A, Line_B, and Line_E; is shifted to be 87.01 • E for Line_C and Line_D; is shifted to be 26.99 • N for Line_F. Because no evaluation point in Line_G to Line_I is located on the tesseroid surface, the three test lines are not shifted for the computation of V zz

Line Parameters
Line_A V zz are only given at exterior points. The results are illustrated in Fig. 17g-i, from which it can be seen that the DEQ method can provide a better approximation of V. When computing V z and V zz , the DEQ method provides better results at interior points, while both methods give similar results at exterior points. Finally, Table 7 shows the relative errors at four vertices of the tesseroid. We omit the other four vertices because they are symmetric to the selected four and no differences are observed. Obviously, the DEQ method performs better than the GLQ2DADE method except for computing V z .
On the basis of the above tests, we find that the DEQ method can compute V more precisely than the GLQ2DADE method at an arbitrary point. When the computation points are on the surfaces of, along the edges of, on the vertices of the tesseroid, and even inside it, the DEQ method usually outperforms the new method. However, the computational accuracy of the new method is still high in these cases. If the computation points are outside the tesseroid, the performance of the two methods is similar. In conclusion, the new method is capable of precisely approximating the gravitational effects of the tesseroid at an arbitrary point, except for the computation of GGT inside the tesseroid.

Tests Based on a Spherical Shell
The parameters for defining the spherical shell are almost the same as those used in Fig. 13. Only one difference is that the shell thickness H s is chosen to be 8 km here. Concerning the high computational burden for the DEQ method, the spherical shell is decomposed into 180 × 360 = 64800 tesseroids with h = 1 • and r = 8 km . The relative error is computed by |b model ∕b analytical − 1| where b model means the modeled gravity field of the spherical shell computed by either the DEQ or GLQ2DADE method and b analytical denotes the analytical value. All computations are conducted in the double-precision environment. Three sets of computation points are distributed along the meridian 0.05 • E from the north pole to the equator with an interval of 1 • , at the height h = 8 km , 10 km, and 100 km. Figure 18 presents the relative errors as a function of the computation point's latitude. In the cases of h = 10 km and 100 km, the relative errors V and V z computed by the DEQ method are smaller than those obtained from the GLQ2DADE method; however, the precision of the computed V zz is slightly better for the GLQ2DADE method. In the case of Table 7 Relative errors V , V z , and V zz at four vertices of the tesseroid It should be noted that, when computing V zz , the site is slightly moved up or down by 10 m along the radial direction, namely h = +10.01 km or h = −40.01 km . For each vertex, the first row shows the relative errors computed by the DEQ method, and the second row corresponds to the GLQ2DADE method h = 8 km , the GLQ2DADE method outperforms the DEQ method in the computation of V z , while the latter method provides more accurate approximations of V and V zz . An interesting fact we noticed is that the accuracy of V zz computed by the DEQ method is evidently affected by high latitudes of the computation points (Fig. 18a). Therefore, the DEQ method should be used with caution when applying it to compute GGT at (near) ground points in the polar region.

Application
To further validate the GLQ2DADE method, we carried out two numerical applications for the computation of topographic effects on a regional scale and atmospheric effects on a global scale in this section.

Topographic Effects
The numerical tests are based on the ETOPO1 DTM dataset (Amante and Eakins, 2009) with 1 ′ resolution. The test region covers the Himalaya region with the most extreme topographic relief existing on the Earth (Fig. 19). The computations have been performed at the Earth's surface, 3 km above the Earth's surface, and at a constant height of 250 km. The three computation height levels roughly simulate terrestrial, airborne, and satellite applications although the aeroplane never flies parallel to the topographic surface in practice (the second height level). At each computation height level, V, g , and V zz due to the given topographic masses are computed. In the case of computing V zz on the Earth's surface, we shift the computation point 1 m above the surface. To reduce the edge effects, the computation area is chosen to be smaller than the DTM area along the latitude and longitude (Fig. 19), where the computation points are regularly distributed with a spacing of 5 ′ . The density is assumed to be constant with a value of 2670 kg/m 3 . Besides using the GLQ-2DADE method, we also implemented the program TC, which uses rectangular prisms, to compute the same results for comparison. It should be noted that all DTM data are used for the computation at each point. Since the results computed by the two methods look very similar, we only show those computed by the GLQ2DADE method as well as the differences between the two methods in Figs. 20,21,and 22. The statistics for all results are summarized in Table 8. From Fig. 20, the computed V over the test region only represents long-wavelength signals, and no detailed information can be seen. This is because V is mainly in the longwavelength domain. The magnitudes of V are quite close for the first two computation height levels, but are clearly reduced for the third height level. The differences between the two methods are significant but look quite similar for different height levels (Table 8), representing a slope descending from the north to south. According to previous numerical tests in Sect. 3, we consider these differences as improvements when using the GLQ-2DADE method instead of using the program TC.
In comparison with Fig. 20, the computed g shown in Fig. 21 contains more detailed gravity field features and more clearly reveals the topographic relief. This is because g has more power at high frequencies than V. The magnitudes of g evidently decrease as long as the computation point moves away from the Earth's surface. As a consequence, the detailed information in g is lost. For instance, g at the height of 250 km represents longwavelength signals. The differences between the two methods also decrease if the computation point moves away from the surface. For the first two height levels, the differences are mainly in the mountainous area. When the computation point is sufficiently high, the differences look like a slope descending from the north to south again, but with rather small magnitudes.
Finally, let us analyze Fig. 22 that presents the computed V zz . In comparison with V and g , V zz has the most power at high frequencies. As a result, their magnitudes dramatically Fig. 19 Topography of the Himalaya test area, where the computation area is marked by the red line. The axes are longitude (degrees East) and latitude (degrees North)

Fig. 20
Top: topographic effects in terms of V evaluated on the Earth's surface (left), at the height 3 km above the Earth's surface (middle), and at a constant height of 250 km above the sea surface (right) by using the GLQ2DADE method in the Himalaya test area. Bottom: the differences between the results computed by the GLQ2DADE method and the program TC Fig. 21 The same as in Fig. 20, but for g decrease with the height increasing for the computation point. Comparing V zz on the Earth's surface to the topography in Fig. 19, the ground V zz shows extreme local gravity field signals over the test region, mainly representing topography variations in the mountainous area. In contrast, V zz at airborne level reveals less local signals from which we can clearly infer the valleys and mountains. Furthermore, V zz at satellite level more precisely depicts the shape of the Tibet plateau than V and g at the same height level. Large differences between the two methods are observed on ground, where most exist in the area with strong topography variations. These differences are greatly reduced in the airborne case. In the satellite case, the differences act as a slope ascending from the north to south. However, their magnitudes are so small that they can be ignored.

Model Setup
The GLQ2DADE method is finally applied to compute global gravitational effects on the Earth's surface caused by atmospheric masses. Evaluation points are on a 45 � × 45 � global grid derived from the ETOPO1 DTM dataset. Both the density and physical bounds of the atmosphere model are needed for the computation. The lower bound of the atmosphere coincides with the Earth's surface that is represented by a global 5 � × 5 � DTM sampled from the ETOPO1. The upper bound is chosen to be 60 km above the sea level, containing more than 99% of the atmospheric mass.
For the sake of simplicity, a static atmosphere model, only considering the vertical density variation, is employed in this application. Thus, we choose the US Standard Atmosphere 1976 (USSA1976). The atmospheric densities derived from the USSA1976 Fig. 22 The same as in Fig. 20, but for V zz model from the sea surface up to 60 km are illustrated in Fig. 23a. In our implementation, we first evenly divide the USSA1976 model into six layers. We then apply a least-squares adjustment to fit the density model h ′ to the densities derived from the USSA1976 model in each layer. Finally, we transform the density model h ′ into the one r ′ that is required for tesseroid modeling. Three density models, namely the linear, quadratic, and cubic density models, are employed here, to examine the impact of using different density models on the computation of atmospheric effects. Table 9 presents the estimated coefficients a n for each density model h ′ in each layer. The relative differences ( | m − us |∕| us | ) between the densities derived from the density model Table 8 Statistics of the topographic effects in terms of V ( m 2 ∕s 2 ), g (mGal), and V zz (E) computed by the GLQ2DADE method (index GLQ) and the program TC (index TC), as well as their differences (index GLQ-TC) The evaluation is performed on the Earth's surface (index S), 3 km above the Earth's surface (index 3 km), and at a constant height of 250 km (index 250 km). Relative difference represents the ratio between the range in the differences and the range in the signals ( m ) and those from the USSA1976 model ( us ) in each layer are shown in Fig. 23b. The corresponding root mean square (RMS) of the relative differences is also given in Table 9. It can be seen that the relative differences are generally below 10%, 1% and 0.1% for using the linear, quadratic and cubic density models. It indicates that the polynomial density model up to a higher order can better fit the USSA1976 model. It should also be noted that the relative differences for the first layer are smaller than those for the other layers. Regarding the fact that the first layer has the largest densities and directly connects the computation point, the good fit can ensure high accuracy for the computed atmospheric effects. The influence of larger relative differences for the other Table 9 Estimated atmospheric density coefficients a n (see Eq. 1) up to cubic order RMS rel represents the root mean square of the relative differences Layer no. Altitude ( km) a 0 (kg/m 3 ) a 1 (kg/m 4 ) a 2 (kg/m 5 ) a 3 (kg/m 6 ) RMS rel (%)  . 23 Densities derived from the USSA1976 model as a function of the altitude (a) and the relative differences between the densities derived from the polynomial density model and the USSA1976 model at different altitudes (b). The RMS of the relative differences for each layer is given in Table 9 1 3 layers should be much less significant because they are far from the computation point and have rather small densities.
Each atmospheric layer is decomposed into a set of tesseroids with the same horizontal dimension h . In the first layer, h is set as 5 ′ to meet the resolution of the used DTM, resulting in 9331200 tesseroids. h is selected as 15 ′ for the remaining layers (i.e., 1036800 tesseroids in each layer) because on the one hand the tesseroids are far from the computation point on ground and on the other hand the computational burden can be considerably reduced. Besides the variable vertical dimensions for the tesseroids in the first layer due to the topographic relief, the other tesseroids have the same vertical dimension as the layer thickness, i.e., r = 10 km.

Results and Discussion
According to the model settings described in Sect. 4.2.1, we compute the global atmospheric effects on the Earth's surface in terms of V, g , and V zz . Notice that the lower surface of the tesseroid that directly contacts the computation point is slightly shifted 1 cm upward when computing V zz . This is because V zz cannot be calculated at the point on the tesseroid surface. As the results computed by using different density models look quite similar, only the results based on the cubic density model are shown in Fig. 24. It clearly shows high correlation between the computed atmospheric effects and the topography. In addition, we give the statistics of the atmospheric effects corresponding to the three density models as well as their differences in Tables 10 and 11 . Table 10 gives that the variation of the atmospheric effects over the Earth is about 2.4 m 2 ∕s 2 for V, 0.23 mGal for g , and 320 mE for V zz . The variations are mainly due to the consideration of the topographic relief in the computation. Table 11 summarizes the influence of using different density models. In general, the results computed by the three density models are quite close. This may be attributed to the fact that the atmospheric densities are very small, thus making the solutions not sensitive to the used density models. More specifically, the quadratic and cubic density models provide almost the same solutions, while both of them are slightly different from the linear model. The RMS of the differences between the quadratic/cubic and linear density models is about 0.1 m 2 ∕s 2 for V, 0.001 mGal for g , and 0.2 mE for V zz . Obviously, the differences are very small for both g and V zz and can be neglected. However, the differences for V cannot be ignored in precise geoid modeling as the corresponding geoid height differences are at the level of about 1 cm. Therefore, the polynomial density model with order higher than 1 is recommended in the computation of atmospheric effects.

Conclusions
In the present work, on the basis of the previous work concerning precise computation of GP and GV due to the tesseroids having constant or linear density in spherical coordinates (Lin and Denker, 2019), we derived a numerical method to evaluate accurately the GP, GV, and GGT of the tesseroids with density varying nonlinearly along the vertical direction. The density model can be a polynomial function up to cubic order which is sufficient for many geodetic and geophysical applications. However, it can further be extended to a higher order if necessary, with the cost of more computational time and lower accuracy. The volume integral in the case of tesseroids is numerically solved by the GLQ rule.
The key point of the new method is the use of a 2DAD technique as well as an extension technique. The latter technique is only applied for the tesseroid whose vertical dimension is smaller than 1 km and for the computation of GV and GGT.
From the numerical tests based on a spherical shell model having a polynomial density model, the new method enables us to compute the gravitational effects of a tesseroid where the evaluation point can be on the surface of, near the surface of, and outside the tesseroid. When computing GP and GV, the evaluation point can also be inside the tesseroid. In this sense, the new method will serve as a reliable tool to evaluate surface, external and/or even internal gravity field. Also, the method is accurate: For approximating a spherical shell with its thickness varying from 1 m to 10 km and the order of the polynomial density model from Fig. 24 The global topography and bathymetry (top left), as well as the global atmospheric effects in terms of gravitational potential V (top right), gravity disturbance g (bottom left), and radial-radial component V zz of GGT (bottom right) evaluated on a 45 � × 45 � grid on the Earth's surface. Notice that the tesseroids with a cubic density model are used for the computation Table 10 Statistics of the atmospheric effects in terms of V ( m 2 ∕s 2 ), g (mGal), and V zz ( mE ) computed by the linear (index l), quadratic (index q), and cubic (index c) density model 0 to 3, the achieved relative errors of the gravity field computed on the shell surface are about 10 −14 −10 −8 for GP, 10 −10 −10 −5 for GV, and 10 −5 −10 −1 for GGT in the double-precision environment. The comparison between the new method and the other already published tesseroid methods also demonstrates its superiority in terms of accuracy, in particular for evaluating a tesseroid with a vertical dimension less than 1 km. Only one method, namely the DEQ method, gives a better performance than the new method when the evaluation is performed on the surfaces of, along the edges of, on the vertices of the tesseroid, and even inside it. For the evaluation at the point exterior to the tesseroid, both methods give similar results. Regarding the lower computational efficiency for the DEQ method, a combined use of the two methods seems to be attractive for practical applications. This topic will be addressed in a future study. The application of the new method for the calculation of topographic effects in the Himalaya region and of atmospheric effects on a global scale verifies its feasibility in real applications. In the first test case (i.e., for topographic effects), the method is also compared to the program TC at three computation height levels. Significant differences between the two methods are obtained, which can be regarded as improvements of using the new method instead of the program TC. In the second test case, the linear, quadratic, and cubic density models are used to describe the atmospheric densities derived from the USSA1976 model. The ground atmospheric effects computed by using different density models do not differ from each other too much, especially for the case of using the quadratic and cubic density models. The RMS of the differences between the quadratic/cubic and linear density models are about 0.1 m 2 ∕s 2 for V. These differences cannot be ignored in precise geoid modeling as the corresponding geoid height differences are at the level of about 1 cm. The relevant differences in terms of g and V zz are very small and negligible. Here, we recommend a polynomial density model with order higher than 1 for precise computation of atmospheric effects.