## 2 Materials and Methods  
### Theoretical Background  
#### Potential Field Theory  
The \(\alpha\)-th component of the gravitational field vector \(\mathcal{G}\) measured in an observation point \(P\) at a position \(r\) due to a source body with mass density \(\rho\) is given by  
\[\mathcal{G}_{\alpha}(\mathbf{r})=\gamma\iiint_{\text{source}}\rho(\mathbf{r}^{\prime}) \left(\frac{\alpha-\alpha^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}dv^{\prime} \right)\, \tag{1}\]  
where \(\alpha=x,y,z\), \(\mathbf{r}^{\prime}\) is a vector pointing to the differential volume element \(dv^{\prime}\) of the source body and \(\gamma=6,67428\times 10^{-11}\text{m}^{3}\text{s}^{-2}\text{kg}^{-1}\) is the gravitational constant. In geophysical surveys, the commonly used units for gravitational field are mGals (1 Gal = \(1\times 10^{-5}\) m/s\({}^{2}\)). As an evolution of the gravity method, FTG method (full tensor gradient method) directly measures the derivatives of the gravitational vector field, which lead to a second-rank (covariant) symmetric tensor \(\mathfrak{G}\), whose \(\alpha\), \(\beta\)-th component is given by  
\[\mathfrak{G}_{\alpha\beta}(r)=\frac{\partial\mathcal{G}_{\beta}(\mathbf{r})}{ \partial\alpha}=\gamma\iiint_{\text{source}}\frac{3\left(\alpha-\alpha^{ \prime}\right)\left(\beta-\beta^{\prime}\right)-\left|\mathbf{r}-\mathbf{r}^{\prime} \right|^{2}\delta_{\alpha\beta}}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|^{5}}\rho \left(\mathbf{r}^{\prime}\right)dv^{\prime}. \tag{2}\]  
The gravity field gradient tensor also has the property that the sum of its trace is zero, satisfying Laplace equation \(\mathfrak{G}_{xx}+\mathfrak{G}_{yy}+\mathfrak{G}_{zz}=0\). In practice, some combinations of the components of the gradiometric tensor are also used, particularly \(\mathfrak{G}_{uv}\), which aids to determine the limits of the source body in N-S and E-W directions [20; 46],  
\[\mathfrak{G}_{uv}=\frac{\mathfrak{G}_{xx}-\mathfrak{G}_{yy}}{2}. \tag{3}\]  
In the case of gravimetry, from Equation (1), each prism contribution to the observed value of the vertical component of \(\mathcal{G}\) in the observation point \(r\) can be numerically approximated by [47]  
\[\mathcal{G}_{z}=\gamma\rho\sum_{i,j,k=1}^{2}\mu_{ijk}\left[z_{k}\arctan\frac{x_ {i}y_{j}}{z_{k}r_{ijk}}-x_{i}ln(r_{ijk}+y_{j})-y_{j}ln(r_{ijk}+x_{i})\right]\, \tag{6}\]  
where \(x_{i}=x_{0}-\zeta_{i},y_{i}=y_{0}-\eta_{i},z_{i}=z_{0}-\zeta_{i}\), \(r_{ijk}=\sqrt{x_{i}^{2}+y_{i}^{2}+z_{i}^{2}}\) and \(\mu_{ijk}=(-1)^{i}(-1)^{i}(-1)^{k}\)[48]. On the other side, the independent components of the gradiometric tensor (Equation (2)) can be numerically computed by [48; 49]  
\[ \begin{split}\mathfrak{G}_{xx}&=\gamma\rho\sum_{i, j,k=1}^{2}\mu_{ijk}\arctan\frac{y_{j}z_{k}}{x_{i}r_{ijk}}\,\qquad\qquad\mathfrak{G}_{xy}&=\gamma\rho\sum_{i,j,k=1}^{2} \mu_{ijk}\ln\left(z_{k}-r_{ijk}\right)\,\\ \mathfrak{G}_{yy}&=\gamma\rho\sum_{i,j,k=1}^{2}\mu_{ ijk}\arctan\frac{x_{i}z_{k}}{y_{j}r_{ijk}}\,\qquad\qquad\mathfrak{G}_{xz}&=\gamma\rho\sum_{i,j,k=1}^{2} \mu_{ijk}\ln\left(y_{j}-r_{ijk}\right)\,\\ \mathfrak{G}_{zz}&=\gamma\rho\sum_{i,j,k=1}^{2}\mu_{ ijk}\arctan\frac{x_{i}y_{j}}{z_{k}r_{ijk}}\,\qquad\qquad\mathfrak{G}_{yz}&=\gamma\rho\sum_{i,j,k=1}^{2} \mu_{ijk}\ln\left(x_{i}-r_{ijk}\right)\,\end{split}  \tag{7}\]  
where \(x_{i}=x_{0}-\xi_{i},y_{i}=y_{0}-\eta_{i},z_{i}=z_{0}-\zeta_{i}\), \(r_{ijk}=\sqrt{x_{i}^{2}+y_{i}^{2}+z_{i}^{2}}\) and \(\mu_{ijk}=(-1)^{i}(-1)^{i}(-1)^{k}\). Finally, in order to numerically compute the total magnetic field anomaly produced in an observation point \(P=(x,y,0)\) at the surface \(z=0\) by a prism with arbitrary constant magnetization, we follow the algorithm presented in [50; 51] which consists in computing the magnetic anomaly due to an infinitely extended in-depth prism (Figure 2). Equation (5) in the schema proposed by [50] takes the form  
\[