Gravity Gradient Tensor of Arbitrary 3D Polyhedral Bodies with up to Third-Order Polynomial Horizontal and Vertical Mass Contrasts

During the last 20 years, geophysicists have developed great interest in using gravity gradient tensor signals to study bodies of anomalous density in the Earth. Deriving exact solutions of the gravity gradient tensor signals has become a dominating task in exploration geophysics or geodetic fields. In this study, we developed a compact and simple framework to derive exact solutions of gravity gradient tensor measurements for polyhedral bodies, in which the density contrast is represented by a general polynomial function. The polynomial mass contrast can continuously vary in both horizontal and vertical directions. In our framework, the original three-dimensional volume integral of gravity gradient tensor signals is transformed into a set of one-dimensional line integrals along edges of the polyhedral body by sequentially invoking the volume and surface gradient (divergence) theorems. In terms of an orthogonal local coordinate system defined on these edges, exact solutions are derived for these line integrals. We successfully derived a set of unified exact solutions of gravity gradient tensors for constant, linear, quadratic and cubic polynomial orders. The exact solutions for constant and linear cases cover all previously published vertex-type exact solutions of the gravity gradient tensor for a polygonal body, though the associated algorithms may differ in numerical stability. In addition, to our best knowledge, it is the first time that exact solutions of gravity gradient tensor signals are derived for a polyhedral body with a polynomial mass contrast of order higher than one (that is quadratic and cubic orders). Three synthetic models (a prismatic body with depth-dependent density contrasts, an irregular polyhedron with linear density contrast and a tetrahedral body with horizontally and vertically varying density contrasts) are used to verify the correctness and the efficiency of our newly developed closed-form solutions. Excellent agreements are obtained between our solutions and other published exact solutions. In addition, stability tests are performed to demonstrate that our exact solutions can safely be used to detect shallow subsurface targets.


Introduction
Gravity exploration methods try to identify the anomalous mass bodies in the Earth (Blakely 1996). Gravity signals such as gravity fields and gravity gradient tensors are measured using gravimeters and gravity gradiometers, which can be located on the air-Earth interface, in boreholes, or be carried by marine ships and aircraft, and even by satellites (Nabighian et al. 2005). The amplitudes of the gravity field are inversely proportional to the square of the distance (R) between the observation site and the causative body, that is O(R −2 ) . As for the gravity gradient tensor, its amplitude is inversely proportional to the third power of distance, that is O(R −3 ) . Therefore, compared to the gravity field, gravity gradient tensor signals are more sensitive to the shallow anomalous structures in the Earth (Droujinine et al. 2007;Beiki and Pedersen 2010;Martinez et al. 2013;Beiki et al. 2014;Gutknecht et al. 2014;Li 2015;Ramillien 2017). Mathematically, both gravity field and gravity gradient tensor signals can be formulated as volume integrals over the causative body. When the observation site is located inside the mass body, a local spherical coordinate system can be introduced at the observation site, so that a factor of distance squared ( R 2 ) would be introduced (such as dv = R 2 sin drd d , Blakely 1996), which can further weaken the singularity appearing in the integrands of the gravity signals (Jin 2002). Therefore, from the mathematical point of view, the gravity field can be evaluated without mathematical singularities, as O(R −2 × R 2 ) = O(1) , but the mathematical singularities always remain in the gravity gradient tensor formulation, when observation sites approach the causative body, that is O(R −3 × R 2 ) = O(R −1 ).
Designing gravity gradiometers for measuring gravity gradient tensor signals is a difficult task. The gravity gradient tensor is a 3 × 3 symmetrical tensor with each entry being the second derivatives of the gravitational potential (or the first derivatives of the gravity field or the gravity acceleration). In terms of differential measurements, accelerations of at least two spatially separated masses were generally used to estimate the components of the gravity gradient tensor. Depending on the accuracies of the available gradiometers, gravity gradient tensor signals can be applied in geodesy (0.01 Eötvös to 0.1 Eötvös), autonomous navigation (0.1 Eötvös to 1 Eötvös) and oilfield and mineral exploration geophysics (1 Eötvös to 10 Eötvös) (Evstifeev 2017). In exploration geophysics, the cause of the increased interest in gravity gradiometers is that, compared to gravity anomalies measured by gravimeters, gravity gradient tensor signals contain more detailed information about the subsurface structures (Chapin 1998;Nabighian et al. 2005). Furthermore, compared to gravity field data, gravity gradient tensor anomaly maps generally provide more contrasting 1 3 and clearer edge delineations, such as those arising from salt domes in oil and gas prospecting (Pedersen and Rasmussen 1990). Up to now, several gravity gradiometers have been adopted in realistic exploration geophysical problems, such as the 3D FTG (full tensor gradiometer) system by Bell Geospace (Bell and Hansen 1998;Bell et al. 1997;Brewster 2016;Abtahi et al. 2016) and the Falcon AGG airborne gravity gradiometer by BHP Billiton (Australia) (Lee 2001).
To invert measured gravity gradient tensor (GGT) signals, we need an accurate forward solver. In general, GGT forward modelling can be performed either numerically or analytically. Using numerical methods such as Gaussian quadrature approaches (Talwani and Ewing 1960), Fourier domain methods (Parker 1973;Wu and Chen 2016), finite-element methods and finite-difference methods (Cai and Wang 2005;Farquharson and Mosher 2009;Jahandari and Farquharson 2013), the GGT signals caused by a mass body can be straightforwardly evaluated. However, the accuracies of numerical GGT signals can be seriously reduced by these numerical methods. For instance, using finite-element and finite-difference methods, improper translation of the computed gravitational potential to its second-order derivatives can lead to serious numerical errors. These undesired precision loss issues can be completely avoided by analytic approaches, which offer highly accurate GGT signals in terms of exact solutions.
Due to early limitations in the considered model geometries, the simple rectangular prismatic element has been widely accepted as a basic element in gravity forward modelling. Using these prismatic elements, the mass distributions in the Earth were simply approximated conceding a certain loss of geometrical accuracy. As for this basic element, several exact solutions of the gravity gradient tensor signals were developed (Forsberg 1984;Li and Chouteau 1998;Montana et al. 1992;Nagy and Papp 2000;Holstein et al. 2013;De Stefano and Panepinto 2016). However, it is difficult to accurately approximate a complex body using prisms. Polyhedral elements are endowed with great flexibility in presentation of 3D mass sources with complex geometries. Compared to prismatic elements, discretisations in terms of polyhedrons need lower numbers of elements to represent complicated mass bodies (Petrović 1996). Several analytical solutions have been successfully derived for the gravity gradient tensor signals of a homogeneous polyhedral body (Okabe 1979;Götze and Lahmeyer 1988;Kwok 1991;Petrović 1996;Werner and Scheeres 1996;Tsoulis and Petrović 2001;Holstein 2002;Tsoulis 2012;Holstein et al. 2013;D'Urso 2014a). Closed-form solutions of the gravity gradient tensor were also derived for other homogeneous simple geometries, such as pyramids (Sastry and Gokula 2016) and cylinders (Rim and Li 2016).
Seeking simplicity, geophysicists often assume that the Earth is composed of 3D anomalies in a layered medium or a succession of strata with horizontally undulating interfaces (e.g., sedimentary basins and underlying bedrock). In each layer, the rock mass density predominantly exhibits depth-dependent variations. To approximate the density variations in sedimentary basins, several closed-form solutions of gravity signals were derived for depth-dependent mass contrast functions, such as exponential functions (Cordell 1973;Chai and Hinze 1988;Chappell and Kusznir 2008), hyperbolic functions (Litinsky 1989;Rao et al. 1995), parabolic functions (Chakravarthi et al. 2002), quadratic polynomials (Rao 1985(Rao , 1990Gallardo-Delgado et al. 2003) and cubic polynomials (García-Abdeslem 2005). Compared to these depthdependent functions, polynomials offer a more flexible way to approximate arbitrarily variable density distributions. Very recently, Jiang et al. (2017) have given a detailed performance comparison between polynomial mass contrast functions and depth-dependent mass contrast functions for the capability of approximating complicated density distributions in the Earth. The result clearly demonstrates the superiority of polynomial mass contrast functions over the 1 3 depth-dependent mass contrast functions. As the real Earth has complicated mass density distributions, it is important to deal with a general polynomial mass contrast function, which not only varies in depth (vertical direction), but also varies in the horizontal direction.
For general polynomial functions, Pohánka (1998), Holstein (2003, Hansen (1999), D'Urso (2014b) and Ren et al. (2017a) have successfully derived closed-form solutions for the gravity field of a polyhedral body, in which the mass contrast linearly varies in both horizontal and vertical directions. Quite recently, analytical expressions for the gravity field of a polyhedral body with cubic polynomial density contrast in both horizontal and vertical directions were derived by D'Urso and Trotta (2017) and Ren et al. (2018). Nevertheless, only Holstein (2003) and D'Urso (2014b) have successfully derived closed-form solutions for gravity gradient tensor of a polyhedral body, in which the mass contrast varies linearly in both horizontal and vertical directions. To the best of our knowledge, closed-form solutions for gravity gradient tensor signals of a polygonal body with high-order polynomial mass contrasts (such as quadratic and cubic orders) varying in both horizontal and vertical directions have not been previously reported.
To overcome this both theoretically and practically important issue, we have successfully derived closed-form solutions for modelling the gravity gradient tensor of the above cases in this study. Inside the polyhedral body, its density contrast is approximated by a general polynomial function (currently, a maximum order up to and including three is considered). The polynomial mass contrast function can vary in both horizontal and vertical directions. Therefore, it has the capability to approximate complicated mass contrasts in realistic Earth models. To begin with, we apply the Gaussian gradient theorem to transform the volume integrals of GGT signals into surface integrals over the faces of polyhedra. Then, we apply the surface divergence theorem to further transform these surface integrals into a sequence of line integrals along the edges of the polyhedra. In the process of reducing the integral dimension, several vector or dyadic identities are employed to reduce the orders of the density polynomials and of the singularities in R in the integrands. The line integrals along the edges of the polyhedra are evaluated in terms of known analytical expressions (e.g., Gradshteyn and Ryzhik 2007).
To verify the accuracies of our new analytical solutions, a right rectangular prism with depth-dependent density contrasts, an irregular polyhedron with linear density contrast and a tetrahedral element with horizontally and vertically varying density contrast are tested. For the prism model, its reference solution for the gravity gradient tensor is computed from the derivatives of García-Abdeslem's (2005) analytical solution for the gravity field. The analytical solution derived by Holstein (2003) is used as reference for the irregular polyhedral model. Results from high-order Gaussian quadrature are used as reference solutions for the tetrahedron model.

Theory
For an arbitrary interior point of the polyhedral body (see Fig. 1, Appendix A and Table 5), let ( ) denote the density contrast in the polyhedral mass target at that point. The difference between two points denotes a vector, for example, the vector from a source point to the observation point ′ is denoted by − � . The gravitational acceleration field at the observation point ′ is (Ren et al. 2018) (1) where G = 6.673 × 10 −11 m 3 kg −1 s −2 is Newton's gravitational constant, ∇ is the gradient operator on point , R = | − � | is the distance from the observation site to the running integral point in the source polyhedral body H, and the gradient theorem (Tai 1997) is applied. Using Eq. (1) for the gravity field, the gravity gradient tensor can be derived as follows: where ∇ � denotes the gradient operator on point ′ and the relation ∇ 1 R = −∇ � 1 R is used. Using the following vector identity (Tai 1997): where can be a scalar or a vector, setting = 1 R and = ∇ and using the gradient theorem (Tai 1997), the volume integral term ∭ H (∇ 1 R )(∇ )dv in Eq.
(2) can be further transformed as  Table 5 r (2) can be transformed as where the surface gradient theorem (Tai 1997) is used. The symbol ∇ s denotes the surface gradient operator, where the subscript s means it is an operator on a 2D surface. In the above, the definition (A.3) of the number h i , constant over facet plane i, as the projection onto the facet normal i of the position vector of a planar point relative to the observation point ′ and the following vector identity were used: where f is an arbitrary scalar function and ̂ is the outward pointing normal unit vector of the considered surface. In Eq. (5), we set f = R . Substituting Eqs. (4) and (5) into Eq. (2), we obtain the following expression for the gravity gradient tensor of a polyhedral body H with arbitrary mass contrast function ( ): The density function is defined as a general polynomial function. It allows for the density variations in both horizontal and vertical directions, which is defined as: where P is the maximum polynomial order, and d are the d-th order density terms. For instance, the zero (constant), first (linear), quadratic and cubic order density terms are given as: where a klm are known coefficients which can be estimated by fitting measured gravity field data (Blakely 1996). The total gravity signal is the sum of the individual contributions from different mass densities, such as 0 , 1 , 2 and 3 . Substituting the density contrast in Eqs. (8) into (7), we get the final gravity gradient tensor: where d denotes the individual contribution from the d-th order density contrast. In this study, we only consider the cases 0 ≤ P ≤ 3 . For simplicity, we denote x, y and z by x p (p = 1, 2, 3) in the following sections, that is

Constant Density Contrast
For the constant term, substituting 0 = a 000 and ∇ 0 = 0 into Eq. (7), we get 0 as follows: Detailed derivations of expressions for the line integral ∫ C ij 1 R dl and the surface integral ∬ H i 1 R 3 ds are given in Appendix B.

Linear Density Contrast
As for the linear terms in the polynomial density contrast, we set 1 = ⋅ , where = (a 100 , a 010 , a 001 ) and = (x, y, z) . Substituting = 1 = ⋅ and ∇ = ∇ 1 = and ∇ ∇ 1 = into Eq. (7), we get The above equation consists of one kind of line integral, i.e., ∫ C ij R dl and two kinds of surface integrals, i.e., ∬ H i R 3 ds and ∬ H i 1 R ds , with their detailed derivations given in Appendix C.

Quadratic Density Contrast
As for the quadratic term in the polynomial density contrast in Eq. (11), its gradient is The second derivative of 2 is Substituting 2 and its derivatives ∇ 2 and ∇ ∇ 2 into the general expression of the gravity gradient tensor in Eq. (7), four kinds of integrals need to be considered (their detailed derivations are given in Appendix D), which are the line integral term ∫ C ij where, using p, q = 1, 2, 3 , the notations in Eq. (14) are adopted, and the product x p x q represents an arbitrary quadratic monomial.

Cubic Density Contrast
Finally, we discuss the gravity gradient tensor due to the cubic density terms given in Eq. (12). For the first-and second-order derivatives in Eq. (7), we have and Substituting 3 , ∇ 3 and ∇ ∇ 3 into the expression of the gravity gradient tensor in Eq. (7), we find that we only need to deal with four types of integrals (their detailed derivations are given in Appendix E) to get the final closed-form solutions of 3 in Eq. (13). These four types of integrals are the line integral ∫ C ij Here, the notations in Eq. (14) are used, such that x p x q and x p x q x t represent arbitrary quadratic and cubic monomials, respectively.

Comparison with Other Solutions
The previously published analytical solutions of gravity gradient tensors for different geometries and mass density contrasts are listed in Table 1. Most of these solutions were designed for rectangular prisms with constant density contrasts (Forsberg 1984;Montana et al. 1992;Li and Chouteau 1998;Nagy and Papp 2000;Rim and Li 2016) (19) ∇ 3 ( ) = a 102 z 2 + a 111 yz + a 120 y 2 + 2a 201 xz + 2a 210 xy + 3a 300 x 2 ̂ + a 012 z 2 + 2a 021 yz + 3a 030 y 2 + a 111 xz + 2a 120 xy + a 210 x 2 ̂ + 3a 003 z 2 + 2a 012 yz + a 021 y 2 + 2a 102 xz + a 111 xy + a 201 x 2 ̂ (20)  (Okabe 1979;Götze and Lahmeyer 1988;Kwok 1991;Werner and Scheeres 1996;Holstein 2002;Holstein et al. 2007a, b;Tsoulis and Petrović 2001;Tsoulis 2012;D'Urso 2014a). Only a few studies have been carried out for linear mass density contrasts (Holstein 2003; D'Urso 2014b; Sastry and Gokula 2016). Using our new findings, we successfully derive a set of analytical expressions of the gravity gradient tensor signals for a general polyhedron with constant, linear, quadratic and cubic polynomial mass functions. Our closed-form solutions for constant and linear polynomial mass functions may have different forms from the previously published solutions, but they should produce the same results. In addition, we were not only the first to find the exact solutions for the quadratic and cubic cases, but also our solutions allow for mass contrasts varying simultaneously in both horizontal and vertical directions. However, it is difficult to derive closed-form solutions of the GGT for higher-order polynomial density contrasts, simply because analytical expressions cannot be found for 1D edge integrals in terms of the existing integral tables (Gradshteyn and Ryzhik 2007).

Verification
Three models were used to validate our closed-form solutions. The first one is a rectangular prism model (Fig. 2), for which the derivatives of closed-form solutions for the gravity field in García-Abdeslem (2005) were taken as references. The second model is an irregular polyhedron with linear density contrast, taken from Holstein (2003), for which a comparison with Holstein's (2003) solution is given. The third model is a relatively complicated tetrahedral body (Fig. 5). As there are no reference solutions for the GGT of polyhedral mass bodies with quadratic and cubic density contrasts, we have used high-order Gaussian quadrature rules (such as 512 × 512 × 512 = 124,217,728 quadrature points) to calculate GGT reference solutions. Additionally, we allowed the density in the tetrahedral body to vary in both horizontal and vertical directions. We should mention that because our formulae require the observation site to be located at the origin of the Cartesian coordinate system, for each observation site, a coordinate translation must be performed to move the observation site to the origin.

A Prismatic Body with Depth-Dependent Density Contrast
The prism is located in a Cartesian coordinate system with the z-axis downward. In Fig. 2a, the coordinate ranges of the prismatic body are x = [10, 20] km, y = [10, 20] km and z = [0, 8] km. The density function is a depth-dependent cubic polynomial which is taken from García-Abdeslem (2005): where the density contrast is in kg∕m 3 and z is in km.
In the García-Abdeslem's (2005) work, an analytic formula for the vertical gravity field caused by a rectangular prism with density contrast varying as a depth-dependent cubic polynomial was derived. Therefore, we took the derivatives of García-Abdeslem's solution (i.e., T xz , T yz , T zz ) as references to verify our solution for the gravity gradient tensor. Furthermore, to verify our solution for an arbitrary polyhedral body, we decomposed the original rectangular prism into two polyhedra (as shown in Fig. 2b). These two polyhedra share the plane with vertices 1, 3, 8 and 9. Note, this test requires points 1, 3, 8 and 9 to be co-planar. If they were not co-planar, the surface (1, 3, 8, 9) would need to be divided into two triangles (1, 3, 8) and (3,8,9). The values of T xz , T yz and T zz were computed at two observation sites. One is located at point (12, 12, − 0.001) km with a 1-m offset right above the top face of the prismatic body, whereas the other one is located at point (20, 10, − 0.001) km, that is near a corner of the rectangular prism.
The results computed by our closed-form solution and by derivatives of García-Abdeslem's (2005) gravity field solution are compared in Table 2. We compare GGT (21) ( ) = − 747.7 + 203.435z − 26.764z 2 + 1.4247z 3 , Table 2 Comparison of gravity gradient tensors calculated by our new closed-form solution and derivatives of García-Abdeslem's (2005) gravity field solution for the prismatic body given in Fig. 2 The relative errors of the tensor trace are calculated as |T xx + T yy + T zz |∕(|T xx | + |T yy | + |T zz |) . Symbol (−) indicates no solution available

Observation sites (km) Component
Our solution (s −2 ) Derivatives of García-Abdeslem's solution (s −2 ) Above the top face (12, 12, -0.001) Above a corner(20, 10, − 0.001) 3.77463390064588E-13 -elements T xz , T yz , T zz , because only g z is given in García-Abdeslem (2005). Excellent agreement is obtained between these two approaches, with relative errors on the order of 10 −13 % at site (12, 12, − 0.001) km and relative errors 10 −6 -10 −7 % at site (20, 10, − 0.001) km. The tensor trace T xx + T yy + T zz = ∇ 2 � U is a useful tool for indicating numerical error. Using the Poisson equation ∇ 2 � U = −4 G for the gravitational potential U, the tensor trace vanishes outside the source body and is −4 G inside the source body. Since our observation points are outside the source body, we have calculated the relative errors shown in Table 2 as |T xx + T yy + T zz |∕(|T xx | + |T yy | + |T zz |) . The relative error of the tensor trace is 2.74 × 10 −13 % at the location above the top face and 3.77 × 10 −11 % above the corner. Both our solution and García-Abdeslem's (2005) solution for the GGT are singular, when observation sites are located on edges and corners. When the observation site gets close to the corner, the mathematical singularity becomes stronger, and this is, why the relative error is larger at site (20, 10, − 0.001) km.
In summary, the excellent agreement of the different solutions has successfully verified the accuracy of our new closed-form solutions for a polyhedral body with depth-dependent polynomials up to and including cubic order.

An Irregular Polyhedron with Linear Density Contrast
An irregular polyhedron is tested to compare our solution with Holstein's (2003) solution for linear media. The target model is composed of 8 faces and 10 points (Fig. 3), originally designed by Holstein et al. (1999, Appendix A). The linear density contrast inside the target body is taken from Holstein (2003, Appendix C), that is where the density is in kg∕m 3 , coordinates are given in km. First, we compute the gravity gradient tensor at two test points, an outer point with coordinates of (0, 0, 0) km and an interior point at the centroid ̄ = 1 88 (40, 250, −1541) km. The results are shown in Table 3. An excellent agreement is obtained between our closed-form solution and Holstein (2003)'s solution.
Second, the gravity gradient anomalies are evaluated on a vertical profile passing through the target body, in order to validate our formulas by the well-known Poisson equation ∇ 2 � U = −4 G , where U is the gravitational potential. The vertical profile is from a point with coordinates of (0, 0, − 34) km to a point with coordinates of (0, 0, 0) km, with observation points placed at a uniform vertical (z) spacing of 1 km. Due to the discontinuity of the gradient tensor signals on the boundary, we only perform tests at two points which are very close to the boundary (with a small distance of ± 10 −10 km). The value of  Fig. 4 a Comparison of the Laplacian term T xx + T yy + T zz calculated by our new closed-form solution to the reference value of −4 G for the causative body in Fig. 3. b Relative errors calculated as −4 G is used as the reference. The relative errors are calculated using the formula (T xx +T yy +T zz )+4 G |T xx |+|T yy |+|T zz | and shown in Fig. 4. The maximum absolute relative error is 3.08389 × 10 −13 %.

A Tetrahedral Mass Body with Horizontal and Vertical Density Contrasts
In practical gravity exploration, the underground mass bodies can have complicated shapes. To calculate their GGT signals, we generally need to discretise the underground mass bodies into sets of disjoint elements with different shapes, such as structured hexahedral or prismatic elements and unstructured tetrahedral elements. Compared to regular prismatic elements, unstructured tetrahedral elements can well approximate arbitrarily complicated anomalies. Using recent Delaunay triangulation techniques (e.g., Si 2015), geophysicists can easily set up discretised triangulated grids to represent complicated anomalous targets. Therefore, unstructured grid techniques have been widely used in the geophysical community, not only in gravity exploration (Jahandari and Farquharson 2013;Ren et al. 2017c), but also in the electromagnetic induction community (Li and Key 2007;Schwarzbach et al. 2011;Ren et al. 2013) and the seismic imaging field (Lelièvre et al. 2011). Testing a single tetrahedral mass body not only aims to test the performance of our new closed-from solutions for complicated density contrasts (horizontal and vertical density contrasts), but also demonstrates its compatibility with other codes developed for unstructured grids. This verification is an important step towards joint inversion for multiple parameters, such as magnetisation vector, conductivity, velocity and mass density on the same unstructured gird. The geometry of the tetrahedron is shown in Fig. 5. We used the density contrast as given in D'Urso and Trotta (2017, equation 166), which includes horizontal variation of the density contrast where the unit of the density is kg /m 3 and the units of the coordinates are km. A measuring plane is located above the tetrahedron in a range of x = [− 0.16, 0.16] km, y = [− 0.16, 0.16] km and z = 0km. On the measuring plane, the gravity gradient tensors were computed using our closed-form solution and a Gaussian quadrature rule with 512 × 512 × 512 points. The high-order Gaussian quadrature rule was built by projecting the tetrahedral element into a hexahedral element (Rathod et al. 2006) and then applying the standard Gauss Legendre Quadrature rule (Golub and Welsch 1969). Since there are no published results of GGTs for high-order polynomial density contrasts, we compare to results computed by high-order Gaussian quadrature as the reference solutions. The six components of the gravity gradient tensor are shown in Fig. 6 for our analytical solutions and in Fig. 7 for Gaussian quadrature solutions. The relative errors with regard to the Gaussian quadrature solution are shown in Fig. 8. Clearly, the results from our analytical solution agree quite well with those computed by the Gaussian quadrature rule, with absolute relative errors less than 2 × 10 −7 % . The computation time for the analytical solution is about 0.9 s, and the computation time for 512 × 512 × 512 points Gaussian quadrature is about 3.7194 h.
Furthermore, the sums of the diagonal entries T xx + T yy + T zz were computed using the analytical solution. When the observation site is outside the source region, as in our example, these three diagonal entries satisfy the Laplace equation, that is ∇ 2 U = T xx + T yy + T zz = 0 , thus providing an independent means of verification of our closed-form solutions. As shown in Fig. 9, the relative error |T xx + T yy + T zz |∕(|T xx | + |T yy | + |T zz |) has a maximum value of 3.26 × 10 −12 % . Therefore, the condition of ∇ 2 U = T xx + T yy + T zz = 0 is satisfied by our solutions. Next, with the aim of testing the case of cross-varying mass contrasts, the following density contrast and the same tetrahedral body as shown in Fig. 5 were used: where the unit of density is kg /m 3 and the units of the coordinates are km. This density contrast function simultaneously varies in both horizontal and vertical directions, ranging from about −300 to 400 kg/m 3 . A measurement profile is positioned at x = [−0.2, 0.2] km, y = 0 km and z = 0 km. The six components of the gravity gradient tensor were calculated by our analytical formulae and Gaussian quadrature with 512 × 512 × 512 points as shown in Fig. 10. The relative errors of our solutions with respect to the Gaussian quadrature solutions are shown in Fig. 11. Note, that the maximum absolute relative errors are less than 5 × 10 −8 %. Therefore, the results computed by our analytical solution have excellent agreement with those computed by high-order Gaussian quadrature, for this case where the density contrast simultaneously varies in both horizontal and vertical directions. Furthermore, a test was conducted with observations sites located on a path at x = [−0.05, 0.05] km, y = 0 km, z = 0.04 km crossing the interior. For this setup, the tensor elements and traces T xx + T yy + T zz were calculated. An excellent agreement between the computed traces of the gravity gradient tensor and the reference value of −4 G (see Fig. 12) is obtained, which clearly verifies our closed-form solution.

Numerical Stability
Theoretically, closed-form solutions should be accurate. However, due to the limitations of floating point arithmetics, these solutions contain inevitable rounding errors, when the amplitudes of the gravity signals approach zero. This phenomenon has been observed in previous studies (Holstein and Ketteridge 1996;Holstein 2003;Holstein et al. 2007a), which shows that the relative error of the gravity field will grow, when the distance between the mass target and the observation site approaches infinity. We use relative error here to avoid discussion of anomaly magnitude and instrument-dependent measurement accuracy. Beyond a certain distance, gravity anomalies can even be totally corrupted by floating point errors. Assuming that the mass body has a size of and the site-to-target distance is denoted by , the growth of the relative error satisfies the following formula (Holstein and Ketteridge 1996;Holstein 2003;Holstein et al. 2007a): where is defined as = ∕ , −1 is referred to as the dimensionless target distance, the exponent indicates the speed of error growth, and is the relative error when 1 = 1 , identified in these papers as the floating point machine precision constant. On a log-log scale plot, Eq. (25) is a linear error growth curve lg( ) = lg 1 + lg( ). We present an experiment, in which the tetrahedron shown in Fig. 5 was used to recover the above linear error growth curve. The size of the tetrahedron was set to the diameter of its circumscribed sphere, that is = 0.134 km. The starting point for the profile is located at the centre of the circumscribed sphere of the tetrahedral model. The measuring profile is formulated as x = y, z = 0 km.
Following Holstein et al. (2007a), the equivalent mass point approach was adopted as the reference solution, by which the tetrahedral body is approximated by a point source located at the centre of the circumscribed sphere. The anomalies calculated by our new closed-form expressions and by the equivalent point-mass approach are denoted by T and T point , respectively. In addition, the relative errors were computed. Four polynomial density contrasts of different orders were used for the stability test: where the unit of the density contrast is kg/m 3 and the units of the coordinates are km. The ranges of density contrasts are about 200-1800 kg/m 3 for the linear case, about 40-2200 kg/m 3 for the quadratic case and about 10-2250 kg/m 3 for the cubic case.
(26) 0 = 1000, (27) 1 = 10,000(x + y + z), (28) 2 = 100,000(x 2 + y 2 + z 2 + xz + yz + xy), (29) 3 = 1,000,000(z 3 + yz 2 + y 2 z + y 3 + xz 2 + xyz + xy 2 + x 2 z + x 2 y + x 3 ), The relative error curves are shown in Fig. 13 for the above four density contrasts. In Fig. 13, at low dimensionless target distances −1 , the relative errors decrease with increasing −1 , because the tetrahedron is better approximated as a mass point with increasing distance. However, the gravity anomalies are progressively corrupted with increasing target distance, and when the dimensionless target distance ( −1 ) is beyond a critical value, the rounding error in the unstable calculations exceeds the solution difference, resulting in a rising trend of the relative errors (cf. Holstein 2003). We reconstructed the entire error growth curves by extrapolating the ascending parts back to ( −1 = 1, = ) . In Fig. 13, the reconstructed curves are denoted by blue dashed lines. We observe that there are two linear error growth trends for each polynomial density contrast, that is T zx and T zy share the same error growth behaviour, and T xx , T yy , T zz and T yx share the same error growth behaviour. The estimated values of are shown in Table 4. Using the estimated , the critical dimensionless target distances were calculated from that value of −1 where the rising part of the error curve indicates a relative error of 1 (Holstein and Ketteridge 1996;Holstein 2003;Zhou 2010) in Eq. (25), which are also shown in Table 4.
As shown in Table 4, for the linear density case, is about 4.31 for T xx , T yy , T zz and T xy and about 5.28 for T xz and T yz , which agrees with previously estimated values of 4 and 5 given by Holstein (2003), respectively. As indicates the error growth speed, Table 4 and Fig. 13 also show that the higher the order of the density contrast polynomial is, the faster the errors accumulate. At the critical dimensionless target distance, the accumulated rounding error has the same amplitude as the true anomaly. Therefore, to safely use our exact solution, we should require our observation sites to be located in a range, which is less than the critical value 1 crit . Finally, interested readers are referred to Table 4 for the specific dimensionless site-target distances for the constant, linear, quadratic and cubic density cases.

Conclusions
In the presented work, we have derived a set of closed-form solutions for the gravity gradient tensor of an arbitrary polyhedral body. The density contrast of the polyhedral body is represented as a polynomial function of up to and including third order. The polynomial function allows density contrasts to simultaneously vary in both horizontal and vertical directions. To our best knowledge, this is the first time that analytical solutions of the gravity gradient tensor are derived for an arbitrary polyhedral body with polynomial orders up to three. Three synthetic models (a prismatic body, an irregular polyhedron and a tetrahedral body) were used to test the correctness and the efficiency of our newly developed closed-form solutions. By comparing to published closed-form solutions and high-order Gaussian quadrature solutions, the high accuracies of our solutions with deviations of less than 2 × 10 −7 % from the Gaussian quadrature solutions were demonstrated. The computation time used by our analytical solution is significantly less than that of the high-order Gaussian quadrature. The numerical stability tests show that, when dealing with cubic density contrast, our closed-form solutions would generate inaccurate solutions, if the dimensionless target distance, which is defined as the ratio of the distance between the observation site and the causative body to the dimension of the causative body, is larger than about 282 for T xx , T yy , T zz , T yx components, and 137 for T zx , T zy components, on a profile with x = y and z = 0 m . This inaccuracy problem is caused by the limited precision of floating point operations, when the amplitudes of the gravity gradient tensors are very close to zero. However, since most applications of gravity gradient tensors aim at detecting anomalous bodies in the very shallow subsurface of the Earth, this problem of precision loss should not cause a serious problem in practical situations, and our closed-form solutions are very safe to be applied in exploration geophysics.
For the cases of constant, linear, second and cubic polynomial order, there are systematic recurrences of previously evaluated integrals as well as occurrences of new types of integrals, when going to higher polynomial order. We assume that this is a trend even for higher polynomial orders ( P > 3 ). However, extending our work to higher polynomial orders will be a nontrivial task as more and more complicated integrals need to be considered. By contrast, using the analytical expressions presented in this study, closed-form solutions for the magnetic field can be derived easily.
Furthermore, the outward pointing unit vector ij of edge C ij and face H i is calculated as: Let us establish a local Cartesian coordinate system, the origin of which is located at the observation point ′ , that is � = (0, 0, 0) . Symbol = (x, y, z) represents an arbitrary source point in the polyhedron, and R = | − � | is the distance from a source point to the observation point. The three unit vectors ̂ i , ̂ ij and ̂ ij form a natural orthonormal basis on edge C ij so that projections of vector ( − � ) along these three unit vectors yield a set of local coordinates (m ij , s ij , h i ) on edge C ij , which are calculated as: In Fig. 1, point ij is the projection of ′ onto edge C ij . Vector ij = ij − � is pointing from the observation site ′ to point ij with its magnitude being denoted by L 0ij . Since − � = ij + ( − ij ) and ij ⋅̂ ij = 0 , the 1D parametrised local coordinate s ij along edge Using the fact that � = (0, 0, 0) , the above equation becomes: Meanwhile, the distance R = | − � | , from any source point on edge C ij to the observation site, is parametrised as The distances from ′ to the vertices 0ij and 1ij are defined as where s 0ij and s 1ij are the parametrised coordinates of the vertices 0ij and 1ij , respectively. In addition, in Fig. 1, i is the projection of the observation site ′ on the i-th face H i , and ̂⊥ ij = ij − i | ij − i | is the unit vector which points from point i to point ij . The direction of ̂⊥ ij can be either identical or opposite to the direction of the outward normal vector ij on edge C ij . i = | − i | is the distance between the projection centre i and a source point ∈ H i .

Appendix B: Gravity Gradient Tensor due to Constant Density Contrast
Using the local coordinate s ij = ( − � ) ⋅̂ ij along edge C ij and using Eq. (A.6), the line integral in (15) can be converted into a definite integral: where ds ij is a differential of variable s ij . The values of s ij at the vertices 0ij and 1ij are s 0ij and s 1ij , respectively. The analytical solution of the above definite integral can be looked up in integral tables (equation 2.261 in Gradshteyn and Ryzhik 2007): When the observation site lies on the extension of edge C ij ( L 0ij = 0 and s 1ij ⋅ s 0ij > 0 ), the integrand function is reduced to 1 |s ij | which leads to the second formula of Eq. (B.2). We should note that when the observation site ′ is located within the edge C ij ( L 0ij = 0 and s 1ij ⋅ s 0ij ≤ 0 ), the linear integral term ∫ C ij 1 R dl has a weak logarithmic singularity. Therefore, we cannot compute gravity gradient tensors, when the observation site is located on an edge of the polyhedral body.
The analytic expression for the surface integral term ∬ H i 1 R 3 ds in Eq. (15) can be derived using the result from Ylä-Oijala and Taskinen (2003) and Ren et al. (2017b) where where the variables m ij , s 1ij , s 0ij , L 0ij , h i , R 1ij and R 0ij are shown in Fig. 1 and all can be calculated in terms of the coordinates of the vertices of the polyhedral body. In Eq. (B.3), the observation site cannot be located on the surface H i of the polyhedral body. In the case that the observation site is located on the surface of the polyhedral body, � ∈ H i , the gravity gradient tensor is singular, which is an inherent feature of gravity gradient tensor (Li and Chouteau 1998;Holstein 2003).

Appendix C: Gravity Gradient Tensor due to Linear Density Contrast
First, we consider the line integral term ∫ C ij R dl in Eq. (16). Using the parametrisation in Second, we deal with the surface integral term ∬ H i R 3 ds . Using the identity ∇ 1 R = � − R 3 , our previous assumption � = (0, 0, 0) , and Eq. (6) (setting f = 1 R ), we have where ∇ s denotes the surface gradient operator, and the surface gradient theorem is used (Tai 1997 where as shown in Fig. 1, i is the projection point of the observation site onto the plane containing the face H i , i = | | − i | | , is a running integral point on edge C ij , i.e., ∈ C ij . Using the geometrical variables given in Fig. 1, we have (Ren et al. 2017a): 14 Illustration of the angular extent ij ( i ) subtended by edge C ij of polygon H i and the total angular extent ( i ) = 2 , when i ∈ H i When m ij → 0 , the limit of Eq. (C.5) exists and is equal to zero. Therefore, there is no singularity in Eq. (C.5). As demonstrated in Fig. 14, ( i ) is the angular extent of the arc region lying within a plane containing surface H i and centred at i . Note that the projection point i can be located inside, outside or on edges of the surface H i . The value of the extent angle ( i ) depends on the geometrical relation between i and the polygon H i . When point i is located inside H i , ( i ) = 2 ; when point i is located on an edge of the polygon H i but not at a vertex, ( i ) = ; when point i is located at a vertex, ( i ) is the angle enclosed by its two adjacent edges. To avoid judgement of the geometrical relation between point i and the polygon H i and facilitate programming, the angle ( i ) can be expressed as a sum of the angles subtended by each edge of the polygon H i (Wilton et al. 1984, page 278): where the contribution of each edge is calculated as to the angular extent ( i ) vanishes, when point i is on the edge C ij or on its extension, and otherwise the sign of ij depends on whether unit vector ̂⊥ ij is identical or opposite to unit normal vector ̂ ij . Finally, the closed-form solution for gravity gradient tensor 1 in Eq. (13) caused by a linear density contrast can be obtained by substituting Eqs. (C.4), (C.3) and (C.1) into Eq. (16).

Appendix D: Gravity Gradient Tensor due to Quadratic Density Contrast
First, using Eqs. (A.5) and (14), the line integral ∫ C ij Using the integral x p x q R 3 ds, (p, q = 1, 2, 3 ). In view of the assumption that � = (0, 0, 0) , we have Using the above equation and the following divergence vector identity: x p x q R 3 can be transformed as (setting = x p̂ q and = 1 R ): Integrating the above equation over the face H i and applying the surface divergence theorem (Tai 1997 where ∫ C ij x p R dl =̂ p ⋅ ∫ C ij R dl and ∬ H i x p R 3 ds =̂ p ⋅ ∬ H i R 3 ds . The analytic solutions for the linear integral ∫ C ij R dl and for the surface integral ∬ H i R 3 ds have already been given in Eqs. (C.1) and in (C.3), respectively. Furthermore, the surface integral term ∬ H i 1 R ds has been derived in Eq. (C.4).
Third, the surface integrals ∬ H i x p R ds for (p = 1, 2, 3) are the three components of the vector surface integral ∬ H i R ds , that is ∬ H i x p R ds =̂ p ⋅ ∬ H i R ds . Using our assumption � = (0, 0, 0) , Eq. (6) (setting f = R ) and the surface gradient theorem (Tai 1997), the integral ∬ H i R ds can be calculated as (Ren et al. (2017a) where the analytic solution for the integral term ∬ H i 1 R ds has been given in Eq. (C.4).

Appendix E: Gravity Gradient Tensor due to Cubic Density Contrast
First, we deal with the line integral ∫ C ij x p x q x t R dl . In terms of the definition of the local coordinate s ij , which is given in Eq. (A.5), we obtain Integrating the above equation over the face H i and applying the surface divergence theorem (Tai 1997) to the first term on the right-hand side, we get where the line integral ∫ C ij x p x q R dl can be analytically evaluated using Eq. (D. x p x q x t The surface integral ∬ H i Rds in Eq. (E.5E.6) can be calculated using the result of Ren et al. (2017a, equation 33). This means where the calculation of the angular extent ( i ) has been presented in Eq. (C.6), and When m ij approaches zero, the last two terms in Eq. (E.9) approach zero. Fourth, volume integral ∭ H x p R dv can be calculated from Ren et al. (2018, equation 24) where the surface integral ∬ H i Rds is given in Eq. (E.8).