## 2 Materials and Methods
### Theoretical Background
#### 2.1.4 Bodies Growth Inversion Algorithm  
The bodies growth algorithm is a non-linear inversion method which estimates the geometry of the anomalous source by the addition of previously known density/magnetic susceptibility contrast prisms to the model. The main advantages of this algorithm are that: (1) it admits irregular prisms ensembles, (2) it allows determination of the regional tendency of the observed data, (3) it does not limit the growth to determined region of the domain, (4) it allows fitting of the positive and negative contrasts, (5) it allows variations on the contrast's values during the process, and (6) it takes into account the measurement's errors. The initial model \(\mathbf{p}_{0}\) can be a set of empty (null contrasts) prisms, or some prisms can be initialized with non-null contrasts, so the anomalous body tends to grow in the periphery of the initial model [40, 44]. The inversion process consists in an iterative procedure to find the adequate prisms so to assign to it a value of the physical property. In the \(k\)-th iteration of the algorithm, \(k\) prisms have been adjusted and \(M-k\) prisms remain unchanged. If the forward response of the model in the \(k\)-th iteration is given by \(\mathbf{m}^{[k]}\), and if in the next iteration the \(j=M-k+1\)-th prism is updated, the forward response of the model in the \(k+1\)-th can be easily updated as  
\[\mathbf{m}^{[k]}=\mathfrak{K}\mathbf{p}_{\lambda}\, \tag{11}\]  
where \(\mathbf{p}_{\lambda}\) is a zero vector except in its \(\lambda\)-th position, where it has the fitted contrast value of the added prism. For an isolated geological structure, the difference between the computed forward response of the model in the \(k+1\)-th iteration and the observed (measured) values \(\mathbf{m}^{\rm(O)}\) is given by  
\[\Delta\mathbf{m}^{[k+1]}=\mathbf{m}^{\rm(O)}-f^{[k+1]}\mathbf{m}^{[k+1]}-d_{\rm reg}^{[k+1 ]}\, \tag{12}\]  
where \(f^{[k+1]}\geq 1\) is a scale factor used to fit \(\mathbf{m}^{\rm(O)}\) and \(\mathbf{m}^{[k+1]}\), and \(d_{\rm reg}^{[k+1]}\) is the regional field at the \(k+1\)-th iteration, linearly approximated by  
\[d_{\rm reg}^{[k+1]}=c_{0}^{[k+1]}+c_{x}^{[k+1]}(x-x_{M})+c_{y}^{[k+1]}(y-y_{M} )\, \tag{13}\]where \((x_{M},y_{M})\) are the coordinates of the mid (central) point of the domain, and \(c_{0},c_{x},c_{y}\) s are coefficients to fit the regional tendency. Thus, the coefficients to be optimized are \(f_{k}\), \(c_{0}^{k}\), \(c_{x}^{k}\), \(c_{y}^{k}\). For gravimetric and magnetic data, the cost or objective function to fit such coefficients is given by  
\[\phi\left(\mathbf{p}^{[k+1]}\right)=e\left(\mathbf{p}^{[k+1]},\lambda\right)^{2}=\Delta \mathbf{m}^{[k+1]^{\mathsf{T}}}\Delta_{m}^{-1}\Delta\mathbf{m}^{[k+1]}+\lambda f^{[k+ 1]^{2}}\mathbf{p}^{[k+1]^{\mathsf{T}}}\Delta_{p}^{-1}\mathbf{p}^{[k+1]}\, \tag{14}\]  
where T stands for transpose, \(\Delta_{m}\) is the measurements (response) covariance matrix (\(N\times N\)) and \(\Delta_{p}=\bar{\kappa}^{\mathsf{T}}\Delta_{d}^{-1}\bar{\kappa}\) is the covariance matrix (\(M\times M\)) parameter. For gradiometric data, the total cost function is the sum of the cost functions of each component of the tensor being inverted. The first term in Equation (14) represents the misfit of the response normalized by \(\Delta_{m}\), the second term is the sum of the parameter vector weighted by the \(\lambda\) parameter (positive), which value should be chosen carefully. Low values of \(\lambda\) provide good fitting of data, but the model may grow excessively, creating fictitious structures, while high values of \(\lambda\) yield simple and shallow models that poorly adjust measurements.