## 1 Introduction  
Mexico's population demands a water volume of 88,839.7 hm\({}^{3}\) annually for several needs [1]. Despite the wrongful belief that industry demands most of it, agriculture use 75.7% of the total water, where 63.4% of it is obtained on the surface and the other 36.6% from the subsoil. However, the major source of water for semi-arid regions corresponds to groundwater systems. Zacatecas state is one of such regions with large periods of drought, with a concession of 1688.4 hm\({}^{3}\) of water [2], where 85.3% is destined for agricultural activities to overcome the low precipitation problems in the region. Historically, the groundwater budget has been the principal tool to determine the availability of water within the aquifers in Mexico; however, several authors have demonstrated the need for the construction of groundwater flow models to support sustainable water management policies [3, 4]. To achieve this goal, the understanding of groundwater flow systems needs to be increased, which requires a better knowledge of the hydrogeological properties of the subsoil, where data are usually scarce and sometimes outdated [5, 6]. A fundamental task for further studies of groundwater functioning is the determination of basement geometry [7], which is necessary for the construction of conceptual and numerical groundwater models to achieve the sustainable management of water resources [8, 9]. This task can be achieved using remote sensing, such as unmanned aerial vehicles [10, 11, 12, 13] for near-shore studies. On the other hand, satellite data have become the most popular solution when equipment is not available, providing high resolution and faster results on a larger scale for multispectral imagery [14] or satellite gravity data. Gravimetric studies have been widely used for depth to the basement estimation for geological mapping, exploring oil reservoirs, and modeling groundwater systems. The basement topography study is important for the understanding of the groundwater flow, which happens within sedimentary deposits above igneous or metamorphic rocks, having notorious density contrasts that produce gravity anomalies that provide structural information about the boundary between those layers. Such gravity interpretations can be grouped into two geophysical data inversion categories: density estimation and depth estimation. For density estimation, the geometry of the model is set to be known, and then the density of the cells or grid points is determined by linear data inversion [15; 16], which can also be performed by incorporating a priori information and depth-dependence weights in the inversion to recover a density model with a more proper structure at different depths [17; 18]. Alternatively, for depth estimation, a density contrast is assumed and non-linear data inversion is performed for the geometry of a parallelepiped. Such a gravity contrast can be considered constant [19; 20] or varying with depth using exponential [21], quadratic [22], cubic [23], or hyperbolic [24] functions. These studies are well-tested for near-shore data in sedimentary basins where borehole data for oil exploration are available; however, such gravity constant functions for depth can be difficult to extrapolate for inland basin studies, where the geology is more complex. With respect to the optimization, these works minimize an objective function performing Gauss-Newton optimization, which leads to a square matrix inversion that can be computationally expensive for high-resolution problems. Nowadays, additional optimization alternatives are being continuously developed for other geophysics areas, such as seismic methods, using conjugate gradient techniques, significantly reducing the computation cost of the inversion. In this context, the conjugate gradient least squares (CGLS) method was developed as an alternative for Gauss-Newton methods, reducing the computational time and storage.