1 Introduction

In spherical approximation, a uniform tesseroid is the mass body with the uniform density \(\rho ^{\mathrm{tess}}\) bound by three surface pairs defined in spherical coordinates by spherical longitudes [\(\lambda _1\), \(\lambda _2\)], latitudes [\(\varphi _1\), \(\varphi _2\)] and radii [\(r_1\), \(r_2\)] (see Fig. 1). Among the common mass bodies used in spatial domain gravity forward modeling (e.g., rectangular prism, vertical cylindrical prism, mass-lines, mass-layers, point mass and polyhedron), the tesseroid is a basic mass body well-suited for global calculations of gravitational terrain and crustal effects at satellite height because of its computational efficiency (Heck and Seitz 2007; Wild-Pfeiffer 2008; Tsoulis et al. 2009; Grombein et al. 2013; Uieda and Barbosa 2017). Furthermore, due to its characterization by spherical coordinates, the tesseroid can be considered as a commonly used mass body when dealing with global digital elevation or crustal models registered with respect to spherical latitude and longitude [e.g., ETOPO1 (Amante and Eakins 2009), Earth2014 (Hirt and Rexer 2015), CRUST1.0 (Laske et al. 2013)].

Fig. 1
figure 1

Illustration of a uniformly magnetized tesseroid with uniform density \(\rho ^{\mathrm{tess}}\) and magnetization \(M^{\mathrm{tess}}\) in spherical coordinates [modified from Kuhn (2003)]. [\(\lambda _1\), \(\lambda _2\)], [\(\varphi _1\), \(\varphi _2\)] and [\(r_1\), \(r_2\)] define the bounding surfaces of the tesseroid in the directions of spherical longitude, latitude and radius. \(P (\lambda , \varphi , r)\) and \(S (\lambda ', \varphi ', r')\) are computation and source points, respectively. Derivatives of the gravitational and magnetic potential are given in a local north-orientated reference frame (LNOF, i.e., x, y, z)

Full size image

Using tesseroids in magnetic forward modeling, Asgharzadeh et al. (2008) derived the magnetic potential (MP), its first-order derivatives (i.e., magnetic vectors MV) and second-order derivatives (i.e., magnetic gradient tensor, MGT) of a tesseroid in terms of spherical integral kernels and evaluated these magnetic effects with the 3D Gauss–Legendre Quadrature (GLQ) method. They investigated crustal magnetic effects using tesseroids over the Iran region within the Middle East at a satellite height of 500 km. Du et al. (2015) compared different numerical approaches [e.g., 3D GLQ and second-order Taylor Series Expansion (TSE)] to compute the MP, MV and MGT of a tesseroid in terms of accuracy and computational efficiency. Their numerical experiments revealed that the accuracy of the MP is higher than for the MV and MGT, and errors of the MGT are largest but decrease more quickly, which is due to faster attenuation of higher-order gravity quantity. Baykiev et al. (2016) derived the MV of a tesseroid applying Poisson’s relation (Poisson 1826; Baranov 1957; Blakely 1995; Li and Oldenburg 1998; Roy 2008; Hinze et al. 2013; Baykiev et al. 2016) based on the gravity gradient tensor (GGT) using Cartesian kernels (Grombein et al. 2013). They presented the magnetic modeling software (Magnetic tesseroids) to calculate the Earth’s induced and remanent MV effects at global and regional scales. Based on derivations of gravitational curvatures [GC, the third-order derivatives of the gravitational potential (GP)] of a tesseroid, Deng and Shen (2018a) proposed the concept of magnetic curvatures (MC), the third-order derivatives of the MP, i.e., the change rates of the MGT with respect to the directions of the local north-orientated reference frame (LNOF, see Fig. 1) in magnetic field modeling. In the above part, a comprehensive overview of studies related to the MP of a uniform tesseroid and their derivatives up to the second order is provided, whereas the related formulas of the MC are not yet introduced. Thus, the existing body of literature is extended here by the MC, which is derived and numerically tested.

A complete summary of all spherical shell formulas in magnetic field was not yet available. Sections 3.2.1 and 3.2.2 in Blakely (1995, pp. 49–54) presented the potential and gravitational attraction of the spherical shell and solid sphere in gravity field, respectively. Section 5.7 put the problem 1 “Use Poisson’s relation to find the magnetic field of a uniformly magnetized spherical shell”, but did not offer the answer clearly in Blakely (1995, p. 98). Appendix B Subroutines in Blakely (1995, p. 370) presented the computation routine named sphere to calculate the gravitational attraction of a sphere in Table B.1, which was also in gravity field. Actually Table B.1 did not give the codes for the spherical shell directly especially for the magnetic field. The GP, GV and GGT of the spherical shell for the gravity field were applied as the reference values with respect to the discretized tesseroids in Grombein et al. (2013), the same as the GC of the spherical shell in Deng and Shen (2018a). The shell formulas of the MP, MV, MGT and MC can be references for the discretized tesseroids in magnetic field, which were not clearly presented in the related studies of Blakely (1995), Du et al. (2015) and Baykiev et al. (2016). In this contribution, the general nth order derivatives of the gravitational potential of the spherical shell are presented. The magnetic field (including MP, MV, MGT and MC) of a uniformly magnetized spherical shell is also presented for the sake of completeness, which answers the problem 1 in Blakely (1995, p. 98).

In recent years, satellite missions are increasingly used to measure magnetic quantities in space to investigate the Earth’s magnetic field. For example, the MGT was obtained by numerical differentiation of the measured MV by Swarm satellite mission to investigate Earth’s magnetic field (Friis-Christensen et al. 2006, 2008; Olsen et al. 2010; Sabaka et al. 2018), and the CHAMP satellite measured the MV components to investigate magnetic effects of the Earth’s lithosphere (Maus et al. 2006, 2007, 2008; Sebera et al. 2019). Further studies combined Swarm and CHAMP magnetic measurements to model the Earth’s lithospheric magnetic field (Olsen et al. 2017; Vervelidou et al. 2018). Recently, Bandyopadhyay et al. (2019) reported the first measurements of magnetic field curvature by using in situ data from the Magnetospheric Multi-Scale mission (MMS). Although there is still a very long way to go to directly measure the MGT and MC in space, this study indicates the possible application of the MC for improved modeling of geophysical phenomena that leave a signal in the Earth’s magnetic field.

Poisson’s relation (Blakely 1995; Roy 2008; Hinze et al. 2013) provides a link between the nth-order derivatives of the MP and the \((n+1)\)th-order derivatives of the GP. Therefore, if expressions for higher-order derivatives of the GP are known, they can be used to derive expressions for the MP and its derivatives. For instance, by using Poisson’s relation, Blakely (1995) derived expressions for the MP of a solid sphere, infinite slab and horizontal cylinder from known expressions of gravity vectors (GV). Baykiev et al. (2016) derived the MV of a uniform tesseroid based on expressions of the GGT given by Grombein et al. (2013). This demonstrates the usefulness of Poisson’s relation for deriving expressions for magnetic quantities, where expressions of higher-order derivatives of the GP are already known.

Various expressions for the GP, GV, GGT and GC of a uniform tesseroid have been derived and applied to evaluate gravitational effects of the Earth’s topographic and crustal masses. Kuhn (2003) applied zero-order TSE approach to calculate the GP and radial GV components of a tesseroid. Heck and Seitz (2007) expanded the 3D TSE approach up to second-order formula to evaluate the GP and radial GV of a tesseroid. Asgharzadeh et al. (2007) calculated the terrain effects of the GP, GV and GGT of a tesseroid in spherical integral kernels with the 3D GLQ approach. Grombein et al. (2013) applied the 3D TSE second-order formula to evaluate the GP, GV and GGT of a tesseroid in terms of the Cartesian kernel functions. Deng et al. (2016) corrected the related formulas in Heck and Seitz (2007) and Grombein et al. (2013) to the correct form (i.e., from denominator \((i+j+k)!\) to i!j!k!) of the general 3D TSE formula for the GP and radial GV components of a tesseroid. Deng and Shen (2018a) derived the formulas of the GC of a tesseroid using spherical integral kernels, and expanded the 3D TSE approach up to sixth-order formula to evaluate the GP, GV, GGT and GC of a tesseroid. Deng and Shen (2018b) derived the formulas of the GC of a tesseroid using the Cartesian kernels, and compared the 2D/3D TSE, GLQ and Newton–Cotes Quadrature (NCQ) methods to validate the GP, GV, GGT and GC of a tesseroid with respect to the reference values of the spherical shell. The advantage of using Cartesian kernels instead of spherical kernels is that the former avoids the polar singularity for the GC formulas (Deng and Shen 2018b).

The near zone (or very-near-area) problem was widely investigated for gravitational functional of a tesseroid. It occurs when the computation point is close to the attracting mass in arbitrary direction using the mathematical explanation 1/l with \(l\rightarrow 0\). The contents of near zone problem for the GP, GV and GGT of the relevant literature (Heck and Seitz 2007; Tsoulis et al. 2009; Li et al. 2011; Grombein et al. 2013; Shen and Deng 2016; Uieda et al. 2016) were reviewed in Deng and Shen (2018b). The near zone problem appears for the GC both using the spherical kernels and Cartesian kernels (Deng and Shen 2018a, b). In magnetic field modeling, Du et al. (2015) demonstrated that the near zone problem also occurred for the MP, MV and MGT of a tesseroid.

Based on the proposed concept of the MC (Deng and Shen 2018a), this contribution derives and analyses the MC expressions of a uniformly magnetized tesseroid and spherical shell. The latter has been included here for numerical validation purposes. The fourth-order derivatives of the GP of the uniform tesseroid and spherical shell are derived first and subsequently used to derive expressions for the MC using Poisson’s relation. The spherical shell is used for validation of the newly derived expressions for the MC, then these expressions have been analyzed numerically with particular focus on the near zone and polar singularity problems.

This paper is organized as follows. Section 2 provides theoretical aspects, where the general principle of Poisson’s relation is presented in Sect. 2.1; expressions for the nth-order derivatives of the GP of the uniform tesseroid using the Cartesian kernels and spherical shell are derived in Sects. 2.2 and 2.3, respectively; Sect. 2.4 provides the general expressions of the MC of a tesseroid in terms of the Cartesian kernel functions; and the magnetic effects up to MC of a uniformly magnetized spherical shell are derived in Sect. 2.5. In Sect. 3 expressions for the MP, MV, MGT and the newly derived expressions for the MC are analyzed numerically. Section 3.1 provides a description of the numerical experiment where the spherical shell is discretized by tesseroids. In Sect. 3.2 the computational precision and computation time of the MC component are investigated based on the number of nodes used in the 3D GLQ. The near zone and polar singularity problems of the MP, MV, MGT and MC have been investigated in Sects. 3.3 and 3.4, respectively. Finally, the main conclusions and an outlook on future studies are provided in Sect. 4.

2 Theoretical Aspects

2.1 Poisson’s Relation Between Magnetic and Gravimetric Functional

Poisson’s relation provides a link between functional of the magnetic and gravity fields (Blakely 1995; Roy 2008; Hinze et al. 2013). According to Poisson’s relation, nth-order derivatives of magnetic effects can be obtained from \((n+1)\)th-order derivatives of gravity effects, and vice versa. The general expression of Poisson’s relation from gravity to magnetic functional is given as (Blakely 1995, pp. 91–93):

$$\begin{aligned} \frac{\partial ^{n}\ {\varPsi }}{\partial ^{n} \alpha } =- \frac{\mu _0\ M}{4\pi G \ \rho }\left[ \sum \frac{\partial ^{n+1} V}{\partial ^{n+1} \alpha } \right] \end{aligned}$$
(1)

where \({\varPsi }\) and V are the magnetic and gravitational potentials, respectively. \(({\partial ^{n} (.)} / {\partial ^{n} \alpha })\) denotes the nth-order derivatives of functional (.) (i.e., \({\varPsi }\) or V), where \(\alpha \in \{x, y, z\}\) indicates the axes (i.e., directions) of the local north-orientated reference frame (LNOF). G is the gravitational constant (i.e., \(G = 6.673 \times 10^{-11}\,\hbox {m}^{3}\, \hbox {kg}^{-1}\, \hbox {s}^{-2}\)), and \(\mu _0\) is the magnetic permeability of free space (i.e., \(\mu _0 = 4 \pi \times 10^{-7} \, \hbox {N} \, \hbox {A}^{-2}\)). \(\rho \) and M are the uniform density and magnetization of the uniform mass body, respectively. \(\sum ({\partial ^{n+1} V} / {\partial ^{n+1} \alpha })\) represents the sum of all \((n+1)\)th-order derivatives of gravitational potential in the three directions indicated by \(\alpha \). It should be noted that Eq. (1) can be used for any order (\(n\ge 0\)), but the condition of uniform magnetization and uniform density of the mass body is required. While for the real Earth’s field, sources of the magnetic field never have magnetization distributions in exact proportion to their density distributions (Blakely 1995, pp. 92–93), the condition of Poisson’s relation can be applied in local applications (Hotine 1969, p. 184). In this case density distributions can be assumed to be in direct proportion to their magnetization distributions. This condition can be used to derive pseudo-gravity from magnetic anomalies, which can aid in the interpretation of magnetic data (Baranov 1957; Blakely 1995, p. 93).

Using Eq. (1), the nth-order derivatives of magnetic effects can be obtained from the \((n+1)\)th-order derivatives of the corresponding gravity effects. For instance, the MP can be obtained from the GV [e.g., see MP expressions of a solid sphere, infinite slab and horizontal cylinder in Blakely (1995)], the MV can be obtained from the GGT [e.g., see MV expressions of a tesseroid in Baykiev et al. (2016)]. Following Poisson’s relation, the derivation of expressions for the MC requires the fourth-order derivatives of the GP, which will be provided in Sect. 2.2 for a uniform tesseroid and for completeness in Sect. 2.3 for a uniform spherical shell.

2.2 nth-Order Derivatives of the GP of a Tesseroid Using the Cartesian Kernels

The first- and second-order derivatives of the GP (i.e., GV and GGT) of a tesseroid in terms of the Cartesian kernel functions were derived in Grombein et al. (2013), and third-order derivatives (i.e., GC) were derived in Deng and Shen (2018b). Based on the same principle, the general expression of the nth-order derivatives of the GP (\(V_{(i, j, k)}\) with \(i+j+k=n\)) of a tesseroid using the Cartesian kernels with respect to the LNOF is expressed as (Grombein et al. 2013; Deng and Shen 2018b):

$$\begin{aligned} V_{(i, j, k)}= & {} G \rho ^{\mathrm{tess}} \int _{r_1}^{r_2} \int _{\varphi _1}^{\varphi _2} \int _{\lambda _1}^{\lambda _2} \kappa \frac{\partial ^{i+j+k} (\frac{1}{\ell }) }{\partial {\varDelta _{x}}^{i}\ \partial {\varDelta _{y}}^{j}\ \partial {\varDelta _{z}}^{k}} \hbox {d}\lambda ' \hbox {d}\varphi ' \hbox {d}r' \end{aligned}$$
(2)
$$\begin{aligned} \ell= & {} \sqrt{r'^2+r^2-2r'r\cos \psi } = \sqrt{\varDelta _x^2 + \varDelta _y^2 + \varDelta _z^2} \end{aligned}$$
(3)
$$\begin{aligned} \varDelta _x= & {} r'\left( \cos \varphi \sin \varphi '-\sin \varphi \cos \varphi ' \cos (\lambda '-\lambda ) \right) \end{aligned}$$
(4)
$$\begin{aligned} \varDelta _y= & {} r' \cos \varphi '\sin (\lambda '-\lambda ) \end{aligned}$$
(5)
$$\begin{aligned} \varDelta _z= & {} r' \cos \psi -r \end{aligned}$$
(6)
$$\begin{aligned} \cos \psi= & {} \sin \varphi \sin \varphi '+\cos \varphi \cos \varphi ' \cos \left( \lambda '-\lambda \right) \end{aligned}$$
(7)
$$\begin{aligned} \kappa= & {} r'^2 \cos \varphi ' \end{aligned}$$
(8)

where \(\rho ^{\mathrm{tess}}\) is the constant density of the uniform tesseroid (see Fig. 1), and \(i,\ j,\ k \in \{0, 1, 2, 3, \cdots ,n\}\) (with \(i+j+k=n\)) indicate the number of derivatives in x, y and z directions, respectively. For the LNOF, the x-axis, y-axis and z-axis point to the north, the east and the zenith direction, respectively. \(P (\lambda ,\varphi ,r )\) and \(S (\lambda ',\varphi ',r' )\) are the computation (or field) point and integration (or running) point, respectively. For detailed derivations of the partial derivatives in Eq. (2) refer to Grombein et al. (2013, Eqs. 15–20).

Based on Eqs. (2)–(8), expressions for the fourth-order gravitational components of a tesseroid using the Cartesian kernels are provided in Appendix 1. Moreover, the same principle can be used for any higher derivatives, i.e., \(n \ge 5\).

2.3 nth-Order Derivatives of the GP of a Uniform Spherical Shell

In this study the analytical formulas for the gravitational potential and its higher-order derivatives induced by a uniform spherical shell are used to validate the gravitational functional of a uniform tesseroid (cf. Sect. 3). For the sake of completeness the nth-order derivatives of the GP of a uniform spherical shell are provided here. The gravitational potential of a uniform spherical shell with constant density \(\rho ^{\mathrm{sh}}\) can be obtained by the differences of the gravitational potential of two solid spheres with the same center O and different radii \(R_1\) and \(R_2\) (\(R_2\) > \(R_1\)).

For a computation point located in the outer space of the shell (e.g., \(r \ge R_2\)) the gravitational potential of the spherical shell (\(V^{\mathrm{sh}}\)) is given as (Blakely 1995; Vanícĕk et al. 2001):

$$\begin{aligned} V^{\mathrm{sh}}=\frac{4}{3}\pi G \rho ^{\mathrm{sh}} \left( R_2^3-R_1^3\right) \frac{1}{r} \end{aligned}$$
(9)

where \(r=\sqrt{x^2+y^2+z^2}\) is the geocentric distance of the computation point P and the thickness of the sphere is given as \(h' = R_2 - R_1\) in Fig. 2.

Based on Eq. (9), expressions for the nth-order gravitational components of a uniform spherical shell are obtained as:

$$\begin{aligned} \begin{aligned} V^{\mathrm{sh}}_{(i,j,k)}= \frac{\partial ^{i+j+k} V^{\mathrm{sh}}}{\partial x^{i} \ \partial y^{j} \ \partial z^{k} } = \frac{4}{3}\pi G \rho ^{\mathrm{sh}} \left( R_2^3-R_1^3\right) \frac{\partial ^{i+j+k} (\frac{1}{r})}{\partial x^{i} \ \partial y^{j} \ \partial z^{k} } \end{aligned} \end{aligned}$$
(10)

where \(i, j, k \in \{0, 1, 2, 3, \cdots , n\}\) with \(i+j+k=n\).

Grombein et al. (2013) provided the first- and second-order gravitational components of the uniform spherical shell as the elements of the gravitational acceleration and the Marussi tensor. Deng and Shen (2018a) further presented the third-order gravitational components of the uniform spherical shell as the GC. Furthermore, the fourth-order gravitational components of the uniform spherical shell are provided in Appendix 2. Moreover, one can derive expressions for the nth-order (\(n\ge 5\)) gravitational components of the uniform spherical shell based on Eq.  (10).

Fig. 2
figure 2

Illustration of a uniformly magnetized spherical shell with uniform density \(\rho ^{\mathrm{sh}}\) and magnetization \(M^{\mathrm{sh}}\). O is the same center of the two spheres with radii \(R_1\) and \(R_2\), indicating the lower and upper bounds of the spherical shell. P is the computation point with radius r from the center O. The orientation of the lines in the hatching shell indicate the direction of the uniform magnetization

Full size image

2.4 Magnetic Curvatures of a Uniformly Magnetized Tesseroid

Du et al. (2015) derived the MV and MGT of a tesseroid from general integral equations of the magnetic potential given by Blakely (1995) instead of using Poisson’s relation. Baykiev et al. (2016) presented the MV of a tesseroid by applying Poisson’s relation using the GGT of Grombein et al. (2013). In this part, the MC of a tesseroid using the Cartesian kernels based on Poisson’s relation with the fourth-order derivatives of the gravitational potential of a tesseroid are derived.

Based on Poisson’s relation using the fourth-order gravitational components of a tesseroid (\(V_{(i,j,k)}\) with \(i,\ j,\ k \in \{0, 1, 2, 3, 4\}\) and \(i+j+k=4\) ) provided in Appendix 1, the general MC of a uniformly magnetized tesseroid (\(B_{\alpha \beta \gamma }\) with \(\alpha ,\ \beta ,\ \gamma \in \{x,\ y, \ z\}\)) are expressed in the LNOF as:

$$\begin{aligned} \begin{aligned} B_{\alpha \beta \gamma }=- \frac{\mu _0\ M^{\mathrm{tess}}}{4\pi G\ \rho ^{\mathrm{tess}}} \Big ( \sum V_{(i,j,k)}\Big ) \end{aligned} \end{aligned}$$
(11)

The summation (\(\sum \)) in Eq. (11) indicates the sum over the fourth-order gravitational derivatives of the GP of the tesseroid in the three directions x, y and z. The detailed Poisson’s relation between MC and fourth-order gravitational derivatives of the GP is provided in Table 6 in Appendix 3. After substituting the expressions of \(V_{(i,j,k)}\) from Appendix 1 in Table 4 into the expressions in Table 6, the MC of a uniformly magnetized tesseroid are provided in Table 7 in Appendix 3.

2.5 Magnetic Effects up to Magnetic Curvatures of a Uniformly Magnetized Spherical Shell

Using Poisson’s relation, the MP (\({\varPsi }^{\mathrm{sh}}\)), MV (\(B^{\mathrm{sh}}_{\alpha }\) with \(\alpha \in \{x,\ y,\ z\}\)) and MGT (\(B^{\mathrm{sh}}_{\alpha \beta }\) with \(\alpha ,\ \beta \in \{x,\ y,\ z\}\)) of a uniformly magnetized spherical shell are obtained from the expressions for the GV (\(V^{\mathrm{sh}}_{\alpha }\) with \(\alpha \in \{x,\ y,\ z\}\)), GGT (\(V^{\mathrm{sh}}_{\alpha \beta }\) with \(\alpha ,\ \beta \in \{x,\ y,\ z\}\)) and GC (\(V^{\mathrm{sh}}_{\alpha \beta \gamma }\) with \(\alpha ,\ \beta ,\ \gamma \in \{x,\ y,\ z\}\)) of the uniform spherical shell as:

$$\begin{aligned} {\varPsi }^{\mathrm{sh}}= & {} - \frac{\mu _0\ M^{\mathrm{sh}}}{4\pi G\ \rho ^{\mathrm{sh}}} \left( \sum V_{\alpha }^{\mathrm{sh}} \right) \end{aligned}$$
(12)
$$\begin{aligned} B^{\mathrm{sh}}_\alpha= & {} - \frac{\mu _0\ M^{\mathrm{sh}}}{4\pi G\ \rho ^{\mathrm{sh}}} \left( \sum V_{\alpha \beta }^{\mathrm{sh}} \right) \end{aligned}$$
(13)
$$\begin{aligned} B^{\mathrm{sh}}_{\alpha \beta }= & {} - \frac{\mu _0\ M^{\mathrm{sh}}}{4\pi G\ \rho ^{\mathrm{sh}}} \left( \sum V_{\alpha \beta \gamma }^{\mathrm{sh}} \right) \end{aligned}$$
(14)

Furthermore, the MC (\(B^{\mathrm{sh}}_{\alpha \beta \gamma }\) with \(\alpha ,\ \beta ,\ \gamma \in \{x,\ y,\ z\}\)) of a uniformly magnetized spherical shell are obtained from the expressions for the fourth-order gravitational components of the uniform spherical shell (\(V^{\mathrm{sh}}_{(i,j,k)}\) with \(i, j, k \in \{0, 1, 2, 3, 4\}\) and \(i+j+k=4\)) provided in Appendix 2 as:

$$\begin{aligned} \begin{aligned} B^{\mathrm{sh}}_{\alpha \beta \gamma }= - \frac{\mu _0\ M^{\mathrm{sh}}}{4\pi G\ \rho ^{\mathrm{sh}}} \left( \sum V_{(i,j,k)}^{\mathrm{sh}} \right) \end{aligned} \end{aligned}$$
(15)

The summation (\(\sum \)) in Eqs.  (12)–(15) indicate the sum of the first-order, second-order, third-order and fourth-order gravitational derivatives of the GP of the spherical shell in three directions x, y and z, respectively.

For the sake of completeness, the detailed expressions for Poisson’s relations for the MP, MV, MGT and MC are given in Table 8 in Appendix 4. After substituting the expressions of \(V^{\mathrm{sh}}_{\alpha }\) and \(V^{\mathrm{sh}}_{\alpha \beta }\) from Grombein et al. (2013), \(V^{\mathrm{sh}}_{\alpha \beta \gamma }\) from Deng and Shen (2018a), and \(V^{\mathrm{sh}}_{(i,j,k)}\) (\(i, j, k \in \{0, 1, 2, 3, 4\}\) and \(i+j+k=4\)) from Appendix 2 into expressions in Table 8, the detailed expressions for the MP, MV, MGT and MC of a uniformly magnetized spherical shell are derived and provided in Table 9 in Appendix 4. It should be noted that the closed analytical expressions for the MP, MV, MGT and MC in Table 9 are used to provide reference values for all numerical tests.

3 Numerical Investigations

3.1 Set-up of the Numerical Experiments

Following the test layout used by Grombein et al. (2013), Shen and Deng (2016), Deng and Shen (2018a, b), and Zhong et al. (2019), expressions for the MC given in Table 7 are used to derive magnetic effects of a uniform spherical shell discretized by uniform tesseroids. For completeness the same test layout is applied to the MP, MV and MGT. Reference values (\(X^{\mathrm{ref}}\)) for all magnetic quantities (i.e., MP, MV, MGT and MC) are obtained from the analytical solutions provided in Table 9 in Appendix 4. Absolute and relative approximation errors \(\delta X^{\mathrm{Abs}}\) and \(\delta X^{\mathrm{Rel}}\) are derived as:

$$\begin{aligned} \delta X^{\mathrm{Abs}}= & {} |\big (\sum X^{\mathrm{tess}}\big ) - X^{\mathrm{ref}}| \end{aligned}$$
(16)
$$\begin{aligned} \delta X^{\mathrm{Rel}}= & {} \frac{|\big (\sum X^{\mathrm{tess}}\big ) - X^{\mathrm{ref}}|}{|X^{\mathrm{ref}}|} \end{aligned}$$
(17)

where X represents any of the magnetic quantities MP, MV, MGT or MC. \(\sum X^{\mathrm{tess}}\) stands for the magnetic effects obtained through discretization of the spherical shell by tesseroids. Hereby, \(\sum X^{\mathrm{tess}}\) expresses the finite element nature of approximation, which is obtained by the sum of every individual effect of all tesseroids forming the entire spherical shell (see parameterization of the numerical experiments below).

The numerical parameters used for the uniformly magnetized spherical shell are listed in Table 1. The value for \(R_2\) (i.e., \(R_2=6{,}371{,}200\,\hbox {m}\)) is chosen as the radius of Earth’s magnetic reference sphere adopted by the International Geomagnetic Reference Field (IGRF) (Finlay et al. 2010; Thébault et al. 2015). For all numerical studies, the thickness of the uniformly magnetized spherical shell is \(h'=R_2-R_1=1, 000\,\hbox {m}\) to be consistent with the reference values used in previous studies (Grombein et al. 2013; Deng and Shen 2018a, b) and for the links between the GC (Deng and Shen 2018a, b) and MC in this paper.

The reference values of the MP, MV, MGT and MC of the uniformly magnetized spherical shell at satellite height \(h_P = 260{,} 000\,\hbox {m}\) are listed in Table 2. For the numerical calculations the computation point P is located on the polar axis. This assumption means that the magnetic field of a uniformly magnetized spherical shell has a magnetic monopole. It should be noted that the real magnetic field does not have monopoles. However, this is of no concern here as this layout is able to provide precise reference values for validation purposes.

Table 1 Numerical parameters of the uniformly magnetized spherical shell used for validation purpose
Full size table
Table 2 Reference values of the MP, MV, MGT and MC of a uniformly magnetized spherical shell characterized in Table 1 at satellite height \(h_P = 260{,} 000\,\hbox {m}\) (e.g., \(r_P=R_2+h_P\))
Full size table

In the following numerical experiments, the entire spherical shell is discretized by tesseroids with a size of \(1^{\circ }\times 1^{\circ }\times 1\, \hbox {km}\), i.e., the horizontal dimensions of the tesseroid are \(\varDelta \lambda =\lambda _2-\lambda _1=1^{\circ }\) and \(\varDelta \varphi =\varphi _2-\varphi _1=1^{\circ }\). The vertical dimension of the tesseroid is \(\varDelta r=r_2-r_1=1000\, \hbox {m}\).

Magnetic effects are obtained by the sum of individual effects from all tesseroids forming the entire spherical shell. The 3D GLQ approach has been used to numerically solve the triple integral in the MC components of a tesseroid (see Table 7 in Appendix 3), which was widely applied for the numerical evaluations in gravity field modeling (Ku 1977; Blakely 1995; Asgharzadeh et al. 2007; Wild-Pfeiffer 2008; Hirt et al. 2011; Hinze et al. 2013; Rexer and Hirt 2015; Roussel et al. 2015; Uieda et al. 2016; Deng and Shen 2017, 2018b, 2019; Asgharzadeh et al. 2018; Zhong et al. 2019) and magnetic field modeling (Ku 1977; Blakely 1995; Asgharzadeh et al. 2008; Hinze et al. 2013; Du et al. 2015; Baykiev et al. 2016; Deng et al. 2019). Detailed expressions for the 3D GLQ approach can be found in Deng and Shen (2018b).

To show the numerical error level of some MGT and MC components, the following Laplace parameters (or Laplacians) can be used:

$$\begin{aligned} \delta \varDelta L_1= & {} \delta B_{xx} + \delta B_{yy} + \delta B_{zz} \end{aligned}$$
(18)
$$\begin{aligned} \delta \varDelta L_2= & {} \delta B_{xxx} + \delta B_{yyx} + \delta B_{zzx} \end{aligned}$$
(19)
$$\begin{aligned} \delta \varDelta L_3= & {} \delta B_{xxy} + \delta B_{yyy} + \delta B_{zzy} \end{aligned}$$
(20)
$$\begin{aligned} \delta \varDelta L_4 = \delta B_{xxz} + \delta B_{yyz} + \delta B_{zzz} \end{aligned}$$
(21)

Theoretically, all Laplacians should be equal to zero. The effect of magnetization on the Laplacian is taken as \(M^{\mathrm{tess}} = 1\, \hbox {A} \, \hbox {m}^{-1}\) from Table 1. The reductions (\(\sim \) 15 to 16 orders of magnitude) indicate the precision level of double precision (cf. Sects. 3.3, 3.4), thus the conclusion can be made that the Laplace equations are practically zero considering a numerical precision of double precision. In the following numerical process, the magnitude of the Laplacian reveals the error level due to the discretization of the spherical shell by tesseroids.

3.2 Precision and Computation Time for the MC Component \(B_{zzz}\) Evaluated by 3D GLQ

Fig. 3
figure 3

Illustration of a relative approximation errors in \(\hbox {Log}_{10}\) scale and b computation times (in seconds) for the MC component \(B_{zzz}\) dependent on the number of 3D GLQ nodes. The x-axes in both panels indicate the value for n where the number of nodes is \(n\times n\times n\)

Full size image

In order to test the effectiveness of the 3D GLQ approach, the MC component \(B_{zzz}\) is chosen to analyze the precision and computation time in relation to the number of nodes (here \(2\times 2\times 2\) to \(8\times 8\times 8\)) used by the 3D GLQ algorithm. For this test the computation point P is set on the polar axis with a height of \(260{,}000\, \hbox {m}\) (i.e., satellite height) above the spherical shell. The height has been chosen to practically eliminate numerical problems due to the near zone problem (cf. Sect. 3.3). Computations have been performed for a uniformly magnetized spherical shell introduced in Sect. 3.1 (see Table 2 for the reference value). All calculations have been performed on a standard PC with a 2.7 GHz Intel Core i5 processor and 8 GB RAM. However, it should be noted that the actual computation time estimates highly depend on the actual hardware and software, thus estimates should largely be interpreted in a relative sense.

Relative approximation errors based on Eq. (17) and computation time using different 3D GLQ nodes are shown in Fig. 3a, b, respectively. As expected, Fig. 3a shows that the relative approximation error decreases with increased number of nodes used in the 3D GLQ algorithm, i.e., the computational precision increases with increased number of nodes. Furthermore, Fig. 3b shows an increase in computation time correlates with the number of nodes.

In order to balance computational precision with the required computation time, the number 3D GLQ nodes are chosen as (\(3\times 3\times 3\)) for all subsequent numerical calculations. According to Fig. 3, this choice results in a relative error of − 1.307 in \(\hbox {Log}_{10}\) scale, which is about 4.9% (i.e., this approximation is sufficient here to test near zone and singularity problems having relative errors with much larger magnitudes.) and a computation time of 8.3 s required to calculate \(B_{zzz}\) at a single computation point.

3.3 Study of the Near Zone Problem of the MP, MV, MGT and MC

Fig. 4
figure 4

Visualization of a relative approximation errors in \(\hbox {Log}_{10}\) scale for the MP (\(\delta {\varPsi }\)), MV (\(\delta B_x\), \(\delta B_y\) and \(\delta B_z\)), MGT (\(\delta B_{xx}\), \(\delta B_{xz}\), \(\delta B_{yy}\), \(\delta B_{yz}\) and \(\delta B_{zz}\)) and MC (\(\delta B_{xxx}\), \(\delta B_{xxy}\), \(\delta B_{xxz}\), \(\delta B_{yyx}\), \(\delta B_{yyy}\), \(\delta B_{yyz}\), \(\delta B_{zzx}\), \(\delta B_{zzy}\) and \(\delta B_{zzz}\)) dependent on the computation point height (in km); b absolute approximation errors for the Laplace parameters of the MGT (\(\delta \varDelta L_1\)) and MC (\(\delta \varDelta L_2\), \(\delta \varDelta L_3\) and \(\delta \varDelta L_4\)) dependent on the computation point height. Based on Table 2, the order of magnitudes of the MGT components are about − 17 or − 16. The order of magnitudes of the MC components are about − 23. In both cases the reduction is about 15–16 orders of magnitude (e.g., \(32-17 = 15\) and \(39-23 = 16\)), which indicates the level of numerical precision (double precision)

Full size image

The following test has been set-up to study the near zone problem for magnetic quantities. The near zone problem, which is also known as “very-near-area” or “innermost zone” problem, refers to the singularity of the integral kernels that cause numerical problems when evaluating very close to the source masses, which was mentioned for the GP and GV by Heck and Seitz (2007), GGT by Grombein et al. (2013) and GC by Deng and Shen (2018a, b). In order to investigate whether evident errors exist in the near zone for the MP, MV, MGT and MC, i.e., when the computation point P is near the source masses, absolute and relative approximation errors [cf. Eqs.  (16), (17)] are analyzed for various heights of the computation point P above the spherical shell. This test layout follows that of Deng and Shen (2018a, b) used to study the GP, GV, GGT and GC. Here we use different heights above the spherical shell to model different geocentric distances of the computation point, i.e., \(r_P = R_2 + h_P\). In this experiment, the range of the height \(h_P\) is \([0,\ 2000\, \hbox {km}]\) with an interval of 5 km. Without loss of generality as stated in Deng and Shen (2018a), the spherical latitude and longitude of the computation point P are set to \(\varphi =0^{\circ }\) and \(\lambda =0^{\circ }\), respectively, because of spherical symmetry.

The relative approximation errors for the MP (\(\delta {\varPsi }\)), MV (\(\delta B_x\), \(\delta B_y\) and \(\delta B_z\)), MGT (\(\delta B_{xx}\), \(\delta B_{xz}\), \(\delta B_{yy}\), \(\delta B_{yz}\) and \(\delta B_{zz}\)) and MC (\(\delta B_{xxx}\), \(\delta B_{xxy}\), \(\delta B_{xxz}\), \(\delta B_{yyx}\), \(\delta B_{yyy}\), \(\delta B_{yyz}\), \(\delta B_{zzx}\), \(\delta B_{zzy}\) and \(\delta B_{zzz}\)) in relation to the height \(h_P\) are shown in Fig. 4a. And the absolute approximation errors for the Laplace parameters for the MGT (\(\delta \varDelta L_1\)) and MC (\(\delta \varDelta L_2\), \(\delta \varDelta L_3\) and \(\delta \varDelta L_4\)) are illustrated in Fig. 4b.

The relative approximation errors for all magnetic quantities (i.e., MP, MV, MGT and MC) show a very similar behavior (cf. Fig. 4a). As the height \(h_P\) increases the relative approximation errors decline rapidly, which is referred here to the “rapid drop zone”. The break point near 600–700 km could be linked to the used discretization \(1^{\circ }\times 1^{\circ }\) for tesseroids. The almost flat behavior is an indication of the error level introduced by the chosen discretization, e.g., at 600–700 km, the discretization error is larger than the error introduced by the near zone problem. At a height between 600 and 800 km (depending on the magnetic quantity), the relative approximation errors become almost constant referred here to the “stable zone”. The magnitude in the “stable zone” seems to represent the numerical precision of the calculations used, i.e., double precision in this case. Apart from the drop off, it is mostly the magnitude of the relative approximation error at or close to \(h_P\) = 0 km that indicates the near zone problem, i.e., the computation point is located either on or very close to the source masses. These behaviors of the MP, MV, MGT and MC are expected as they “inherit” the same properties as the gravitational functional as shown in Deng and Shen (2018a, b), though not including the fourth-order derivatives of the GP required for the MC quantity.

The overall minimum and maximum values (i.e., range) of the relative approximation errors in \(\hbox {Log}_{10}\) scale are about \([-\,16, 9]\) for the MC, \([-\,15, 5]\) for the MGT, \([-\,15, 3]\) for the MV and \([-\,14, 0]\) for the MP. In the rapid drop zone, the curves of the MC share the same properties expressed by two curves that overlap several MC components, i.e., one curve contains the 7 MC components (\(\delta B_{xxx}\), \(\delta B_{xxz}\), \(\delta B_{yyy}\), \(\delta B_{yyz}\), \(\delta B_{zzx}\), \(\delta B_{zzy}\) and \(\delta B_{zzz}\)); while the other overlapped curve contains the 2 MC components (\(\delta B_{xxy}\) and \(\delta B_{yyx}\)). Similarly, the curves of the 5 MGT components (\(\delta B_{xx}\), \(\delta B_{xz}\), \(\delta B_{yy}\), \(\delta B_{yz}\) and \(\delta B_{zz}\)) overlap together as one curve in the rapid drop zone, the same for the 3 MV components (\(\delta B_x\), \(\delta B_y\) and \(\delta B_z\)). Among the magnetic quantities in Fig. 4a, the relative approximation errors of the MC are largest in the rapid drop zone under the same condition. This resembles the same behavior for gravitational functional, i.e., the near zone problem worsens for higher derivatives.

Figure 4b shows the variations of absolute approximation errors in \(\hbox {Log}_{10}\) scale for the MGT and MC Laplace parameters in relation to the height. With increased height, the absolute approximation errors decline steadily. Specifically, the minimum and maximum values in \(\hbox {Log}_{10}\) scale of the MGT Laplace parameter (\(\delta \varDelta L_1\)) and MC Laplace parameters (\(\delta \varDelta L_2\), \(\delta \varDelta L_3\) and \(\delta \varDelta L_4\)) are about \([-\,32, -\,25]\) and \([-\,39, -\,30]\), respectively, which are very close to zero, where all Laplace parameters should theoretically be equal to zero. Based on Table 2, the order of magnitudes of the MGT components in Laplacian are about − 17 or − 16. The magnitudes of the MC components in Laplacians are about − 23. In both cases the magnitude reduces by about 15–16 orders of magnitude (e.g., \(32-17 = 15\) and \(39-23 = 16\)), which indicates the level of numerical precision (double precision). This also indicates that large part of the near zone problem cancels out. Thus, the precision level of the MGT and MC Laplace parameters shows that the sum of the MC (\(\delta B_{xxx}\), \(\delta B_{yyx}\) and \(\delta B_{zzx}\); \(\delta B_{xxy}\), \(\delta B_{yyy}\) and \(\delta B_{zzy}\); \(\delta B_{xxz}\), \(\delta B_{yyz}\) and \(\delta B_{zzz}\)) in Eqs.  (19)–(21) satisfies the Laplace equations at proper precision, the same for the sum of the MGT (\(\delta B_{xx}\), \(\delta B_{yy}\) and \(\delta B_{zz}\)) in Eq. (18). Effectively, these also provide the same information as the relative approximation errors, i.e., large errors in the near zone and numerical precision (i.e., magnitude 14–15 digits) level further away. While having a larger magnitude in the near zone, the absolute errors in the Laplace parameters are still negligible (very close to zero) while errors in the magnetic quantities are very large.

3.4 Study of the Polar Singularity Problem of the MP, MV, MGT and MC

Fig. 5
figure 5

For computation points at satellite height (\(h_P\) = 260 km), visualization of a relative approximation errors in \(\hbox {Log}_{10}\) scale for the MP (\(\delta \varPsi \)), MV (\(\delta B_x\), \(\delta B_y\) and \(\delta B_z\)), MGT (\(\delta B_{xx}\), \(\delta B_{xz}\), \(\delta B_{yy}\), \(\delta B_{yz}\) and \(\delta B_{zz}\)) and MC (\(\delta B_{xxx}\), \(\delta B_{xxy}\), \(\delta B_{xxz}\), \(\delta B_{yyx}\), \(\delta B_{yyy}\), \(\delta B_{yyz}\), \(\delta B_{zzx}\), \(\delta B_{zzy}\) and \(\delta B_{zzz}\)) in relation to the latitude of the computation point; b absolute approximation errors for the Laplace parameters for the MGT (\(\delta \varDelta L_1\)) and MC (\(\delta \varDelta L_2\), \(\delta \varDelta L_3\) and \(\delta \varDelta L_4\)) in relation to the latitude of the computation point. Based on Table 2, the order of magnitudes of the MGT components of the Laplacians are about − 17 or − 16. The order of magnitudes of the MC components of the Laplacians are about − 23. The colors and styles of curves for different parameters are the same as used in Fig. 4

Full size image

The following test has been set-up to study the polar-singularity problem for magnetic quantities. The polar singularity problem was found for the GC component (\(\delta V_{xxz}\)) in spherical integral kernels in Deng and Shen (2018a). In Deng and Shen (2018b), the polar singularity problem could be avoided for the GC component (\(\delta V_{xxz}\)) using the Cartesian kernels. In this experiment, in order to study whether the polar singularity problem exists for the MC, the latitude of the computation point is gradually changed from the equator to the pole while the computation point is at satellite height \(h_P = 260\) km. Here only the satellite height is considered because this kind of effect cannot be studied at lower altitudes (i.e., \(h_P = 0\) or 1 km) due to a bias from the near zone problem as revealed in Sect. 3.3. Importantly, using other satellite heights (e.g., 300 km, 400 km or 500 km) have no effect on the main conclusions of the paper. For completeness, the polar singularity problem is also studied for the remaining MP, MV and MGT.

Following the same layout as in Grombein et al. (2013), Shen and Deng (2016) and Deng and Shen (2018a, b), without loss of generality (i.e., spherical symmetry) the longitude of the computation point P is set to \(\lambda =0^{\circ }\) and the geocentric radius according to the computation point height is selected (see Table 2). Because of the symmetry between the northern and southern hemispheres, the spherical latitude \(\varphi \) of the computation point varies from the equator to the north pole in \(1^{\circ }\) intervals, which is the same as used by Deng and Shen (2018a, b) to study the polar singularity problem of the GC.

The relative approximation errors of the MP, MV, MGT and MC in relation to the latitude are illustrated in Fig. 5a. In addition, the absolute approximation errors of the MGT and MC Laplace parameters for the satellite height are also shown in Fig. 5b.

Figure 5a reveals the variations of the relative approximation errors for the MP, MV, MGT and MC in relation to the latitude of the computation point at satellite height \(h_P = 260\,\hbox {km}\). It can be used to study the polar singularity problem as it is not or only minimally impacted by the near zone problem. At the polar region (\(\varphi \ge 85^{\circ }\)), relative approximation errors rapidly increase for all magnetic quantities. Again relative approximation errors of the MC are generally larger than those of the MP, MV and MGT for all latitudes. With increased latitude, some curves (i.e., \(\delta {\varPsi }\), \(\delta B_{y}\), \(\delta B_{z}\), \(\delta B_{xx}\), \(\delta B_{yy}\), \(\delta B_{yz}\), \(\delta B_{zz}\), \(\delta B_{xxy}\), \(\delta B_{xxz}\), \(\delta B_{yyx}\), \(\delta B_{yyy}\), \(\delta B_{yyz}\), \(\delta B_{zzy}\) and \(\delta B_{zzz}\)) decline gradually at first, then rise quickly at the pole, whereas other curves (i.e., \(\delta B_{x}\), \(\delta B_{xz}\), \(\delta B_{xxx}\) and \(\delta B_{zzx}\)) increase slowly from \(0^{\circ }\) to \(85^{\circ }\). The minimum and maximum values in \(\hbox {Log}_{10}\) scale of the relative approximation errors for the MP, MV, MGT and MC are approximately \([-\,9, -6]\) for the MP, \([-\,8, -\,4]\) for the MV, \([-\,6, -\,2]\) for the MGT and \([-\,6, -\,1]\) for the MC, respectively, where larger errors are present at the pole.

The minimum and maximum values (i.e., magnitude) of the approximation errors of the MGT and MC Laplace parameters illustrated in Fig. 5b are \([-\,15, -\,11]\) for the MC Laplace parameters (\(\delta \varDelta L_2\), \(\delta \varDelta L_3\) and \(\delta \varDelta L_4\)) and \([-\,15, -\,12]\) for the MGT Laplace parameter (\(\delta \varDelta L_1\)) compared to the actual magnitude of the quantities. For the satellite case, the ranges of the Laplace parameters validate the reliability of the calculations as near zone problems are minimized. It is likely that a canceling effect is present at lower magnitude for the satellite height case. In addition, the ranges are consistent with the findings in Fig. 4b for the rapid drop zone, being less affected by the near zone problems.

As for the polar singularity problem of the MC, the curves of the MC in Fig. 5a have no numerical overflow at the latitude \(90^{\circ }\). In other words, the polar singularity does not occur in Fig. 5a for the MC at the North Pole. It should be noted that the errors increase significantly towards the pole, i.e., some numerical problems may occur probably in part related to the “deformation” of the tesseroids (due to meridian convergence) towards the pole. The same for this character can be found for the MP, MV and MGT in Fig. 5a. The reason can be inferred from the use of Cartesian integral kernels in the tesseroid formulas for the MC with similar behavior experienced for the GC using the Cartesian kernels as shown in Deng and Shen (2018b). Thus, it can be expected that the same properties are inherited by the magnetic quantities.

4 Conclusions

In recent years, the concept of gravitational curvatures has been proposed in gravity field modeling, and the theoretical and practical aspects of the GC have been investigated. Analogously, the general concept of magnetic curvatures was outlined in the “Conclusions and outlook” part of Deng and Shen (2018a), while the formulas of the MC using Cartesian integral kernels were presented in this contribution. In order to derive expressions for the MC, the general expressions for the nth-order derivatives of the GP not only for the uniform tesseroid using the Cartesian kernels but also for the uniform spherical shell have been derived and provided in Appendices 1 and 2. Using Poisson’s relation together with the components of the fourth-order derivatives of the GP of the uniform tesseroid and spherical shell, expressions for the MC of a uniformly magnetized tesseroid and spherical shell are derived in Appendices 3 and 4, respectively. The correctness of the newly derived expressions has indirectly been confirmed through the numerical satisfaction at the double precision level of the Laplace equations for the MC, where the order of magnitudes of the MC are about − 23 from Table 2, and the magnitudes of the absolute approximation errors of the MC Laplacians are about − 35 from Fig. 5b.

The near zone problem has been studied for the MP, MV, MGT and MC with the 3D GLQ approach. Numerical results indicate that when the computation point P approaches the surface of the spherical shell, i.e., gets increasingly closer to the source masses, largely increased relative approximation errors are present not only for the MP, MV and MGT, but also for the MC. As demonstrated in Shen and Deng (2016) and Deng and Shen (2018a, b), the near zone problem also exists for the GP, GV, GGT and GC using the TSE and GLQ approaches. Thus, the weak singularity of the kernel functions is a main limitation for the given discretization problems both in gravity and magnetic fields modeling. For this conclusion refer to the fact that formulas of magnetic quantities are based on those of the gravitational functional, thus it is expected that they inherit similar properties. Moreover, when the computation point P is in the near zone, the MC are impacted to a larger extent than the MP, MV and MGT, i.e.,the MC have larger relative approximation errors than the MP, MV and MGT.

In addition, we performed numerical experiments at satellite height \(h_P=260\, \hbox {km}\) to study the polar singularity for computation points with latitude from \(0^{\circ }\) to \(90^{\circ }\). For the satellite height application at the polar region (\(\varphi \ge 85^{\circ }\)), all relative approximation errors increase with larger magnitudes for the MC than those of the MP, MV and MGT. Numerical experiments revealed that polar singularity problems as manifested by a lack of numerical overflow do not occur for the MC at satellite height, which is the same behavior as for the GC using the Cartesian kernels as shown in Deng and Shen (2018b). Moreover, the relative accuracy range of the MC in \(\hbox {Log}_{10}\) form could be achieved as about [− 6, − 1] with tesseroid size \(1^{\circ }\times 1^{\circ }\) at satellite height. The larger discretizations lead to larger errors for the MV in Baykiev et al. (2016) and the GP in Shen and Deng (2016).

As the concept of magnetic curvatures was proposed in magnetic field, other MC formulas for different mass elements (e.g., rectangular prism, vertical cylindrical prism, lines, layers, polyhedron, tetrahedron and spherical cap) both in the spatial domain and in spectral domain (e.g., based on the spherical harmonic (SH) expansion), are left for future studies. Analogous to the first laboratory measurement of the GC (Rosi et al. 2015), direct measurements of the MC could be expected in near future (Qi et al. 2019; Bandyopadhyay et al. 2019). Applications of the MC components are expected to provide supplementary information for geomagnetic field modeling. For instance, the global magnetic models (e.g., SH models to degree/order 720 or higher) could be used for a follow-up study of how the MC look like and what geophysical information they convey.