Comparison of Different Methods for a Moho Modeling Under Oceans and Marginal Seas: A Case Study for the Indian Ocean

Abstract

Since marine seismic studies are relatively sparse and unevenly distributed, detailed tomographic images of the Moho geometry under large parts of the world’s oceans and marginal seas are not yet available. Marine gravity data is, therefore, often used to detect the Moho depth in these regions. Alternatively, Airy’s isostatic theory can be applied for this purpose. In this study, we compare different isostatic and gravimetric methods for a Moho recovery under the oceanic crust and continental margins, particularly focusing on a numerical performance of Airy, Vening Meinesz–Moritz (VMM), direct gravity inversion, and generalized (for the Earth’s spherical approximation) Parker–Oldenburg methods. Numerical experiments are conducted to estimate the Moho depth beneath the Indian Ocean. Results reveal that, among these investigated methods, the VMM model is probably the most suitable for a gravimetric Moho recovery beneath the oceanic crust and continental margins, when taking into consideration the lithospheric mantle density information. This method could to some extent model realistically a Moho geometry beneath mid-oceanic spreading ridges, oceanic subductions, most of oceanic volcanic formations, and marine sediment deposits. Nonetheless, this model still cannot fully reproduce a gradual Moho deepening caused by a conductive cooling and a subsequent isostatic rebalance of the oceanic lithosphere, which can functionally be described by a Moho deepening with the increasing ocean-floor age. Results also indicate that the Airy method typically overestimates the Moho depth under oceanic volcanic formations, while the direct gravity inversion and generalized Parker–Oldenburg methods could not reproduce more detailed features in the Moho geometry. Since Pratt’s theory better describes a large-scale isostatic mechanism of the oceanic lithosphere by means of compensation density variations, but does not account for additional changes in compensation depth (i.e., Moho depth) that are caused by these density changes, we tested a possibility of combining Pratt and Airy’s isostatic theories in order to estimate the Moho depth under the oceanic crust. Even this combined model cannot fully reproduce a gradual Moho deepening with the increasing ocean-floor age.

Appendices

Appendix A: Direct Gravity Inversion Method—Functional Relation Between the Moho Depth and Gravity Data

The gravitational contribution generated by a Moho geometry can be expressed as follows (Tenzer and Chen 2014a; 2014b):

$$ \delta g^{B} \left( {r,\Omega } \right) = - {\text{G}}\;\iint\limits_{\Phi } {\Delta \rho^{{\text{c/m}}} \left( {\Omega^{^{\prime}} } \right)\;\int_{{\;r^{^{\prime}} = {\text{R}} - {\text{D}}_{{0}} }}^{{\;{\text{R - }}D_{M} \left( {\Omega^{^{\prime}} } \right)}} {\frac{{\partial \,l^{ - 1} \left( {r,\psi ,\,r^{^{\prime}} } \right)}}{\partial \;r}\;} r^{^{\prime}2} \;dr^{^{\prime}} d}\Omega^{^{\prime}} , $$

(A.1)

where the Moho depths $D_{M}$ are measured relative to the geoid surface and ${\text{D}}_{{0}}$ is the mean Moho depth.

The spectral representation of the radial derivative of the reciprocal spatial distance $\partial \,l^{ - 1} /\partial \,r$ for the external convergence domain $r \ge r^{^{\prime}}$ reads (e.g., Heiskanen and Moritz 1967)

$$ \frac{{\partial \,l^{ - 1} \left( {r,\psi ,\,r^{^{\prime}} } \right)}}{\partial \,r} = - \frac{1}{{r^{^{\prime}2} }}\sum\limits_{n = 0}^{\infty } {\left( {\frac{{r^{^{\prime}} }}{r}} \right)^{n + 2} \left( {n + 1} \right)\;{\text{P}}_{{\text{n}}} \left( t \right)} . $$

(A.2)

Substitution from Eq. (A.2) to Eq. (A.1) yields

$$ \delta g^{B} \left( {r,\Omega } \right) = {\text{G}}\;\iint\limits_{\Phi } {\Delta \rho^{{\text{c/m}}} \left( {\Omega^{^{\prime}} } \right)\int_{{\;r^{^{\prime}} = {\text{R}} - {\text{D}}_{{0}} }}^{{\;{\text{R - }}D_{M} \left( {\Omega^{^{\prime}} } \right)}} {\sum\limits_{n = 0}^{\infty } {\left( {\frac{{r^{^{\prime}} }}{r}} \right)^{n + 2} \left( {n + 1} \right)\;{\text{P}}_{{\text{n}}} \left( t \right)\;dr^{^{\prime}} } } \,d}\Omega^{^{\prime}}. $$

(A.3)

Defining the radial integral on the right-hand side of Eq. (A.3) as

$$ \begin{aligned} K\left( {r,t,D_{M} } \right) & = \int_{{\;r^{\prime} = {\text{R}} - {\text{D}}_{{0}} }}^{{\;{\text{R}} - D_{M} \left( {\Omega^{^{\prime}} } \right)}} {\sum\limits_{n = 0}^{\infty } {\left( {\frac{{r^{^{\prime}} }}{r}} \right)^{n + 2} \left( {n + 1} \right)\;{\text{P}}_{{\text{n}}} \left( t \right)\;dr^{^{\prime}} } } \\ \, & = - r\sum\limits_{n = 0}^{\infty } {{\text{P}}_{{\text{n}}} \left( t \right)} \left( {\frac{{{\text{R}} - {\text{D}}_{{0}} }}{r}} \right)^{n + 3} \frac{n + 1}{{n + 3}}\left[ {1 - \left( {\frac{{{\text{R}} - D_{M} }}{{{\text{R}} - {\text{D}}_{{0}} }}} \right)^{n + 3} } \right]\,, \\ \end{aligned} $$

(A.4)

and combing Eqs. (A.4) and (A.3), we arrive at

$$ \delta g^{B} \left( {r,\Omega } \right) = {\text{G}}\;\iint\limits_{\Phi } {\Delta \rho^{{\text{c/m}}} \left( {\Omega^{^{\prime}} } \right)\;K\left( {r,t,D_{M} } \right)\;d}\Omega^{^{\prime}} . $$

(A.5)

The expression in Eq. (A.5) is a nonlinear Fredholm integral equation of the first kind (e.g., Sjöberg 2009). Its linearization is done by applying a Taylor series with respect to a mean Moho depth ${\text{D}}_{{0}}$, while disregarding the higher than first-order terms. We then write

$$ K\left( {r,t,D_{M} } \right) = K\left( {r,t,{\text{D}}_{{0}} } \right) + \left. {\frac{\partial }{{\partial D_{M} }}{\text{K}}(r,t,D_{M} )} \right|_{{D_{M} = {\text{D}}_{{0}} }} (D_{M} - {\text{D}}_{{0}} ). $$

(A.6)

It can readily be shown that $K\left( {r,t,{\text{D}}_{{0}} } \right) = 0$. Hence, after inserting from Eq. (A.6) back to Eq. (A.5), we get

$$ \delta g^{B} \left( {r,\Omega } \right) \cong {\text{G}}\iint\limits_{\Phi } {\Delta \rho^{{\text{c/m}}} \left( {\Omega^{^{\prime}} } \right)\;\left. {\frac{\partial }{{\partial D_{M} }}{\text{K}}(r,t,D_{M} )} \right|_{{D_{M} = {\text{D}}_{{0}} }} \,dD^{^{\prime}} \;d\Omega^{^{\prime}} }, $$

(A.7)

where $dD = D_{M} - {\text{D}}_{{0}}$.

We further define

$$ \begin{aligned} T(r,t,{\text{D}}_{{0}} ) & = \left. {\frac{\partial }{{\partial D_{M} }}{\text{K}}(r,t,D_{M} )} \right|_{{D_{M} = {\text{D}}_{{0}} }} \\ \, & { = } - r\sum\limits_{n = 0}^{\infty } {\left( {\frac{{{\text{R}} - {\text{D}}_{{0}} }}{r}} \right)}^{n + 3} \frac{n + 1}{{n + 3}}{\text{P}}_{{\text{n}}} (t)\left. {\left[ { - (n + 3)\left( {\frac{{{\text{R}} - D_{M} }}{{{\text{R}} - {\text{D}}_{{0}} }}} \right)^{n + 2} \frac{ - 1}{{{\text{R}} - {\text{D}}_{{0}} }}} \right]} \right|_{{D_{M} = {\text{D}}_{{0}} }} \\ \, & = - \sum\limits_{n = 0}^{\infty } {\left( {\frac{{{\text{R}} - D_{M} }}{r}} \right)}^{n + 2} (n + 1){\text{P}}_{{\text{n}}} (t)|_{{D_{M} = {\text{D}}_{{0}} }} \\ \, & = - \sum\limits_{n = 0}^{\infty } {\left( {\frac{{{\text{R}} - {\text{D}}_{{0}} }}{r}} \right)}^{n + 2} (n + 1){\text{P}}_{{\text{n}}} (t). \\ \end{aligned} $$

(A.8)

Denoting $\tau = \left( {{\text{R}} - {\text{D}}_{{0}} } \right)/r$, the spectral representation of the integral kernel ${\text{T}}(r,t,{\text{D}}_{{0}} )$ is rewritten as

$$ {\text{T}}\left( {t,\tau } \right) = - \sum\limits_{n = 0}^{\infty } {\tau^{n + 2} \,\left( {n + 1} \right)\,{\text{P}}_{{\text{n}}} \left( t \right)} . $$

(A.9)

To find the spatial form of ${\text{T}}\left( {t,\tau } \right)$, we first separate its spectral representation into zero- and higher-order terms (with respect to ${\text{P}}_{{\text{n}}}$) as follows:

$$ - {\text{T}}\left( {t,\tau } \right) = \tau^{2} \;{\text{P}}_{{0}} \left( t \right) + \sum\limits_{n = 1}^{\infty } {\tau^{n + 2} \,\left( {n + 1} \right)\,{\text{P}}_{{\text{n}}} \left( t \right)} , $$

(A.10)

and further rearrange Eq. (A.10) into the following form

$$ - {\text{T}}\left( {t,\tau } \right) = \tau^{2} \;\left[ {\,1 + \sum\limits_{n = 1}^{\infty } {n\tau^{n} \,{\text{P}}_{{\text{n}}} \left( t \right)} + \sum\limits_{n = 1}^{\infty } {\tau^{n} \,{\text{P}}_{{\text{n}}} \left( t \right)} } \right], $$

(A.11)

with

$$ \sum\limits_{n = 1}^{\infty } {\tau^{n} \,{\text{P}}_{{\text{n}}} \left( t \right)} = \frac{1}{\zeta } - 1, $$

(A.12)

and

$$ \sum\limits_{n = 1}^{\infty } {n\tau^{n} \,{\text{P}}_{{\text{n}}} \left( t \right)} = - \frac{{\tau \left( {\tau - t} \right)}}{{\zeta^{3} }}, $$

(A.13)

where

$$ \zeta = \sqrt {1 + \tau^{2} - 2\tau \,t} . $$

(A.14)

Combining Eqs. (A.9)–(A.14), the closed form of ${\text{T}}\left( {t,\tau } \right)$ is found to be

$$ \begin{gathered} {\text{T}}\left( {t,\tau } \right) = - \,\tau^{2} \;\left[ {\,1 - \frac{{\tau \left( {\tau - t} \right)}}{{\zeta^{3} }} + \frac{1}{\zeta } - 1} \right] \hfill \\ \, = - \tau^{2} \;\left( {1 - \tau \,t} \right)\;\left( {1 + \tau^{2} - 2\tau \,t} \right)^{{ - {\kern 1pt} \frac{3}{2}}} . \hfill \\ \end{gathered} $$

(A.15)

The functional relation between the (complete) Bouguer gravity disturbances $\delta g^{B}$ and the unknown Moho depth corrections $dD^{^{\prime}} = D_{M} - {\text{D}}_{{0}}$ in Eq. (A.7) is defined by means of the integral convolution of ${\text{T}}\left( {t,\tau } \right)$ and $\Delta \rho^{{\text{c/m}}}$. The discretization of the integral in Eq. (A.7) is done by means of dividing the entire surface integration domain into a finite number of surface elements. The most convenient way is to apply discretization on a regular grid of geographical coordinates. Introducing the surface elements $\Delta \Omega_{{\text{j}}}^{^{\prime}} = \;\;\cos \varphi_{{\text{j}}}^{^{\prime}} \,\Delta \varphi_{{\text{j}}}^{^{\prime}} \;\Delta \lambda_{{\text{j}}}^{^{\prime}}$ (where $\Delta \varphi^{\prime}$ and $\Delta \lambda^{\prime}$ are discretization steps in latitude and longitude, respectively), the integral equation in Eq. (A.7) becomes

$$ \delta g_{i}^{B} \cong {\text{G}}\,\sum\limits_{j = 1}^{J} {\Delta {\uprho }_{j}^{{\text{c/m}}} \;\int\limits_{{\Delta \varphi^{\prime}_{{\text{j}}} }} {\int\limits_{{\Delta \lambda^{\prime}_{{\text{j}}} }} {{\text{T}}\left( {t,\tau } \right)} } \;\cos \varphi^{^{\prime}} \;d\varphi^{^{\prime}} d\lambda^{^{\prime}} dD_{j}^{^{\prime}} } \;\,\,\,\,\left( {i = 1,\,2,\,...,I} \right), $$

(A.16)

where $I$ is the total number of input gravity data and $J$ is the total number of discretization elements. In our discrete model, we assume that $I = J$. When $I > J$, the Moho depths can be determined by applying a least-squares analysis. The Moho density contrast values $\left\{ {\,\Delta {\uprho }_{j}^{{\text{c/m}}} :\;j = 1,\,2,\,...,\,J} \right\}\;$ are considered within each of the discretized surface elements $\left\{ {\,\Delta \Omega^{\prime}_{{\text{j}}} :\;j = 1,\,2,\,...,\,J} \right\}^{{}}$.

Since the surface integral in Eq. (A.16) does not have a closed analytical form in terms of geographical coordinates, it is solved (for each surface element) numerically. Applying the simplest numerical scheme, Eq. (A.16) becomes

$$ \delta g_{i}^{m} \cong {\text{G}}\,\sum\limits_{j = 1}^{J} {\Delta {\uprho }_{j}^{{\text{c/m}}} \;{\text{T}}_{{\text{i,j}}} \;\Delta \Omega_{j}^{^{\prime}} \;dD_{j}^{^{\prime}} } \;\,\,\,\,\left( {i = 1,\,2,\,...,I} \right), $$

(A.17)

where a surface integral for each individual surface element is evaluated from

$$ \int\limits_{{\Delta \varphi_{{\text{j}}} }} {\int\limits_{{\Delta \lambda_{{\text{j}}} }} {{\text{T}}\left( {t,\tau } \right)} } \;\cos \varphi^{^{\prime}} \;d\varphi^{^{\prime}} \;d\lambda^{^{\prime}} \approx {\text{T}}_{{\text{i,j}}} \;\Delta \Omega_{j}^{^{\prime}} . $$

(A.18)

For $i \ne j$, the kernel ${\text{T}}_{{\text{i,j}}}$ in Eq. (A.18) is computed as

$$ {\text{T}}_{{\text{i,j}}} = \, - \left( {\frac{{{\text{R}} - D_{0} }}{{{\text{R}} + H_{i} }}} \right)^{2} \;\left( {1 - \frac{{{\text{R}} - D_{0} }}{{{\text{R}} + H_{i} }}\,t_{i,j} } \right)\;\left[ {1 + \left( {\frac{{{\text{R}} - D_{0} }}{{{\text{R}} + H_{i} }}} \right)^{2} - 2t_{i,j} \left( {\frac{{{\text{R}} - D_{0} }}{{{\text{R}} + H_{i} }}} \right)} \right]^{{ - \frac{3}{2}}} . $$

(A.19)

For $i = j$, $t_{i,j} = 1$, after some algebra, the kernel ${\text{T}}_{{\text{i,i}}}$ simplifies to

$$ {\text{T}}_{{\text{i,i}}} = \, - \left( {\frac{{{\text{R}} - D_{0} }}{{H_{i} + D_{0} }}} \right)^{2} \;. $$

(A.20)

The integral kernel ${\text{T}}$ has a singularity for $\tau \to 1 \wedge t \to 1$. However, this singularity does not occur in the numerical solution because a minimum crustal thickness is everywhere at least a few kilometers below the geoid surface.

Appendix B: Generalized Parker–Oldenburg Method—Functional Relation Between the Moho Depth and Gravity Data

The disturbing potential generated by a Moho geometry can be expressed as follows (Chen and Tenzer, 2017b):

$$ T^{B} \left( {r,\Omega } \right) = {\text{G}}\;\iint\limits_{\Phi } {\Delta \rho^{{\text{c/m}}} \left( {\Omega^{^{\prime}} } \right)\;\int_{{\;r^{\prime} = {\text{R}} - {\text{D}}_{0} }}^{{\;{\text{R}} - D_{M} \left( {\Omega^{\prime}} \right)}} {l^{ - 1} \left( {r,\psi ,\,r^{^{\prime}} } \right)} \;r^{^{\prime}2} \;dr^{^{\prime}} d}\Omega^{^{\prime}}. $$

(B.1)

The spectral representation of a reciprocal spatial distance $l^{ - 1}$ for the external convergence domain $r \ge r^{^{\prime}}$ reads (e.g., Heiskanen and Moritz 1967),

$$ l^{ - 1} \left( {r,\psi ,\,r^{^{\prime}} } \right) = \frac{1}{{r^{^{\prime}} }}\sum\limits_{n = 0}^{\infty } {\left( {\frac{{r^{^{\prime}} }}{r}} \right)^{n + 1} \;{\text{P}}_{{\text{n}}} \left( t \right)} . $$

(B.2)

Using the addition theorem, the Legendre polynomial ${\text{P}}_{{\text{n}}}$ can be expressed as follows:

$$ {\text{P}}_{{\text{n}}} (t) = \frac{1}{2n + 1}\sum\limits_{m = - n}^{n} {{\text{Y}}_{{{\text{nm}}}} (\Omega ){\text{Y}}_{{{\text{nm}}}} (\Omega^{^{\prime}} )} . $$

(B.3)

Combining Eqs. (B.1)–(B.3), we arrive at

$$ \begin{aligned} T^{B} \left( {r,\Omega } \right) & = \;\frac{{\text{G}}}{r}\iint\limits_{\Phi } {\Delta \rho^{{\text{c/m}}} \left( {\Omega^{^{\prime}} } \right)\;\int_{{\;r^{^{\prime}} {\text{R}} - {\text{D}}_{0} }}^{{\;{\text{R}} - D_{M} \left( {\Omega^{^{\prime}} } \right)}} {\sum\limits_{n = 0}^{\infty } {\frac{1}{2n + 1}\left( {\frac{{r^{^{\prime}} }}{r}} \right)^{n} \;\sum\limits_{m = - n}^{n} {{\text{Y}}_{{{\text{nm}}}} (\Omega ){\text{Y}}_{{{\text{nm}}}} (\Omega^{^{\prime}} )} } } \;dr^{^{\prime}} d}\Omega^{^{\prime}} \\ & = \;\frac{{{\text{GM}}}}{{\text{R}}}\sum\limits_{n = 0}^{\infty } {\frac{1}{2n + 1}\left( {\frac{{\text{R}}}{r}} \right)^{n + 1} \;\sum\limits_{m = - n}^{n} {{\text{T}}_{{{\text{nm}}}}^{B,0} \;{\text{Y}}_{{{\text{nm}}}} (\Omega )} } . \\ \end{aligned} $$

(B.4)

The disturbing potential coefficients ${\text{T}}_{{{\text{nm}}}}^{B}$ (corrected for gravitational contributions of topography and crustal density heterogeneities) in Eq. (B.4) are given by

$$ \begin{aligned} {\text {T}}_{{\text nm}}^{B} & = \;\frac{1}{\text M}\frac{1}{{2n + 1}}\iint\limits_{\Phi } {\Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;\int_{{\;r^{'} = R - D_{0} }}^{{\; {\text R} - {\text D}_{\text M} \left( {\Omega ^{'} } \right)}} {\left( {\frac{{r^{'} }}{\text R}} \right)^{n} } Y_{{\text nm}} (\Omega ^{'} )\,dr^{'} d}\Omega ^{'} \\ & = \;\frac{1}{{4\pi \bar{\rho }^{{\text Earth}} {\text R}^{3} }}\frac{3}{{2n + 1}}\iint\limits_{\Phi } {\Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;\int_{{\;r^{'} = {\text R} - {\text D}_{0} }}^{{\; {\text R} - {\text D}_{\text M} \left( {\Omega ^{'} } \right)}} {\left( {\frac{{r^{'} }}{\text R}} \right)^{n} r^{{'2}} \;Y_{{\text nm}} (\Omega ^{'} )} \;dr^{'} d}\Omega ^{'} , \\ \end{aligned} $$

(B.5)

where ${\text{M}} = \;\frac{4}{3}{{\uppi \overline{\uprho }}}^{{{\text{Earth}}}} {\text{R}}^{3}$ is the Earth’s total mass.

The radial integral in Eq. (B.5) is formally divided into two parts (e.g., Sjöberg 2009)

$$ \begin{aligned} {\text{T}}_{{nm}}^{B} + \Delta T_{{nm}}^{{B,0}} & = \;\frac{3}{{4\pi \bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\sum\limits_{{k = 1}}^{{n + 3}} {C_{{n + 3}}^{k} } \iint\limits_{\Phi } {\tau ^{k} \Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;Y_{{nm}} (\Omega ^{'} )\;d}\Omega ^{'} \\ & = \;\frac{3}{{\bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\sum\limits_{{k = 1}}^{{n + 3}} {C_{{n + 3}}^{k} (\Delta \rho \tau ^{k} )_{{nm}} } \\ & = \;\frac{3}{{\bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\left[ {\left( {n + 3} \right)(\Delta \rho \tau )_{{nm}} + \sum\limits_{{k = 2}}^{{n + 3}} {C_{{n + 3}}^{k} (\Delta \rho \tau ^{k} )_{{nm}} } } \right]. \\ \end{aligned} $$

(B.6)

The radial integral kernels of two constituents on the right-hand side of Eq. (B.6) read

$$ {\text{K}}_{1} = - \frac{{{\text{R}}^{3} }}{n + 3}\left[ {\sum\limits_{k = 1}^{n + 3} {\tau_{0}^{k} {\text{C}}_{n + 3}^{k} } } \right],\,\,\,\,{\text{K}}_{2} = - \frac{{{\text{R}}^{3} }}{n + 3}\left[ {\sum\limits_{k = 1}^{n + 3} {\tau^{k} {\text{C}}_{n + 3}^{k} } } \right], $$

(B.7)

where $\tau_{0} = - {\text{D}}_{{0}} /{\text{R}}$ and $\tau = - D_{M} \left( {\Omega^{^{\prime}} } \right)/{\text{R}}.$

Inserting from Eqs. (B.6) and (B.7) back to Eq. (B.5), we arrive at

$$ \begin{aligned} {\text{T}}_{{nm}}^{B} & = \;\frac{1}{{4\pi \bar{\rho }^{{{\text{Earth}}}} R^{3} }}\frac{3}{{2n + 1}}\iint\limits_{\Phi } {\left( {K_{1} + K_{2} } \right)\Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;Y_{{nm}} (\Omega ^{'} )\;d}\Omega ^{'} , \\ & = \;\frac{3}{{4\pi \bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\sum\limits_{{k = 1}}^{{n + 3}} {C_{{n + 3}}^{k} } \left[ {\iint\limits_{\Phi } {\tau ^{k} \Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;Y_{{nm}} (\Omega ^{'} )\;d}\Omega ^{'} } \right. \\ & \,\,\,\,\,\, - \left. {\iint\limits_{\Phi } {\tau _{0}^{k} \Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;Y_{{nm}} (\Omega ^{'} )\;d}\Omega ^{'} } \right]. \\ \end{aligned} $$

(B.8)

Now, we define

$$ \Delta T_{{nm}}^{{B,0}} = \;\frac{3}{{4\pi \bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\sum\limits_{{k = 1}}^{{n + 3}} {C_{{n + 3}}^{k} } \iint\limits_{\Phi } {\tau _{0}^{k} \Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;Y_{{nm}} (\Omega ^{'} )\;d}\Omega ^{'} . $$

(B.9)

Moreover, we introduce the following relation:

$$ \begin{aligned} T_{{nm}}^{B} + \Delta T_{{nm}}^{{B,0}} & = \;\frac{3}{{4\pi \bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\sum\limits_{{k = 1}}^{{n + 3}} {C_{{n + 3}}^{k} } \iint\limits_{\Phi } {\tau ^{k} \Delta \rho ^{{c/m}} \left( {\Omega ^{'} } \right)\;Y_{{nm}} (\Omega ^{'} )\;d}\Omega ^{'} \\ & = \;\frac{3}{{\bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\sum\limits_{{k = 1}}^{{n + 3}} {C_{{n + 3}}^{k} (\Delta \rho \tau ^{k} )_{{nm}} } \\ & = \;\frac{3}{{\bar{\rho }^{{{\text{Earth}}}} \left( {2n + 1} \right)\left( {n + 3} \right)}}\left[ {\left( {n + 3} \right)(\Delta \rho \tau )_{{nm}} + \sum\limits_{{k = 2}}^{{n + 3}} {C_{{n + 3}}^{k} (\Delta \rho \tau ^{k} )_{{nm}} } } \right]. \\ \end{aligned} $$

(B.10)

The expression in Eq. (B.10) is further rewritten as

$$ (\Delta \rho \tau )_{nm} = \frac{{\left( {{\text{T}}_{{{\text{nm}}}}^{B} + {\Delta T}_{{{\text{nm}}}}^{B,0} } \right){\overline{\rho }}^{{{\text{Earth}}}} \left( {2n + 1} \right)}}{3} - \frac{1}{n + 3}\sum\limits_{k = 2}^{n + 3} {{\text{C}}_{n + 3}^{k} (\Delta \rho \tau^{k} )_{nm} } . $$

(B.11)

Since both sides of Eq. (B.11) comprise the term $(\Delta \rho \tau )_{nm}$, the solution is carried out iteratively.