## 3 Methodology  
[PERSON]'s law of universal gravitation provides the vertical component of the gravity acceleration at a point \(\vec{r_{0}}=(x_{0},y_{0},z_{0})\) given by  
\[g_{z}(x,y,z)=\gamma\int_{\Omega}\rho(x,y,z)\frac{z_{0}-z}{\left[(x_{0}-x)^{2}+ (y_{0}-y)^{2}+(z_{0}-z)^{2}\right]^{3/2}}\mathrm{d}x\mathrm{d}y\mathrm{d}z, \tag{1}\]  
where \(\gamma=6.672\times 10^{-11}\,\mathrm{m}^{3}\mathrm{kg}^{-1}\mathrm{s}^{-2}\) is the gravitational constant, \(\rho\) the density, and \(\Omega\) the volume for the Cartesian coordinates of the body. We considered the gravitational attraction due to a rectangular prism of constant density [32; 33], where the analytic solution of Equation (1) will be  
\[g_{z}=\rho\gamma\ \left[z\tan^{-1}\Bigl{(}\frac{xy}{z\tau}\Bigr{)}-x\ln(y+r)-y \ln(x+r)\Bigr{]}\Bigr{|}^{\Delta x_{2}}_{\Delta x_{1}}\Bigr{|}^{\Delta y_{2}} _{\Delta y_{1}}\Bigr{|}^{\Delta z_{2}}_{\Delta z_{1}}, \tag{2}\]  
where \(r=\sqrt{x^{2}+y^{2}+z^{2}}\), \(\Delta x_{k}=x_{0}-x_{k}\), \(\Delta y_{k}=y_{0}-y_{k}\), and \(\Delta z=z_{0}-z_{k}\) for \(k=1,2\) for the corners of the prism (Figure 2a). The gravity observation vector \(g_{i}\) (\(i=1,N\)) was computed by adding the contribution for all of the \(j\)-th prisms (\(j=1,M\)) of the density model, where \(N\) is the total of gravimetric stations measured along the surface and \(M=n_{x}\times n_{y}\times n_{z}\) is the number prism of the model. Typically, the 3D gravimetric forward modeling is accomplished by computing Equation (2) for a considerable number of \(n_{z}\) grid points. For the depth of the basement estimation, only two prisms on the \(z\)-direction are considered (Figure 2b) for a large number of \(x,y\) grid points. The upper \(j\)-th prism will contain the density contrast (\(\Delta\rho_{1}\)) between the surface (\(z_{\mathrm{min}}\)) and depth \(z_{j}\) and the lower \(j\)-th prism will contain the density (\(\Delta\rho_{2}\)) for the bottom part of the model up to \(z_{\mathrm{max}}\) for the maximum depth of study. In this case, \(z_{j}\) aimed for the depth of the basement, and the density contrast was considered to be known. Figure 1: Location of Zacatecas State in Mexico (**a**) and the Aguanaval Aquifer polygon (**b**) for the study. For the inverse problem, the objective function \(Q\) proposed seeks to minimize the discrepancies between the observed gravity data \(\mathbf{d}^{\text{obs}}\) and the synthetic \(\mathbf{d}^{\text{cal}}\) using the L-2 norm as follows:  
\[Q(\mathbf{m})=||\mathbf{C}_{d}(\mathbf{d}^{\text{obs}}-\mathbf{d}^{\text{cal}} )||^{2}+\alpha_{\text{reg}}^{2}||\
abla\mathbf{m}||^{2}+||\mathbf{C}_{m}( \mathbf{m}-\mathbf{m}_{\text{apr}})||^{2}, \tag{3}\]  
where \(\mathbf{m}\) is the model vector that minimizes \(Q\) and \(\mathbf{C}_{d}\) is the diagonal matrix containing the covariance of the data. Here, a Tikhonov regularization constraint on the model was implemented, pondered by the parameter \(\alpha_{\text{reg}}\). The right term corresponds to the a priori information of the model, being \(\mathbf{C}_{m}=\alpha_{\text{apr}}\mathbf{I}\) for a whole constraint on the model at the inversion. If borehole data are available, the matrix \(\mathbf{C}_{m}\) will act as a weight on the regions of the model where the well is located, which will significantly reduce the non-uniqueness of the problem and will improve the accuracy of the inversion. For density estimation of the model, the model parameter \(\mathbf{m}\equiv\mathbf{m}^{\rho}\). It is well known that the inverse problem is linear for density estimation, as seen in Equation (2).