Exact solutions of the vertical gravitational anomaly for a polyhedral prism with vertical polynomial density contrast of arbitrary orders

SUMMARY

We present general closed-form solutions for the vertical gravitational anomaly caused by a polyhedral prism with mass density contrast varying with depth. Our equations are the first ones to implement a polynomial vertical mass density contrast of arbitrary order. Singularities in the gravity field which arise when the observation site is close to or in the anomalous polyhedral prism are removed in our analytic expressions. Therefore, the observation site can be located outside, on the faces of or inside the anomalous mass bodies. A simple prismatic body of anomalous density is adopted to test the accuracy of our newly developed closed-form solution. Cases of constant, linear, quadratic, cubic and quartic polynomial orders of mass density contrast are tested. For cases of constant, linear, quadratic and cubic polynomial orders, the relative errors between our results and other published exact solutions are less than 10⁻¹¹%. For the case of quartic polynomial order, relative errors less than 10⁻¹⁰% are obtained between our solutions and those computed by a high-order Gaussian quadrature rule (512 × 512 |$\times 512=134\, 217\, 728$| quadrature points), where our new analytic solution needs significantly less computational time (0.0009  versus 31.106 s). These numerical experiments not only verified the accuracy of our new formula but also demonstrated their potential in computing exact gravity anomalies for complicated mass density distributions in the Earth.

Geopotential theory, Gravity anomalies and Earth structure, Numerical approximations and analysis, Numerical modelling, Numerical solutions

1 INTRODUCTION

Given a mass distribution, the calculation of its gravitational signals such as the gravity potential and gravity field is a core task in the study of the Earth’s gravity field (Blakely 1996; Li & Oldenburg 1998; Hofmann-Wellenhof & Moritz 2006). The gravitational signals of a mass body can be calculated by either solving Poisson-type partial differential equations or evaluating volume integral expressions. The former approach generally can employ the well-established finite-difference and finite-element methods (Zhang et al.2004; Cai & Wang 2005; Farquharson & Mosher 2009). These Poisson-equation approaches first calculate the gravity potential distribution over the entire domain which is composed of both mass bodies and free space. Then, the gravity fields at the observation sites are obtained by taking derivatives of the numerical gravity potential. Using standard linear shape functions for the gravity potential, the gravity fields obtained by finite-element methods are constant in each finite element and discontinuous across the element faces. Averaging these discontinuous gravity fields surrounding observation sites leads to large numerical errors in the final numerical solutions of the gravity field. Similar issues exist in finite-difference type methods. Therefore, elaborate strategies are required to obtain highly accurate gravity fields for these Poisson-equation approaches, such as p-type (adjusting orders of shape functions) adaptive refinement (Grayver & Kolev 2015) or h-type (adjusting size of elements) adaptive mesh refinement (Key & Weiss 2006; Franke et al.2007; Li & Key 2007; Schwarzbach et al.2011; Ren et al.2013; Ren & Tang 2014; Ren et al.2018a). In recent times, researchers have focused more on the latter approach of evaluating volume integrals. Different closed-form solutions could be developed for these volume integrals for the gravity potential and gravity field. These closed-form solutions offer high accuracy for numerical gravity potentials and their derivatives. They effectively avoid the accuracy loss encountered in the Poisson-equation approaches.

To derive closed-form solutions of gravity signals, only the underground mass anomaly needs to be discretized into disjoint elements. These elements can have different geometric shapes, such as circular disks and cylinders (Singh 1977; Krogh et al.1982; Damiata & Lee 2002; Rim & Li 2016; Asgharzadeh et al.2018), prisms and polyhedra (Li & Chouteau 1998; Wu & Chen 2016; Wu 2016, 2018; Jiang et al.2018; Ren et al.2018b). The prismatic body has been the focus of interest for an extended period of time due to its simplicity and high effectiveness. A large number of analytic expressions were developed for the prismatic mass body with different mass density contrasts such as polynomial density contrasts of constant, linear and quadratic orders or parabolic polynomial density contrast functions (Nagy 1966; Goodacre 1973; Banerjee & Das Gupta 1977; Rao 1990; Kwok 1991; Garcia-Abdeslem 1992; Rao & Babu 1993; Murthy & Swamy 1996; Li & Chouteau 1998; Chakravarthi et al.2001; Gallardo-Delgado et al.2003; Garcia-Abdeslem 2005; Swank 2006; Zhou 2009). Similarly, intensive studies have been conducted to find closed-form solutions for constant and linear density contrasts for the polyhedral body (Talwani & Ewing 1960; Paul 1974; Barnett 1976; Okabe 1979; Pohanka 1988, 1998; Holstein & Ketteridge 1996; Petrović 1996; Hansen 1999; Singh & Guptasarma 2001; Holstein 2002, 2003; Hamayun et al.2009; Tsoulis 2012; D’Urso 2013, 2014; Sudhakar et al.2014; Kim & Wessel 2016; Ren et al.2017a). In gravity problems, consideration of a polynomial variation of density in the vertical direction is particularly important, because the density of many geological bodies increases continuously with depth, for instance, in sedimentary basins. This is a direct result of lithostatic pressure increasing with depth.

Prisms and polyhedra are the fundamental mass elements used in the computation of the Earth’s gravity field. Exact solutions for prisms and polyhedra with low-order polynomial variations of mass density contrasts with depth (such as constant, linear and quadratic orders) have been widely used in both geophysics and geodesy (Blakely 1996; Li & Chouteau 1998; Hofmann-Wellenhof & Moritz 2006). Their success leads to an important but still unresolved theoretical question, which is to find exact solutions of the gravity field for a general density function that is continuous with depth. As demonstrated by Zhang & Jiang’s (2017) recent work for complicated underground mass anomalies, a general continuous density function can be better approximated by a polynomial density function (of arbitrarily high order) than by other depth-dependent density functions. The necessity to consider general density functions is also evident by the complicated physical and chemical processes of the Earth’s evolution (such as plate tectonics, volcano eruptions and earthquakes), which formed complicated mass density contrasts in the crust and mantle. In addition, from the gravity inversion point of view, we might prefer large discretized mass elements to represent the anomalous mass bodies, because this can significantly reduce the computational cost. Large elements can also be used to reduce the non-uniqueness of the inversion process (Li & Oldenburg 1998). However, for these cases with large discretized elements, low-order polynomial functions might be inappropriate. Therefore, high-order polynomial density contrasts might be more practical to approximate the complex mass distribution in the Earth.

Until now, only a limited number of studies has been carried out to find closed-form solutions of the gravity fields for prisms and polyhedra with general polynomial density contrasts. For instance, Zhang et al. (2001) derived analytical formulae for gravity anomalies caused by 2-D bodies with density contrast varying as polynomial functions of arbitrary orders. Garcia-Abdeslem (2005) formulated analytic expressions for gravity fields of a prismatic body with constant, linear, quadratic and cubic orders of polynomial density contrast. D’Urso (2015) derived exact closed-form equations for 2-D polygonal bodies with polynomial density contrasts of constant, linear, quadratic and cubic orders. More recently, similarly to the previously presented 2-D results (D’Urso 2015), D’Urso & Trotta (2017) and Ren et al. (2018c) independently considered the more difficult 3-D case, where exact solutions of gravity fields for polyhedral bodies with constant, linear, quadratic and cubic orders of polynomial density contrasts were found. In two recent contributions, Jiang et al. (2017) and Zhang & Jiang (2017) describe exact solutions for gravity fields of a prismatic mass body with arbitrary orders of polynomial density functions.

Now, the remaining natural task is to find exact gravity field solutions of a general polyhedral mass body with polynomial density contrasts of arbitrary orders. This task is a difficult and hard one to be accomplished. In this study, instead of considering the general polyhedral mass body, we chose the polyhedral prism which can be considered as a body between the simple prism and the complicated polyhedron. The polyhedral prism is a general prism with polygonal cross-section, which can be applied in gravity terrain corrections, gravity inversion and estimating the geometry of an isolated 3-D source (Plouff 1976; Kwok 1991; Murthy & Swamy 1996; Tsoulis 2003; Oliveira et al.2011; Oliveira & Barbosa 2013). After careful treatment, here we report our new findings about the exact solutions for both the gravity potential and gravity field of a polyhedral prism with vertical polynomial density contrast of arbitrary order. To derive the exact solution, the divergence theorem and gradient theorem were used leading to a new set of recursive equations. In these new recursive equations, the initial terms are simple surface integrals. Singularity-free closed-form solutions were derived for these initial terms. Therefore, the final recursive formulae for the gravity field are completely singularity-free. This means that the observation sites can be placed outside, inside and on the surface of the polyhedral prism.

The structure of this paper is as follows. First, the geometry of the polyhedral prism and its mass density contrast are defined. Second, the detailed procedures are given to derive the formulae for the gravitational signals, and its singularity-free properties are discussed. Then, prisms with constant, linear, quadratic and cubic orders of polynomial density contrast are used to verify the accuracy of our newly developed closed-form solution. For the case of quartic polynomial density contrast, high-order Gaussian quadrature rules are used to verify our solution. Excellent agreements between these solutions were obtained. Finally, the numerical stability of our new formula is tested by observing its accuracy reduction when the distance from the observation site to the mass target gradually raises. The stability tests show that our new algorithms can be applied safely in realistic exploration problems.

2 THEORY

2.1 Expressions for vertical gravitational anomaly

The prismatic mass body with polyhedral-cross section is defined in Fig. 1. The observation site |$\mathbf {r}^{\prime }$| is located at the origin of the right-handed Cartesian coordinate system, |$\mathbf {r}^{\prime }=(0,0,0)\, \text{m}$|⁠, where the xy plane is horizontal and the z-axis is vertically downward. For clarity, we present all symbols used in both Fig. 1 and the remaining part of this study in Table 1. The density contrast in the polyhedral prism is expressed as a polynomial function of arbitrary order with respect to depth, that is

\begin{eqnarray*}\lambda (z)=\sum _{n=0}^{Q}a_{n}z^{n}, \end{eqnarray*}

(1)

where n is the power index, and Q is a non-negative integer to represent the density contrast function, which can be an arbitrary non-negative integer. The constant coefficients a_n can be estimated by fitting the measured data sets such as borehole data (Blakely 1996). Please be aware that once the observation site moves to another position, coefficients a_n need to be re-evaluated, as we assume |$\mathbf {r}^{\prime }=(0,0,0)\, \text{m}$|⁠. The gravitational potential ϕ caused by this polyhedral prism mass body is (Blakely 1996)

\begin{eqnarray*}\phi =G \iiint _{H}\frac{\lambda (z)}{R} \text{d}v, \end{eqnarray*}

(2)

where λ(z) is the density contrast in eq. (1), H denotes the volume of the polyhedral prism mass body, R is the distance between the observation site |$\mathbf {r}^{\prime }=(0,0,0)\, \text{m}$| and an arbitrary running integral point in the polyhedral prism mass body |$\mathbf {r}=(x,y,z)\in H$|⁠, |$R=|\mathbf {r}^{\prime }-\mathbf {r}|$|⁠, and G is the gravitational attraction constant, |$G=6.673\times 10^{-11}\, \text{m}^3\text{kg}^{-1}\text{s}^{-2}$|⁠.

Figure 1.

Illustration of a vertical polyhedral prism H of polygonal cross-section with polynomial density contrast λ(z) of arbitrary order in depth (z) in a right-handed Cartesian coordinate system, where the xy plane is horizontal and the z-axis is vertically downward. The ith planar facet is represented by ∂H_i, and |$\partial H=\bigcup _{i=1}^{N}\partial H_{i}$| is the union of all the facets, representing the prism H, with N ≥ 5. For simplicity, we include the bottom facet ∂H₊ and top facet ∂H₋ among the N surfaces. The meanings of all symbols are given in Table 1. The observation site |$\mathbf {r}^{\prime }$| is located at the origin of the Cartesian coordinate system.

Table 1.

List of symbols for the polyhedral prism as shown in Fig. 1.

H	a polyhedral prism or its volume
λ(z)	density contrast along depth
N	number of faces on the polyhedral prism, N ≥ 5
∂H−	top surface of the polyhedral prism
∂H+	bottom surface of the polyhedral prism
∂Hi	ith surface of the polyhedral prism, 1 ≤ i ≤ N
z+, z−	z-coordinates of ∂H+ and ∂H−, respectively
|$\mathbf {r}^{\prime }$|	observation site with coordinates |$(0,0,0)\, \text{m}$|
|$\mathbf {r}$|	arbitrary source point with coordinates (x, y, z) in H
R	distance from |$\mathbf {r}^{\prime }$| to |$\mathbf {r}$|
|$\mathbf {o}_i$|	projection of |$\mathbf {r}^{\prime }$| on the ith surface ∂Hi, that is, ∂H− in the example in Fig. 1
Mi	number of edges on the ith surface
Cij	jth edge of the ith surface
|$\mathbf {\hat{n}}_i$|	outward-pointing normal unit vector of the ith surface
|$\mathbf {\hat{m}}_{ij}$|	outward-pointing normal unit vector of edge Cij lying within its defining surface ∂Hi
|$\mathbf {\hat{e}}_{ij}$|	tangential unit vector of edge Cij
|$\mathbf {v}_{0ij}, \mathbf {v}_{1ij}$|	start and end vertices of edge Cij, respectively
R0ij, R1ij	distances from |$\mathbf {r}^{\prime }$| to |$\mathbf {v}_{0ij}$| and |$\mathbf {v}_{1ij}$|⁠, respectively
sij	projected coordinate of vector |$(\mathbf {r}-\mathbf {r}^{\prime })$| along |$\hat{\mathbf {e}}_{ij}$|
s0ij, s1ij	parameterised coordinates sij for |$\mathbf {v}_{0ij}$| and |$\mathbf {v}_{1ij}$|⁠, respectively
hi	projection of (⁠|$\mathbf {r}^{\prime } - \mathbf {r}$|⁠) onto |$\hat{\mathbf {n}}_i$|
mij	projection of |$(\mathbf {r}-\mathbf {o}_i)$| onto |$\mathbf {\hat{m}}_{ij}$|⁠, |$\mathbf {r}\in C_{ij}$|
|$\mathbf {p}_{ij}^{\perp }$|	vertical projection of |$\mathbf {r}^{\prime }$| onto edge Cij
|$\hat{\boldsymbol {\rho }}_{ij}^{\perp }$|	unit vector from point |$\mathbf {o}_i$| to point |$\mathbf {p}_{ij}^{\perp }$|
τij	distance between |$\mathbf {o}_i$| and |$\mathbf {r}$|⁠, |$\mathbf {r}\in C_{ij}$|
|$\beta (\mathbf {o}_i)$|	angle of the circular region centred at |$\mathbf {o}_i$| on ∂Hi
βij	angle between the points |$\mathbf {v}_{0ij}$| and |$\mathbf {v}_{1ij}$| as measured from point |$\mathbf {o}_{i}$|

Substituting eq. (1) into eq. (2), the gravity potential is expressed as

\begin{eqnarray*}\phi &=&G\iiint _{H}\sum _{n=0}^{Q}a_{n}\frac{z^{n}}{R}\text{d}v =G\sum _{n=0}^{Q}a_{n}\phi ^{n}, \end{eqnarray*}

(3)

where we define |$\phi ^{n} = \iiint _{H}z^{n}\frac{1}{R}\text{d}v$|⁠.

Applying the gradient operator with respect to the coordinates |$\mathbf {r}$| of a source point and using the fact |$\mathbf {r}^{\prime }=(0,0,0)\, \text{m}$|⁠, we have

\begin{eqnarray*}z^{n}\frac{1}{R}&=&z^{n-1}\hat{\mathbf {z}}\cdot \nabla {R}. \end{eqnarray*}

(4)

Now, using eq. (4), the volume integral for the summand ϕⁿ of the gravity potential in eq. (3) can be further written as

\begin{eqnarray*}\phi ^{n}=\iiint _{H}z^{n}\frac{1}{R}\text{d}v &=&\iiint _{H}z^{n-1}\hat{\mathbf {z}}\cdot \nabla {R}\text{d}v\nonumber \\&=&\iiint _{H}\left[ \nabla \cdot \left(Rz^{n-1}\hat{\mathbf {z}})-R\nabla \cdot (z^{n-1}\hat{\mathbf {z}}\right) \right]\text{d}v \nonumber \\&=&\subset\!\!\!\!\!\!\!\!{\iint}\!\!\!\!\!\!\! \supset _{\partial H}Rz^{n-1}\left(\hat{\mathbf {z}}\cdot \hat{\mathbf {n}}\right)\text{d}s-\iiint _{H}(n-1)Rz^{n-2}\text{d}v, \end{eqnarray*}

(5)

for n > 1. Here, ∂H represents the closed surface of the polyhedral prism. In the above, as shown in Fig. 1, the vertical polygonal faces of the polyhedral prism are perpendicular to the horizontal plane.

Assuming that |$\hat{\mathbf {n}}$| is the outward normal unit vector to the surface of the polyhedral prism body H, |$\hat{\mathbf {z}}$| is the unit vector of the z-axis, we have |$\hat{\mathbf {z}}\cdot \hat{\mathbf {n}}=0$| for the vertical polygonal faces of the polyhedral prism, and |$\hat{\mathbf {z}}\cdot \hat{\mathbf {n}}=\pm 1$| for its two horizontal (bottom and top) surfaces. Hence, ϕⁿ in eq. (5) becomes:

\begin{eqnarray*}\phi ^{n}=z_{+}^{n-1}\iint _{\partial H_+}R\text{d}s- z_{-}^{n-1}\iint _{\partial H_-}R\text{d}s-(n-1) \iiint _{H}Rz^{n-2}\text{d}v, \end{eqnarray*}

(6)

where, |$\iint _{\partial H_\pm }$| (⁠|$\iint _{\partial H_+}$| and |$\iint _{\partial H_-})$| represents the surface integrals on bottom (∂H₊) and top (∂H₋) surfaces, and symbols z₊ and z₋ are the z-coordinates of the bottom surface and the top surface, respectively.

If n = 0, the gravity potential becomes

\begin{eqnarray*}\phi ^{0}=\iiint _{H}\frac{1}{R}\text{d}v . \end{eqnarray*}

(7)

If n = 1, eq. (6) becomes

\begin{eqnarray*}\phi ^{1}=\iiint _{H}z\frac{1}{R}\text{d}v &=&\iiint _{H}\hat{\mathbf {z}}\cdot \nabla {R}\text{d}v=\iiint _{H}\nabla \cdot (\hat{\mathbf {z}}{R})\text{d}v\nonumber \\&=&\subset\!\!\!\!\!\!\!\!{\iint}\!\!\!\!\!\!\! \supset _{\partial H}R\left(\hat{\mathbf {z}}\cdot \hat{\mathbf {n}}\right)\text{d}s\nonumber \\&=& \iint _{\partial H_+}R\text{d}s - \iint _{\partial H_-}R\text{d}s. \end{eqnarray*}

(8)

Applying the gradient operation to the gravitational potential in eq. (2) with respect to the coordinates |$\mathbf {r}^{\prime }$| of the observation point , we obtain the gravity or gravitational field vector, |$\mathbf {g}=\nabla ^{\prime } \phi$|⁠. In most geophysical problems, we are only interested in the gravity acceleration along the vertical (or the z-axis) direction. This is simply so, because most instruments are designed to measure only the vertical component of the gravity anomaly or gravity acceleration. Therefore, in this study, only the scalar vertical (z) component of the gravity anomaly or gravity field is considered. The horizontal components of the gravity field can be derived using procedures similar to the ones introduced in the following. The vertical component of the gravity anomaly (g, subscript is dropped for simplicity) can be expressed as

\begin{eqnarray*}g&=&\nabla ^{\prime } \phi \cdot \hat{\mathbf {z}} = G\iiint _{H}\lambda (z)\frac{z}{R^{3}}\text{d}v\nonumber \\&=& G\iiint _{H}\sum _{n=0}^{Q}a_{n}z^{n}\frac{z}{R^{3}}\text{d}v \nonumber \\&=&G\sum _{n=0}^{Q}a_{n}\iiint _{H}z^{n}\frac{z}{R^{3}}\text{d}v \nonumber \\&=&G\sum _{n=0}^{Q}a_{n}g^{n}, \end{eqnarray*}

(9)

where |$\hat{\mathbf {z}}$| is the unit vector along the z axis, and |$g^{n}=\iiint _{H}z^{n}\frac{z}{R^{3}}\text{d}v$|⁠.

We observe that the nth contribution of the summation in eq. (9) can be further represented as

\begin{eqnarray*}g^{n}&=&\iiint _{H}z^{n}\frac{z}{R^{3}}\text{d}v={\color {black}-}\iiint _{H}z^{n}\frac{\partial }{\partial {}z}\frac{1}{R}\text{d}v \nonumber \\&=&-\iiint _{H}z^{n}\hat{\mathbf {z}}\cdot \nabla \frac{1}{R}\text{d}v, \end{eqnarray*}

(10)

where the gradient operator is on the source point |$\mathbf {r}$|⁠. Using the following vector identity (Jin 2002)

\begin{eqnarray*}z^{n}\hat{\mathbf {z}}\cdot \nabla \frac{1}{R}+\frac{1}{R}\nabla \cdot \left(z^{n}\hat{\mathbf {z}}\right)=\nabla \cdot \left(\frac{1}{R}z^{n}\hat{\mathbf {z}}\right), \end{eqnarray*}

(11)

eq. (10) is rewritten as

\begin{eqnarray*}g^{n} & = & -\iiint _{H}\left[\nabla \cdot \left(\frac{1}{R}z^{n}\hat{\mathbf {z}}\right)-\frac{1}{R}\nabla \cdot \left(z^{n}\hat{\mathbf {z}}\right)\right]\text{d}v. \end{eqnarray*}

(12)

In terms of the divergence theorem (Jin 2002), for n > 0, eq. (12) becomes

\begin{eqnarray*}g^{n} & = & -\subset\!\!\!\!\!\!\!\!{\iint}\!\!\!\!\!\!\! \supset _{\partial H}\frac{1}{R}z^{n}\left(\hat{\mathbf {z}}\cdot \hat{\mathbf {n}}\right)\text{d}s+n\iiint _{H}\frac{z^{n-1}}{R}\text{d}v, \end{eqnarray*}

(13)

where |$\hat{\mathbf {n}}$| is the outward unit normal vector on the closed surface ∂H of the polyhedral prism H.

Using the geometry of the polyhedral prism as shown in Fig. 1, eq. (13) becomes

\begin{eqnarray*}g^{n} & = & z_{-}^{n}\iint _{\partial H_-}\frac{1}{R}\text{d}s- z_{+}^{n}\iint _{\partial H_+}\frac{1}{R}\text{d}s+n\iiint _{H}\frac{z^{n-1}}{R}\text{d}v \nonumber \\& = & z_{-}^{n}\iint _{\partial H_-}\frac{1}{R}\text{d}s- z_{+}^{n}\iint _{\partial H_+}\frac{1}{R}\text{d}s+n\phi ^{n-1}, \end{eqnarray*}

(14)

for n > 0.

For n = 0, eq. (13) becomes

\begin{eqnarray*}g^{0} & = & \iint _{\partial H_-}\frac{1}{R}\text{d}s-\iint _{\partial H_+}\frac{1}{R}\text{d}s. \end{eqnarray*}

(15)

Known from eq. (14), once closed-form solutions of ϕ^{n − 1} are given, only the surface integrals |$\iint _{\partial H_\pm }R^{-1}\text{d}s$| remain to be evaluated to derive the closed-form solutions for gⁿ.

2.2 Closed-form equations

To deal with ϕⁿ in eq. (6), we first need to evaluate volume integrals of type |$\iiint _{H}z^{n-2} R\text{d}v$| for n ≥ 2. For ϕ⁰ and ϕ¹, please refer to eqs (7) and (8), respectively. The above volume integral can be expressed in a compact form,

\begin{eqnarray*}K(t,q)=\iiint _{H}z^{t}R^{q}\text{d}v. \end{eqnarray*}

(16)

Here, as n ≥ 2 in |$\iiint _{H}z^{n-2} R\text{d}v$|⁠, the integer number t satisfies t ≥ 0. Furthermore, as shown below, we will transform the integrals iteratively such that the order q of R increases by 2 in each step. Thus, since the initial value of q is 1, the subsequent values of q are odd integers q ≥ 1, q = 1, 3, 5, …, ∞.

Using the fact |$\mathbf {r}^{\prime }=(0,0,0)\, \text{m}$| and applying the gradient operator on the source point, we get the following relationship

\begin{eqnarray*}z^{t}R^{q}&=& \frac{1}{q+2}z^{t-1}\hat{\mathbf {z}}\cdot \nabla R^{q+2}, \end{eqnarray*}

(17)

where t ≥ 0 and q ≥ 1. Thus, the general volume integral in eq. (16) is written as

\begin{eqnarray*}K(t,q)=\iiint _{H}z^{t}R^{q}\text{d}v &=& \iiint _{H}\frac{1}{q+2}z^{t-1}\hat{\mathbf {z}}\cdot \nabla R^{q+2}\text{d}v \nonumber \\&=& \frac{1}{q+2}\iiint _{H}\left[ \nabla \cdot (R^{q+2}z^{t-1}\hat{\mathbf {z}})-R^{q+2}\nabla \cdot (z^{t-1}\hat{\mathbf {z}}) \right]\text{d}v \nonumber \\&=& \frac{1}{q+2}\subset\!\!\!\!\!\!\!\!{\iint}\!\!\!\!\!\!\! \supset _{\partial H}R^{q+2}z^{t-1}\left(\hat{\mathbf {z}}\cdot \hat{\mathbf {n}}\right)\text{d}s-\frac{t-1}{q+2}\iiint _{H}R^{q+2}z^{t-2}\text{d}v\nonumber \\&=& \frac{z_{+}^{t-1}}{q+2}\iint _{\partial H_+}R^{q+2}\text{d}s-\frac{z_{-}^{t-1}}{q+2}\iint _{\partial H_-}R^{q+2}\text{d}s-\frac{t-1}{q+2}\iiint _{H}z^{t-2}R^{q+2}\text{d}v\nonumber \\&=& \frac{z_{+}^{t-1}}{q+2}\iint _{\partial H_+}R^{q+2}\text{d}s-\frac{z_{-}^{t-1}}{q+2}\iint _{\partial H_-}R^{q+2}\text{d}s-\frac{t-1}{q+2}K(t-2,q+2), \end{eqnarray*}

(18)

for t ≥ 2.

If t = 0, we get

\begin{eqnarray*}K(0,q) & = & \iiint _{H}R^{q} \text{d}v. \end{eqnarray*}

(19)

If t = 1, we get

\begin{eqnarray*}K(1,q) & = & \iiint _{H}zR^{q} \text{d}v = \frac{1}{q+2}\iint _{\partial H_+}R^{q+2}\text{d}s-\frac{1}{q+2}\iint _{\partial H_-}R^{q+2}\text{d}s, \end{eqnarray*}

(20)

that is, there is not any volume integral to be evaluated.

Eq. (18) shows that the volume integral K(t, q) is transformed into a volume integral term K(t′, q′) plus surface terms, such that t′ = t − 2 and q′ = q + 2. Note that t + q = t′ + q′. It follows that when t is an even integer greater than 0, K(t, q) is transformed to K(0, t + q) plus a sum of surface integrals in t/2 steps, and when t is an odd integer greater than 0, K(t, q) is transformed to K(1, q + t − 1) plus a sum of surface integrals in (t − 1)/2 steps. Now, we further observe that when the order of z becomes 1 (order of R changes to q + t − 1), eq. (20) is straight-forward to apply and, consequently, only the surface integrals |$\iint _{\partial H_\pm }R^{q+t+1}\text{d}s$| remain to be evaluated.

Thus, according to eqs (18)–(20), for any number t ≥ 0, the volume integral |$\iiint _{H}z^{t}R^{q}\text{d}v$| is eventually transformed into a simple volume integral |$\iiint _{H}R^{k}\text{d}v$| and two surface integrals |$\iint _{\partial H_\pm }R^{p}\text{d}s$| over the bottom and top surfaces of the polyhedral prism mass body. Note that both k and p are odd integers, because q is an odd integer (cf. discussion after eq. 16). The ranges of these odd integers and their final values after an iterative application of eq. (18) are listed in the following for the different cases:

1. When t = 2, 4, …, ∞, k = q + t, |$p=q+2*\frac{t}{2}$|⁠, and, as q ≥ 1, it follows that k ≥ 3 and p ≥ 3. Both volume and surface integrals need to be evaluated.

2. When t = 3, 5, …, ∞, |$p=q+2*\frac{t+1}{2}$|⁠, and, as q ≥ 1, it follows p ≥ 5. No volume integral needs to be evaluated.

3. When t = 1, p = q + 2, and, as q ≥ 1, it follows p ≥ 3. No volume integral needs to be evaluated.

4. When t = 0, k = q, and, as q ≥ 1, it follows k ≥ 1. No surface integrals need to be evaluated.

Finally, by setting t = n − 2, the task of evaluating the volume integral |$\iint _{H}z^{n-2}R\text{d}v$| (i.e. K(n − 2, 1) in eq. 16) for the gravitational potential ϕⁿ (n ≥ 2) is transformed into the simple volume integral |$\iiint _{H}R^{k}\text{d}v$| and the simple surface integrals |$\iint _{\partial H_\pm }R^{p}\text{d}s$|⁠. These integrals are also representative of the remaining integrals |$\iiint _{H}\frac{1}{R}\text{d}v$| and |$\iint _{\partial H_\pm }R\text{d}s$| in eqs (6)–(8) for the gravitational potentials ϕ⁰, ϕ¹ and ϕⁿ and |$\iint _{\partial H_\pm }\frac{1}{R}\text{d}s$| for the gravity field gⁿ in eq. (14). Therefore, our final task of evaluating ϕⁿ (n ≥ 0, also gⁿ) becomes to evaluate:

\begin{eqnarray*}K(0,k) &=& \iiint _{H}R^{k}\text{d}v, \quad k\ge -1 , \end{eqnarray*}

(21)

\begin{eqnarray*}S(\partial H_\pm ,p) &=& \iint _{\partial H_\pm }R^{p}\text{d}s, \quad p\ge -1, \end{eqnarray*}

(22)

where k and p are two odd integers. Eqs (21) and (22) act as fundamental building blocks for deriving closed-form solutions in our recursive formulae.

First, we proceed with the general volume integral K(0, k) in eq. (21). Using the following vector identity:

\begin{eqnarray*}R^{k} &=& -\frac{1}{k+3}\nabla \cdot \left[ \left(\mathbf {r}^{\prime }-\mathbf {r} \right)R^k \right], \end{eqnarray*}

(23)

and the divergence theorem (Jin 2002), K(0, k) becomes

\begin{eqnarray*}K(0,k)=\iiint _{H}R^{k}\text{d}v &=& \iiint _{H}-\frac{1}{k+3}\nabla \cdot \left[ \left(\mathbf {r}^{\prime }-\mathbf {r} \right)R^k \right] \text{d}v , \quad k\ge -1 \nonumber \\&=& -\frac{1}{k+3}\subset\!\!\!\!\!\!\!\!{\iint}\!\!\!\!\!\!\! \supset _{\partial H}\hat{\mathbf {n}}\cdot \left(\mathbf {r}^{\prime }-\mathbf {r} \right)R^k \text{d}s\nonumber \\&=& -\frac{1}{k+3}\sum _{i=1}^{N}h_i\iint _{\partial H_{i}}R^k \text{d}s\nonumber \\&=& -\frac{1}{k+3}\sum _{i=1}^{N}h_i S(\partial H_{i},k), \end{eqnarray*}

(24)

where |$S(\partial H_{i},k)=\iint _{\partial H_{i}}R^k\text{d}s$|⁠, and |$h_i = \hat{\mathbf {n}}_i\cdot \left(\mathbf {r}^{\prime }-\mathbf {r} \right)$| is the projection of the vector |$\left(\mathbf {r}^{\prime }-\mathbf {r} \right)$| onto the normal vector |$\hat{\mathbf {n}}_i$| of the ith surface of the polyhedral body H. |$\partial H = \bigcup _{i=1}^{N}\partial H_{i}$| is the whole surface of the polyhedral prism, including the top surface, the bottom surface and the vertical surfaces. If the ith surface is the bottom or top surface, ∂H_i = ∂H₊ or ∂H_i = ∂H₋. Using eq. (24), calculation of K(0, k) has been transformed to the task of calculating a number of surface integrals similar to the one in eq. (22). According to the summation in eq. (24), these surface integrals need to be evaluated over all faces of the polyhedral prism. To proceed, we have to find the closed-form solutions for surface integrals with integrand R^p over the surfaces of the polyhedral prism.

Now, we have to derive the exact solution for the following surface integrals:

\begin{eqnarray*}S(\partial H_{i},p) &=& \iint _{\partial H_{i}}R^{p}\text{d}s, \quad p\ge -1. \end{eqnarray*}

(25)

If p = −1, the integrand of integral |$\iint _{\partial H_{i}}R^{-1}\text{d}s$| is singular when R → 0. Here, we use the singularity-free method to analytically calculate the integral (Yla-Oijala & Taskinen 2003). We use the following identity (Ren et al.2017a):

\begin{eqnarray*}R^{p} &=& \nabla _{s}\cdot \left[ \frac{1}{p+2}\frac{R^{p+2}}{\tau _{ij}^2}\left(\mathbf {r}-\mathbf {o}_i\right) \right], \end{eqnarray*}

(26)

where ∇_s · is the surface divergence operator. As shown in Fig. 1, taking the top surface as the surface ∂H_i for demonstration, the observation point |$\mathbf {r}^{\prime }$| is projected on the top surface with point |$\mathbf {o}_i$| as the projection point, which can be evaluated by |$\mathbf {r}^{\prime }$| and the vertices of facet ∂H_i (see Appendix  A). Point |$\mathbf {r}$| is an arbitrary point on the surface ∂H_i, |$\mathbf {r}\in \partial H_{i}$|⁠. The variable |$\tau _{ij}=|\mathbf {r}-\mathbf {o}_i|$| (cf. Fig. 1) is the distance from the point |$\mathbf {r}$| to point |$\mathbf {o}_i$| and can also be considered as the polar radial coordinate on the top surface with the origin at the projection point |$\mathbf {o}_i$|⁠.

We assume that the ith polygonal surface ∂H_i of the polyhedral prism has M_i edges (each side face has 4 edges), that is |$\partial (\partial {H}_{i}) = \bigcup _{j=1}^{M_i}C_{ij}$|⁠, C_ij is the jth edge of the ith surface. Then, using eq. (26), the surface divergence theorem (Jin 2002) and results from our previous analysis (Ren et al.2017a,b), S(∂H_i, p) becomes

\begin{eqnarray*}S(\partial {H}_{i},p)=\iint _{\partial {H}_{i}}R^{p}\text{d}s&=&\iint _{\partial {H}_{i}}\nabla _{s}\cdot \left[\frac{1}{p+2}\frac{R^{p+2}}{\tau _{ij}^2}\left(\mathbf {r}-\mathbf {o}_i\right)\right]\text{d}s \nonumber \\&=&\frac{1}{p+2}\sum _{j=1}^{M_i}m_{ij}\int _{C_{ij}}\frac{R^{p+2}}{\tau _{ij}^2}\text{d}l-\frac{\beta (\mathbf {o}_i)}{p+2}|h_i|^{p+2}\nonumber \\&=&\frac{1}{p+2}\sum _{j=1}^{M_i}B_{ij}^{p+2}-\frac{\beta (\mathbf {o}_i)}{p+2}|h_i|^{p+2}, \end{eqnarray*}

(27)

where p = −1, 1, 3, ......, ∞, |$h_i=\hat{\mathbf {n}}_i \cdot \left(\mathbf {r}^{\prime }-\mathbf {r}\right)$|⁠, |$m_{ij}=\left( \mathbf {r}-\mathbf {o}_i\right)\cdot \hat{\mathbf {m}}_{ij}$| is constant over edge C_ij, and |$\hat{\mathbf {m}}_{ij}$| is the outgoing normal unit vector of edge C_ij. The detailed procedure to derive eq. (27) can be found in Ren et al. (2017a). It is worth noting that the formula for the surface integral in eq. (27) is still valid, when the integral surface is one of the vertical surfaces of the prismatic mass body H.

The variable |$\beta (\mathbf {o}_i)$| in eq. (27) is the angle that facet ∂H_i subtends at the point |$\mathbf {o}_i$|⁠. As known from geometric analysis, this angle has the following values: if the centre point |$\mathbf {o}_i$| is inside the polygon, |$\beta (\mathbf {o}_i)=2\pi$|⁠; if the point |$\mathbf {o}_i$| is on the edge of the polygon, |$\beta (\mathbf {o}_i)=\pi$|⁠; |$\beta (\mathbf {o}_i)=\theta$| if point |$\mathbf {o}_i$| is at a corner of the polygon with the corner angle being θ; and if point |$\mathbf {o}_i$| is outside of the polygon, |$\beta (\mathbf {o}_i)=0$|⁠. Although values of |$\beta (\mathbf {o}_i)$| can be determined by the geometrical relationship between the observation site and the mass body, we prefer its alternative stable presentation:

\begin{eqnarray*}\beta (\mathbf {o}_i) = \sum _{j=1}^{M_i}\beta _{ij}= \sum _{j=1}^{M_i}\hat{\mathbf {m}}_{ij}\cdot \hat{\boldsymbol {\rho }}_{ij}^{\perp }\left(\arctan {\frac{s_{1ij}}{|m_{ij}|}}-\arctan {\frac{s_{0ij}}{|m_{ij}|}}\right), \end{eqnarray*}

(28)

where β_ij is the contribution from edge C_ij. The in-plane normal unit vector to edge C_ij is |$\hat{\mathbf {m}}_{ij}$|⁠, and |$\hat{\boldsymbol {\rho }}_{ij}^{\perp }$| is the unit vector from point |$\mathbf {o}_i$| to point |$\mathbf {p}_{ij}^{\perp }$|⁠, which is the projection point of point |$\mathbf {o}_i$| onto edge C_ij (Fig. 1). The position of source point |$\mathbf {r}$| along edge C_ij (⁠|$\mathbf {r}\in C_{ij}$|⁠) is parametrized by a 1-D variable s_ij,

\begin{eqnarray*}s_{ij}=(\mathbf {r}-\mathbf {o}_i)\cdot \hat{\mathbf {e}}_{ij}, \end{eqnarray*}

(29)

where |$\hat{\mathbf {e}}_{ij}$| is the tangential unit vector of edge C_ij. The values of the parametrized 1-D variable at the two vertices |$\mathbf {v}_{0ij}$| and |$\mathbf {v}_{1ij}$| of edge C_ij are s_0ij and s_1ij, respectively. When the observation site |$\mathbf {r}^{\prime }$| is located on edge C_ij or at one of its vertices, we have m_ij = 0 and h_i = 0, and one may suspect a singularity to occur in the arctan functions in eq. (28). However, the value range of the arctan function is bounded to the interval |$]-\frac{\pi }{2},\frac{\pi }{2}[$|⁠. Hence, β_ij is finite and no singularities exist.

The line integral in eq. (27) can be calculated by another recursive formula:

\begin{eqnarray*}B_{ij}^{p+2} = m_{ij}\int _{C_{ij}}\frac{R^{p+2}}{\tau _{ij}^2}\text{d}l&=& m_{ij}\int _{C_{ij}}R^{p+1}\frac{R}{\tau _{ij}^2}\text{d}l\nonumber \\&=& m_{ij}\int _{C_{ij}}R^{p+1}(\frac{1}{R}+\frac{h_i^2}{R\tau _{ij}^2})\text{d}l\nonumber \\&=&m_{ij}\int _{C_{ij}}R^{p}\text{d}l+h_i^2m_{ij}\int _{C_{ij}}\frac{R^{p}}{\tau _{ij}^2}\text{d}l\nonumber \\&=&m_{ij}\int _{C_{ij}}R^{p}\text{d}l+h_i^2 B_{ij}^{p}, \end{eqnarray*}

(30)

where we used |$R=\sqrt{h_i^2+m_{ij}^2+s_{ij}^2}$| and |$\tau _{ij}=\sqrt{m_{ij}^2+s_{ij}^2}$|⁠. In the second term on the right hand side of the above formula, the initial line integral is as follows (Ren et al.2017a, eq. 51):

\begin{eqnarray*}B_{ij}^{1} &=&m_{ij}\int _{s_{0ij}}^{s_{1ij}}\frac{\sqrt{h_i^2+m_{ij}^2+s_{ij}^2}}{m_{ij}^2+s_{ij}^2}\text{d}s_{ij} \nonumber \\&=&|h_i|\left( \arctan {\frac{|h_i|s_{1ij}}{m_{ij} R_{1ij}}}-\arctan {\frac{|h_i|s_{0ij}}{m_{ij} R_{0ij}}} \right)+m_{ij}\ln {\frac{s_{1ij}+R_{1ij}}{s_{0ij}+R_{0ij}}}, \end{eqnarray*}

(31)

where R_0ij and R_1ij are the distances from the observation site to the two vertices (⁠|$\mathbf {v}_{0ij}$| and |$\mathbf {v}_{1ij}$|⁠) of edge C_ij, respectively. Note that the index p in eq. (27) is an odd integer p = −1, 1, 3, 5, …, ∞. For p = 1, 3, 5, …, ∞, the computation of |$B_{ij}^{1}$| is sufficient to initialise the resulting recursive formula in eq. (30). In the case p = −1, eq. (27) leads to the computation of |$B_{ij}^{1}$| in eq. (31) directly without any recursion using eq. (30). When the observation site |$\mathbf {r}^{\prime }$| is located on edge C_ij (including its two vertices), where m_ij = 0 and h_i = 0, we simply have |$B_{ij}^{1}=0$| so that the line integral in eq. (31) is singularity-free.

Additionally, a closed-form solution for the line integral |$\int _{C_{ij}}R^{p}\text{d}l$| in the first term on the right-hand side of eq. (30) is needed. Using integration by parts yields (Gradshteyn & Ryzhik 2007, eq. (2.260)):

\begin{eqnarray*}L_{ij}^p=\int _{C_{ij}}R^{p}\text{d}l &=& \int _{s_{0ij}}^{s_{1ij}}\left(\sqrt{h_i^2+m_{ij}^2+s_{ij}^2}\right)^{p}\text{d}s_{ij}\nonumber \\&=&\left.\frac{s_{ij}R^{p}}{p+1}\right|_{s_{0ij}}^{s_{1ij}}+\frac{p}{p+1}\left(h_i^2+m_{ij}^2\right)\int _{s_{0ij}}^{s_{1ij}}R^{p-2}\text{d}s_{ij}\nonumber \\&=&\left.\frac{s_{ij}R^{p}}{p+1}\right|_{s_{0ij}}^{s_{1ij}}+\frac{p}{p+1}\left(h_i^2+m_{ij}^2\right)L_{ij}^{p-2}, \end{eqnarray*}

(32)

and the initial integral is

\begin{eqnarray*}L_{ij}^1=\int _{C_{ij}}R \text{d}l = \frac{1}{2}\left[ (h_i^2+m_{ij}^2)\ln {\frac{s_{1ij}+R_{1ij}}{s_{0ij}+R_{0ij}}}+s_{1ij}R_{1ij}-s_{0ij}R_{0ij} \right]. \end{eqnarray*}

(33)

Note that eq. (27) leads to eq. (30) and, with that, to the integral |$L_{ij}^p$| only for indices p = 1, 3, 5, …, ∞, but not for p = −1. Hence, eq. (33) is sufficient to initialize the recursive formula in eq. (32). When the observation site |$\mathbf {r}^{\prime }$| is located on edge C_ij or its vertices (m_ij = 0 and h_i = 0), we set the first term of |$L_{ij}^1$| in eq. (33) to zero. By contrast, the second and third terms may have finite values, and as the observation site |$\mathbf {r}^{\prime }$| approaches to one of the vertices, the corresponding term s_0ijR_0ij or s_1ijR_1ij will tend to zero. For locations on the edge, the first term of |$L_{ij}^p$| in eq. (32) is finite, and the second term of |$L_{ij}^p$| is zero. Hence, when the observation site approaches to the edges of the polyhedral prism, the line integral |$L_{ij}^p$| in eq. (32) is singularity-free.

2.3 Algorithm and singularity-free property

In summary, by using the recursive formulae for S(∂H_i, p) in eq. (27), |$B_{ij}^{p}$| in eq. (30) and |$L_{ij}^p$| in eq. (32), we can find the exact solutions of the gravity potential and the vertical gravity field of a polyhedral prism with vertical polynomial mass density contrast of arbitrary orders.

First, we consider ϕⁿ. According to eqs (6)–(8), (18) and (24), we have the following formulae:

For n = 0, we have

\begin{eqnarray*}\phi ^{0} = -\frac{1}{2}\sum _{i=1}^{N}h_i S(\partial H_i,-1). \end{eqnarray*}

(34)

For n = 1, we have

\begin{eqnarray*}\phi ^1=S(\partial H_+,1)-S(\partial H_-,1). \end{eqnarray*}

(35)

For n = 2, we have

\begin{eqnarray*}\phi ^2=z_+S(\partial H_+,1)-z_-S(\partial H_-,1)+\frac{1}{4}\sum _{i=1}^{N}h_i S(\partial H_i,1). \end{eqnarray*}

(36)

For n = 3, 5, 7, …, ∞, we have:

\begin{eqnarray*}\phi ^{n} &=& z_+^{n-1}S(\partial H_+,1)-z_-^{n-1}S(\partial H_-,1)\nonumber \\&+&\sum _{j=2}^{\frac{n+1}{2}}\frac{(n-1)\cdot \cdot \cdot (n-2j+3)}{(-1)^{j-1}\cdot 1\cdot \cdot \cdot (2j-3)}\left[ \frac{z_+^{n-2j+1}}{2j-1}S(\partial H_+,2j-1)-\frac{z_-^{n-2j+1}}{2j-1}S(\partial H_-,2j-1) \right]. \end{eqnarray*}

(37)

For n = 4, 6, …, ∞, we have

\begin{eqnarray*}\phi ^{n} &=& z_+^{n-1}S(\partial H_+,1)-z_-^{n-1}S(\partial H_-,1)\nonumber \\&+&\sum _{j=2}^{\frac{n}{2}}\frac{(n-1)\cdot \cdot \cdot (n-2j+3)}{(-1)^{j-1}\cdot 1\cdot \cdot \cdot (2j-3)}\left[ \frac{z_+^{n-2j+1}}{2j-1}S(\partial H_+,2j-1)-\frac{z_-^{n-2j+1}}{2j-1}S(\partial H_-,2j-1) \right]+(-1)^{\frac{n}{2}+1}\cdot \frac{1}{n+2}\sum _{i=1}^{N}h_i S(\partial H_i,n-1). \end{eqnarray*}

(38)

Second, we consider |$\mathbf {g}^{n}$|⁠. According to eqs (14) and (15), we have

\begin{eqnarray*}g^{n} = \left\lbrace {\begin{array}{*{20}{c}}S(\partial H_-,-1)-S(\partial H_+,-1), \qquad \text{for } n=0\\z_-^nS(\partial H_-,-1)-z_+^nS(\partial H_+,-1)+n\phi ^{n-1}, \qquad \text{for } n\ge 1 \end{array}} \right.. \end{eqnarray*}

(39)

From eqs (35), (37) and (39), we know that when n is an odd number, only the contributions from the top and bottom surfaces are needed to evaluate ϕⁿ and g^{n − 1}. Thus, we have to evaluate the surface integrals S(∂H_i, t), t = 1, 3, 5, .., ∞ and S(∂H_±, t), t = −1, 1, 3, ..., ∞. The latter surface integral S(∂H_±, t) is the special case of the former surface integral S(∂H_i, t), so that here we only describe the solution to calculate S(∂H_i, t). By setting p = t in eq. (27), we need to compute the angle β using eq. (28), the line integral |$B_{ij}^{p+2}$| using eq. (30) and the initial line integral in eq. (31).

As for the gravity potential ϕⁿ, the integrand of the volume integral in eq. (3), which contains a factor 1/R, is mathematically singular, when the observation point |$\mathbf {r}^{\prime }$| is close to the polyhedral prism body H (R → 0). Using the divergence theorem, the weak singularity of the gravity potential is eliminated except for term ϕ⁰, as shown in eqs (34) to (39) where only S(∂H_i, −1) is singular. Fortunately, we have derived the singularity-free exact solution for S(∂H_i, −1) in eqs (27) and (31). Therefore, our solutions for the gravity potential ϕⁿ are completely singularity-free.

The gravity field gⁿ has a stronger mathematical singularity of 1/R² in its integrand (R → 0). As ϕⁿ and S(∂H_i, −1) in eq. (39) are now both singularity-free, it is easy to show that gⁿ is also singularity-free.

2.4 Closed-form solutions for rectangular prism with low-order polynomial density contrast functions

The above singularity-free formulae are developed for a general prism with polyhedral cross-section. Since the rectangular prism is a commonly used element in the computation of the Earth’s gravity field, it is quite interesting to investigate the explicit closed-form solutions for this standard element. Thus far, closed-form solutions for the gravity field of the prism with low-order polynomial density contrasts for n ≤ 4 were widely considered (Garcia-Abdeslem 2005; Ren et al.2017a, 2018c; D’Urso & Trotta 2017). So, here, we only investigate this low-order case of n ≤ 4.

2.4.1 Constant order—Case n = 0

The zero-order density contrast is λ(z) = a₀. Using eqs (7), (24) and (15), the gravity potential and the vertical gravity field are

\begin{eqnarray*}\phi & = & Ga_{0}\phi ^{0}=Ga_{0}\iiint _{H}\frac{1}{R}\text{d}v=-\frac{Ga_{0}}{2}\sum _{i=1}^{N}h_i\iint _{\partial H_{i}}\frac{1}{R}\text{d}s,\end{eqnarray*}

(40)

\begin{eqnarray*}g & = & Ga_{0}g^{0}=Ga_{0}\left[ \iint _{\partial H_-}\frac{1}{R}\text{d}s-\iint _{\partial H_+}\frac{1}{R}\text{d}s, \right]. \end{eqnarray*}

(41)

By setting p = −1 in eq. (27), exact solutions for ϕ and g are obtained.

2.4.2 Linear order—Case n = 1

The density contrast is λ(z) = a₀ + a₁z. Using eqs (8), (14), (35) and (39), the gravity potential and the vertical gravity field are

\begin{eqnarray*}\phi & = & Ga_{0}\phi ^{0}+Ga_{1}\phi ^{1}=Ga_{0}\phi ^{0} + Ga_{1}\iint _{\partial H_+}R\text{d}s - Ga_{1}\iint _{\partial H_-}R\text{d}s,\end{eqnarray*}

(42)

\begin{eqnarray*}g & = & Ga_{0}g^{0}+Ga_{1}\left[ z_{-}\iint _{\partial H_-}\frac{1}{R}\text{d}s- z_{+}\iint _{\partial H_+}\frac{1}{R}\text{d}s+\iiint _{H}\frac{1}{R}\text{d}v \right] \nonumber \\& = & Ga_{0}g^{0}+Ga_{1}\left[ z_{-}\iint _{\partial H_-}\frac{1}{R}\text{d}s- z_{+}\iint _{\partial H_+}\frac{1}{R}\text{d}s-\sum _{i=1}^{N}\frac{h_i}{2}\iint _{\partial H_i}\frac{1}{R}\text{d}s \right]. \end{eqnarray*}

(43)

By setting p = 1 and p = −1 for the gravity potential and gravity field, respectively, in eq. (27) and using exact solutions for the terms of constant density contrast, exact solutions for ϕ and g are obtained.

2.4.3 Quadratic order—Case n = 2

The density contrast is λ(z) = a₀ + a₁z + a₂z². Using eqs (6), (14) and (36) and by setting n = 2 in eq. (39), the gravity potential and the vertical gravity field are

\begin{eqnarray*}\phi & = & \sum _{i=0}^{1}Ga_{i}\phi ^{i}+Ga_{2}\left[z_{+}\iint _{\partial H_+}R\text{d}s - z_{-}\iint _{\partial H_-}R\text{d}s-\iiint _{H}R\text{d}v\right] \nonumber \\& = & \sum _{i=0}^{1}Ga_{i}\phi ^{i}+Ga_{2}\left[z_{+}\iint _{\partial H_+}R\text{d}s - z_{-}\iint _{\partial H_-}R\text{d}s+\frac{1}{4}\sum _{i=1}^{N}h_i\iint _{\partial H_i}R\text{d}s\right],\end{eqnarray*}

(44)

\begin{eqnarray*}g & = & \sum _{i=0}^{1}Ga_{i}g^{i}+Ga_{2}\left[z_{-}^{2}\iint _{\partial H_-}\frac{1}{R}\text{d}s - z_{+}^{2}\iint _{\partial H_+}\frac{1}{R}\text{d}s+2\iiint _{H}\frac{z}{R}\text{d}v \right] \nonumber \\& = & \sum _{i=0}^{1}Ga_{i}g^{i}+Ga_{2}\left[z_{-}^{2}\iint _{\partial H_-}\frac{1}{R}\text{d}s - z_{+}^{2}\iint _{\partial H_+}\frac{1}{R}\text{d}s+ 2\iint _{\partial H_+}{R}\text{d}s - 2\iint _{\partial H_-}{R}\text{d}s \right]. \end{eqnarray*}

(45)

By setting p = 1 and p = −1 in eq. (27) and using exact solutions for the terms of constant and linear density contrasts, exact solutions for ϕ and g are obtained.

2.4.4 Cubic order—Case n = 3

The density contrast is λ(z) = a₀ + a₁z + a₂z² + a₃z³. In this case, using eqs (6), (14) and (20) and by setting n = 3 in eqs (37) and (39), the gravity potential and the vertical gravity field are

\begin{eqnarray*}\phi & = & \sum _{i=0}^{2}Ga_{i}\phi ^{i}+Ga_{3}\left[z_{+}^{2}\iint _{\partial H_+}R\text{d}s - z_{-}^{2}\iint _{\partial H_-}R\text{d}s-2\iiint _{H}zR \text{d}v \right] \nonumber \\& = & \sum _{i=0}^{2}Ga_{i}\phi ^{i}+ Ga_{3}\left[z_{+}^{2}\iint _{\partial H_+}R\text{d}s - z_{-}^{2}\iint _{\partial H_-}R\text{d}s-\frac{2}{3}\iint _{\partial H_+}{R^3}\text{d}s+\frac{2}{3}\iint _{\partial H_-}{R^3}\text{d}s \right],\end{eqnarray*}

(46)

\begin{eqnarray*}g & = & \sum _{i=0}^{2}Ga_{i}g^{i}+Ga_{3}\left[z_{-}^{3}\iint _{\partial H_-}\frac{1}{R}\text{d}s - z_{+}^{3}\iint _{\partial H_+}\frac{1}{R}\text{d}s+3\iiint _{H}\frac{z^2}{R}\text{d}v \right]\nonumber \\& = & \sum _{i=0}^{2}Ga_{i}g^{i}+Ga_{3}\left[\begin{array}{l}z_{-}^{3}\iint _{\partial H_-}\frac{1}{R}\text{d}s - z_{+}^{3}\iint _{\partial H_+}\frac{1}{R}\text{d}s +3z_{+}\iint _{\partial H_+}{R}\text{d}s\\-3z_{-}\iint _{\partial H_-}{R}\text{d}s+\frac{3}{4}\sum _{i=1}^{N}h_i\iint _{\partial H_i}R\text{d}s \end{array} \right]. \end{eqnarray*}

(47)

By setting p = 1, p = 3 and p = −1 in eq. (27), and using exact solutions from the constant, linear and quadratic cases, exact solutions for ϕ and g are obtained.

2.4.5 Quartic order—Case n = 4

The density contrast is λ(z) = a₀ + a₁z + a₂z² + a₃z³ + a₄z⁴. In this case, using eqs (6), (14) and (18) and by setting n = 4 in eqs (38) and (39), the gravity potential and the vertical gravity field are

\begin{eqnarray*}\phi & = & \sum _{i=0}^{3}Ga_{i}\phi ^{i}+Ga_{4}\left[z_{+}^{3}\iint _{\partial H_+}R\text{d}s- z_{-}^{3}\iint _{\partial H_-}R\text{d}s-3 \iiint _{H}z^{2}R\text{d}v \right]\nonumber \\& = & \sum _{i=0}^{3}Ga_{i}\phi ^{i}+Ga_{4}\left[\begin{array}{l}z_{+}^{3}\iint _{\partial H_+}R\text{d}s- z_{-}^{3}\iint _{\partial H_-}R\text{d}s - z_{+}\iint _{\partial H_+}R^3\text{d}s \\+ z_{-}\iint _{\partial H_-}R^3ds - \frac{1}{6}\sum _{i=1}^{N}h_i\iint _{\partial H_i}R^3ds \end{array} \right],\end{eqnarray*}

(48)

\begin{eqnarray*}g & = & \sum _{i=0}^{3}Ga_{i}g^{i}+Ga_{4}\left[ z_{-}^{4}\iint _{\partial H_-}\frac{1}{R}\text{d}s- z_{+}^{4}\iint _{\partial H_+}\frac{1}{R}\text{d}s+4\iiint _{H}\frac{z^{3}}{R}\text{d}v \right]\nonumber \\& = & \sum _{i=0}^{3}Ga_{i}g^{i}+Ga_{4}\left[\begin{array}{l}z_{-}^{4}\iint _{\partial H_-}\frac{1}{R}\text{d}s- z_{+}^{4}\iint _{\partial H_+}\frac{1}{R}\text{d}s + 4z_{+}^{2}\iint _{\partial H_+}R \text{d}s\\- 4z_{-}^{2}\iint _{\partial H_-}R \text{d}s - \frac{8}{3}\iint _{\partial H_+}R^3 \text{d}s + \frac{8}{3}\iint _{\partial H_-}R^3 \text{d}s \end{array}\right]. \end{eqnarray*}

(49)

By setting p = 1, p = 3 and p = −1 in eq. (27), and using exact solutions from the constant, linear, quadratic and cubic cases, singularity-free exact solutions for ϕ and g are obtained.

3 VERIFICATION

Following numerical experiments were carried out in a right-handed Cartesian coordinate system in which the positive z-axis is downward. As already explained, the observation site is required to be located at the origin of the Cartesian coordinate system. The simple rectangular prismatic body as shown in Fig. 2 was used to verify the accuracy of our newly developed formula. The dimension of the prismatic body is |$x=[10\, \text{km}, 20\, \text{km}]$|⁠, |$y=[10\, \text{km}, 20\, \text{km}]$| and |$z=[0\, \text{km}, 8\, \text{km}]$|⁠. It has a depth-dependent residual density contrast which is a polynomial function of cubic order (Garcia-Abdeslem 2005):

\begin{eqnarray*}\lambda (z) = -747.7 + 203.435 z - 26.764 z^{2} + 1.4247 z^{3}, \end{eqnarray*}

(50)

where λ(z) is in kgm⁻³ and z is given in a unit of km. The observation profile is the same as that of Garcia-Abdeslem’s (2005).

Figure 2.

Garcia-Abdeslem’s (2005) rectangular prism model. Dashed lines in red colour represent two observation profiles. Profile B is located right on the top surface of the prism, while profile A has a vertical offset of 15 cm from this top surface. Two observation sites (S₁ and S₂) of profile B are on edges of the top surface (with z = 0) where the gravity fields in Garcia-Abdeslem’s (2005) formulae are singular.

We computed the vertical gravity fields at 61 equally spaced observation sites along profile A from |$x=0\, \text{km}$| to |$x=30\, \text{km}$| (⁠|$y=15\, \text{km}$|⁠, |$z=-0.15\, \text{m}$|⁠) (as shown in Fig. 2). Since there is a small vertical offset between profile A and the top surface of the prism, there is no singularity for the computed gravity fields. The results from our singularity-free method, Garcia-Abdeslem (2005)’s method and the equivalent point mass approach (Holstein et al.2007) are compared in Table (2) for the vertical gravity field. The contributions to the vertical gravity field by each term of the polynomial density contrast in eq. (50) are shown in Fig. 3. Excellent agreement is obtained between our closed-form solution and Garcia-Abdeslem (2005)’s method for all four orders (constant, linear, quadratic and cubic orders). Due to the close distance from the observation site to the prism, the equivalent point mass approach cannot produce good results. Note that the coefficients of density contrasts for constant and quadratic orders are negative, so that the vertical gravity fields of these two terms are negative, as shown in Figs 3(a) and (c). As shown in Fig. 3e, the relative errors referring to Garcia-Abdeslem’s (2005) solutions are less than 10⁻¹¹% which verifies the accuracy of our newly developed algorithm.

Figure 3.

Comparisons of gravity anomalies computed by Garcia-Abdeslem’s (2005) analytical expression and our formula for a prismatic body (Fig. 2) with polynomial density contrasts of (a) constant order, (b) linear order, (c) quadratic order and (d) cubic order; (e) relative errors with regard to Garcia-Abdeslem’s (2005) solutions. The observation sites are located at |$y=15\, \text{km}$| and |$z=-0.15\, \text{m}$|⁠, thus avoiding singularities. Excellent agreements are obtained for all orders.

Table 2.

Comparison of solutions for the vertical gravity field computed by our analytical formula, Garcia-Abdeslem’s (2005) formula and the equivalent point mass approach for the prismatic body as shown in Fig. 2. The profile is located above the prism, that is, at |$y=15\, \text{km}$| and |$z=-0.15\, \text{m}$|⁠, thus avoiding singularities. Excellent agreement between the results by our analytical formula and Garcia-Abdeslem’s (2005) formula is obtained.

	g (mGal) along profile A (⁠|$y=15\, \text{km}$| and |$z=-0.15\, \text{m}$|⁠)
|$x\, \text{(km)}$|	Our new solutions	Garcia-Abdeslem’s solutions	Equivalent point mass solutions
0	−1.41666286151468E+00	−1.41666286151481E+00	−1.84108618104241E+00
1	−1.73422227639846E+00	−1.73422227639855E+00	−2.23149257955557E+00
2	−2.15234264546948E+00	−2.15234264546958E+00	−2.73741793052409E+00
3	−2.71326520931830E+00	−2.71326520931837E+00	−3.40343974343350E+00
4	−3.48203673411649E+00	−3.48203673411646E+00	−4.29552616723902E+00
5	−4.56231001247872E+00	−4.56231001247878E+00	−5.51326925463090E+00
6	−6.12675013291898E+00	−6.12675013291993E+00	−7.21003321285478E+00
7	−8.48173961731087E+00	−8.48173961731099E+00	−9.62627299077403E+00
8	−1.22299031940987E+01	−1.22299031940998E+01	−1.31437976407972E+01
9	−1.88269449325808E+01	−1.88269449325800E+01	−1.83688884601663E+01
10	−3.62664287162128E+01	−3.62664287162135E+01	−2.62366314127874E+01
11	−5.36259783186966E+01	−5.36259783186970E+01	−3.80499008405665E+01
12	−5.99739916027339E+01	−5.99739916027357E+01	−5.51012744441785E+01
13	−6.32743074931516E+01	−6.32743074931500E+01	−7.70049860116706E+01
14	−6.49254770325312E+01	−6.49254770325319E+01	−9.82615386701048E+01
15	−6.54308299900759E+01	−6.54308299900765E+01	−1.07615318326458E+02

Profile B with a range of |$x=[0\, \text{km},30\, \text{km}]$| is located on the top surface of the prism at |$y=15\, \text{km}$| and |$z=0\, \text{km}$|⁠, as shown in Fig. (2). The formula presented by Garcia-Abdeslem (2005) is singular, if the observation point is located at an edge of the anomalous prism (i.e. at sites S₁ and S₂ in the case of profile B). However, our formulae are singularity-free and show a regular behaviour as demonstrated in Table (3). We observe that solutions computed by our analytical formula, Garcia-Abdeslem’s (2005) formula and a high-order Gaussian quadrature method (using 512 × 512 × 512 quadrature points) (Davis & Rabinowitz 1984) show excellent agreement with each other. Remarkably, our solution shows smooth behaviour even at the edges where the gravity signals are mathematically singular in the analytical expression presented by Garcia-Abdeslem (2005). If the observation site is close to the rectangular prism, compared to closed-form solutions, even the high-order Gaussian quadrature method cannot achieve high-precision results, demonstrating the superiority of our proposed closed-form solution.

Table 3.

Comparison of solutions for the vertical gravity field computed by our analytical formula, Garcia-Abdeslem’s (2005) formula and a high-order Gaussian quadrature method (using 512 × 512 × 512 quadrature points) (Davis & Rabinowitz 1984) for the prismatic body as shown in Fig. 2. Similar to Fig. 3, but for the singular case where profile B is located on the top of the prism, that is, at |$y=15\, \text{km}$| and |$z=0\, \text{km}$|⁠. At the edges of the prism, that is, at |$x=10\, \text{km}$| (site S₁), the analytical solution presented by Garcia-Abdeslem (2005) is singular, but our solution is completely singularity-free. Due to symmetry, solutions for only one half of profile B are shown (results for |$x=(15\, \text{km},30\, \text{km}]$| are omitted).

	g along profile B (⁠|$y=15\, \text{km}$| and |$z=0\, \text{km}$|⁠)
|$x\, \text{(km)}$|	Our new solutions (mGal)	Garcia-Abdeslem’s solutions (mGal)	Gauss quadrature solutions (mGal)
0	−1.41659381299933E+00	−1.41659381299899E+00	−1.41659381299908E+00
1	−1.73413869984550E+00	−1.73413869984593E+00	−1.73413869984550E+00
2	−2.15224028284275E+00	−2.15224028284243E+00	−2.15224028284317E+00
3	−2.71313815047598E+00	−2.71313815047617E+00	−2.71313815047556E+00
4	−3.48187657349074E+00	−3.48187657349082E+00	−3.48187657349141E+00
5	−4.56210442191832E+00	−4.56210442191851E+00	−4.56210442191737E+00
6	−6.12648027897631E+00	−6.12648027897630E+00	−6.12648027897633E+00
7	−8.48137503186591E+00	−8.48137503186615E+00	−8.48137503186620E+00
8	−1.22293900434146E+01	−1.22293900434145E+01	−1.22293900434128E+01
9	−1.88261712992561E+01	−1.88261712992562E+01	−1.88261712992559E+01
10	−3.62673071958274E+01	−	−3.61614752500613E+01
11	−5.36285124167034E+01	−5.36285124167031E+01	−5.34172694362954E+01
12	−5.99762760875470E+01	−5.99762760875471E+01	−5.97502011647203E+01
13	−6.32764627789341E+01	−6.32764627789341E+01	−6.30471992541654E+01
14	−6.49275676133833E+01	−6.49275676133832E+01	−6.46962805858505E+01
15	−6.54329007321985E+01	−6.54329007321983E+01	−6.51862869517059E+01

In order to test the accuracy of our formula for high-order cases, the following mass density contrast is tested:

\begin{eqnarray*}\lambda (z) = { z^{4}}, \end{eqnarray*}

(51)

where λ(z) is in |$\, \text{kg} \, \text{m}^{-3}$| and z is in |$\, \text{km}$|⁠. Right above the top surface of the prismatic body (as shown in Fig. 2), a measuring grid is arranged with coordinates |$x=[0\, \text{km},30\, \text{km}]$| and |$y=[0\, \text{km},30\, \text{km}]$| and |$z=0\, \text{km}$|⁠. This model is tested by our new formula and the high-order Gaussian quadrature rule presented in Davis & Rabinowitz (1984). The comparison is shown in Fig. 4. Since a large number of 512 × 512 × 512 quadrature points is used, the Gaussian quadrature rule can offer a reasonably reliable gravity anomaly in this singular case. The relative absolute errors between these two results are less than 10⁻¹⁰%, which is shown in Fig. 4(c). As suggested by Fig. 4, the excellent agreement between these two approaches verifies the accuracy of our new solution in cases of high-order polynomial density contrasts.

Figure 4.

Comparison of gravity anomalies computed by (a) our formula and (b) the high-order Gauss quadrature rule (512 × 512 × 512 quadrature points) for the prismatic body in Fig. 2 but with a quartic polynomial density contrast. The observation sites are uniformly distributed on the top surface of the prism at |$z=0\, \text{km}$|⁠, |$x=[0\, \text{km},30\, \text{km}]$| and |$y=[0\, \text{km},30\, \text{km}]$|⁠. The maximum absolute relative error (c) between these two approaches is less than |$10^{-10}\,\%$|⁠.

4 ASSESSMENT OF NUMERICAL STABILITY

To test the stability of our newly derived formulae, we performed a third test using the same prismatic mass body with the residual density contrast given by eq. (50). Here, the approach to test numerical stability follows Holstein’s (2003) idea. The prism has a size of |$\alpha =4\cdot\sqrt{66}\, \text{km}$| (the length of a diagonal of the prism). We define δ as the site-to-target distance and a dimensionless target distance 1/γ as (Holstein 2003; Holstein et al.2013)

\begin{eqnarray*}1/\gamma = \frac{\delta }{\alpha }. \end{eqnarray*}

(52)

Here, we vary 1/γ from 0.86 to 48 759, which means δ becomes progressively larger. Due to the limitations of machine precision, an analytical formula may lead to invalid gravity anomalies when 1/γ is beyond a certain critical value. Therefore, we intend to determine this critical value 1/γ, below which our proposed formula can be safely applied to calculate high-precision vertical gravity fields.

The vertical gravity fields for different 1/γ are calculated by our formula, the high-order Gaussian quadrature rule with 512 × 512 × 512 quadrature points (Davis & Rabinowitz 1984), and the equivalent point mass approach (Holstein 2003). The results of these three approaches are denoted by g, g_gauss and g_point and shown in Fig. 5a. The computing times of our method and the high-order Gaussian quadrature method were 0.0018 s and 85.469 s for 25 observation sites, respectively. These computations were run on a computer with a Intel core i5 processor at 3.30 GHz using four cores. Two kinds of numerical relative errors (Fig.5b) are calculated. One is |g − g_gauss|/|g_gauss|, which measures the differences between our closed-form solutions and the high-order Gaussian quadrature solutions. Another is |g − g_point|/|g_point|, which measures the differences between our results and those of the equivalent point mass approach. In addition, absolute errors between our results and those by the high-order Gaussian quadrature rule are shown in Fig. 5(c).

Figure 5.

Comparisons of results computed by our formula, the high-order Gaussian quadrature rule (512 × 512 × 512 quadrature points) and the equivalent point mass approach when the dimensionless target distance 1/γ increases. Panel (a) shows the computed values of the vertical gravity fields, panel (b) is for their relative errors and panel (c) shows the absolute error between our results and those calculated by the high-order Gaussian quadrature rule.

In Fig. 5(b), when the observation site is located close to the prism, that is, at 1/γ ≈ 1, |g − g_gauss|/|g_gauss| is less than |$10^{-10}\,\%$|⁠. But |g − g_point|/|g_point| shows values larger than |$10^{0}\,\%$| as the prism cannot be approximated as an equivalent point mass in this case. When 1/γ increases, the relative error |g − g_gauss|/|g_gauss| increases over the whole range of 1/γ. From 1/γ ≈ 1 to 1/γ ≈ 80, |g − g_point|/|g_point| shows a decreasing behaviour which means that the prism can be gradually well approximated as an equivalent point mass when 1/γ increases. When 1/γ ≈ 80, |g − g_gauss|/|g_gauss| and |g − g_point|/|g_point| finally converge (see Fig. 5a). This is because the gravity signals are close to zero at this critically far distance. Due to the limited machine precision of the computer (Intel core i5 CPU 3.3GHz and 16GB RAM), our formula has difficulty to produce accurate gravity signals near zero. It is interesting to observe that beyond this critical distance, solutions of the Gaussian quadrature rules and the equivalent point mass still match well with each other (see Fig. 5a). This is because solutions of Gaussian quadrature rules can be considered as simple summations of individual contributions from a set of point masses. The amplitudes of these point-mass contributions are generally at the same level. Therefore, there are no additions of small numbers to large numbers which leads to the precision loss in our analytical formulae. Factually, the numerical quantities generated in our solution grow in magnitude with increasing 1/γ, while the signal is decreasing. Our solution requires these large quantities to cancel under summation. With limited machine precision ε (approximately 10⁻¹⁶ in our computer), however, the representation error of a number W is of order Wε. When W is sufficiently large, Wε can exceed the signal that is sought, in which case no significance is left in the summed result.

This precision issue exists for all analytical solutions of gravity anomalies (Holstein 2003; Holstein et al.2013). Readers should pay attention to this limit of our new formula, when they deal with large values of 1/γ. However, we should not worry about this issue in exploration gravity problems. As shown in Fig. 5(b), when 1/γ reaches a critical value of 180 (the intersection between the 100% error line and the fitted error growth line), the residuals between our solutions and those calculated by the Gaussian quadrature rule are around 10⁻⁵ mGal (Fig. 5c) which is smaller than the amplitude of the minimum detectable gravity anomaly in exploration geophysics (about |$5\cdot 10^{-3}\, \text{mGal}$| for the typical CG-5 AutoGrav gravity meter, Reudink et al.2014). As |$\alpha =4\cdot \sqrt{66}\, \text{km}$|⁠, using eq. (52), the observation site can be located at a distance of around |$5849\, \text{km}$| to the target mass. If higher order polynomial formulae are used for the density contrast, then the ”horizon” will reduce compared to that found here.

5 CONCLUSIONS

New formulae were developed for the gravity potential and the vertical gravity field of a polyhedral prism. The polyhedral prism is a vertical prism with a polygonal cross-section. In this polyhedral prism, the density contrast is assumed to vary with depth as a polynomial function. The orders of the polynomial density contrast functions can be arbitrary, such as constant, linear, quadratic, cubic, etc. Using our newly developed formulae, exact solutions are obtained for arbitrary order.

No new singularities exist in our formulation of polynomial density contrast anomaly, over those that already exist in the constant density case. These are known to be removable. Thus, our algorithm can be safely used for all observation points placed in, on or outside the prism. Considering a rectangular prism, we have demonstrated the accuracy of our new approach by comparing our solution to those of previously published analytic solutions and of Gaussian quadrature rules. Compared to our new approach, these previously published methods have the drawback that they do not have singularity removal implemented and, hence, may deliver less accurate solutions, if the observation sites are located in or on the border of an anomalous mass body.

Using our new formulae, interested readers may be able to derive closed-form solutions for the gravity gradient tensor, the magnetic field and even the magnetic gradient tensor of the polyhedral prism. Although we have considered only depth-dependent density functions, the identities and formulae presented in this study may be further used to derive exact solutions for a polyhedral body with both horizontal and vertical polynomial density contrasts of arbitrary orders.

ACKNOWLEDGEMENTS

This work was financially supported by the National Science Foundation of China (41574120, 41474103 and 41704100), the National Basic Research Program of China (973 - 2015CB060200), the Natural Science Foundation of Hunan Province of China (2018JJ3643), the Innovation-Driven Project of Central South University (2018CX042), an award for outstanding young scientists by Central South University (Lieying program 2013) and the Pioneer Hundred Talents Program, Chinese Academy of Sciences. The authors thank the editor Bert Vermeersen and the reviewer Jianzhong Zhang for their comments to improve the draft. Authors want to give deep thanks to the reviewer H. Holstein for his word-by-word and symbol-by-symbol checking which significantly improved the clarity and quality of the manuscript.

APPENDIX A

FORMULA FOR o_i

Here, we present the formula to evaluate point |$\mathbf {o}_{i}$|⁠. We assume that two connected edges of facet ∂H_i are defined by C_i1 and C_i2. Two vertices of C_i1 are denoted by |$\mathbf {v}_{0i1}$| and |$\mathbf {v}_{1i1}$|⁠, and vertices of C_i2 are represented by |$\mathbf {v}_{0i2}$| and |$\mathbf {v}_{1i2}$|⁠, with |$\mathbf {v}_{1i1}$| and |$\mathbf {v}_{0i2}$| being the same point, as shown in Fig. A1. Therefore, the outward normal unit vector |$\hat{\mathbf {n}}_{i}$| of surface ∂H_i can be expressed as

\begin{eqnarray*}\hat{\mathbf {n}}_{i}=\frac{\overrightarrow{\mathbf {v}_{0i1}\mathbf {v}_{1i1}}\times \overrightarrow{\mathbf {v}_{0i2}\mathbf {v}_{1i2}}}{\mid \overrightarrow{\mathbf {v}_{0i1}\mathbf {v}_{1i1}}\times \overrightarrow{\mathbf {v}_{0i2}\mathbf {v}_{1i2}}\mid }. \end{eqnarray*}

(A1)

Figure A1.

Illustration of the relationship between |$\mathbf {r}^{\prime }$| and facet ∂H_i.

As |$\mathbf {o}_i$| is the projection point of |$\mathbf {r}^{\prime }$| on facet ∂H_i, the vector pointing from |$\mathbf {r}^{\prime }$| to |$\mathbf {o}_{i}$| is perpendicular to ∂H_i, so that |$\overrightarrow{\mathbf {r}^{\prime }\mathbf {o}_{i}}\bot \overrightarrow{\mathbf {v}_{0i1}\mathbf {v}_{1i1}}$| and |$\overrightarrow{\mathbf {r}^{\prime }\mathbf {o}_{i}}\bot \overrightarrow{\mathbf {v}_{0i2}\mathbf {v}_{1i2}}$|⁠. Then, considering |$\mathbf {o}_i$| is on the plane ∂H_i, we have three independently nonlinearly related equations:

\begin{eqnarray*}\overrightarrow{\mathbf {r}^{\prime }\mathbf {o}_{i}}\cdot \overrightarrow{\mathbf {v}_{0i1}\mathbf {v}_{1i1}}=(\mathbf {r}^{\prime }-\mathbf {o}_{i})\cdot (\mathbf {v}_{0i1}-\mathbf {v}_{1i1})=0,\end{eqnarray*}

(A2)

\begin{eqnarray*}\overrightarrow{\mathbf {r}^{\prime }\mathbf {o}_{i}}\cdot \overrightarrow{\mathbf {v}_{0i2}\mathbf {v}_{1i2}}=(\mathbf {r}^{\prime }-\mathbf {o}_{i})\cdot (\mathbf {v}_{0i2}-\mathbf {v}_{1i2})=0,\end{eqnarray*}

(A3)

\begin{eqnarray*}\hat{\mathbf {n}}_{i}\cdot \left(\mathbf {o}_i-\mathbf {v}_{0i1}\right)=0. \end{eqnarray*}

(A4)

With the fact that |$\mathbf {r}^{\prime }=(0,0,0)$|⁠, assuming |$\mathbf {v}_{0i1}=(x_{0i1},y_{0i1},z_{0i1})$|⁠, |$\mathbf {v}_{1i1}=(x_{1i1},y_{1i1},z_{1i1})$|⁠, |$\mathbf {v}_{1i2}=(x_{1i2},y_{1i2},z_{1i2})$|⁠, by solving eqs (A2)–(A4), we obtain the coordinates of |$\mathbf {o}_{i}=(x,y,z)$|⁠:

\begin{eqnarray*}\left( \begin{array}{c}x \\y \\z \\\end{array} \right)= \left( \begin{array}{ccc}A & B & C \\x_{0i1}-x_{1i1} & y_{0i1}-y_{1i1} & z_{0i1}-z_{1i1} \\x_{1i2}-x_{1i1} & y_{1i2}-y_{1i1} & z_{1i2}-z_{1i1} \\\end{array} \right)^{-1} \left( \begin{array}{c}D \\0 \\0 \\\end{array} \right) , \end{eqnarray*}

(A5)

where

\begin{eqnarray*}A&=&(y_{0i1}-y_{1i1})*(z_{1i2}-z_{1i1})-(z_{0i1}-z_{1i1})*(y_{1i2}-y_{1i1}),\nonumber \\B&=&-(x_{0i1}-x_{1i1})*(z_{1i2}-z_{1i1})+(z_{0i1}-z_{1i1})*(x_{1i2}-x_{1i1}),\nonumber \\C&=&(x_{0i1}-x_{1i1})*(y_{1i2}-y_{1i1})-(y_{0i1}-y_{1i1})*(x_{1i2}-x_{1i1}),\nonumber \\D&=&A*x_{0i1}+B*y_{0i1}+C*z_{0i1}.\nonumber\end{eqnarray*}