# A special case of the Poisson PDE formulated for Earth's surface and its capability to approximate the terrain mass density employing land-based gravity data, a case study in the south of Iran

Yahya Allahtavakoli, Abdolreza Safari, Petr Vaníček

► To cite this version:

Yahya Allahtavakoli, Abdolreza Safari, Petr Vaníček. A special case of the Poisson PDE formulated for Earth's surface and its capability to approximate the terrain mass density employing land-based gravity data, a case study in the south of Iran. Geophysical Journal International, 2016, 207 (3), pp.1529-1553. 10.1093/gji/ggw340. hal-04925096

HAL Id: hal-04925096

https://hal.science/hal-04925096v1

Submitted on 1 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A special case of the Poisson PDE formulated for Earth's surface and its capability to approximate the terrain mass density employing land-based gravity data, a case study in the south of Iran

Yahya AllahTavakoli,<fnref n="1" /> Abdolreza Safari<fnref n="1" /> and Petr Vaníček<fnref n="2" />

Accepted 2016 September 12. Received 2016 August 2; in original form 2015 July 21

## SUMMARY

This paper resurrects a version of Poisson's Partial Differential Equation (PDE) associated with the gravitational field at the Earth's surface and illustrates how the PDE possesses a capability to extract the mass density of Earth's topography from land-based gravity data. Herein, first we propound a theorem which mathematically introduces this version of Poisson's PDE adapted for the Earth's surface and then we use this PDE to develop a method of approximating the terrain mass density. Also, we carry out a real case study showing how the proposed approach is able to be applied to a set of land-based gravity data. In the case study, the method is summarized by an algorithm and applied to a set of gravity stations located along a part of the north coast of the Persian Gulf in the south of Iran. The results were numerically validated via rock-samplings as well as a geological map. Also, the method was compared with two conventional methods of mass density reduction. The numerical experiments indicate that the Poisson PDE at the Earth's surface has the capability to extract the mass density from land-based gravity data and is able to provide an alternative and somewhat more precise method of estimating the terrain mass density.

**Key words:** Inverse theory; Gravity anomalies and Earth structure; Geopotential theory.

## 1 INTRODUCTION

Variations in the Earth's gravitational field are a vital clue to a better understanding of the Earth's mass density distribution. Nowadays, various land-based gravimetric techniques straightforwardly enable Earth's gravitational field to be observable at the Earth's surface, and based on Newtonian gravitational potential theory, the land-based gravitational observations can lead geoscientists to beneficial knowledge about the density distribution of Earth's masses. Newtonian gravitational potential theory indicates that Earth's gravitational field is fundamentally entangled with the concept of the Earth's mass density, and it is inherently based on two famous equations: (1) an Integral Equation (IE) called Newton's IE, and (2) a Partial Differential Equation (PDE) called the Poisson PDE. The Poisson PDE as the inverse form of Newton's IE is a well-known mathematical model in Newtonian gravitational potential theory. Newtonian potential theory straightforwardly shows that the Earth's gravitational field satisfies Poisson's PDE within Earth's interior space if the mass-density functional satisfies a continuity condition called the Holder condition (Kellogg 1929; p. 156). Nevertheless, because of the discontinuity in the mass-density function

when crossing from the Earth's interior space to its exterior space, the theory of Newtonian potential cannot specifically mention the situation of Poisson's PDE at the Earth's boundary. However, there exists an unfamiliar case of this PDE which may be probable at the Earth's surface as follows (e.g. see Vaníček & Krakiwsky 1986, p. 462):

$$ \nabla^2 V = -k\pi G\rho \quad (1) $$

where the functional $V: \mathbb{R}^3 \to \mathbb{R}$ is the Earth's gravitational potential, $\nabla^2 = \nabla^T \nabla$ is a Laplace operator, $G$ is the universal gravitation constant, $\rho$ is the Earth's mass-density function, and the topographic coefficient $k$ may differ from zero to four depending on the topographic form of the Earth's surface (i.e. $k \in (0, 4)$), see Fig. 1. Eq. (1) can be easily reshaped into the following ordinary form of the Poisson PDE associated with the Earth's interior space (e.g. see Moritz 1980; Hofmann-Wellenhof & Moritz 2006) which can be considered as a rough justification for satisfying the Poisson PDE at the Earth's surface.

$$ \nabla^2 V = -4\pi G \bar{\rho} \quad (2) $$

<footnote n="1">¹College of Engineering, School of Surveying and Geospatial Engineering, University of Tehran, Tehran, Iran. E-mail: asafari@ut.ac.ir</footnote>

<footnote n="2">²Department of Geodesy and Geomatics Engineering, University of New Brunswick, P.O. Box 4400, Fredericton, New Brunswick E3B 5A3, Canada</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic representation of the topographic coefficient as well as the differential neighbourhood schematically illustrated by a circle around a point on the Earth' surface for various conditions (see Vaniček & Krakiwsky 1986, pp. 512, 462).</figcaption></figure>

where $\bar{\rho} = \bar{k}\rho$ and $\bar{k} = k/4 \in (0, 1)$. The above equation indicates that a point located on the Earth's surface may have a density ($\bar{\rho}$) equal to a fraction ($\bar{k}$) of the certain density ($\rho$) corresponding to a point located just beneath the Earth's surface entirely within the Earth (Vaniček & Krakiwsky 1986, p. 462). Any differential neighbourhood of the point on the Earth's surface is partially empty and partially dense; hence, the density $\bar{\rho}$ is defined as the ratio of the mass enclosed in the dense part of the differential neighbourhood to the total volume of the differential neighbourhood, see Fig. 1.

Although there exist certain basic justifications (such as eq. 2) for the Poisson PDE at the Earth's surface (i.e. eq. 1), it hitherto has been left mathematically unconfirmed. Also, it has not still found significant applications in any geo-scientific subjects. Herein, we mathematically introduce the above version of the Poisson PDE at the Earth's surface. For this purpose, some definitions are presented and a theorem guaranteeing the above Poisson PDE is proved. Then, we use the PDE to propose a method for determining the terrain mass-density distribution from land-based gravity data without using any prior knowledge about the mass density.

There exist certain methods of estimating the lateral mass-density distribution of the terrain, using only land-based gravity data, such as the profile method of Nettleton (1939), the F-H method (Parasnis 1952, 1979), the F-H method of Fukao et al. (1981), the Akaike's Bayesian Information Criterion (ABIC) method of Murata (1993) and the method of Nawa et al. (1997). First, Nettleton (1939) proposed an idea of determining the terrain density by means of minimizing the correlation between profiles of the complete Bouguer anomaly and topography. Then Parasnis (1952) used the Nettleton's idea for proposing the F-H method that extracted the mass density from the correlation between the free-air anomaly and topographic corrections. Then Fukao et al. (1981) developed the F-H method of Parasnis (1952) based on reducing the correlation of complete Bouguer anomaly with topography via stepwise functions. Murata (1993) also developed the F-H method by means of using a Bayesian approach to obtain optimal mass density. Then, Nawa et al. (1997) also developed the Bayesian method of Murata (1993) for proposing a method of inverting grav-

ity data for the lateral density variation of the terrain above sea level.

All the aforementioned methods, which are so adequate for estimating the unknown mass density required for the gravity reduction, have the benefit of approximating the mass density via the land-based gravity data without using any additional knowledge of the mass density. Hence, they are usually categorized as the methods of estimating the mass density reduction (there exist other names referring to these methods, such as methods of estimating Bouguer density, surficial density, and terrain density). Similarly, the proposed method of this paper, provided by the Poisson PDE, is also able to determine the lateral mass-density variation of the terrain only via land-based gravity data, without using any prior information about the mass density. Herein, the theories, as well as the satisfactory numerical experiments, have shown that the proposed approach can be considered as an appropriate alternative method of determining terrain mass density. In the current research, we see that the Poisson PDE adapted for the Earth's surface provides the proposed method with certain theoretical advantages over its conventional predecessors. The above-mentioned conventional methods employ only the vertical gradient of land-based gravity data and they do not use the horizontal derivatives for evaluating the mass density. Moreover, these methods have no capability of extracting any knowledge of the mass density from the deflection of the vertical and they neglect the effects of the mass density on the deflection of the vertical. However, the proposed method of the current research is also enabled by the Poisson PDE to use the horizontal derivatives as well as any knowledge of the deflection of the vertical for assessing the terrain mass density. The proposed method needs the deflection of the vertical to change the scalar land-based gravity data into the Earth's gravitational vector field which is used by the Poisson PDE for approximating the terrain mass density.

Also, the numerical experiments of the current research show that the proposed approach has the ability to integrate land-based gravity data with global gravitational field models (GGMs) for inverting the data into the Earth's mass density. Furthermore, the numerical experiments also indicate that the approach is able to provide the terrain mass density with a more accurate estimation in comparison with two conventional ones, the F-H method of Parasnis (1952) and the Bayesian method of Murata (1993).

In this paper, the plan of the presentation will be as follows: in Section 2, we introduce basic definitions and theorems required for the explanation of Poisson's PDE at the Earth's surface. This is followed in Section 3 by certain formulae originating a method of approximating the terrain mass density. Section 4 shows how the proposed approach is applied to a real case study. In this section, the case study numerically proves capabilities of the Poisson PDE as well, where we performed numerical validations for the proposed method. Finally, Section 5 ends the paper with the conclusions and the discussions.

## 2 THE POISSON PDE AT THE EARTH'S SURFACE

In this section, we bring some basic definitions and crucial theorems in order to introduce a version of Poisson's PDE to the Earth's surface. For this purpose, first it is necessary to present readers with the following definition of the Earth's gravitational gradient tensor.
**Definition 1** For a differentiable gravitational potential $V(\mathbf{r})$: $\mathbb{R}^3 \to \mathbb{R}$, a matrix such as $\Gamma(\mathbf{r})$: $\mathbb{R}^3 \to \mathbb{R}^{3 \times 3}$, called *gradiometric or gravitational gradient tensor*, is defined as follows:

$$ \Gamma(\mathbf{r}) = \begin{bmatrix} \frac{\partial^2 V(\mathbf{r})}{\partial x^2} & \frac{\partial^2 V(\mathbf{r})}{\partial y \partial x} & \frac{\partial^2 V(\mathbf{r})}{\partial z \partial x} \\ \frac{\partial^2 V(\mathbf{r})}{\partial x \partial y} & \frac{\partial^2 V(\mathbf{r})}{\partial y \partial z} & \frac{\partial^2 V(\mathbf{r})}{\partial z \partial y} \\ \frac{\partial^2 V(\mathbf{r})}{\partial x \partial z} & \frac{\partial^2 V(\mathbf{r})}{\partial y \partial z} & \frac{\partial^2 V(\mathbf{r})}{\partial z^2} \end{bmatrix} \quad (3) $$

where $\mathbf{r} = (x, y, z) \in \mathbb{R}^3$ denotes a position vector in a Cartesian coordinate system.

Nowadays the gradient tensor plays important roles in many geoscientific issues and it has found various applications to studying certain physical characteristics of the Earth. Herein, the gradient tensor is also an important quantity that simply enables the current research to implement its ideas.

Then, it is also necessary to recall a basic definition from Calculus, which is the definition of 'Divergence'. We know that the concept of the Poisson PDE is involved with the following definition corresponding to the divergence of Earth's gravitational field.

**Definition 2** The divergence of Earth's gravitational vector field $\mathbf{f}_g = \nabla V: \mathbb{R}^3 \to \mathbb{R}^3$ is defined as follows:

$$ \mathrm{div}(\mathbf{f}_g) = \nabla^2 V = \mathrm{trace}(\Gamma) \quad (4) $$

where $\nabla = [\partial/\partial x \ \partial/\partial y \ \partial/\partial z]^T$, $\nabla^2 = \partial^2/\partial x^2 + \partial^2/\partial y^2 + \partial^2/\partial z^2$ and trace are a gradient operator, a Laplace operator and an operator of the trace of matrix, respectively.

Eq. (4) implies the divergence to be an invariant value of the gradient tensor since trace is invariant with respect to a change of the basis of the coordinate system. Hence, the gradient tensor as well as its trace can always be used for investigating the divergence of Earth's gravitational field regardless of the Coordinate system.

According to the following theorem (Kellogg 1929, p. 156), the gradient tensor as well as the divergence will exist and satisfy the Poisson PDE within the Earth's interior space, if the Earth's mass density ($\rho$) is considered as a piecewise-continuous function.

**Theorem 1** Let $V$ be the potential of a distribution with piecewise continuous density $\rho$ in a regular region $\Theta$. Then at any interior point of $\Theta$, at which $\rho$ satisfies a Holder continuity condition, the derivatives of the second order of $V$ (i.e. the gradient tensor) exist and satisfy Poisson's PDE as follows:

$$ \nabla^2 V = -4\pi G\rho \quad (5) $$

where the regular region $\Theta$ is a differentiable surface encompassing all masses of the Earth, such as the Brillouin sphere, for example, see Moritz (1980) and Hofmann-Wellenhof & Moritz (2006). In Theorem 1, in fact, the mere continuity of the density does not suffice to ensure the above Poisson's PDE; therefore a continuity condition called the Holder condition is imposed on the density function $\rho$ (Kellogg 1929, p. 156). The density function $\rho$ is said to satisfy the Holder continuity condition at position $\mathbf{r}$ if there exist three positive constants $c$, $A$ and $\alpha$ such that

$$ |\rho(\mathbf{r}) - \rho(\mathbf{r}')| \le A \|\mathbf{r} - \mathbf{r}'\|^\alpha \quad (6) $$

for all points $\mathbf{r}'$ for which $\|\mathbf{r} - \mathbf{r}'\| \le c$ (Kellogg 1929, p. 156).

Accordingly, Theorem 1 leaves unmentioned the situation at discontinuities of $\rho$ such as boundary points of the Earth's surface. Hence, we cannot say that in general the gradient tensor exists at

points of the Earth's surface. Therefore, instead we here define the gradient tensor at the Earth's surface in the following way. For this purpose, consider a point of the Earth's surface with a known position such as $\mathbf{r} \in \mathbb{R}^3$; here we call this point the computation point. For the computation point, consider a neighbourhood such as an open-ball of radius $\delta$ centred at $\mathbf{r}$, denoted by $B_\delta(\mathbf{r})$. The average of the gradient tensor in the neighbourhood $B_\delta(\mathbf{r})$ is defined as follows:

**Definition 3** Let $\Omega_0$ be the set of points of the Earth's surface and let $\mathbf{f}_g$ have continuous derivatives (i.e. a continuous gradient tensor) at each point of $B_\delta(\mathbf{r}) - \Omega_0$. Then, the average of the gradient tensor in the neighbourhood $B_\delta(\mathbf{r})$ is defined by

$$ \bar{\Gamma}_{\delta}(\mathbf{r}) = \frac{1}{v} \oint_{B_{\delta}(\mathbf{r}) - \Omega_0} \Gamma dv \quad (7) $$

where $v$ is the volume of the neighbourhood $B_\delta(\mathbf{r})$.

So, according to the above definition, the gradient tensor at the Earth's surface may be defined as follows:

**Definition 4** Let $\Omega_0$ be the set of points of the Earth's surface and let $\lim_{\delta \to 0^+} \bar{\Gamma}_\delta(\mathbf{r})$ exist at $\mathbf{r} \in \Omega_0$. Then, the gravitational gradient tensor at $\mathbf{r}$ is defined as follows:

$$ \Gamma(\mathbf{r}) = \Gamma_0(\mathbf{r}) = \lim_{\delta \to 0^+} \bar{\Gamma}_\delta(\mathbf{r}) \quad (8) $$

Moreover, according to the piecewise-continuity of $\rho$, the average of the Earth's mass densities enclosed within $B_\delta(\mathbf{r})$ can be defined by the following definition:

**Definition 5** Let $\Omega_{in}$ be the Earth's interior space. The average of the Earth's mass density in $B_\delta(\mathbf{r}) \cap \Omega_{in}$, denoted by $\bar{\rho}_\delta(\mathbf{r})$, is defined by

$$ \bar{\rho}_{\delta}(\mathbf{r}) = \frac{1}{v_{\text{Mass}}} \oint_{B_{\delta}(\mathbf{r}) \cap \Omega_{\text{in}}} \rho \, dv \quad (9) $$

where $v_{\text{Mass}}$ is the volume of the Earth's masses enclosed in $B_\delta(\mathbf{r}) \cap \Omega_{in}$ (i.e. $v_{\text{Mass}} = \iiint_{B_\delta(\mathbf{r}) \cap \Omega_{in}} dv$; and so, we have $0 \le v_{\text{Mass}} \le v$). Hence, it is clear that if the mass-density function $\rho$ is continuous in $B_\delta(\mathbf{r}) \cap (\Omega_{in} \cup \Omega_0)$, we will have $\rho(\mathbf{r}) = \lim_{\delta \to 0^+} \bar{\rho}_\delta(\mathbf{r})$.

Now, according to the above definitions, Theorem 1 can be developed at the Earth's surface by the following theorem.

**Theorem 2** Let $\mathbf{r}$ belong to the Earth's surface $\Omega_0$, and suppose there exists $\delta > 0$ which allows $\rho$ to satisfy the Holder condition at every point of $B_\delta(\mathbf{r}) - \Omega_0$. Then, the average of the gradient tensor satisfies the following equation:

$$ \mathrm{trace}(\bar{\Gamma}_\delta) = -k_\delta \pi G \bar{\rho}_\delta \quad (10) $$

where $k_\delta$ is called the topographic coefficient and it is defined by

$$ k_\delta = 4 \times \frac{v_{\text{Mass}}}{v} \quad (11) $$

and the topographic coefficient $k_\delta$ as a function of $\mathbf{r}$ and $\delta$, depends on both the shape of the Earth's surface and the shape of $B_\delta(\mathbf{r})$. Moreover, since $0 < v_{\text{Mass}} < v$, we have $0 < k_\delta < 4$. Furthermore, if $\rho$ is also continuous in $B_\delta(\mathbf{r}) \cap (\Omega_{in} \cup \Omega_0)$ and if both $\lim_{\delta \to 0^+} k_\delta$ and $\lim_{\delta \to 0^+} \bar{\Gamma}_\delta(\mathbf{r})$ exists, the Poisson PDE at the Earth's surface (i.e. eq. 1)

will be finally obtained as follows:

$$ \mathrm{trace}(\Gamma) = -k\pi G\rho \quad (12) $$

where $k := \lim_{\delta \to 0^+} k_\delta$.
Proof: according to Definitions 2 and 3, we have

$$
\begin{align}
\operatorname{trace}(\bar{\Gamma}_{\delta}) &= \frac{1}{v} \oint_{B_{\delta}(\mathbf{r}) \cap \Omega_{0}} \operatorname{trace}(\Gamma) \, dv \nonumber \\
&= \frac{1}{v} \oint_{B_{\delta}(\mathbf{r}) \cap \Omega_{0}} \nabla^2 V \, dv \tag{13}
\end{align}
$$

which can be split into two integrals as follows:

$$
\operatorname{trace}(\bar{\Gamma}_{\delta}) = \frac{1}{v} \oint_{B_{\delta}(\mathbf{r}) \cap \Omega_{\text{in}}} \overbrace{\nabla^2 V}^{-4\pi G\rho} \, dv + \frac{1}{v} \oint_{B_{\delta}(\mathbf{r}) \cap \Omega_{\text{ex}}} \overbrace{\nabla^2 V}^{0} \, dv \quad (14)
$$

where $\Omega_{\text{ex}} = \mathbb{R}^3 - \Omega_{\text{in}} - \Omega_0$ is the exterior space of the Earth.

In eq. (14), we see that the right-hand side of eq. (13) has been split into two integrals: the first one is an integration over the intersection of the neighbourhood and the Earth's interior space, that is, $B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}$; and the second is an integration over the intersection of the neighbourhood and the Earth's exterior space, that is, $B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}$.

According to Theorem 1, in the Earth's exterior space where the Earth's mass density is zero, the right-hand side of the Poisson PDE (i.e. $-4\pi G\rho$) is trivial. It means that in the Earth's exterior space, the mass density of the atmosphere is supposed to be equal to zero in comparison with the mass density of the Earth's rocks. So, in eq. (14), the second integration over $B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}$ belonging to the Earth's exterior space will be equal to zero.

Moreover, according to the assumptions as well as Theorem 1, since $\rho$ is supposed to satisfy the Holder condition at each point of $B_\delta(\mathbf{r}) - \Omega_0$, the gradient tensor exists at each point of $B_\delta(\mathbf{r}) \cap \Omega_{\text{in}} \subset B_\delta(\mathbf{r}) - \Omega_0$ and also satisfies the Poisson PDE. Hence, according to both Theorem 1 and Definition 5, eq. (14) changes into eq. (10) as follows:

$$
\begin{align}
\operatorname{trace}(\bar{\Gamma}_\delta) &= \frac{-4\pi G}{v} \underbrace{\oint \rho \, dv}_{\substack{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}} \\ v_{\text{Mass}} \bar{\rho}_\delta}} \nonumber \\
&= -\underbrace{\left(4 \times \frac{v_{\text{Mass}}}{v}\right)}_{k_\delta} \pi G \bar{\rho}_\delta \tag{15}
\end{align}
$$

Now, according to the assumptions and Definition 4, if $\delta \to 0^+$, we will have $\operatorname{trace}(\bar{\Gamma}_\delta) \to \operatorname{trace}(\Gamma)$, and so, the proof of Theorem 2 will be completed as follows:

$$
\operatorname{trace}(\Gamma) = -4 \lim_{\substack{\delta \to 0^+ \\ \underbrace{k}_{\rho}}} \left( \frac{v_{\text{Mass}}}{v} \right) (\pi G) \lim_{\substack{\delta \to 0^+ \\ \overbrace{\rho}^{\bar{\rho}_{\delta}}}} \bar{\rho}_{\delta} \quad (16)
$$

3 APPROXIMATION OF THE MASS
DENSITY VIA THE POISSON PDE AT
THE EARTH'S SURFACE, USING THE
LAND-BASED GRAVITY DATA

In this section, we propose a method for extracting the lateral mass-
density variation of the terrain from land-based gravity data, without
using any prior knowledge about the mass density. The interesting
feature of the method is that it has directly employed the Poisson
PDE for the purpose of the mass density estimation at the Earth's
surface. For this purpose, certain required formulae are proposed to
convert land-based gravity data into the mass density. The formulae
are interestingly based on the aforementioned version of the Poisson

<figure><img src="image_3.png" /><figcaption>**Figure 2.** A schematic figure representing the neighbourhood, the computation point and its neighbour points.</figcaption></figure>

PDE at the Earth's surface; and they are developed for areas with
a sufficient rough topography. So, herein it is supposed that the
Earth's gravitational vector is observed at a rough part of the Earth's
surface.

Hence, first consider a computation point with a known posi-
tion vector such as **r** belonging to the Earth's surface, at which
the Earth's mass density is required to be estimated. Then as-
sume a neighbourhood for the computation point such as $B_\delta(\mathbf{r})$
which is defined here by a $2\delta \times 2\delta \times 2\delta$ - cube such as Cube(**r**,**δ**)
$= \{\mathbf{r} \in \mathbb{R}^3 \mid \|\mathbf{r} - \mathbf{r}_\infty\| < \delta\}$. Now, consider *N* points of the Earth's
surface, located in this neighbourhood, as the neighbours of the
computation point. According to the presumed roughness of the
Earth's surface, the neighbour points are expected to be nearly
3D-scattered around the computation point, see Fig. 2 illustrating
a schematic representation of the neighbourhood, the computation
point and its neighbour points.

Now, assume that **f**<sub>g</sub> is an observed value of the Earth's gravita-
tional vector at the known position of the computation point **r**. In
addition, suppose that the position vectors of all the *N* neighbour
points are known and their gravitational vectors are also observed.
The observed gravitational vectors and the known positions of the
*N* neighbour points are considered as follows (N ≥ 3):

$$
\begin{align}
\mathbf{f}_g(\mathbf{r}_n) &= \mathbf{f}_{gn} \nonumber \\
\mathbf{r}_n &\in \text{Cube} (\mathbf{r}, \delta) \cap \Omega_0 \nonumber \\
n &= 1, 2, \dots, N \tag{17}
\end{align}
$$

According to the aforementioned definition of the gradient tensor at the Earth's surface (i.e. Definition 4), the fact that the Earth's gravitational field is continuous at the Earth's surface implies the following lemma for the points of the Earth's surface.

**Lemma 1** Let $\Omega_0$ be the set of points of the Earth's surface and let $\Gamma$ be the gradient tensor defined via Definition 4. Also, let $\mathbf{f}_g$ and $\mathbf{f}_{gn}$ be the gravitational vectors observed on the Earth's surface at the computation point $\mathbf{r} \in \Omega_0$ and its neighbour point $\mathbf{r}_n \in \Omega_0$, respectively. Then, the gradient tensor satisfies the following equation for the computation point and its neighbour at the Earth's surface.

$$
\Delta \mathbf{f}_{gn} = \Gamma \Delta \mathbf{r}_n + \mathbf{e}_{\text{MDL}} (\Delta \mathbf{r}_n) \quad (18)
$$

where $n = 1, 2, \dots, N$, $\Delta \mathbf{r}_n = \mathbf{r}_n - \mathbf{r}$ and $\Delta \mathbf{f}_{gn} = \mathbf{f}_{gn} - \mathbf{f}_g$; the vector $\mathbf{e}_{\text{MDL}}$ is called ‘error of the model’ or ‘truncation error’. If $\Delta \mathbf{r}_n \to \mathbf{0}$, we will have $\frac{\mathbf{e}_{\text{MDL}}}{\|\Delta \mathbf{r}_n\|} \to \mathbf{0}$.
*Proof:* Although this Lemma seems so simple, its proof is slightly long; hence, the proof is presented by Appendix A.

Now, we are going to use the above linear model for approximating the gradient tensor (i.e. $\Gamma$) at the computation point $\mathbf{r}$. Hence, defining the following matrices is a requirement:

$$ 
\begin{aligned}
\Delta \mathbf{F}_g &= [\Delta \mathbf{f}_g^1 \ \Delta \mathbf{f}_g^2 \ \cdots \ \Delta \mathbf{f}_g^N] \in \mathbb{R}^{3 \times N} \\
\Delta \mathbf{R} &= [\Delta \mathbf{r}_1 \ \Delta \mathbf{r}_2 \ \cdots \ \Delta \mathbf{r}_N] \in \mathbb{R}^{3 \times N} \\
\mathbf{e}_{\text{MDL}}(\Delta \mathbf{R}) &= [\mathbf{e}_{\text{MDL}}(\Delta \mathbf{r}_1) \ \mathbf{e}_{\text{MDL}}(\Delta \mathbf{r}_2) \ \cdots \ \mathbf{e}_{\text{MDL}}(\Delta \mathbf{r}_N)] \in \mathbb{R}^{3 \times N}
\end{aligned}
\quad (19) $$

It means that $\{\Delta \mathbf{f}_g^n\}$, $\{\Delta \mathbf{r}_n\}$ and $\{\mathbf{e}_{\text{MDL}}(\Delta \mathbf{r}_n)\}$ are columns of $\Delta \mathbf{F}_g$, $\Delta \mathbf{R}$ and $\mathbf{e}_{\text{MDL}}(\Delta \mathbf{R})$, respectively. So, we will have

$$ \Delta \mathbf{F}_g = \Gamma \Delta \mathbf{r} + \mathbf{e}_{\text{MDL}}(\Delta \mathbf{R}) \quad (20) $$

where we will also have $\frac{\mathbf{e}_{\text{MDL}}(\Delta \mathbf{R})}{\|\Delta \mathbf{R}\|} \to \mathbf{0}$ if $\Delta \mathbf{R} \to \mathbf{0}$. By disregarding the error of the model in the above equation, we will arrive at the following approximation

$$ \Delta \mathbf{F}_g \approx \Gamma \Delta \mathbf{R} \quad (21) $$

In the following, the above approximation is used for approximating both the gradient tensor and the mass density at the Earth's surface.

As it was mentioned before, herein the Earth's surface has been supposed to be sufficiently rough or mountainous. So, the neighbour points can be supposed to be 3-D-scattered around the computation point, which implies the matrix $\Delta \mathbf{R} \in \mathbb{R}^{3 \times N}$ to be a full-rank matrix. Hence, its Moore-Penrose Pseudo Inverse is equal to $\Delta \mathbf{R}^{+} = \Delta \mathbf{R}^{\mathrm{T}}(\Delta \mathbf{R}\Delta \mathbf{R}^{\mathrm{T}})^{-1}$ satisfying $\mathbf{I} = \Delta \mathbf{R}\Delta \mathbf{R}^{+}$, for example, see Bjork (1988). Besides, it is clear that $\Delta \mathbf{R}^{+}$ (i.e. $\Delta \mathbf{R}^{\mathrm{T}}(\Delta \mathbf{R}\Delta \mathbf{R}^{\mathrm{T}})^{-1}$) is also a right inverse of $\Delta \mathbf{R}$. Hence, we will have the following approximation of the gradient tensor via right multiplying $\Delta \mathbf{R}^{+}$ by eq. (21):

$$ \Gamma \approx \Delta \mathbf{F}_g \Delta \mathbf{R}^{+}. \quad (22) $$

The above formula can be considered as a finite difference method for estimating the gradient tensor. Hence, the error of the model $\mathbf{e}_{\text{MDL}}$ can be considered as the truncation error of the finite differences. According to eqs (20) and (22), the effect of the truncation error on the approximated gradient tensor will be equal to $\mathbf{e}_{\text{MDL}}(\Delta \mathbf{R})\Delta \mathbf{R}^{+}$.

Consequently, the Earth's mass density at the computation point is finally approximated via a combination of eqs (22) and (12) as follows:

$$ \rho = \frac{1}{-k\pi G} \mathrm{trace}(\Gamma) \approx \frac{1}{-k\pi G} \mathrm{trace}(\Delta \mathbf{F}_g \Delta \mathbf{R}^+) \quad (23) $$

where $k \in (0, 4)$.

Also, according to eqs (20) and (23), the effect of the truncation error (i.e. $\mathbf{e}_{\text{MDL}}$) on the mass density obtained by eq. (23), will be equal to

$$ \frac{1}{-k\pi G} \mathrm{trace}(\mathbf{e}_{\text{MDL}}(\Delta \mathbf{R}) \Delta \mathbf{R}^+) \quad (24) $$

The above formula needs an adequate prediction about the truncation error $\mathbf{e}_{\text{MDL}}$ to present accurate knowledge about the effect of the truncation error on the estimated mass density. An accurate prediction about the effect of the truncation error can lead us to a suitable value of the neighbourhood dimension (i.e. $\delta$) as the scale of the finite difference method.

In a case study discussed in the next section, the above procedure of the mass density approximation is numerically examined to investigate the main idea of the paper.

<figure><img src="image_4.png" /><figcaption>Figure 3. The topography of the case study area, the black dots signify the locations of the gravity stations and the geodetic datum of this map is WGS-84.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. The spatial distribution of the observed gravity over the area of the case study.</figcaption></figure>

## 4 NUMERICAL EXPERIMENTS

According to the previous sections, the main idea of the paper is that the aforementioned version of the Poisson PDE is able to extract Earth's mass density from land-based gravity data, which was numerically examined here in a real case study. For this purpose, we applied the formulae, proposed by Section 3, to a set of real gravity data comprising 6350 gravity stations located along the north coast of the Persian Gulf in the south of the geographic region of Iran.

The spatial resolution of the data was approximately 0.5–1 km and 2–5 km in the directions of south–north and west–east, respectively. These observations were acquired by National Iranian Oil Company (NIOC) within a gravimetry project for certain petroleum exploration studies. Fig. 3 shows the locations of the gravity stations on a 90 m × 90 m-DEM over the area of the case study (this figure displays only 20 per cent of the locations to avoid mixing up the locations). It should be mentioned that the DEM was extracted from a version of the SRTM's DEM prepared by Reuter et al. (2007), in which topographic elevations are referenced to the WGS84/EGM96 geoid, see Farr et al. (2007).

Besides, Fig. 4 presents the spatial distribution of the observed gravity over the area of the case study.

According to the methodology described in Section 3, first we needed to have the Earth's gravitational vector at each gravity station. Nevertheless, the gravitational observation associated with each gravity station was only the gravity (i.e. only the norm of the Earth' gravity field). Hence, we had to accomplish the following steps to approximate the Earth's gravitational vector, that is, $\mathbf{f}_g = \nabla V$ at each station:

(1) Considering the observed gravity as the norm of the Earth's gravity vector at each station.

(2) Computing a reference gravitational vector at each station by means of both a set of fully normalized spherical harmonic coefficients of EIGEN-6C4 (Förste et al. 2014) and the standard series expansion of spherical harmonics with $n$EIGEN = 2190 as the maximum degree of the series expansion.
(3) Computing the Earth's centrifugal acceleration vector at each station by means of its geodetic position and the Earth's rotational velocity, for example, see Hofmann-Wellenhof & Moritz (2006, p. 43).

(4) Computing a reference gravity vector by adding the Earth's centrifugal acceleration vector (obtained via the third step) to the reference gravitational vector (obtained via second step).

(5) Considering the unit vector of the reference gravity vector (obtained via the fourth step) as the unit vector of the real gravity vector.

(6) Computing the Earth's gravity vector at each station, via multiplying the observed value of gravity by its unit vector obtained via the 5th step.

(7) Finally, obtaining the Earth's gravitational vector at each station by means of subtracting the centrifugal acceleration vector from the gravity vector computed in the sixth step.

In brief, the above procedures find the norm of $\mathbf{f}_g$ via the observation and obtain its direction via EIGEN-6C4. We can abridge the above procedures via the following formula:

$$ \mathbf{f}_g = g\mathbf{u}_g - \mathbf{f}_c \quad (25) $$

where g is the observed value of the Earth's gravity acquired at each gravity station; $\mathbf{u}_g$ is the unit vector of the Earth's gravity field, which is approximated via $\mathbf{u}_g \doteq \mathbf{g}^{\text{EIGEN-6C4}}/\|\mathbf{g}^{\text{EIGEN-6C4}}\|$; $\mathbf{g}^{\text{EIGEN-6C4}}$ is the Earth's reference gravity field obtained via $\mathbf{g}^{\text{EIGEN-6C4}} = \mathbf{f}_g^{\text{EIGEN-6C4}} + \mathbf{f}_c$; $\mathbf{f}_c$ is the Earth's centrifugal acceleration vector; and $\mathbf{f}_g^{\text{EIGEN-6C4}}$ is the Earth's gravitational field obtained by the set of fully normalized spherical harmonic coefficients $\{\tilde{C}_{nm}, \tilde{S}_{nm}\}$ of EIGEN-6C4, as well as standard series expansion of spherical harmonics, as follows:

$$ \mathbf{f}_g^{\text{EIGEN-6C4}} = \nabla \left( \frac{GM}{r} \sum_{n=0}^{n\text{EIGEN}} \left(\frac{a}{r}\right)^n \times \sum_{m=0}^{n} \left( \tilde{C}_{nm} \cos m\lambda + \tilde{S}_{nm} \sin m\lambda \right) \bar{P}_{nm}(\cos\theta) \right) \quad (26) $$

where n and m are degree and order of the harmonic coefficients, respectively; nEIGEN = 2190 indicating the maximum degree of the series expansion; GM (geocentric gravitational constant) and a (semi major axis) are the EIGEN-6C4 scaling parameters; $\bar{P}_{nm}$ is the fully normalized associated Legendre function; $(r, \theta, \lambda)$—(radius, geocentric co-latitude and geocentric longitude)—indicates the spherical coordinates; and $\nabla$ denotes the gradient differential operator in the Cartesian coordinate system.

The direction of the gravitational vector can be estimated also by knowledge about the deflections of the vertical. The components of the deflection of the vertical (i.e. the meridian component $\xi$ and the prime vertical component $\eta$), as well as the geodetic coordinate $(\varphi, \lambda)$, can immediately lead us to the astronomical coordinates $(\Phi, \Lambda)$ which are necessary for estimating the unit vector of the Earth's gravity field (i.e. $\mathbf{u}_g$) as follows, for example, see Hofmann-Wellenhof & Moritz (2006, p. 55):

$$ \mathbf{u}_g = \begin{bmatrix} \cos \Phi \cos \Lambda \\ \cos \Phi \sin \Lambda \\ \sin \Phi \end{bmatrix} \quad (27) $$

where $\Phi = \varphi + \xi$ and $\Lambda = \lambda + \frac{1}{\cos\varphi}\xi$ for all $\varphi \neq \pm\frac{\pi}{2}$. Herein, the spatial distribution of the deflection of the vertical in the area

<figure><img src="image_6.png" /><figcaption>Figure 5. The spatial distribution of the deflection of the vertical over the area of the case study, estimated via EIGEN-6C4: the deflection of the vertical is required for estimating the direction of the Earth's gravity vector at each station.</figcaption></figure>

of the case study, estimated via EIGEN-6C4, is illustrated by Fig. 5.

Obviously, the use of a GGM such as EIGEN-6C4 for estimating functionals of Earth's gravity field, such as the components of the deflection of the vertical, implies propagated and omission errors. The propagated error is because of any errors occurring at observations used for developing the GGM and also, the omission (or truncation) error is due to the fact that the resolution of the GGM is finite and the GGM lost a part of Earth's gravity field spectrum which extends to infinity (Pavlis et al. 2012). Hence, herein, we have investigated the accuracy of the deflection of the vertical evaluated by EIGEN-6C4. Table 1 has numerically examined EIGEN-6C4 on the accuracy of the deflection of the vertical by means of 9 stations located in the region of the Persian country. At these stations, the deflection of the vertical is known and they were derived from the geodetic and astronomical latitudes and longitudes, where the astronomical coordinates were observed by National Cartographic Center of Iran (NCC) for certain geoid computation purposes. Table 1 indicates that the RMS of the errors associated with the meridian and prime vertical deflection components (i.e. $\xi$ and $\eta$) was assessed at 1.2 and 2.7 arc second, respectively. The effect of these errors on the estimated mass density is investigated by Section 4.5.

According to all the formulae introduced by Section 3, first all the vectors necessarily needed to be presented in a Cartesian coordinate system. Herein, we investigated all the vectors in the geocentric Cartesian system corresponding to WGS-84 since all the positions were observed via GNSS-positioning in this coordinate system. However, we also used a local Cartesian coordinate system at which the neighbourhood Cube$(\mathbf{r}, \delta)$ is defined and used for seeking neighbours around each gravity station. At this local Cartesian coordinate system, the x-axis is directed toward the east, the z-axis is a tangent on the plumb line and outward from the Earth and the y-axis is chosen so that the system is right-handed.

Then, according to the methodology described by Section 3, each gravity station as a computation point needed a neighbourhood to choose as its neighbours from other gravity stations. Here, we defined the neighbourhood for each gravity station by means of Cube$(\mathbf{r}, \delta)$ introduced in Section 3. We selected the neighbours for each gravity station by various values of the neighbourhood dimension $\delta$. We had to take into consideration that despite the fact that the growth of $\delta$ increases the number of the neighbours around each gravity station, and consequently, despite the fact that the growth of $\delta$ increases the redundancy (i.e. degrees of freedom) corresponding to the problem, it decreases the certainty of the approximation. It is clear that the growth of $\delta$ causes the growth of the
<table><caption>**Table 1.** Evaluation of the errors corresponding to the deflection of the vertical obtained by EIGEN-6C4, via the astro-geodetic data at nine stations in the region of Iran.</caption><thead><tr><th></th><th>Station</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>RMS</th></tr></thead><tbody><tr><td rowspan="9">Deflection of the vertical [arcsec]</td><td>Astronomical coordinates [dms]</td><td>Φ</td><td>37.22489</td><td>36.09328</td><td>35.54490</td><td>36.39233</td><td>37.49447</td><td>35.45025</td><td>33.51365</td><td>31.57070</td><td>29.35576</td></tr><tr><td></td><td>Λ</td><td>48.01376</td><td>49.33235</td><td>53.44567</td><td>55.08258</td><td>56.43133</td><td>59.21166</td><td>59.37185</td><td>59.45344</td><td>60.10597</td></tr><tr><td>Geodetic coordinates [dms]</td><td>φ</td><td>37.22483</td><td>36.09364</td><td>35.54539</td><td>36.39250</td><td>37.49418</td><td>35.44564</td><td>33.51300</td><td>31.57038</td><td>29.35579</td></tr><tr><td></td><td>λ</td><td>48.01349</td><td>49.33090</td><td>53.44374</td><td>55.08147</td><td>56.43111</td><td>59.21158</td><td>59.37098</td><td>59.45297</td><td>60.10567</td></tr><tr><td rowspan="3">Astro-geodetic</td><td>ξ</td><td>0.6</td><td>-3.7</td><td>-5.0</td><td>-1.7</td><td>2.9</td><td>6.1</td><td>6.5</td><td>3.2</td><td>-0.3</td></tr><tr><td>η</td><td>2.1</td><td>11.7</td><td>15.6</td><td>8.9</td><td>1.7</td><td>0.6</td><td>7.2</td><td>4.0</td><td>2.6</td></tr><tr><td><b>EIGEN-6C4</b></td><td><b>ξ</b></td><td><b>0.7</b></td><td><b>-2.3</b></td><td><b>-2.6</b></td><td><b>-2.0</b></td><td><b>4.0</b></td><td><b>7.0</b></td><td><b>8.2</b></td><td><b>3.5</b></td><td><b>-0.2</b></td></tr><tr><td rowspan="3">Error</td><td>η</td><td><b>3.6</b></td><td><b>10.8</b></td><td><b>13.5</b></td><td><b>10.8</b></td><td><b>1.5</b></td><td><b>-0.6</b></td><td><b>6.0</b></td><td><b>10.6</b></td><td><b>0.2</b></td></tr><tr><td>ξ</td><td>0.1</td><td>1.3</td><td>2.3</td><td>-0.3</td><td>1.0</td><td>1.0</td><td>1.7</td><td>0.3</td><td>0.1</td><td><b>1.2</b></td></tr><tr><td>η</td><td>1.5</td><td>-1.0</td><td>-2.1</td><td>1.9</td><td>-0.2</td><td>-1.2</td><td>-1.2</td><td>6.6</td><td>-2.4</td><td><b>2.7</b></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The RMS of the effect of the truncation error against the neighbourhood dimension.</figcaption></figure>

error of model (i.e. the truncation error) in eq. (20). Fig. 6 illustrates a prediction about the effect of the truncation error on the estimated mass density. This prediction was obtained via eq. (24) as well as EIGEN-6C4. For this purpose, first the truncation error (i.e. $e_{MDL}$ introduced in eq. (20)) was predicted by EIGEN-6C4 and then it was applied to eq. (24). In Fig. 6, the root mean square (RMS) of the effect of the truncation error is plotted against the neighbourhood dimension $\delta$; this figure also shows that the growth of $\delta$ causes the growth of the effect of the truncation error on the estimated mass density.

By finding the neighbours of each gravity station, the two matrices $\Delta F_g$ and $\Delta R$ defined by eq. (19) could be straightforwardly made at each gravity station. Therefore, the mass density determination could be easily performed by means of eq. (23) for each station.

According to eq. (23), in order to avoid ill-conditioning of $\Delta R$, the Earth's surface around the computation point should be sufficiently rough or mountainous. Otherwise, the ill-conditioning of $\Delta R$ may cause an instability in the problem and then the instability also causes the errors to rise sharply in the solutions. Hence, here we decided to omit some gravity stations for which the condition num-

ber of $\Delta R$ was higher than a threshold such as $\tau$. High values of the threshold $\tau$ may strengthen the ill-conditioning and moreover, its low values neglect many gravity stations. Hence, it was required to take an optimum choice of the threshold $\tau$ that weakens the effect of the ill-conditioning on the solutions, and also minimizes the number of neglected gravity stations. For this purpose, we used the L-Curve parameter-choice method (Lawson & Hanson 1995) for simultaneously decreasing the effect of the ill-conditioning and minimizing the number of the neglected stations. Appendix B describes how the threshold $\tau = 86.7$ as the optimum choice has been obtained by the L-Curve method. Since the current numerical experiments were executed just to examine the capability of the proposed approach and the Poisson PDE at the Earth's surface adequately, the authors wittingly used neither a prior knowledge about the mass density nor additional constraints on the vertical gradients of the gravity field, except the abovementioned gravity-based data, that is, the land-based gravity observations and EIGEN-6C4. Hence, it was required to apply the approach to the area with a rough topography which enabled the approach to determine the vertical gradients of the gravity. However, the approach can also be practical for flat areas if it uses prior knowledge about either the vertical components of the gradient tensor or the mass density. Moreover, it is worth remembering that similarly, all the previous methods of estimating the mass density reduction are also developed for areas with rough topography to extract the mass density from the vertical gradient of land-based gravity data, see Nettleton (1939), Parasnis (1952, 1979), Fukao et al. (1981), Murata (1993) and Nawa et al. (1997).

Furthermore, according to the definition of $k$, this coefficient is expected to be typically around 2 at the Earth's surface. Hence, herein, it was roughly considered that $k \approx 2$. However, based on the definition of the topographic coefficient $k$ mentioned by Theorem 2, Appendix C discusses a method of assessing this coefficient by means of a terrain model such as the SRTM's DTM. In this appendix, the true values of $k$ corresponding to the region of the case study were evaluated and applied to the proposed approach. Appendix C indicates that the rough approximation $k \approx 2$ was enough for the current research to show the ability of the Poisson PDE to estimate the mass density. However, it shows that the precise values of $k$, which are straightforwardly estimated by the terrain model, are able to improve the precision of the estimated mass density.
**Algorithm 1.** A summary of the proposed approach

**Input:** Land-based gravity observations at *n* gravity stations ($g_{i=1,2,...,n}$); geodetic coordinates of the gravity stations (($\phi, \lambda, h)_{i=1,2,...,n}$); neighborhood dimension ($\delta$); threshold of the condition number ($\tau$); Earth's rotational velocity ($\omega$); universal gravitation constant ($G$); a GGM such as EIGEN-6C4;

**Output:** Mass density and gradient tensor at all the stations (i.e. $\rho_{i=1,2,...,n}$ and $\Gamma_{i=1,2,...,n}$);

**Step 1: Preparing the requirements:**

(Preparing the Earth's gravitational vector and the position vector at each station, i.e. $\mathbf{f}_{g,i=1,2,...,n}$ and $\mathbf{r}_i=1,2,...,n$)

For $i=1,2,...,n$ do:

**Step 1.1:** Change the geodetic coordinates ($\phi, \lambda, h)_i$ into the Cartesian coordinates $\mathbf{r}_i = (x_i, y_i, z_i)$ in WGS-84;

**Step 1.2:** Use $\mathbf{r}_i$, $\omega$ and the following equation for assessing the centrifugal acceleration vector at the $i^{th}$ station ($\mathbf{f}_{c,i}$), see (Hofmann-Wellenhof and Moritz, 2006, p. 43);
$$ \mathbf{f}_{c,i} = \omega^2 \mathbf{r}'_i, \quad \mathbf{r}'_i = (x_i, y_i, 0) $$

**Step 1.3:** Apply $g_i$, $\mathbf{f}_{c,i}$ and the GGM to Eq.(25) for assessing the Earth's gravitational vector at the $i^{th}$ gravity station (i.e. $\mathbf{f}_{g,i}$);

End for;

**Step 2: Estimating the unknowns:**

(Estimating the mass density and gradient tensor at each station, i.e. $\rho_{i=1,2,...,n}$ and $\Gamma_{i=1,2,...,n}$)

For $i=1,2,...,n$ do:

**Step 2.1:** Transform all the position and gravitational vectors from WGS-84 into the local Cartesian coordinate system at the $i^{th}$ station;

**Step 2.2:** Choose neighbors of the $i^{th}$ station from other gravity stations by means of the neighborhood Cube ($\mathbf{r}_i, \delta$);

**Step 2.3:** Construct the following matrices via the differences between the gravitational vectors and position vectors of the neighbors (denoted by $\mathbf{f}_{g,i}^j=1,2,...,N_i$ and $\mathbf{r}_j^j=1,2,...,N_i$, respectively) and the ones of the $i^{th}$ station ($\mathbf{f}_{g,i}, \mathbf{r}_i$) (see Section 3);

$$ \Delta \mathbf{F}_{g,i} = [ \Delta \mathbf{f}_{g,i}^1 \quad \Delta \mathbf{f}_{g,i}^2 \quad \dots \quad \Delta \mathbf{f}_{g,i}^{N_i} ] \in \mathbb{R}^{3 \times N_i} $$

$$ \Delta \mathbf{R}_i = [ \Delta \mathbf{r}_i^1 \quad \Delta \mathbf{r}_i^2 \quad \dots \quad \Delta \mathbf{r}_i^{N_i} ] \in \mathbb{R}^{3 \times N_i} $$

$$ \Delta \mathbf{f}_{g,i}^j = \mathbf{f}_{g,i}^j - \mathbf{f}_{g,i}, \quad \Delta \mathbf{r}_i^j = \mathbf{r}_i^j - \mathbf{r}_i, \quad j = 1, 2, \dots, N_i $$

**Step 2.4:** Estimate the condition number of $\Delta \mathbf{R}_i$ (i.e. $\text{cond}(\Delta \mathbf{R}_i)$);

**Step 2.5:** If $\text{cond}(\Delta \mathbf{R}_i) < \tau$ then

**Step 2.5.1:** Estimate the gradient tensor corresponding to the $i^{th}$ station as follows (see Eq.(21)):

$$ \Gamma_i \approx \Delta \mathbf{F}_{g,i} \Delta \mathbf{R}_i^+; $$

**Step 2.5.2:** Set $k_i = 2$ or estimate it via the method discussed in Appendix C;

**Step 2.5.3:** Estimate the mass density corresponding to the $i^{th}$ station as follows (see Eq.(23)):

$$ \rho_i = \frac{1}{-k_i \pi G} \text{trace}(\Gamma_i); $$

End If;

End for;
<figure><img src="image_8.png" /><figcaption>Figure 7. The histograms of the estimated mass density obtained by various values of the neighbourhood dimension at all the gravity stations.</figcaption></figure>

Accordingly, the mass density was estimated at all the gravity stations by means of eq. (23) as well as the above-mentioned procedures. Algorithm 1 summarizes the above-mentioned procedures via two main steps. The first step prepares the requirements including the Earth's gravitational vector field at each gravity station, as well as the corresponding position vectors and then the second step estimates the mass density and gradient tensor at each gravity station, using the Earth's gravitational vectors and the position vectors as well as formulae proposed by Section 3.

Fig. 7 illustrates the histograms of the estimated mass density for various values of $\delta$. Also, the spatial distribution of the estimated mass density over the region of the case study is illustrated by Fig. 8. It is interesting to remember that these results have been obtained only via the land-based gravity data as well as the GGM without using any prior knowledge about the mass density.

## 4.1 Validation of the results via geological rock-samplings

In the area of the case study, certain rock-samplings were done to directly determine the mass densities corresponding to geological formations. 86 rock-samplings were gathered in 13 rock-sampling stations while the gravity data was acquiring at the 6350 gravity stations. Immediately the samples were collected from out-cropped geological formations around each rock-sampling station, their corresponding mass densities were directly measured by a digital solid mass density meter with resolution 0.01 g cm$^{-3}$ (NIOC 2004). The sample values of the mass density corresponding to the 13 sampling stations are included in Table 2. Herein, we utilized these sample

values as a criterion for directly assessing the results while we knew the facts that this criterion might have certain drawbacks because the outcrops might not be truly representative of the formations below, see Parasnis (1952) and Murata (1993).

These 13 rock-sampling stations did not exactly coincide with any gravity stations. Hence, for comparison of the above sample values with the estimated mass density obtained via the methodology of this paper, we had to interpolate the estimated mass density at each rock-sampling station. For this purpose, each rock-sampling station needed to find its neighbours among the gravity stations by means of defining the neighbourhood Cube$(r, \delta)$. Then the neighbours of each rock-sampling station were used to interpolate the mass density via a simple weighted arithmetic mean. Herein, the Gaussian kernel was used to assign a weight to each neighbour, where the bandwidth of the Gaussian kernel was supposed to be equal to the neighbourhood dimension $\delta$. Therefore, we could interpolate the estimated mass density in the 13 rock-sampling stations via the arithmetic mean. Then we compared them with the average values of the mass density obtained via the rock samplings, see Table 3.

Table 3 includes the mass densities estimated by various values of the neighbourhood dimension at all the sampling stations. Also, Table 3 shows the RMS of the errors corresponding to the mass densities obtained by various values of $\delta$; this table shows that the minimum value of the RMS is 157 kg m$^{-3}$ that occurred at $\delta = 4.5$ km.

Moreover, Table 4 presents all the 13 sampling stations with the error of the estimated mass density obtained by $\delta = 4.5$ km. In Table 4, we see that the minimum value (the best value) of the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The spatial distribution of the estimated mass density in the area of the case study, obtained by various values of the neighbourhood dimension.</figcaption></figure>

<table><caption>**Table 2.** 86 sample values of the Earth's mass density in 13 rock-sampling stations.</caption><thead><tr><th></th><th>S1</th><th>S2</th><th>S3</th><th>S4</th><th>S5</th><th>S6</th><th>S7</th><th>S8</th><th>S9</th><th>S10</th><th>S11</th><th>S12</th><th>S13</th></tr></thead><tbody><tr><td>Rock-sampling stations</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>φ (wgs84)</td><td>27.0641</td><td>27.0718</td><td>27.0670</td><td>27.0839</td><td>27.0826</td><td>26.8531</td><td>26.8532</td><td>26.8016</td><td>26.7961</td><td>26.8387</td><td>26.8224</td><td>27.0088</td><td>27.0093</td></tr><tr><td>λ (wgs84)</td><td>53.6216</td><td>53.6287</td><td>53.6323</td><td>53.6383</td><td>53.6385</td><td>53.9111</td><td>53.9114</td><td>54.2826</td><td>54.2870</td><td>54.3889</td><td>54.4026</td><td>54.5821</td><td>54.5821</td></tr><tr><td>Sample values</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>1</td><td>2050</td><td>2360</td><td>2580</td><td>2450</td><td>2440</td><td>1930</td><td>1690</td><td>2030</td><td>2140</td><td>1950</td><td>1940</td><td>2380</td><td>2310</td></tr><tr><td></td><td>2</td><td>2060</td><td>2320</td><td>2260</td><td>2390</td><td>2430</td><td>2040</td><td>2200</td><td>2020</td><td>2120</td><td>1900</td><td>1930</td><td>2380</td><td>2290</td></tr><tr><td></td><td>3</td><td>2090</td><td>2320</td><td>2570</td><td>2410</td><td>2430</td><td>2160</td><td>1500</td><td>2040</td><td>2130</td><td>1800</td><td></td><td>2380</td><td>2290</td></tr><tr><td></td><td>4</td><td>2120</td><td>2420</td><td>2600</td><td>2450</td><td>2450</td><td>1960</td><td>2110</td><td></td><td>2070</td><td>2010</td><td></td><td>2350</td><td>2280</td></tr><tr><td></td><td>5</td><td>2060</td><td>2370</td><td>2600</td><td>2390</td><td>2440</td><td>1830</td><td>1650</td><td></td><td>2140</td><td>1970</td><td></td><td>2390</td><td>2280</td></tr><tr><td></td><td>6</td><td></td><td>2560</td><td>2270</td><td></td><td></td><td>1930</td><td>1670</td><td></td><td></td><td>1950</td><td></td><td>2350</td><td>2300</td></tr><tr><td></td><td>7</td><td></td><td></td><td>2610</td><td></td><td></td><td></td><td>2210</td><td></td><td></td><td>1950</td><td></td><td>2340</td><td>2290</td></tr><tr><td></td><td>8</td><td></td><td></td><td>2610</td><td></td><td></td><td></td><td>2270</td><td></td><td></td><td>2050</td><td></td><td>2370</td><td>2310</td></tr><tr><td></td><td>9</td><td></td><td></td><td>2520</td><td></td><td></td><td></td><td>2100</td><td></td><td></td><td>2040</td><td></td><td>2340</td><td>2290</td></tr><tr><td></td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td><td>1650</td><td></td><td></td><td>1780</td><td></td><td>2360</td><td>2270</td></tr><tr><td>The average of the sample values</td><td>2076</td><td>2392</td><td>2513</td><td>2418</td><td>2438</td><td>1975</td><td>1905</td><td>2030</td><td>2120</td><td>1940</td><td>1935</td><td>2364</td><td>2291</td></tr><tr><td>The standard deviation</td><td>29</td><td>90</td><td>144</td><td>30</td><td>8</td><td>113</td><td>296</td><td>10</td><td>29</td><td>91</td><td>7</td><td>18</td><td>13</td></tr></tbody></table>

error is equal to $10 \text{ kg m}^{-3}$ occurred at the fifth station. Also, Table 4 shows that the large values of the error occurred at the third, sixth and seventh stations. Comparing Table 2 with Table 3, we can easily realize that the large errors are expected to occur at these three stations because of their large standard deviations.

## 4.2 Comparison of the results with mass density provided by the F-H method

In this step, we also used the F-H conventional method of Parasnis (1952) for estimating mass density at the 13 rock-sampling stations,
<table><caption>**Table 3.** The estimated mass density corresponding to the 13 sampling-stations; in this table, the estimated mass density was compared with the average of the sample values mentioned in Table 2 (unit: kg m⁻³).</caption><thead><tr><th>Station</th><th>ρ<br>The average of<br>the sample values<br>mentioned by Table 2</th><th>ρ<br>(δ = 2 km)</th><th>ρ<br>(δ = 3 km)</th><th>ρ<br>(δ = 4 km)</th><th>ρ<br>(δ = 4.5 km)</th><th>ρ<br>(δ = 5 km)</th><th>ρ<br>(δ = 6 km)</th><th>ρ<br>(δ = 7 km)</th></tr></thead><tbody><tr><td>S1</td><td><b>2076</b></td><td>2935</td><td>2520</td><td>2299</td><td><b>2239</b></td><td>2171</td><td>2127</td><td>2152</td></tr><tr><td>S2</td><td><b>2392</b></td><td>2939</td><td>2512</td><td>2380</td><td><b>2302</b></td><td>2209</td><td>2151</td><td>2159</td></tr><tr><td>S3</td><td><b>2513</b></td><td>3001</td><td>2518</td><td>2372</td><td><b>2292</b></td><td>2204</td><td>2152</td><td>2165</td></tr><tr><td>S4</td><td><b>2418</b></td><td>3165</td><td>2662</td><td>2532</td><td><b>2470</b></td><td>2283</td><td>2184</td><td>2147</td></tr><tr><td>S5</td><td><b>2438</b></td><td>3121</td><td>2639</td><td>2512</td><td><b>2448</b></td><td>2274</td><td>2181</td><td>2151</td></tr><tr><td>S6</td><td><b>1975</b></td><td>2219</td><td>2251</td><td>2272</td><td><b>2259</b></td><td>2272</td><td>2279</td><td>2289</td></tr><tr><td>S7</td><td><b>1905</b></td><td>2219</td><td>2251</td><td>2272</td><td><b>2261</b></td><td>2272</td><td>2279</td><td>2289</td></tr><tr><td>S8</td><td><b>2030</b></td><td>2146</td><td>2080</td><td>2123</td><td><b>2115</b></td><td>2082</td><td>2179</td><td>2110</td></tr><tr><td>S9</td><td><b>2120</b></td><td>1859</td><td>2012</td><td>2062</td><td><b>2020</b></td><td>2023</td><td>2139</td><td>2068</td></tr><tr><td>S10</td><td><b>1940</b></td><td>2619</td><td>1804</td><td>1847</td><td><b>1901</b></td><td>1924</td><td>1964</td><td>2042</td></tr><tr><td>S11</td><td><b>1935</b></td><td>3633</td><td>1958</td><td>2001</td><td><b>1901</b></td><td>1974</td><td>1904</td><td>2032</td></tr><tr><td>S12</td><td><b>2364</b></td><td>2211</td><td>2251</td><td>2273</td><td><b>2275</b></td><td>2338</td><td>2362</td><td>2377</td></tr><tr><td>S13</td><td><b>2291</b></td><td>2210</td><td>2250</td><td>2273</td><td><b>2275</b></td><td>2338</td><td>2361</td><td>2376</td></tr><tr><td>RMS (kg m⁻³)</td><td></td><td>674</td><td>207</td><td>163</td><td><b>157</b></td><td>180</td><td>210</td><td>218</td></tr></tbody></table>

<table><caption>**Table 4.** The comparison of the mass density estimated by eq. (23) with the average of the sample values at the 13 sampling stations, when δ = 4.5 km (unit: kg m⁻³).</caption><thead><tr><th>Rock-sampling</th><th>Station</th><th>S1</th><th>S2</th><th>S3</th><th>S4</th><th>S5</th><th>S6</th><th>S7</th><th>S8</th><th>S9</th><th>S10</th><th>S11</th><th>S12</th><th>S13</th></tr><tr><th></th><th>ρ (average)</th><th>2076</th><th>2392</th><th>2513</th><th>2418</th><th>2438</th><th>1975</th><th>1905</th><th>2030</th><th>2120</th><th>1940</th><th>1935</th><th>2364</th><th>2291</th></tr></thead><tbody><tr><td>ρ ≈ <sup>1</sup>&frasl;<sub>&minus;2&pi;&pi;G</sub> trace(&Delta;F<sub>g</sub>&Delta;R<sup>+</sup>)</td><td></td><td>2239</td><td>2302</td><td>2292</td><td>2470</td><td>2448</td><td>2259</td><td>2261</td><td>2115</td><td>2020</td><td>1901</td><td>1901</td><td>2275</td><td>2275</td></tr><tr><td>Error (kg m<sup>&minus;3</sup>)</td><td></td><td><b>163</b></td><td><b>&minus;90</b></td><td><b>&minus;221</b></td><td><b>52</b></td><td><b>10</b></td><td><b>284</b></td><td><b>356</b></td><td><b>85</b></td><td><b>&minus;100</b></td><td><b>&minus;39</b></td><td><b>&minus;34</b></td><td><b>&minus;89</b></td><td><b>&minus;16</b></td></tr></tbody></table>

and then we compared its results with the above-mentioned results obtained by the proposed method of the current research. The Parasnis method as a version of the Nettleton method (Nettleton 1939) can be considered as the basis of other methods of the mass density reduction, cf. Fukao et al. (1981), Murata (1993) and Nawa et al. (1997). The Parasnis method can extract the mass density from the derivative of the free-air gravity anomaly with respect to the topographic correction coefficient (H), and hence, it is usually called ‘F-H’ method. Similar to the proposed method of our current research, the Parasnis’ method is also suited for areas with a sufficient rough topography and Parasnis (1952) suggested that the simplest way to apply his method is to make the measurements at various points on the slope of a hill.

The Parasnis method (i.e. the F-H method) is based on an assumption that the Bouguer gravity anomaly is a random error, which implies the following minimization estimating the terrain mass density ρ (Murata 1993):

$$ \sum_{n=1}^{N} \| F_n - \rho H_n - \tilde{g}_{\text{Bou}} \|^2 \rightarrow \text{Min} \qquad (28) $$

where $F_n$ and $H_n = c_{Bo}h_n - T_n/\rho_0$ are the free-air anomaly and the topographic correction coefficient at the neighbour stations $\{(x_n, y_n)\}(n = 1, 2, \dots, N)$, respectively; $\tilde{g}_{\text{Bou}}$ is the mean value of the Bouguer anomalies; $c_{Bo}$ is also the so-called ‘simple Bouguer slab’ correction due to nth gravity station; and $T_n$ is the apparent topographic correction. The above optimization indicates that the F-H method of Parasnis (1952) certainly intends to calculate the mass density via the gradient of the Free-air anomaly with respect to the topographic correction coefficient H. Hence, Parasnis (1952) proposed the following simple formula for extracting the terrain mass density from land-based gravity data (Parasnis 1952, 1979):

$$ \Delta g - \Delta \gamma_0 + c_{\text{FA}} \Delta h = \rho \Delta H \qquad (29) $$

where Δg is the difference between the observed values of gravity at two gravity stations; Δγ₀ is the difference between the calculated values of gravity as given by the International Gravity Formula (i.e. the normal gravity field at the reference ellipsoid) at the latitude of the two station; Δh is also the difference between their orthometric heights (elevations) above the sea level; $c_{FA}\Delta h$ is the difference between the so-called ‘free air’ corrections due to the elevations; $\Delta H = c_{Bo}\Delta h - \Delta T/\rho_0$ is the difference between the topographic correction coefficients at the two stations, $c_{Bo}\Delta h$ is also the difference between the so-called ‘simple Bouguer slab’ corrections; and ΔT is the difference between the apparent topographic corrections when it is calculated for an assumed mass density of $\rho_0$ (see Parasnis 1952) (herein, for calculating the apparent topographic corrections, we employed a version of the SRTM’s 90 m × 90 m − DEM provided by Reuter et al. 2007). Accordingly, it is clear that the left-hand side of eq. (29) is the difference of the free-air gravity anomalies at the two stations (i.e. $\Delta F = \Delta g - \Delta \gamma_0 + c_{FA}\Delta h$); and its right-hand side is a linear function of the mass density ρ with respect to ΔH (Parasnis 1952).

Herein, first we made the Parasnis’s equation (i.e. eq. (29)) between each gravity station and its neighbours located within the neighbourhood Cube(r, δ), and then the F-H method could perform the mass density estimation for each gravity station by means of solving all the associated Parasnis equations simultaneously. Hence, the mass density could also be estimated by the F-H method at the 13 rock-sampling stations. Table 5 includes the estimated mass densities obtained by the F-H method at the 13 gravity stations for various values of the neighbourhood dimension δ. This table indicates that the best value of the RMS of the errors of the mass density estimated by the F-H method is equal to 342 kg m⁻³ that occurred at δ = 6 km.

Besides, Table 6 presents all the 13 sampling stations with the error of the mass density estimated by the F-H method when δ = 6 km.
<table><caption>**Table 5.** The mass density, estimated by the F-H method at the 13 sampling-stations for various value of the neighbourhood dimension; in this table, the estimated mass density was compared with the average of the sample values mentioned in Table 2 (unit: kg m⁻³).</caption><thead><tr><th>Station</th><th>ρ<br>The average of<br>the sample values<br>mentioned in Table 2</th><th>ρ<br>(δ = 3 km)</th><th>ρ<br>(δ = 4 km)</th><th>ρ<br>(δ = 5 km)</th><th>ρ<br>(δ = 6 km)</th><th>ρ<br>(δ = 7 km)</th><th>ρ<br>(δ = 8 km)</th><th>ρ<br>(δ = 9 km)</th></tr></thead><tbody><tr><td>S1</td><td><strong>2076</strong></td><td>3033</td><td>2691</td><td>2421</td><td>2314</td><td>2300</td><td>2314</td><td>2342</td></tr><tr><td>S2</td><td><strong>2392</strong></td><td>3129</td><td>2777</td><td>2440</td><td>2315</td><td>2298</td><td>2305</td><td>2334</td></tr><tr><td>S3</td><td><strong>2513</strong></td><td>3083</td><td>2727</td><td>2419</td><td>2314</td><td>2300</td><td>2311</td><td>2337</td></tr><tr><td>S4</td><td><strong>2418</strong></td><td>3356</td><td>2886</td><td>2506</td><td>2317</td><td>2289</td><td>2294</td><td>2322</td></tr><tr><td>S5</td><td><strong>2438</strong></td><td>3327</td><td>2873</td><td>2494</td><td>2316</td><td>2290</td><td>2295</td><td>2322</td></tr><tr><td>S6</td><td><strong>1975</strong></td><td>2395</td><td>2364</td><td>2384</td><td>2400</td><td>2413</td><td>2437</td><td>2446</td></tr><tr><td>S7</td><td><strong>1905</strong></td><td>2395</td><td>2359</td><td>2384</td><td>2401</td><td>2413</td><td>2437</td><td>2446</td></tr><tr><td>S8</td><td><strong>2030</strong></td><td>2267</td><td>2289</td><td>2265</td><td>2314</td><td>2321</td><td>2368</td><td>2390</td></tr><tr><td>S9</td><td><strong>2120</strong></td><td>2219</td><td>2220</td><td>2240</td><td>2268</td><td>2302</td><td>2366</td><td>2382</td></tr><tr><td>S10</td><td><strong>1940</strong></td><td>1948</td><td>2515</td><td>2644</td><td>2589</td><td>2603</td><td>2589</td><td>2603</td></tr><tr><td>S11</td><td><strong>1935</strong></td><td>2205</td><td>2503</td><td>2853</td><td>2566</td><td>2646</td><td>2639</td><td>2628</td></tr><tr><td>S12</td><td><strong>2364</strong></td><td>2343</td><td>2393</td><td>2443</td><td>2470</td><td>2490</td><td>2504</td><td>2516</td></tr><tr><td>S13</td><td><strong>2291</strong></td><td>2344</td><td>2393</td><td>2444</td><td>2470</td><td>2491</td><td>2505</td><td>2516</td></tr><tr><td>RMS (kg m<sup>-3</sup>)</td><td></td><td>556</td><td>400</td><td>390</td><td><strong>342</strong></td><td>363</td><td>372</td><td>376</td></tr></tbody></table>

<table><caption>**Table 6.** The comparison of the mass density, estimated by the F-H method, with the average of the sample values at the 13 sampling stations, when $δ = 6$ km (unit: kg m⁻³).</caption><thead><tr><th>Rock-sampling</th><th>Station</th><th>S1</th><th>S2</th><th>S3</th><th>S4</th><th>S5</th><th>S6</th><th>S7</th><th>S8</th><th>S9</th><th>S10</th><th>S11</th><th>S12</th><th>S13</th></tr><tr><td></td><td>ρ (average)</td><td>2076</td><td>2392</td><td>2513</td><td>2418</td><td>2438</td><td>1975</td><td>1905</td><td>2030</td><td>2120</td><td>1940</td><td>1935</td><td>2364</td><td>2291</td></tr></thead><tbody><tr><td><strong>F-H (Parasnis) method</strong></td><td></td><td>2314</td><td>2315</td><td>2314</td><td>2317</td><td>2316</td><td>2400</td><td>2401</td><td>2314</td><td>2268</td><td>2589</td><td>2566</td><td>2470</td><td>2470</td></tr><tr><td><em>Error</em> (kg m<sup>-3</sup>)</td><td></td><td><strong>238</strong></td><td>-77</td><td>-199</td><td>-101</td><td>-122</td><td><strong>425</strong></td><td><strong>496</strong></td><td><strong>284</strong></td><td><strong>148</strong></td><td><strong>649</strong></td><td><strong>631</strong></td><td><strong>106</strong></td><td><strong>179</strong></td></tr></tbody></table>

By comparing Table 3 with Table 5, as well as comparing Table 4 with Table 6, we can easily understand that the results of formula $ρ \approx (1/−kπG)trace(ΔF_gΔR^+)$ which is based on the Poisson PDE at the Earth’s surface, is sufficiently more accurate than those of the F-H method. The best values of the RMS of the errors, corresponding to the proposed method of the Poisson PDE and the F-H method, are 157 and 342 kg m⁻³, respectively, which certainly shows that the results of the proposed method of the current research is at least two times more accurate than the ones of the F-H method.

## 4.3 Comparison of the results with results provided by the Murata method

Murata (1993) proposed a Bayesian method of the mass density reduction based on an assumption that the Bouguer anomaly varies smoothly as compared to topography. The Murata method finds an optimum value of the mass density reduction due to fitting a two-dimensional differentiable function to the observed Bouguer anomaly. This method is able to accomplish an optimal trade-off between smoothness and goodness of fit to the Bouguer anomalies, while it performs the trade-off via an objective Bayesian approach minimizing ABIC (Murata 1993). Like the F-H method performed in the previous subsection, herein we have applied the Murata method to the gravity stations for estimating the mass density corresponding to each sampling station. For this purpose, the neighbour gravity stations were found around each sampling station by various values of the neighbourhood dimension $δ$, and then the Murata method was applied to the gravity data of the neighbour stations. It should be mentioned that the neighbourhood, which is used in the current research for finding the neighbours, is similar to the spatial window, used by Murata (1993), which specifies the analysed area. Although Murata (1993) supposed that the size of the spatial window (i.e. the neighbourhood), required for finding the

gravity stations, was fixed to 20 km × 20 km throughout, herein we examined various values of the window size. The Murata method is performed in the spatial window (i.e. the neighbourhood) via minimizing the following objective function (Murata 1993; Nawa et al. 1997):

$$ \sum_{n=1}^{N} \|F_n - \rho H_n - f(x_n, y_n|\mathbf{s})\|^2 + \sum_{j=1}^{2} \omega_j \iint \|\nabla^{[j]} f\|^2 dx dy \quad (30) $$

where $F_n$ and $H_n$ are the free-air anomaly and the topographic correction coefficient at the neighbour stations $\{(x_n, y_n)\}(n = 1, 2, \dots, N)$, respectively; $f$ is the cubic B-spline function fitted to the observed Bouguer gravity anomalies $F_n - \rho H_n$; similar to the F-H method done in the previous subsection, we applied the SRTM’s 90 m × 90 m − DEM provided by Reuter et al. (2007) for calculating the apparent topographic corrections as well as the topographic correction coefficient $H_n$; the vector $\mathbf{s} \in \mathbb{R}^M$ includes the spline parameters; and $M$ is the number of the spline parameters, that is, the number of spline knots. It was proposed by Murata (1993) to suppose $M = 20 \times 20$ when the window size is 20 km × 20 km.

There is a trade-off in the above optimization problem between the roughness of the function $f$ fitted to the Bouguer anomaly data and the residual of the Bouguer anomaly from the function $f$. The above minimization can also be considered as a type of Generalized Tikhonov regularization method which regularizes the function $f$ and controls its roughness via a type of Sobolev semi-norm (i.e.

$$ \sum_{j=1}^{2} \omega_j \iint \|\nabla^{[j]} f\|^2 dx dy), \text{ for example, see (Hansen 1998). The two parameters } \omega_1 \text{ and } \omega_2 \text{ are also considered as the trade-off parameters of the optimization, which control both the first-order roughness (i.e. the gradient) and the second-order roughness (i.e. the curvature) of the function } f \text{ (Nawa et al. 1997). Although Murata (1993) preferred to use the cubic B-spline function for modelling the Bouguer}$$
gravity anomaly, several other two-dimensional differentiable functions can be used for this purpose such as polynomials. Murata (1993) believed that the spline function is able to represent the Bouguer gravity anomaly over a wide area with a relatively efficient calculation by comparison to the other basis functions.

For the sake of the simplicity, the Murata optimization, introduced by eq. (30), is converted into the following linear minimization (Murata 1993; Nawa et al. 1997):

$$ \|\mathbf{F} - \rho\mathbf{H} - \mathbf{E}\mathbf{s}\|^2 + \|\mathbf{D}\mathbf{s}\|^2 \to \text{Min} \qquad (31) $$

where $\mathbf{E} \in \mathbb{R}^{N \times M}$ such that $\mathbf{E}\mathbf{s}$ is satisfied for the spline parameters at the $N$ gravity stations. Also, $\mathbf{D} \in \mathbb{R}^{M \times M}$, which includes the trade-off parameters $\omega_1$ and $\omega_2$, is obtained by

$$ \|\mathbf{D}\mathbf{s}\|^2 = \sum_{j=1}^{2} \omega_j \iint \|\nabla^{[j]} f(x_n, y_n|\mathbf{s})\|^2 dx dy \qquad (32) $$

The matrices $\mathbf{E} \in \mathbb{R}^{N \times M}$, $\mathbf{D} \in \mathbb{R}^{M \times M}$, and $\mathbf{s} \in \mathbb{R}^M$ are given explicitly in Murata (1993; Appendix A). Finally, the Murata optimization (i.e. eq. 31) is reduced to solving the following equation (Murata 1993).

$$ \mathbf{a} = \mathbf{Z}\mathbf{v} \qquad (33) $$

where $\mathbf{a} = [F_1 \ F_2 \ \dots \ F_N]^T$ is the observation vector, $\mathbf{v} = [s_1 \ s_2 \ \dots \ s_M \ \rho]^T$ is the unknown vector, and $\mathbf{Z} \in \mathbb{R}^{(N+M) \times (M+1)}$ is defined as follows:

$$ \mathbf{Z} = \begin{bmatrix} E_{11} & \cdots & E_{1M} & H_1 \\ \vdots & \ddots & \vdots & \vdots \\ E_{N1} & \cdots & E_{NM} & H_N \\ D_{11} & \cdots & D_{1M} & 0 \\ \vdots & \ddots & \vdots & \vdots \\ D_{M1} & \cdots & D_{MM} & 0 \end{bmatrix} \qquad (34) $$

Besides, the last step of the Murata method is to determine optimal values of the trade-off parameters. Murata (1993) makes his

method a suitable choice for both parameters $\omega_1$ and $\omega_2$ by means of minimizing ABIC as follows:

$$ \text{ABIC} = N \log(2\pi\hat{\sigma}^2) - \log\phi + \log\lambda + N + 6 \qquad (35) $$

where $\phi$ is the determinant of the cofactor matrix of $\mathbf{D}^T\mathbf{D}$ with respect to the last diagonal element (see Murata 1993; Appendix A), $\lambda$ is the determinant of $\mathbf{Z}^T\mathbf{Z}$, and $\hat{\sigma}^2$ is defined by

$$ \hat{\sigma}^2 = \frac{1}{N} \|\mathbf{a} - \mathbf{Z}\mathbf{v}\|^2. \qquad (36) $$

Accordingly, we were able to perform the Murata method and estimate the mass density at the 13 rock-sampling stations. Table 7 shows the mass density estimated via the Murata method at the sampling stations for various values of the window size. This table indicates that the optimum value of the RMS of errors, associated with the estimated mass density of the sampling stations, was accomplished when the window size and the number of the knots were equal to 12 km × 12 km and $M = 12 \times 12$, respectively.

Besides, Table 8 illustrates the results obtained by the Murata method at the 13 sampling stations when the window size is 12 km × 12 km. Table 8 includes the number of the neighbour gravity stations, the minimum values of the ABIC criterion, the optimum values of the trade-off parameters obtained by the ABIC minimization, the mass density estimated by the Murata method, and the errors associated with the estimated mass densities, at each rock-sampling station. In this table, the errors have been calculated by subtracting the average of the geological sample values, mentioned in Table 2, from the estimated mass density at each sampling station.

By comparing Table 3 with Table 7, as well as comparing Table 4 with Table 8, we can easily realize that the results of the proposed method of the current research and those of the Murata method are mutually comparable. However, it is easily seen that the best values of the RMS of the errors, corresponding to the proposed method and the Murata method, are 157 and 185 kg m⁻³, respectively, which indicates that the proposed method was able to surpass the Murata method with respect to the RMS of the errors.

<table><caption>**Table 7.** The mass density, estimated by the Murata method at the 13 sampling-stations for various values of the window size and the number of spline knots; Also, this table includes the RMS of the error obtained by the comparison between the mass density estimated by the Murata method and the average values of the geological samples mentioned in Table 2.</caption><thead><tr><th></th><th colspan="9">Mass density estimated via the Murata method (kg m<sup>-3</sup>)</th></tr><tr><th></th><th>4×4</th><th>6×6</th><th>8×8</th><th>10×10</th><th>12×12</th><th>14×14</th><th>16×16</th><th>18×18</th><th>20×20</th></tr></thead><tbody><tr><th colspan="10">Window size (km×km)</th></tr><tr><th colspan="10">Number of the spline knots (M)</th></tr><tr><th>S1</th><td>2431</td><td>2350</td><td>2553</td><td>2336</td><td><strong>2095</strong></td><td>2206</td><td>2209</td><td>2162</td><td>2222</td></tr><tr><th>S2</th><td>2673</td><td>3039</td><td>2414</td><td>2572</td><td><strong>2565</strong></td><td>2095</td><td>2165</td><td>2176</td><td>2214</td></tr><tr><th>S3</th><td>2327</td><td>2431</td><td>2225</td><td>2533</td><td><strong>2336</strong></td><td>2160</td><td>2185</td><td>2166</td><td>2206</td></tr><tr><th>S4</th><td>2674</td><td>1915</td><td>2100</td><td>2212</td><td><strong>2581</strong></td><td>2620</td><td>2599</td><td>2243</td><td>2216</td></tr><tr><th>S5</th><td>1322</td><td>1905</td><td>2050</td><td>2258</td><td><strong>2613</strong></td><td>2597</td><td>2603</td><td>2244</td><td>2209</td></tr><tr><th>S6</th><td>2302</td><td>2317</td><td>2321</td><td>2310</td><td><strong>2307</strong></td><td>2325</td><td>2328</td><td>2309</td><td>2324</td></tr><tr><th>S7</th><td>2301</td><td>2317</td><td>2321</td><td>2310</td><td><strong>2307</strong></td><td>2325</td><td>2328</td><td>2309</td><td>2324</td></tr><tr><th>S8</th><td>1934</td><td>1957</td><td>2075</td><td>2251</td><td><strong>2248</strong></td><td>2331</td><td>2393</td><td>2350</td><td>2420</td></tr><tr><th>S9</th><td>1969</td><td>1938</td><td>1981</td><td>2206</td><td><strong>2185</strong></td><td>2259</td><td>2315</td><td>2397</td><td>2354</td></tr><tr><th>S10</th><td>1425</td><td>1498</td><td>1748</td><td>1771</td><td><strong>1903</strong></td><td>2061</td><td>2074</td><td>2123</td><td>2107</td></tr><tr><th>S11</th><td>1706</td><td>1689</td><td>1849</td><td>1995</td><td><strong>1963</strong></td><td>2122</td><td>2135</td><td>2131</td><td>2170</td></tr><tr><th>S12</th><td>2092</td><td>2150</td><td>2182</td><td>2266</td><td><strong>2301</strong></td><td>2321</td><td>2372</td><td>2407</td><td>2439</td></tr><tr><th>S13</th><td>2092</td><td>2161</td><td>2195</td><td>2266</td><td><strong>2301</strong></td><td>2322</td><td>2372</td><td>2407</td><td>2439</td></tr><tr><th>RMS (kg m<sup>-3</sup>)</th><td>419</td><td>360</td><td>272</td><td>206</td><td><strong>185</strong></td><td>241</td><td>244</td><td>245</td><td>256</td></tr></tbody></table>
<table><caption>**Table 8.** The number of the neighbour gravity stations (N), the minimum values of the ABIC criterion, the optimum values of the trade-off parameters $\omega_1$ and $\omega_2$ obtained by the ABIC minimization, the mass density $\rho$ estimated by the Murata method, and the errors associated with the estimated mass densities, at the 13 rock-sampling stations, when the window size is 12 km x 12 km.</caption><thead><tr><th>Rock-sampling stations</th><th>S1</th><th>S2</th><th>S3</th><th>S4</th><th>S5</th><th>S6</th><th>S7</th><th>S8</th><th>S9</th><th>S10</th><th>S11</th><th>S12</th><th>S13</th></tr></thead><tbody><tr><td><i>N</i></td><td>60</td><td>51</td><td>56</td><td>38</td><td>39</td><td>125</td><td>124</td><td>116</td><td>114</td><td>124</td><td>121</td><td>127</td><td>127</td></tr><tr><td>Min(ABIC) x 10<sup>-3</sup></td><td>1.12</td><td>0.98</td><td>1.07</td><td>0.73</td><td>0.75</td><td>2.40</td><td>2.38</td><td>2.04</td><td>2.01</td><td>2.14</td><td>2.08</td><td>2.19</td><td>2.19</td></tr><tr><td>ω<sub>1</sub> x 10<sup>-2</sup></td><td>0.02</td><td>0.01</td><td>0.03</td><td>0.06</td><td>0.02</td><td>0.12</td><td>0.10</td><td>0.05</td><td>0.10</td><td>0.08</td><td>0.08</td><td>0.07</td><td>0.07</td></tr><tr><td>ω<sub>2</sub> x 10<sup>-2</sup></td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.60</td><td>0.85</td><td>0.02</td><td>0.36</td><td>0.17</td><td>0.01</td><td>0.03</td><td>0.06</td></tr><tr><td>ρ (kg m<sup>-3</sup>)</td><td><b>2095</b></td><td><b>2565</b></td><td><b>2336</b></td><td><b>2581</b></td><td><b>2613</b></td><td><b>2307</b></td><td><b>2307</b></td><td><b>2248</b></td><td><b>2185</b></td><td><b>1903</b></td><td><b>1963</b></td><td><b>2301</b></td><td><b>2301</b></td></tr><tr><td>Error (kg m<sup>-3</sup>)</td><td><b>19</b></td><td><b>173</b></td><td>-177</td><td><b>163</b></td><td><b>175</b></td><td><b>332</b></td><td><b>402</b></td><td><b>218</b></td><td><b>65</b></td><td>-37</td><td><b>28</b></td><td>-63</td><td><b>10</b></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The geological map corresponding to the region of the case study as well as the spatial distribution of the gravity stations over the geological formations: in this figure, the blue dots signify the locations of the gravity stations, and the geodetic datum of the map is WGS-84.</figcaption></figure>

## 4.4 Validation of the results via a geological map

In order to have another validation, we compared the results with a geological map of the region of the case study. Fig. 9 illustrates the geological map corresponding to the area of the case study, which was produced by NIOC. In Fig. 9, the situations of the gravity stations over the geological formations can be clearly comprehended, where the areas of the geological formations were separated from each other by various colours. From this geological map, we could relate each gravity station to a geological formation at which the gravity station was located.

Herein, we had 201 sample values of the mass density belonging to the geological formations. 86 of these 201 sample values were used by the previous subsections for the validation of the results. The numbers of the samples, their lithology, their average values, and their standard deviations are included in Table 9 for all the geological formations. Also, Fig. 10(a) graphically gives an overview of these sample values corresponding to the geological formations.

Also, Fig. 10(b) graphically shows both the average value and the standard deviation associated with the sample values of each geological formation.

<table><caption>**Table 9.** The number of the samples, the lithology, the average values of the mass density obtained by the sample values, and the standard deviation of the sample values, corresponding to the geological formations out cropped in the case study region.</caption><thead><tr><th>Formation name<br>Formation<br>abbrev.</th><th>AghaJari<br>Aj</th><th>Bangestan<br>Bgn</th><th>Bakhtiari<br>Bk</th><th>Jahrom<br>Ja</th><th>Pabdeh<br>Pd</th><th>Hormuz<br>Sp</th><th>Asmari-<br>Jahrom<br>As-Ja</th><th>Gachsaran<br>Gs</th><th>Khami<br>Kgp</th><th>Mishan<br>Mn</th><th>Pabdeh-<br>Gurpi<br>Pd-Gu</th><th>Tarbur<br>Tb</th></tr></thead><tbody><tr><th>Number of samples</th><td>9</td><td>19</td><td>8</td><td>4</td><td>4</td><td>49</td><td>25</td><td>31</td><td>20</td><td>15</td><td>12</td><td>5</td></tr><tr><th>Lithology</th><td>sandstone</td><td>limestone</td><td>sandstone<br>conglomerate</td><td>dolomite</td><td>limestone</td><td>salt</td><td>limestone</td><td>gips<br>marl</td><td>limestone</td><td>limestone</td><td>marl</td><td>limestone</td></tr><tr><th>Average value<br>of the mass<br>density<br>[kg m<sup>-3</sup>]</th><td>2283</td><td>2488</td><td>2074</td><td>2593</td><td>2413</td><td>2098</td><td>2525</td><td>2216</td><td>2547</td><td>2313</td><td>2365</td><td>2586</td></tr><tr><th>Standard<br>deviation<br>[kg m<sup>-3</sup>]</th><td>222</td><td>139</td><td>265</td><td>174</td><td>139</td><td>241</td><td>150</td><td>100</td><td>19</td><td>214</td><td>142</td><td>36</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) A graph of the sample values of the geological formations; (b) a graph of the average values of the mass density corresponding to the geological formations, obtained by the sample values, as well as an overview of their associated standard deviations.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** The spatial distribution of the mass density ($\rho_{geo}$) obtained via both the geological map and the rock-samplings.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** The spatial distribution of the standard deviation ($\sigma_{geo}$) associated with the geological mass density illustrated in Fig. 10.</figcaption></figure>

Accordingly, by means of the geological map we could assign to each gravity station a geological formation at which the gravity station was placed. Then, by means of the above-mentioned average value of the mass density, obtained by the sample values associated with each geological formation, the geological map could give us a prediction about the mass density at 3879 of the gravity stations. Fig. 11 illustrates the spatial distribution of the mass density extracted from the geological map. Similarly, by means of the geological map we could assign the standard deviations, mentioned in Table 9, to the 3879 gravity stations. Fig. 12 shows the spatial distribution of the standard deviation of the extracted mass density over the region of the case study. This figure gives us an insight into the uncertainty associated with the assigned mass density extracted from the geological map.

We used the mass density, extracted from the geological map, for testing the mass density estimated via the proposed approach. In the following, the mass density extracted from the geological map is called 'geological mass density', and also denoted by $\rho_{geo}$. Furthermore, the mass density estimated via the proposed method (i.e. via eq. 23) is called 'estimated mass density', and denoted by $\rho$.

In order to validate the main idea of the paper, the difference between the estimated mass density ($\rho$) and the geological mass density ($\rho_{geo}$) was prepared and investigated. In Fig. 13, the histograms of $\rho - \rho_{geo}$ corresponding to various values of the neighbourhood dimension (i.e. $\delta$) is shown. The histograms clearly show that the average value of $\rho - \rho_{geo}$ occurred around zero. Besides, Fig. 14 gives an overview of the spatial distribution of the absolute value of the difference (i.e. $|\rho - \rho_{geo}|$) for various values of the neighbourhood dimension.

It was required to statistically test the difference between the estimated mass density and the geological mass density at each gravity station to find out about the acceptability of the difference. For this purpose, the following hypothesis test was a necessary requirement at each station.

$$ \begin{cases} H_0: E(\rho - \rho_{geo}) = 0 \\ H_a: E(\rho - \rho_{geo}) \neq 0 \end{cases} \qquad (37) $$

where $E$ is a mathematical expectation operator. The above test could be changed into the following hypothesis test at each station

$$ \begin{cases} H_0: E(\rho_{geo}) = E(\rho) \\ H_a: E(\rho_{geo}) \neq E(\rho) \end{cases} \qquad (38) $$

The above test is a hypothesis test for statistically examining the expected value of the estimated mass density (i.e. for examining a value of $E(\rho)$).

Here, it was considered that $\rho_{geo}$ is distributed normally with expectation $E(\rho_{geo})$ and variance $\sigma_{geo}^2$. The standard deviation $\sigma_{geo}$ was previously obtained at each gravity station by means of the geological map as well as the standard deviations of the sample values. The spatial distributions of $\rho_{geo}$ and $\sigma_{geo}$ over the area of the case study were shown in Figs 11 and 12, respectively.

Accordingly, the above hypothesis test yields the following confidence interval of $E(\rho)$ with confidence level $1 - \alpha$

$$ |\rho_{geo} - E(\rho)| < z_{1-\frac{\alpha}{2}} \sigma_{geo}. \qquad (39) $$

In the above confidence interval, $z$ has a standard normal distribution; $\alpha$ is a significance level, and herein, we have conventionally considered that $\alpha = 5\%$ (i.e. $z_{1-\frac{\alpha}{2}} = 1.96$).

Herein, we assumed that also the estimated mass density ($\rho$) satisfies the above confidence interval, the same as its expected value $E(\rho)$. Hence, we supposed that if the estimated mass density (i.e. $\rho$) occurred at the above confidence interval, it would be accepted according to the above hypothesis test. In Fig. 15, the percentage of the accepted mass density was plotted against the neighbourhood
<figure><img src="image_14.png" /><figcaption>Figure 13. The histograms of the difference between the estimated mass density and the geological mass density for various values of the neighbourhood dimension (i.e. δ).</figcaption></figure>

dimension. In this figure, we clearly see that the maximum of the percentage of the accepted mass densities has taken place at δ = 5 km. Fig. 15 indicates that the estimated mass density, obtained by δ = 5 km, has found an acceptable value in more than 75 per cent of the gravity stations. Besides, Fig. 16 shows the RMS of ρ – ρ<sub>geo</sub> with respect to the neighbourhood dimension, which obviously indicates that the corresponding minimum also occurred at δ = 5 km.

Also, in Fig. 17, the above-mentioned confidence interval corresponding to the estimated mass density was plotted vertically against the station number for all the gravity stations when δ = 5 km. This figure notably shows how the estimated mass density frequently occurred among the confidence interval and was accepted at many gravity stations.

## 4.5 The propagation of the errors into the estimated mass density

The mass density estimated by the proposed approach is polluted by the truncation error as well as the propagated errors. The truncation error was discussed by Section 3, and at the beginning of this section, the truncation error could be predicted by both eq. (24) and EIGEN-6C4. However, in this subsection, we are aiming at assessing the effects of the propagated errors, associated with the land-based gravity data as well as the deflection of the vertical, on the mass density evaluated by the proposed approach. For this

purpose, defining the following vectors is a requirement of each computation point:

$$ \Theta = [\xi \ \eta \ \xi_1 \ \eta_1 \ \xi_2 \ \eta_2 \ \dots \ \xi_N \ \eta_N]^T \in \mathbb{R}^{2N+2} \quad (40) $$

$$ G = [g \ g_1 \ g_2 \ \dots \ g_N]^T \in \mathbb{R}^{N+1} $$

where the vectors Θ and G are comprised of the components of the deflection of the vertical and the scalar land-based gravity observations, respectively (ξ, η and g are the meridian component of the deflection of the vertical, the prime vertical component of the deflection of the vertical and the gravity observation corresponding to the computation point, respectively; similarly, ξn, ηn and gn belong to the nth neighbour station of the computation point where n = 1, 2, ..., N). Accordingly, the general law of propagation of variances implies that the variance of the estimated mass density (σρ<sup>2</sup>) is obtained by

$$ \sigma_{\rho}^{2} = \frac{\partial \rho}{\partial \Theta} C_{\Theta} \frac{\partial \rho}{\partial \Theta}^{T} + \frac{\partial \rho}{\partial G} C_{G} \frac{\partial \rho}{\partial G}^{T}. \quad (41) $$

The matrices C<sub>Θ</sub> and C<sub>G</sub> are the covariance matrices corresponding to Θ and G, respectively. Herein, the vectors Θ and G are supposed to be statistically uncorrelated with each other and their covariance matrices C<sub>Θ</sub> and C<sub>G</sub> are diagonal matrices as

$$ \begin{aligned} C_{\Theta} &= \operatorname{diag}(\sigma_{\xi}^{2}, \sigma_{\eta}^{2}, \sigma_{\xi_{1}}^{2}, \sigma_{\eta_{1}}^{2}, \dots, \sigma_{\xi_{N}}^{2}, \sigma_{\eta_{N}}^{2}) \\ C_{G} &= \operatorname{diag}(\sigma_{g}^{2}, \sigma_{g_{1}}^{2}, \dots, \sigma_{g_{N-1}}^{2}, \sigma_{g_{N}}^{2}) \end{aligned} \quad (42) $$
<figure><img src="image_15.png" /><figcaption>**Figure 14.** The spatial distribution of the absolute value of the difference between the estimated mass density and the geological mass density (i.e. |ρ - ρ<sub>geo</sub>|) for various values of the neighbourhood dimension, over the region of the case study.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** The percentage of the accepted mass densities against the neighbourhood dimension.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** The RMS of ρ - ρ<sub>geo</sub> against the neighbourhood dimension δ.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** The graph of the confidence interval vertically plotted against station number: in this figure, the black lines indicate the lower as well as the upper bounds of the confidence intervals, and also, the dots signify the estimated mass densities; it is clear that the estimated values of the mass density frequently occurred within the confidence intervals.</figcaption></figure>

where $\sigma_{\xi}^2$, $\sigma_{\eta}^2$ and $\sigma_g^2$ denote the variances of the components of the deflection of the vertical and the gravity observation at the computation point, respectively. Similarly, $\sigma_{\xi_n}^2$, $\sigma_{\eta_n}^2$ and $\sigma_{g_n}^2$ represent the variances of the components of the deflection of the vertical and the gravity observation at the nth neighbour station, respectively ($n = 1, 2, \dots, N$). Furthermore, we have

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial \Theta} &= \left[ \frac{\partial \rho}{\partial \xi} \quad \frac{\partial \rho}{\partial \eta} \quad \frac{\partial \rho}{\partial \xi_1} \quad \frac{\partial \rho}{\partial \eta_1} \quad \frac{\partial \rho}{\partial \xi_2} \quad \frac{\partial \rho}{\partial \eta_2} \quad \dots \quad \frac{\partial \rho}{\partial \xi_N} \quad \frac{\partial \rho}{\partial \eta_N} \right] \\
\frac{\partial \rho}{\partial \mathbf{G}} &= \left[ \frac{\partial \rho}{\partial g} \quad \frac{\partial \rho}{\partial g_1} \quad \frac{\partial \rho}{\partial g_2} \quad \dots \quad \frac{\partial \rho}{\partial g_N} \right]
\end{aligned}
\qquad (43) $$

According to eq. (23), we can straightforwardly obtain that

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial \xi} &\approx \frac{1}{-k\pi G} \operatorname{trace} \left( \frac{\partial \Delta \mathbf{F}_g}{\partial \xi} \Delta \mathbf{R}^+ \right), \\
\frac{\partial \rho}{\partial \eta} &\approx \frac{1}{-k\pi G} \operatorname{trace} \left( \frac{\partial \Delta \mathbf{F}_g}{\partial \eta} \Delta \mathbf{R}^+ \right), \\
\frac{\partial \rho}{\partial \xi_n} &\approx \frac{1}{-k\pi G} \operatorname{trace} \left( \frac{\partial \Delta \mathbf{F}_g}{\partial \xi_n} \Delta \mathbf{R}^+ \right), \\
\frac{\partial \rho}{\partial \eta_n} &\approx \frac{1}{-k\pi G} \operatorname{trace} \left( \frac{\partial \Delta \mathbf{F}_g}{\partial \eta_n} \Delta \mathbf{R}^+ \right), \\
\frac{\partial \rho}{\partial g} &\approx \frac{1}{-k\pi G} \operatorname{trace} \left( \frac{\partial \Delta \mathbf{F}_g}{\partial g} \Delta \mathbf{R}^+ \right), \\
\frac{\partial \rho}{\partial g_n} &\approx \frac{1}{-k\pi G} \operatorname{trace} \left( \frac{\partial \Delta \mathbf{F}_g}{\partial g_n} \Delta \mathbf{R}^+ \right), \\
n &= 1, 2, \dots, N
\end{aligned}
\qquad (44) $$

where eqs (19), (25) and (27) can provide the following formulae for assessing the derivatives of $\Delta \mathbf{F}_g$ with respect to $\xi$, $\eta$ and $g$:

$$ \frac{\partial \Delta \mathbf{F}_g}{\partial \xi} = \operatorname{diag} \left( g \begin{bmatrix} \sin \Phi \cos \Lambda \\ \sin \Phi \sin \Lambda \\ -\cos \Phi \end{bmatrix} \right) \mathbf{1}_{3 \times N}, $$

$$ 
\begin{aligned}
\frac{\partial \Delta \mathbf{F}_g}{\partial \eta} &= \operatorname{diag} \left( \frac{g}{\cos \varphi} \begin{bmatrix} \cos \Phi \sin \Lambda \\ -\cos \Phi \cos \Lambda \\ 0 \end{bmatrix} \right) \mathbf{1}_{3 \times N} \\
\frac{\partial \Delta \mathbf{F}_g}{\partial g} &= \operatorname{diag} \left( \begin{bmatrix} \cos \Phi \cos \Lambda \\ \cos \Phi \sin \Lambda \\ \sin \Phi \end{bmatrix} \right) \mathbf{1}_{3 \times N}
\end{aligned}
\qquad (45) $$

where $\mathbf{1}_{3 \times N}$ is a $3 \times N$ matrix of ones. Similarly, the derivatives of $\Delta \mathbf{F}_g$ with respect to $\xi_n$, $\eta_n$ and $g_n$ ($n = 1, 2, \dots, N$) are obtained by the following equations:

$$ 
\begin{aligned}
\left[ \frac{\partial \Delta \mathbf{F}_g}{\partial \xi_n} \right]_{i,j \neq n} &= \left[ \frac{\partial \Delta \mathbf{F}_g}{\partial \eta_n} \right]_{i,j \neq n} = \left[ \frac{\partial \Delta \mathbf{F}_g}{\partial g_n} \right]_{i,j \neq n} = 0 \\
\left[ \frac{\partial \Delta \mathbf{F}_g}{\partial \xi_n} \right]_{i=1:3, j=n} &= g_n \begin{bmatrix} -\sin \Phi_n \cos \Lambda_n \\ -\sin \Phi_n \sin \Lambda_n \\ \cos \Phi_n \end{bmatrix}, \\
\left[ \frac{\partial \Delta \mathbf{F}_g}{\partial \eta_n} \right]_{i=1:3, j=n} &= \frac{g_n}{\cos \varphi_n} \begin{bmatrix} -\cos \Phi_n \sin \Lambda_n \\ \cos \Phi_n \cos \Lambda_n \\ 0 \end{bmatrix} \\
\left[ \frac{\partial \Delta \mathbf{F}_g}{\partial g_n} \right]_{i=1:3, j=n} &= \begin{bmatrix} \cos \Phi_n \cos \Lambda_n \\ \cos \Phi_n \sin \Lambda_n \\ \sin \Phi_n \end{bmatrix} \\
&i = 1, 2, 3 \quad j = 1, 2, \dots, N \\
n &= 1, 2, \dots, N \quad \varphi \neq \pm \frac{\pi}{2}
\end{aligned}
\qquad (46) $$

It is clear that the above propagation of variances (i.e. eq. 41) does not take into account errors associated with the positioning and coordinates of the gravity stations. It is because of the fact that the errors associated with the coordinates of the gravity stations can be counted as an error source of the gravity observation. In a general form, random errors corresponding to the coordinates such as the Cartesian coordinates $\mathbf{r} = (x, y, z) \in \mathbb{R}^3$ can be taken into consideration within the variance of the gravity observation (i.e. $\sigma_g^2$) by means of the following propagation of the variances:

$$ \sigma_g^2 = \left[ \frac{\partial g}{\partial x} \quad \frac{\partial g}{\partial y} \quad \frac{\partial g}{\partial z} \right] \mathbf{C}_r \left[ \frac{\partial g}{\partial x} \quad \frac{\partial g}{\partial y} \quad \frac{\partial g}{\partial z} \right], \qquad (47) $$

where the matrix $\mathbf{C}_r$ is the covariance matrix associated with the position vector $\mathbf{r} = (x, y, z) \in \mathbb{R}^3$. It is worth mentioning that the variance $\sigma_g^2$ depends mainly on the errors of the vertical coordinate in comparison with the ones of the horizontal components. If the vertical gradient of gravity at the Earth's surface is considered to be equal to the Bouguer gradient (i.e. $\frac{\partial g}{\partial z} \approx 0.1967 \frac{\mathrm{mGal}}{\mathrm{m}}$ (Vaniček & Krakiwsky 1986, p. 502)), a $\pm 1$ cm vertical position error will cause an approximate $\pm 2$ µGal error in the gravity observation.

The final precision of the land-based gravity data depends mostly on the errors associated with the positioning, instrumentation of the acquisition and gravity base network used by the gravimetry project. According to the internal report (NIOC 2004) prepared by the national Iranian oil company as the organization which acquired the land-based gravity data, the precision of the gravity observation at
<figure><img src="image_19.png" /><figcaption>**Figure 18.** The graph of the average of the standard deviation of the mass density against the neighbourhood dimension for the two cases: the first case supposes that $\sigma_{\epsilon} = 1.2''$, $\sigma_{\eta} = 2.7''$ and $\sigma_{g} = 0$, and the second one supposes that $\sigma_{g} = 70 \mu\text{Gal}$ and $\sigma_{\epsilon} = \sigma_{\eta} = 0$.</figcaption></figure>

each station is claimed to be approximately around $\pm 70 \mu\text{Gal}$ considering the precisions of the coordinates, gravimeters and gravity base network associated with the gravimetry project.

Herein, we considered the RMS of the errors obtained by Table 1 as the estimates of the variances associated with the components of the deflection of the vertical. Thus, we could also evaluate the variance of the estimated mass density at each gravity station according to the aforementioned propagation of variances (i.e. eq. 41) as well as the known precisions of land-based gravity data and the precision of the deflections of the vertical. In brief, it was considered that $\sigma_g = 70 \mu\text{Gal}$, $\sigma_\epsilon = 1.2 \text{ arcsec}$ and $\sigma_\eta = 2.7 \text{ arcsec}$, and this consideration was used for assessing the diagonal covariance matrices $C_\odot$ and $C_G$ at each gravity station. Then, the evaluated matrices $C_\odot$ and $C_G$ were applied to eq. (41) for estimating the variance of the mass density ($\sigma_\rho^2$) at each gravity station. Fig. 18 illustrates the average of the standard deviation of the estimated mass density (i.e. $\sigma_\rho$) with respect to the neighbourhood dimension $\delta$. In this figure, the standard deviation has been shown for two cases to compare the effects of both $\sigma_\epsilon$ and $\sigma_\eta$ on the mass density with the ones of $\sigma_g$. At the first case, the estimated mass density was affected just by the noises belonging to the components of the deflection of the vertical and it supposed that $\sigma_\epsilon = 1.2 \text{ arcsec}$, $\sigma_\eta = 2.7 \text{ arcsec}$ and $\sigma_g = 0$.

<figure><img src="image_20.png" /><figcaption>**Figure 19.** (a) The graph of the average of the standard deviation of estimated mass density against the neighbourhood dimension; (b) the graph of the average of the norm of the residual against the neighbourhood dimension.</figcaption></figure>

On the other hand, the second case supposes that the estimated mass density is polluted just by the noises associated with the land-based gravity observations where $\sigma_g = 70 \mu\text{Gal}$ and $\sigma_\epsilon = \sigma_\eta = 0$. The two cases are compared with each other by Fig. 18 which indicates how the mass density is mainly influenced by the noises corresponding to the components of the deflection of the vertical in comparison with the noises of the land-based gravity observations.

## 4.6 Tuning the neighbourhood dimension

Accordingly, it is clear that the neighbourhood dimension $\delta$ plays an important role in the proposed method of estimating the mass density. According to the methodology described by Section 3, the growth of the neighbourhood dimension causes the number of neighbours $N$ to increase and so, the redundancy of the associated problem, as well as the precision of the estimated mass density, is improved by increasing the neighbourhood dimension. Fig. 19(a) illustrates the precision of the estimated mass density (i.e. $\sigma_\rho$) with respect to the neighbourhood dimension $\delta$. In this figure, first the variance $\sigma_\rho$ at all the gravity stations was assessed by eq. (41), considering $\sigma_g = 70 \mu\text{Gal}$, $\sigma_\epsilon = 1.2 \text{ arcsec}$ and $\sigma_\eta = 2.7 \text{ arcsec}$, and then the average of $\sigma_\rho$ was plotted against the neighbourhood dimension $\delta$. Fig. 19(a) indicates that the average of $\sigma_\rho$ decreases when the neighbourhood dimension $\delta$ becomes greater in amount, that is, the growth of the neighbourhood dimension $\delta$ improves the precision of the estimated mass density. On the other hand, the growth of the neighbourhood dimension makes the effect of the truncation error be greater in size. The truncation error, so-called the errors of the model, $e_{\text{MDL}}(\Delta\mathbf{R})$, defined by eq. (20), will be considerable for great values of the neighbourhood dimension. Hence, the growth of the neighbourhood dimension causes the residual $\Delta\mathbf{F}_g - \Gamma\Delta\mathbf{R}$ to grow in norm and consequently, the certainty of the associated approximation (i.e. eq. 21) used for approximating the gradient tensor as well as the mass density at the earth's surface, will be decreased by the growth of the neighbourhood dimension. In Fig. 19(b), the average of the norm of the residual (i.e. $\|\Delta\mathbf{F}_g - \Gamma\Delta\mathbf{R}\|$) was plotted against the neighbourhood dimension. This figure clearly shows that the average of the norm of the residual is strictly increased by increasing the neighbourhood dimension.

Herein, the authors strongly recommend the geological knowledge such as the mass densities extracted from the geological map to find an optimum value of the neighbourhood dimension. According to the results described by the previous Subsection, the geological
<figure><img src="image_21.png" /><figcaption>Figure 20. The L-Curve corresponding to $\sigma_\rho$ and $\|\Delta\mathbf{F}_g - \Gamma\Delta\mathbf{R}\|$.</figcaption></figure>

knowledge about the mass density indicates that the optimum value of the neighbourhood dimension should be around $\delta = 5$ km which provided the RMS of the corresponding error and the percentage of the accepted mass density with their minimum and maximum values, respectively, see Figs 15 and 16.

Nevertheless, if the geological knowledge about the mass density is not available, the above-mentioned characteristics of the neighbourhood dimension can be so useful for tuning it. As the growth of the neighbourhood dimension increases both the precision of the estimated mass density and the uncertainty of the model, it is required to find an optimum value of the neighbourhood dimension $\delta$ which can simultaneously minimize the standard deviation of the estimated mass density $\sigma_\rho$ and the norm of the residual $\|\Delta\mathbf{F}_g - \Gamma\Delta\mathbf{R}\|$. For this purpose, the authors recommend the L-Curve parameter-choice method (Lawson & Hanson 1995) for finding the optimum value which is able to minimize both $\sigma_\rho$ and $\|\Delta\mathbf{F}_g - \Gamma\Delta\mathbf{R}\|$ simultaneously. In Fig. 20, the average of $\sigma_\rho$ is plotted against the average of $\|\Delta\mathbf{F}_g - \Gamma\Delta\mathbf{R}\|$ where the graph forms the L-Curve. In this figure, it is clear that $\delta = 5$ km as the optimum value of the neighbourhood dimension occurred at the corner of the graph.

## 4.7 The contributions of the diagonal elements of the gradient tensor toward estimating the mass density

According to the proposed approach based on the Poisson PDE introduced by Theorem 2, the estimated mass density was constructed out of the diagonal elements of the gradient tensor. First the gradient tensor was estimated by eq. (22) and then it was applied to the Poisson PDE for assessing the mass density via eq. (23). Fig. 21 illustrates the spatial distribution of the diagonal elements of the gradient tensor estimated over the region of the case study (when $\delta$ is equal to it's the optimum value, i.e. $\delta = 5$ km). This figure shows the diagonal elements that were assessed in the local coordinate system at which (1) the x-axis is directed toward the east, (2) the z-axis is a tangent on the plumb line and outward from the Earth, and (3) the y-axis is chosen so that the system is right-handed.

Also, Table 10 provides the diagonal elements of the gradient tensor, as well as their effects on the estimated mass density, with the summary statistics. Moreover, Fig. 22 illustrates the spatial distribution of the estimated mass density as well as the ones of the corresponding effects of the diagonal elements. In this figure, it is plainly seen how the diagonal elements of the gradient tensor contribute towards assessing the mass density. According to the summary statistics of Table 10 as well as the corresponding spatial distributions illustrated by Fig. 22, it is straightforwardly comprehended that the mass density was mainly influenced by the second vertical derivative of the gravitational potential (i.e. $\partial^2 V / \partial^2 z$). However, it is also realized that the horizontal gradients, especially $\partial^2 V / \partial^2 y$ (i.e. the second derivative in the north direction), also had considerable effects on the evaluated mass density.

## 5 CONCLUSION AND DISCUSSION

The main purpose of this paper was to familiarize readers with a version of Poisson's PDE corresponding to the gravitational field at the Earth's surface and then to show its capability in approximating the mass density at the boundary points of the Earth's surface. According to the concepts of Newtonian potential theory, the Poisson PDE of the gravitational field could not generally exist

<figure><img src="image_22.png" /><figcaption>Figure 21. The spatial distribution of the diagonal elements of the estimated gradient tensor over the region of the case study.</figcaption></figure>

<table><caption>Table 10. The summary statistics of the diagonal elements and their effects on the evaluated mass density.</caption><thead><tr><th rowspan="3"></th><th colspan="3">The diagonal elements<br>of the gradient<br>tensor [Eötvös]</th><th colspan="3">The effects of the diagonal elements on<br>the estimated mass density<br>[kg m<sup>-3</sup>]</th><th rowspan="3">The estimated mass<br>density<br>[kg m<sup>-3</sup>]</th></tr><tr><th>∂<sup>2</sup>V</th><th>∂<sup>2</sup>ν</th><th>∂<sup>2</sup>V</th><th>1 / (-kπG) (∂<sup>2</sup>V / ∂x<sup>2</sup>)</th><th>1 / (-kπG) (∂<sup>2</sup>V / ∂y<sup>2</sup>)</th><th>1 / (-kπG) (∂<sup>2</sup>V / ∂z<sup>2</sup>)</th></tr><tr><th>∂<sup>2</sup>V / ∂x<sup>2</sup></th><th>∂<sup>2</sup>ν / ∂y<sup>2</sup></th><th>∂<sup>2</sup>V / ∂z<sup>2</sup></th><th>1 / (-kπG) (∂<sup>2</sup>V / ∂x<sup>2</sup>)</th><th>1 / (-kπG) (∂<sup>2</sup>V / ∂y<sup>2</sup>)</th><th>1 / (-kπG) (∂<sup>2</sup>V / ∂z<sup>2</sup>)</th></tr></thead><tbody><tr><td>Min</td><td>-1588</td><td>-1641</td><td>895</td><td>3335</td><td>3085</td><td>-6261</td><td>790</td></tr><tr><td>Average</td><td>-1512</td><td>-1506</td><td>2046</td><td>3606</td><td>3592</td><td>-4878</td><td>2319</td></tr><tr><td>Max</td><td>-1399</td><td>-1294</td><td>2626</td><td>3787</td><td>3912</td><td>-2135</td><td>4689</td></tr><tr><td>Range</td><td>189</td><td>347</td><td>1731</td><td>452</td><td>827</td><td>4126</td><td>3899</td></tr><tr><td>Standard deviation</td><td>19</td><td>48</td><td>136</td><td>46</td><td>114</td><td>324</td><td>331</td></tr></tbody></table>
<figure><img src="image_23.png" /><figcaption>**Figure 22.** The spatial distribution of the effects of the diagonal elements of the gradient tensor on the estimated mass density, over the region of the case study.</figcaption></figure>

at the Earth's boundary; hence, we had to mathematically introduce the Poisson PDE to the Earth's surface. For this purpose, first we had to provide certain definitions such as the definition of the gradient tensor at the Earth's surface. Then, we were able to use the definitions, as well as an important theorem of Newtonian potential theory (i.e. Theorem 1 mentioned in Section 2), for proposing Theorem 2 that introduces the Poisson PDE to the Earth's boundary points (see Section 2).

In Section 3, it was shown that the definition of the gradient tensor at the Earth's surface (i.e. Definition 4) is able to relate the land-based gravitational vectors to each other (see Lemma 1), and then the gradient tensor could be formulated by this relation. Consequently, the formulated gradient tensor was applied to the Poisson PDE to propose a formula for estimating the terrain mass density.

In Section 4, we were able to numerically examine the main idea of the research in a case study. In the case study, we had 6350 gravity stations at which the norm of the Earth's gravity field was observed. We combined the observed land-based gravity data with EIGEN-6C4 to gain the Earth's gravitational vector field at each gravity station. Then, we applied the proposed method of approximating the mass density to the gravitational vectors corresponding to the gravity stations. The methodology was able to approximate the mass density at each gravity station. We validated the approximated mass density via certain rock-samplings acquired in 13 stations. The best value of the RMS of the differences at the 13 rock-sampling stations was obtained around $157 \text{ kg m}^{-3}$ when $\delta = 4.5 \text{ km}$. Also, at the 13 rock-sampling stations, the results were compared with results obtained by a conventional method of the mass density reduction called the F-H method (i.e. the Parasnis method). The comparison numerically indicated the superiority of the proposed method over the conventional F-H method with respect to the RMS of the errors. Besides, the estimated mass densities at the 13 rock-sampling stations were compared with values provided by another conventional method of mass density reduction belonging to Murata (1993). The comparison of the results of the proposed method of the current research with the ones of the Murata method indicated that the accuracies of the two methods can be mutually comparable, although the proposed method of the current research was implemented so simply in comparison with the Murata method, and its RMS of the errors at the sampling stations was somewhat superior to that of the Murata method.

Furthermore, the numerical experiments were confirmed by means of a geological map. According to the comparison between the mass density estimated by the proposed method and the geological mass density obtained by the geological map, the estimated mass density was statistically accepted in 75.9 per cent of the gravity stations. It should be taken into account that this results have been obtained just by gravity data without using any prior knowledge about the mass density.

In conclusion, according to the proposed theories and the numerical experiments, we can briefly say that the resurrected version of Poisson's PDE at the Earth' surface has the capability of extracting the mass density from land-based gravity data.

## REFERENCES

Bjork, A., 1988. *Least Squares Methods: Handbook of Numerical Analysis*, Elsevier Science Publisher.

Farr, T.G. et al., 2007. The shuttle radar topography mission, *Rev. Geophys.*, **45**, RG2004, doi:10.1029/2005RG000183.

Förste, C. et al., 2014. EIGEN-6C4-The latest combined global gravity field model including GOCE data up to degree and order 1949 of GFZ Potsdam and GRGS Toulouse, in *EGU General Assembly Conference Abstracts*, pp. 3707, doi:10.5880/icgem.2015.1.

Fukao, Y., Yamamoto, A. & Nozaki, K., 1981. A method of density determination for gravity correction, *J. Phys. Earth*, **29**, 163-166.

Hansen, P.C., 1998. *Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion*, SIAM.

Hofmann-Wellenhof, B. & Moritz, H., 2006. *Physical Geodesy*, Springer Science & Business Media.

Kellogg, O.D., 1929. *Foundations of Potential Theory*, Read Books.

Lawson, C. & Hanson, R., 1995. *Solving Least Squares Problems*, SIAM.

Moritz, H., 1980. *Advanced Physical Geodesy*, Wichmann.

Murata, Y., 1993. Estimation of optimum average surficial density from gravity data: an objective Bayesian Approach, *J. geophys. Res.*, **98**, 12 097-12 109.

Nawa, K., Fukao, Y., Shichi, R. & Murata, Y., 1997. Inversion of gravity data to determine the terrain density distribution in southwest Japan, *J. geophys. Res.*, **102**, 27 703-27 719.

Nettleton, L., 1939. Determination of density for reduction of gravimeter observations, *Geophysics*, **4**, 176-183.

NIOC, 2004. *Gravity and Magnetic Data Acquisition Plan in Fars Area*, National Iranian Oil Company (NIOC).

Parasnis, D., 1952. A study of rock densities in the English Midlands, *Geophys. Suppl. Mon. Not. R. Astron. Soc.*, **6**, 252-271.
Parasnis, D.S., 1979. *Principles of Applied Geophysics*, Chapman and Hall.

Pavlis, N.K., Holmes, S.A., Kenyon, S.C. & Factor, J.K., 2012. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008), *J. geophys. Res.*, **117**, B04406, doi:10.1029/2011JB008916.

Reuter, H.I., Nelson, A. & Jarvis, A., 2007. An evaluation of void-filling interpolation methods for SRTM data, *Int. J. Geogr. Inf. Sci.*, **21**, 983–1008.

Vaniček, P. & Krakiwsky, J.E., 1986. *Geodesy, The Concepts*, North Holland.

## APPENDIX A: THE PROOF OF LEMMA 1

Now, remember the assumptions of Definition 3, that is, assume that $\Omega_0$ be the set of points of Earth's surface and also assume that there exists a positive value such as $\delta$ implying that the gravitational field $\mathbf{f}_g$ has continuous derivatives at each point of the neighbourhood $B_\delta(\mathbf{r})$ except the points belonging to the Earth's surface.

Hence, in the neighbourhood $B_\delta(\mathbf{r})$, the gravitational field $\mathbf{f}_g$ can be represented by a piecewise differentiable function as follows:

$$ \mathbf{f}_g(\mathbf{r}') = \begin{cases} \mathbf{f}_g^{\text{ex}}(\mathbf{r}') & \mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}} \\ \mathbf{f}_g^{\text{in}}(\mathbf{r}') & \mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_{\text{in}} \end{cases} \qquad (\text{A1}) $$

where $\mathbf{f}_g^{\text{ex}}$ and $\mathbf{f}_g^{\text{in}}$ are differentiable vector functions in $B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}$ and $B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}$, respectively. According to the continuity of the gravitational field at the boundary points of Earth's surface (see Kellogg 1929, p. 151), we have

$$ \forall \mathbf{r} \in \Omega_0, \quad \mathbf{f}_g(\mathbf{r}) = \mathbf{f}_g^{\text{ex}}(\mathbf{r}) = \mathbf{f}_g^{\text{in}}(\mathbf{r}) \qquad (\text{A2}) $$

Also, the presumed differentiability of $\mathbf{f}_g^{\text{ex}}$ in $B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}$ implies the following Equation for all $\mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}$

$$ \mathbf{f}_g^{\text{ex}}(\mathbf{r}'+\Delta\mathbf{r}) = \mathbf{f}_g^{\text{ex}}(\mathbf{r}') + \Gamma^{\text{ex}}(\mathbf{r}')\Delta\mathbf{r} + \mathbf{e}^{\text{ex}}(\Delta\mathbf{r}) \qquad (\text{A3}) $$

where we have $\frac{e^{\text{ex}}(\Delta\mathbf{r})}{\|\Delta\mathbf{r}\|} \to 0$ if $\|\Delta\mathbf{r}\| \to 0$. In a same way, the presumed differentiability of $\mathbf{f}_g^{\text{in}}$ implies the following equation for all $\mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}$

$$ \mathbf{f}_g^{\text{in}}(\mathbf{r}'+\Delta\mathbf{r}) = \mathbf{f}_g^{\text{in}}(\mathbf{r}') + \Gamma^{\text{in}}(\mathbf{r}')\Delta\mathbf{r} + \mathbf{e}^{\text{in}}(\Delta\mathbf{r}) \qquad (\text{A4}) $$

where we also have $\frac{e^{\text{in}}(\Delta\mathbf{r})}{\|\Delta\mathbf{r}\|} \to 0$ if $\|\Delta\mathbf{r}\| \to 0$.

According to the above assumptions, the Earth's gravitational gradient tensor in the neighbourhood $B_\delta(\mathbf{r})$ can be also represented by a piecewise continuous function as follows:

$$ \Gamma(\mathbf{r}') = \begin{cases} \Gamma^{\text{ex}}(\mathbf{r}') & \mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}} \\ \Gamma_0(\mathbf{r}') & \mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_0 \\ \Gamma^{\text{in}}(\mathbf{r}') & \mathbf{r}' \in B_\delta(\mathbf{r}) \cap \Omega_{\text{in}} \end{cases} \qquad (\text{A5}) $$

where $\Gamma^{\text{ex}} = \mathbf{f}_g^{\text{ex}}\nabla^T$, $\Gamma^{\text{in}} = \mathbf{f}_g^{\text{in}}\nabla^T$, and $\Gamma_0$ is defined by Definition 4. According to the discontinuity of the mass density at the Earth's surface, it is clear that $\Gamma(\mathbf{r})$ should certainly be discontinuous at the Earth's surface, see Kellogg (1929, p. 156) (i.e. we have $\Gamma^{\text{ex}}(\mathbf{r}) \neq \Gamma_0(\mathbf{r}) \neq \Gamma^{\text{in}}(\mathbf{r})$).

Now, Definition 4 implies the following relations for $\mathbf{r} \in \Omega_0$.

$$ \begin{aligned} \Gamma(\mathbf{r}) &:= \Gamma_0(\mathbf{r}) \\ &= \lim_{\delta \to 0^+} \bar{\Gamma}_\delta(\mathbf{r}) \\ &= \lim_{\delta \to 0^+} \frac{1}{v} \oint_{B_\delta(\mathbf{r}) - \Omega_0} \Gamma dv \end{aligned} $$

$$ \begin{aligned} &= \lim_{\delta \to 0^+} \frac{1}{v} \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}} \Gamma dv + \lim_{\delta \to 0^+} \frac{1}{v} \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}} \Gamma dv \\ &= \lim_{\delta \to 0^+} \left( \frac{1}{v} \bar{\Gamma}^{\text{ex}}(\mathbf{r}) \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}} dv \right) \\ &\quad + \lim_{\delta \to 0^+} \left( \frac{1}{v} \bar{\Gamma}^{\text{in}}(\mathbf{r}) \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}} dv \right) \end{aligned} \qquad (\text{A6}) $$

where $v$ is the volume of the neighbourhood $B_\delta(\mathbf{r})$. Also, $\bar{\Gamma}^{\text{ex}}(\mathbf{r})$ and $\bar{\Gamma}^{\text{in}}(\mathbf{r})$ are the mean of $\Gamma^{\text{ex}}$ and $\Gamma^{\text{in}}$ in $B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}$ and $B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}$, respectively, which are defined as follows:

$$ \begin{aligned} \bar{\Gamma}^{\text{ex}}(\mathbf{r}) &= \frac{\oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}} \Gamma dv}{\oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}} dv} \\ \bar{\Gamma}^{\text{in}}(\mathbf{r}) &= \frac{\oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}} \Gamma dv}{\oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}} dv} \end{aligned} \qquad (\text{A7}) $$

Hence, if there exists $\bar{k} = \lim_{\delta \to 0^+} \frac{1}{v} \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}} dv$ at $\mathbf{r} \in \Omega_0$, the presumed continuity of both $\Gamma^{\text{ex}}$ and $\Gamma^{\text{in}}$ will imply the following relations for all $\mathbf{r} \in \Omega_0$.

$$ \begin{aligned} \Gamma(\mathbf{r}) &:= \Gamma_0(\mathbf{r}) \\ &= \underbrace{\left( \lim_{\delta \to 0^+} \bar{\Gamma}^{\text{ex}}(\mathbf{r}) \right)}_{\Gamma^{\text{ex}}(\mathbf{r})} \underbrace{\left( \lim_{\delta \to 0^+} \frac{1}{v} \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{ex}}} dv \right)}_{1 - \bar{k}(\mathbf{r})} \\ &\quad + \underbrace{\left( \lim_{\delta \to 0^+} \bar{\Gamma}^{\text{in}}(\mathbf{r}) \right)}_{\Gamma^{\text{in}}(\mathbf{r})} \underbrace{\left( \lim_{\delta \to 0^+} \frac{1}{v} \oint_{B_\delta(\mathbf{r}) \cap \Omega_{\text{in}}} dv \right)}_{\bar{k}(\mathbf{r})} \end{aligned} \qquad (\text{A8}) $$

and

$$ \forall \mathbf{r} \in \Omega_0, \quad \Gamma(\mathbf{r}) = (1 - \bar{k}(\mathbf{r}))\Gamma^{\text{ex}}(\mathbf{r}) + \bar{k}(\mathbf{r})\Gamma^{\text{in}}(\mathbf{r}) \qquad (\text{A9}) $$

According to the definition of the topographic coefficient ($k$) in the Theorem 2, it is clear that $\bar{k} = \frac{1}{4}k$ and $0 < \bar{k} < 1$. Interestingly, eq. (A9) indicates that the gradient tensor at the Earth's surface can be considered as a weighted average of the gradient tensors corresponding to the exterior and the interior spaces of the Earth (i.e. $\Gamma^{\text{ex}}$ and $\Gamma^{\text{in}}$).

Now, if $\mathbf{r}' \to \mathbf{r} \in \Omega_0$, eq. (A2) will be able to change eqs (A3) and (A4) into the following equations.

$$ \forall \mathbf{r} \in \Omega_0, \quad \mathbf{f}_g^{\text{ex}}(\mathbf{r} + \Delta\mathbf{r}) = \mathbf{f}_g(\mathbf{r}) + \Gamma^{\text{ex}}(\mathbf{r})\Delta\mathbf{r} + \mathbf{e}^{\text{ex}}(\Delta\mathbf{r}) \qquad (\text{A10}) $$

and

$$ \forall \mathbf{r} \in \Omega_0, \quad \mathbf{f}_g^{\text{in}}(\mathbf{r} + \Delta\mathbf{r}) = \mathbf{f}_g(\mathbf{r}) + \Gamma^{\text{in}}(\mathbf{r})\Delta\mathbf{r} + \mathbf{e}^{\text{in}}(\Delta\mathbf{r}) \qquad (\text{A11}) $$

Then, for $\mathbf{r} \in \Omega_0$ the following equation will be also achieved via a combination of eqs (A9)–(A11).

$$ \begin{aligned} (1 - \bar{k}(\mathbf{r}))\mathbf{f}_g^{\text{ex}}(\mathbf{r} + \Delta\mathbf{r}) + \bar{k}(\mathbf{r})\mathbf{f}_g^{\text{in}}(\mathbf{r} + \Delta\mathbf{r}) \\ = \mathbf{f}_g(\mathbf{r}) + \Gamma(\mathbf{r})\Delta\mathbf{r} + \mathbf{e}_{\text{MDL}}(\Delta\mathbf{r}) \end{aligned} \qquad (\text{A12}) $$
where $\mathbf{e}_{\text{MDL}}(\Delta\mathbf{r})$ is the error of the model and is defined as follows:

$$ \mathbf{e}_{\text{MDL}}(\Delta\mathbf{r}) := (1 - \bar{k}(\mathbf{r}))\mathbf{e}^{\text{ex}}(\Delta\mathbf{r}) + \bar{k}(\mathbf{r})\mathbf{e}^{\text{in}}(\Delta\mathbf{r}) \quad (\text{A13}) $$

Both eqs (A3) and (A4) imply that if $\|\Delta\mathbf{r}\| \to 0$, $\frac{e_{\text{MDL}}(\Delta\mathbf{r})}{\|\Delta\mathbf{r}\|} \to 0$.

Now, considering the continuity of the gravitational vector at the Earth's surface, if $\mathbf{r} + \Delta\mathbf{r}$ also belongs to the Earth's surface (i.e. $\mathbf{r} + \Delta\mathbf{r} \in \Omega_0$), the following equation will exist similar to eq. (A2).

$$ \forall \mathbf{r} + \Delta\mathbf{r} \in \Omega_0, \quad \mathbf{f}_g(\mathbf{r} + \Delta\mathbf{r}) = \mathbf{f}_g^{\text{ex}}(\mathbf{r} + \Delta\mathbf{r}) = \mathbf{f}_g^{\text{in}}(\mathbf{r} + \Delta\mathbf{r}) \quad (\text{A14}) $$

Hence, now the proof of Lemma 1 will be finally completed via the following equations obtained by a combination of eqs (A12) and (A14).

$$ \forall \{\mathbf{r}, \mathbf{r} + \Delta\mathbf{r}\} \subset \Omega_0 $$

$$ (1 - \bar{k}(\mathbf{r})) \underbrace{\mathbf{f}_g^{\text{ex}}(\mathbf{r} + \Delta\mathbf{r})}_{\mathbf{f}_g(\mathbf{r} + \Delta\mathbf{r})} + \bar{k}(\mathbf{r}) \underbrace{\mathbf{f}_g^{\text{in}}(\mathbf{r} + \Delta\mathbf{r})}_{\mathbf{f}_g(\mathbf{r} + \Delta\mathbf{r})} = \mathbf{f}_g(\mathbf{r}) + \Gamma(\mathbf{r})\Delta\mathbf{r} + \mathbf{e}_{\text{MDL}}(\Delta\mathbf{r}) \quad (\text{A15}) $$

and

$$ \mathbf{f}_g(\mathbf{r} + \Delta\mathbf{r}) = \mathbf{f}_g(\mathbf{r}) + \Gamma(\mathbf{r})\Delta\mathbf{r} + \mathbf{e}_{\text{MDL}}(\Delta\mathbf{r}) \quad (\text{A16}) $$

## APPENDIX B: TUNING THE THRESHOLD OF THE CONDITION NUMBER

As it was mentioned by Section 4, in order to avoid the corresponding ill-conditioning of the problem, we chose the gravity stations located at areas with a rough or mountainous topography. For this purpose, we selected gravity stations at which the associated condition number of $\Delta\mathbf{R}$ was less than a threshold such as $\tau$. Hence, finding the optimum value of the threshold $\tau$, which minimized the number of omitted stations and maximized the precision of the solutions, was a requirement. For this purpose, first the mass density at all the gravity stations was estimated by the proposed formula (i.e. eq. 23) for the various values of the neighbourhood dimension $\delta = 2$ km, 3 km, ..., 20 km. Then, we evaluated the standard deviation of all the estimated mass densities corresponding to gravity stations at which the condition number was less than an arbitrary value of the threshold $\tau$. For assessing the standard deviation of the estimated mass density, we

<figure><img src="image_24.png" /><figcaption>Figure B2. The L-Curve formed by the graph of the average of $\sigma_\rho$ against the number of omitted stations.</figcaption></figure>

used the variance propagation explained by Section 4.5, considering $\sigma_g = 70$ µGal, $\sigma_\xi = 1.2$ arcsec and $\sigma_\eta = 2.7$ arcsec. Fig. B1(a) illustrates the average of the evaluated standard deviation of the mass density with respect to the threshold of the condition number $\tau$. This figure obviously shows how the standard deviation of the mass density increases by the growth of $\tau$, which means that the instability of the solutions is amplified by increasing the ill-conditioning.

On the other hand, low values of the threshold $\tau$ neglected many gravity stations. In Fig. B1(b), the number of the omitted gravity stations is plotted against the threshold $\tau$. This figure indicates that the number of the neglected stations is decreasing by the growth of the threshold $\tau$.

It was so important to find an optimum value of the threshold $\tau$ minimizing the standard deviation of the estimated mass density and the number of the neglected stations simultaneously. For this purpose, we used the L-Curve parameter-choice method (Lawson & Hanson 1995) for finding the optimum value which was able to diminish both the standard deviation and the number of the neglected stations. Fig. B2 illustrates the L-Curve formed by plotting a graph of the average of the corresponding standard deviation against the number of the neglected stations. In Fig. B2, the L-Curve indicates that the optimum value of the threshold $\tau$ is equal to 86.7 which occurred at the corresponding corner of the L-Curve graph.

<figure><img src="image_25.png" /><figcaption>Figure B1. (a) The graph of the average of the standard deviation, belonging to the estimated mass density, against the threshold of the condition number; (b) the graph of the number of omitted stations against the threshold of the condition number.</figcaption></figure>
APPENDIX C: APPROXIMATING THE
TOPOGRAPHIC COEFFICIENT *k*

According to the definition of the topographic coefficient *k* mentioned by Theorem 2, this coefficient is the limit of *k*<sub>δ</sub> as the neighbourhood dimension *δ* approaches zero. Hence, in practice the topographic coefficient *k* can be approximated by the coefficient *k*<sub>δ</sub> which is able to be assessed by a Digital Terrain Model (DTM). In the following, a method of evaluating the topographic coefficient *k*<sub>δ</sub> is originated by means of applying a DTM such as the SRTM to its corresponding definition. In brief, the coefficient *k*<sub>δ</sub> is defined by the ratio of the volume of topographic masses, enclosed in the neighbourhood, to the whole volume of the neighbourhood, see Theorem 2.

Suppose that the topography of the earth surface is represented
by a gridded model of a DTM comprised of a set of uniform gridded
prisms discretely representing the earth's shape.

Hence, first, let the neighbourhood Cube(r, δ), introduced by Section 3, be defined in a local geodetic system which (1) x-axis is directed towards the east, (2) z-axis is outwards from (also perpendicular to) the vertical datum and (3) y-axis is chosen so that the system is right-handed. Then, suppose the intersection of the neighbourhood Cube(r, δ) and the DTM as follows:

$$
\text{Cube}(\mathbf{r}, \delta) \cap \text{DTM} = \{ \underline{\mathbf{r}} \in \text{DTM} \mid \|\mathbf{r} - \mathbf{r}'\|_{\infty} < \delta \} \quad (C1)
$$

where this intersection includes prisms of the DTM surrounded by
the neighbourhood Cube(***r***,**δ**). The volume of the above intersection
represents the volume of the Earth’s topographic masses enclosed
in the neighbourhood, which is required for approximating the topo-
graphic coefficient k***<sub>δ</sub>***. Also, the estimation of this volume requires
the intersection of each prism and the neighbourhood to be assessed.
For this purpose, the height of each prism is reduced to the height
of its enclosed part in the neighbourhood as follows:

$$
\underline{h} = \begin{cases} (1) & 0 \qquad h - h_0 < -\delta \\ (2) & h - h_0 + \delta \qquad |h - h_0| \le \delta \\ (3) & 2\delta \qquad h - h_0 > +\delta \end{cases} \tag{C2}
$$

where *h* is the height of each prism belonging to Cube(**r**, **δ**) ∩ DTM,
$h_0$ is the height of the computation point, and the reduced height
$\underline{h} \in [0, 2\delta]$ is the height of the enclosed part of the prism in the
neighbourhood. The above reduction means that (1) if a point of the
DTM is located under the bottom face of Cube(**r**, **δ**) (i.e. $h - h_0 < -\delta$),
its reduce height $\underline{h}$ will be equal to zero, (2) if the point
is within Cube(**r**, **δ**), its reduce height $\underline{h}$ will be referenced to the
bottom of Cube(**r**, **δ**) (i.e. $\underline{h} = h - h_0 + \delta$), and (3) if the point
is located above the top face of Cube(**r**, **δ**) (i.e. $h - h_0 > \delta$),
its reduced height $\underline{h}$ will be equal to $2\delta$. Fig. C1 schematically illustrates
the corresponding reduced height $\underline{h}$, where the intersection of the
neighbourhood Cube(**r**, **δ**) and the DTM (i.e. Cube(**r**, **δ**) ∩ DTM) is
also distinguished from other parts of the DTM by a grey colour.
Therefore, the reduced height $\underline{h}$ enables the volume of the enclosed
part of the topography to be approximated by the following formula
(see Fig. C1)

$$
V_{\text{Mass}}(\mathbf{r}, \delta) \approx \sum_{\text{Cube} \cap \text{DTM}} \underline{h} \times \Delta x \times \Delta y, \quad (C3)
$$

where $V_{\text{Mass}}(\mathbf{r}, \delta)$ is the volume of the topographic masses en-
closed in the neighbourhood; and $\Delta x$ and $\Delta y$ are the pixel-sizes
of the DTM in the directions of the horizontal plane. If the sub-
set Cube($\mathbf{r}, \delta$) $\cap$ DTM is considered as a $n_x \times n_y$ - grid, we will
have $\Delta x = 2\delta/n_x$ and $\Delta y = 2\delta/n_y$. So, $n_{\text{DTM}}$ as the number of the
prisms associated with Cube($\mathbf{r}, \delta$) $\cap$ DTM, is equal to $n_x \times n_y$ (i.e.

Computation Point

<figure><img src="image_26.png" /><figcaption>**Figure C1.** A schematic representation of the intersection of the neighbour-
hood Cube(**r**, *δ*) and the DTM, and the reduced height *h* ∈ [0, 2*δ*].</figcaption></figure>

$$
n_{\text{DTM}} = n_x \times n_y). \text{ Hence, eq. (C3) is changed into}
$$

$$
\begin{align}
V_{\text{Mass}}(\mathbf{r}, \delta) &\approx \sum_{\text{Cube} \cap \text{DTM}} \underline{h} \left( \frac{2\delta}{n_x} \times \frac{2\delta}{n_y} \right) \nonumber \\
&= \left( \frac{4\delta^2}{n_{\text{DTM}}} \right) \sum_{\text{Cube} \cap \text{DTM}} \underline{h} \tag{C4}
\end{align}
$$

Moreover, suppose that $V_{\text{Cube}}(\delta)$ is the whole volume of the neigh-
bourhood. Therefore, we have

$$
V_{\text{Cube}}(\delta) = (2\delta)(2\delta)(2\delta) = 8\delta^3 \quad (\text{C5})
$$

Furthermore, it is clear that the volume of the topographic masses
enclosed in the neighbourhood (i.e. $V_{\text{Mass}}$) is less than or equal to
the whole volume of the neighbourhood (i.e. $V_{\text{Cube}}$). Thus,

$$
0 \le V_{\text{Mass}} \le V_{\text{Cube}}, \quad 0 \le \frac{V_{\text{Mass}}}{V_{\text{Cube}}} \le 1 \qquad (\text{C6})
$$

Besides, the definition of the topographic coefficient *k*<sub>δ</sub> indicates
that this coefficient is 4 times the ratio of the volume of the enclosed
masses to the whole volume of the neighbourhood (see Theorem 2).
Accordingly, the topographic coefficients *k* and *k*<sub>δ</sub> can be approx-
imated by the corresponding definition and the above-mentioned
equations as follows:

$$
\begin{align}
k &\approx k_{\delta} \approx 4 \times \frac{V_{\text{Mass}}(\mathbf{r}, \delta)}{V_{\text{cube}}(\delta)} \nonumber \\
  &\approx \left( \frac{2}{\delta \times n_{\text{DEM}}} \right) \sum_{\text{Cube} \cap \text{DEM}} \underline{h} \tag{C7}
\end{align}
$$

Fig. C2(a) illustrates a graph of the average of the topographic
coefficient *k*<sub>δ</sub> over the region of the case study, against the neigh-
bourhood dimension *δ*. These results were obtained by applying the
90 m × 90 m-DEM of SRTM to eq. (C7). Also, Fig. C2(b) illustrates
a graph of the corresponding standard deviation of the estimated to-
pographic coefficient *k*<sub>δ</sub> against the neighbourhood dimension *δ*. In
these figures, it is clearly seen that the average value of the estimated
topographic coefficient *k*<sub>δ</sub> decreases to 2 and its corresponding stan-
dard deviation approaches zero as the neighbourhood dimension *δ*
grows in size.

Besides, Fig. C3 shows the spatial distribution of the estimated
values of the topographic coefficient *k*<sub>δ</sub> over the region of the case
study. These values of the topographic coefficient were obtained
by *δ* = 4.5 km (this value of the neighbourhood dimension had the
best performance at the 13 rock-sampling stations, see Section 4.1).
<figure><img src="image_27.png" /><figcaption>**Figure C2.** (a) A graph of the average of the topographic coefficient $k_δ$ against the neighbourhood dimension $δ$, (b) a graph of the standard deviation of $k_δ$ against the neighbourhood dimension $δ$.</figcaption></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>**Figure C3.** The spatial distribution corresponding to the estimated values of the topographic coefficient $k_δ$ over the region of the case study.</figcaption></figure>

Comparing Fig. C3 with the Earth's topography shown by Fig. 3 in Section 4, the considerable correlation between the topographic coefficient $k_δ$ and the Earth's shape can straightforwardly be recognized. Also, Table C1 provides the estimated topographic coefficient $k_δ$, corresponding to the gravity stations, with its summary statistics. This table indicates that the correlation between the estimated values of the topographic coefficient $k_δ$ and the elevations of the gravity stations is equal to $-0.64$.

Moreover, Table C1 indicates that the average value of the estimated topographic coefficient over the region of the case study is around 2 and we can suppose that $k = 2 \pm 0.055$.

In Section 4, we supposed that $k \approx 2$. Now, based on the above results, we can assess the effect of the approximation $k \approx 2$ on the estimation of the mass density. According to the formula $\rho \approx \frac{1}{-k\pi G}\text{trace}(\Delta\mathbf{F}_g\Delta\mathbf{R}^+)$ (i.e. eq. 22), the general law of propagation of variances implies that

$$ \sigma_{\rho}^{2} = \left(\frac{\partial \rho}{\partial k}\right)^{2} \sigma_{k}^{2} $$

$$ \begin{aligned} & \approx \left( \frac{1}{k^2 \pi G} \text{trace}(\Delta \mathbf{F}_g \Delta \mathbf{R}^+) \right)^2 \sigma_k^2 \\ &= \frac{\sigma_k^2}{k^2} \rho^2 \end{aligned} \qquad (C8) $$

where $\sigma_{\rho}^2$ is the variance of the mass density, caused by the variance of $k$ (i.e. $\sigma_k^2$). Hence, we have

$$ \sigma_{\rho} \approx \frac{\sigma_k \rho}{k} \qquad (C9) $$

Accordingly, the standard deviation $\sigma_{\rho}$ corresponding to the mass density is assessed by the following relation, considering $k = 2 \pm 0.055$ and $\rho = 2200 \text{ kg m}^{-3}$.

$$ \sigma_{\rho} \approx \frac{\sigma_k \rho}{k} = \frac{0.055 \times 2200}{2} = 60.5 \text{ kg m}^{-3}. \qquad (C10) $$

Comparing the above error with the other error sources such as the truncation error, we realize that the effect of the topographic coefficient on the estimated mass density will not be significantly considerable if $\delta = 4.5$ km. Nevertheless, Fig. C2(b) indicates that the standard deviation of the topographic coefficient will increase and its corresponding effect on the estimated mass density will also amplify if the neighbourhood decreases in size.

Besides, we investigated the results of the 13 rock-sampling stations at two cases. The first case supposed that $k = 2$, and the second case supposed that $k \approx k_δ$ obtained by both the DTM and eq. (C7). The comparison between these two cases indicated that the accurate values of the topographic coefficient, evaluated by the above-mentioned procedures, could improve the results and enhance the accuracy of the estimated mass density to some extent, see Table C2.

<table><caption>**Table C1.** The summery statistics of the topographic coefficient $k_δ$ corresponding to the gravity stations, when $\delta = 4.5$ km.</caption><thead><tr><th></th><th>Min</th><th>Average</th><th>Max</th><th>Standard deviation</th><th>Correlation with elevation</th></tr></thead><tbody><tr><td>$k_δ$</td><td>1.734</td><td>2.003</td><td>2.147</td><td>0.055</td><td>-0.0640</td></tr></tbody></table>

<table><caption>**Table C2.** The comparison of the mass density estimated at the two cases: (1) $k = 2$, and (2) $k \approx k_δ$ obtained by the SRTM's DTM (unit: kg m⁻³).</caption><thead><tr><th></th><th>Station</th><th>S1</th><th>S2</th><th>S3</th><th>S4</th><th>S5</th><th>S6</th><th>S7</th><th>S8</th><th>S9</th><th>S10</th><th>S11</th><th>S12</th><th>S13</th><th>RMS</th></tr><tr><th></th><th>ρ (average)</th><th>2076</th><th>2392</th><th>2513</th><th>2418</th><th>2438</th><th>1975</th><th>1905</th><th>2030</th><th>2120</th><th>1940</th><th>1935</th><th>2364</th><th>2291</th><th></th></tr></thead><tbody><tr><td>$k = 2$</td><td>ρ</td><td>2239</td><td>2302</td><td>2292</td><td>2470</td><td>2448</td><td>2259</td><td>2261</td><td>2115</td><td>2020</td><td>1901</td><td>1901</td><td>2275</td><td>2275</td><td><b>157</b> [kg m<sup>-3</sup>]</td></tr><tr><td></td><td>Error</td><td>163</td><td>-90</td><td>-221</td><td>52</td><td>10</td><td>284</td><td>356</td><td>85</td><td>-100</td><td>-39</td><td>-34</td><td>-89</td><td>-16</td><td></td></tr><tr><td>$k \approx k_δ$</td><td>ρ</td><td>2235</td><td>2377</td><td>2320</td><td>2356</td><td>2356</td><td>2172</td><td>2172</td><td>2033</td><td>2050</td><td>1846</td><td>1912</td><td>2249</td><td>2249</td><td><b>128</b> [kg m<sup>-3</sup>]</td></tr><tr><td></td><td>Error</td><td>159</td><td>-15</td><td>-193</td><td>-62</td><td>-82</td><td>197</td><td>267</td><td>3</td><td>-70</td><td>-94</td><td>-23</td><td>-115</td><td>-42</td><td></td></tr></tbody></table>