Non-Linear 3D Satellite Gravity Inversion for Depth to the Basement Estimation in a Mexican Semi-Arid Agricultural Region

Abstract

In Mexico, agriculture in semi-arid regions is highly dependent on groundwater resources, where most of the aquifers’ characterization is a pending task. In particular, the depth to the basement is unknown for most of the Mexican territory. Hence, the development and performance of new techniques for the basement relief estimation is imperative for further hydrogeological studies. In this paper, we present a depth to the basement estimation using non-linear gravimetric inversion employing satellite data. Gravity forward modeling was implemented using both gravitational attraction due to juxtapositioned blocks and gravimetric non-linear inversion using conjugate gradient least squares to minimize the objective function in terms of a depth model. All of this took place under the sparse system framework. We present a synthetic result using the SEG-Bishop depth model taken for calibration purposes. Then, we recollected gravity data from The Satellite Geodesy group from SCRIPPS for the depth to the basement estimation of an unconfined aquifer in the northern-central semi-arid region of Zacatecas, Mexico. Both synthetic and satellite data were recovered, consistent depth models for both cases were presented, and a comparison with conventional gravimetric linear inversion for density estimation was performed.

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.

We present a non-linear gravimetric inversion for depth to the basement estimation for an unconfined aquifer of the north-central region of Zacatecas, Mexico. The inverse problem was optimized by developing a CGLS algorithm for obtaining a depth model. A constant density contrast was assumed due to the lack of a priori information and the complexity of the region of interest.The algorithm was calibrated by computing synthetic gravity data using the Bishop model, recovering the depth model with the inversion scheme proposed. The procedure was repeated using gravity inversion with satellite data recollected from The Satellite Geodesy group from SCRIPPS. Finally, a comparison between density estimations and depth to the basement estimation was performed.

2. Geology and Previous Studies

The study area was located in the Mexican Altiplano, constituting part of the Central Terrane and most of the Guerrero Composite Terrane, formed between Late Jurassic and early Cretaceous [25]. This terrane is characterized by volcanic–sedimentary successions formed in a back-arc tectonic setting; we address Centeno-Garcia [26] for a complete compilation of works in this area. Locally, the study area was in Zacatecas State (Figure 1a), at the region of the Aguanaval (unconfined) aquifer (Figure 1b), situated in Fresnillo and Rio Grande municipalities. Historically, this area has been an important mine production source, mostly for silver exploitation, since the XVI century. The stratigraphy of this region is well-exposed and made up of shale, quartz-rich sandstone, scarce limestone, and basaltic pillow lavas [27]; however, basement rocks do not crop out in the Fresnillo area [28] and information on its depth remain unknown.

Agricultural and industrial activities predominate within the Aguanaval aquifer, which over-exploit the groundwater that, in itself, has low recoverability due to low precipitations. Other than geological, works in this region involve hydro-geological [29] and technical reports such as Comisión Nacional del Agua (CONAGUA) [2]. Despite the mining, industrial and agricultural importance of this region, there is a lack of geophysical academic works published—only AMT [30] and TEM [31] studies are found—where gravity and seismic interpretations need to be carried out for a better understanding of the basement relief; therefore, the importance of this study for depth to the basement estimation is emphasized.

3. Methodology

Newton’s law of universal gravitation provides the vertical component of the gravity acceleration at a point $\overrightarrow{r_0}=(x_0,y_0,z_0)$ given by

$$g_z(x,y,z)=\gamma {\int }_{\mathsf{\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,$$

(1)

where $\gamma =6.672\times {10}^{-11}$m${}^3$kg${}^{-1}$s${}^{-2}$ is the gravitational constant, $\rho$ the density, and $\mathsf{\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

$$\begin{array}{c}\hfill g_z=\rho \gamma {\left.{\left.{\left.\left[z{tan}^{-1}\left(\frac{xy}{zr}\right)-xln\left(y+r\right)-yln\left(x+r\right)\right]\right|}_{\Delta x_1}^{\Delta x_2}\right|}_{\Delta y_1}^{\Delta y_2}\right|}_{\Delta z_1}^{\Delta z_2},\end{array}$$

(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_{\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_{\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.

For the inverse problem, the objective function Q proposed seeks to minimize the discrepancies between the observed gravity data ${\mathbf{d}}^{\mathrm{obs}}$ and the synthetic ${\mathbf{d}}^{\mathrm{cal}}$ using the L-2 norm as follows:

$$Q\left(\mathbf{m}\right)=\left|\right|{\mathbf{C}}_d({\mathbf{d}}^{\mathrm{obs}}-{\mathbf{d}}^{\mathrm{cal}}){\left|\right|}^2+{\alpha }_{\mathrm{reg}}^2{\left|\right|\nabla \mathbf{m}\left|\right|}^2+\left|\right|{\mathbf{C}}_m(\mathbf{m}-{\mathbf{m}}_{\mathrm{apr}}){\left|\right|}^2,$$

(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 }_{\mathrm{reg}}$. The right term corresponds to the a priori information of the model, being ${\mathbf{C}}_m={\alpha }_{\mathrm{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). The minimization of 3 with respect to the density model will result in the system

$$\left[{\mathbf{A}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{}\mathbf{A}+{\alpha }^2{\mathbf{D}}^{\mathbf{T}}\mathbf{D}+{\mathbf{C}}_{\mathbf{m}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{m}}^{}\right]{\mathbf{m}}^{\rho }={\mathbf{A}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{}{\mathbf{d}}^{\mathrm{obs}}+{\mathbf{C}}_{\mathbf{m}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{m}}^{}{\mathbf{m}}_{\mathrm{apr}}^{\mathbf{\ae }},$$

(4)

where $\mathbf{D}$ is the discrete operator of ∇ and $\mathbf{A}$ is the density kernel equal to the derivative of Equation (2) with respect to $\rho$ (linear). For depth estimation of the model, the model parameter $\mathbf{m}\equiv {\mathbf{m}}^z$, where z is the depth of the prism (Figure 2b) and z does not denote exponent. By inspecting Equation (2), it is evident that the inversion by depth is a non-linear problem. This issue is solved by performing a linearization around a starting model ${\mathbf{m}}_{\mathbf{0}}^{\mathbf{z}}$, resulting in the following system:

$$\left[{\mathbf{J}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{}\mathbf{J}+{\alpha }^2{\mathbf{D}}^{\mathbf{T}}\mathbf{D}+{\mathbf{C}}_{\mathbf{m}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{m}}^{}\right]{\mathbf{m}}^z={\mathbf{J}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{d}}^{}({\mathbf{d}}^{\mathrm{obs}}+{\mathbf{Jm}}_0^z)+{\mathbf{C}}_{\mathbf{m}}^{\mathbf{T}}{\mathbf{C}}_{\mathbf{m}}^{}{\mathbf{m}}_{\mathrm{apr}}^{\mathbf{z}}.$$

(5)

In this case, the matrix $\mathbf{J}$ is the Jacobian corresponding to the derivative of Equation (2) with respect to z evaluated on the starting model ${\mathbf{m}}_0^z$; with this, Frechet derivatives are avoided, e.g, methods such as perturbations and sensitivity approaches, which can be computationally expensive. The solution of Equations (4) and (5) can be achieved by inverting the symmetric square matrix of the left size, commonly called ${\mathbf{A}}^{\mathbf{T}}\mathbf{A}$ inversion. However, this implementation requires the storage and computation of such a matrix, which results in a large computational time, and sometimes the memory is surpassed for large scale problems. One solution for this problem is to perform conjugate gradient least squares to solve a system in the form of $\mathbf{Gm}=\mathbf{d}$. This method internally computes the square matrix of Equation (4), allocating only the matrix $\mathbf{A}$ and the sparse matrix $\mathbf{D}$, ${\mathbf{C}}_m$, and ${\mathbf{C}}_d$. For example, for depth inversion, the system $\mathbf{Gm}=\mathbf{d}$ is solved as

$$\left[\begin{array}{c}{\mathbf{C}}_{\mathbf{d}}\mathbf{J}\\ {\alpha }_{\mathrm{reg}}\mathbf{D}\\ {\mathbf{C}}_{\mathbf{m}}\end{array}\right]{\mathbf{m}}^z=\left[\begin{array}{c}{\mathbf{C}}_{\mathbf{d}}({\mathbf{d}}^{\mathrm{obs}}+{\mathbf{Jm}}_0^z)\\ \mathbf{0}\\ {\mathbf{C}}_m^{}{\mathbf{m}}_{\mathrm{apr}}^z\end{array}\right],$$

(6)

In this case, the augmented matrix $\mathbf{G}$ were large and sparse due to the discrete operations for the Tikhonov regularization, and the model vector ${\mathbf{m}}^{\rho }$ was not modified. Notice that multiplication of Equation (6) by ${\mathbf{G}}^T$ yields to Equation (5). The algorithm to solve the system $\mathbf{Gm}={\mathbf{d}}^{\mathrm{CG}}$ is presented in Algorithm 1, modified from Nocedal & Wright [34]. Notice that an efficient implementation of this algorithm requires that all the matrices be stored in a sparse representation. We used coordinate list (COO) for storing such sparse matrices, given that the matrix and its transpose have to be accessed simultaneously [35].

Algorithm 1 CGLS algorithm to iteratively solve the problem $Gm=d^{\mathrm{CG}}$

Require:$m_0$ (Starting model), G (Augmented kernel), $d^{\mathrm{CG}}$ (Data Vector), n (Iterations).

Ensure:$Gm=d^{\mathrm{CG}}$

1: $d_0\Leftarrow d^{\mathrm{CG}}$

2: $r_0\Leftarrow G^T{d}^{\mathrm{CG}}$

3: $p_0\Leftarrow r_0$

4: $t_0\Leftarrow G{p}_0$

5: while $k<n$ or stop criteria is satisfied do

6:   ${\alpha }_k\Leftarrow \parallel r_{k-1}{\parallel }^2/{\parallel t_{k-1}\parallel }^2$

7:   $m_k\Leftarrow m_{k-1}+{\alpha }_k{p}_{k-1}$

8:   $d_k\Leftarrow d_{k-1}-{\alpha }_k{t}_{k-1}$

9:   $r_k\Leftarrow G^T{d}_k$

10:  ${\beta }_k\Leftarrow \parallel r_k{\parallel }^2/{\parallel r_{k-1}\parallel }^2$

11:  $p_k\Leftarrow r_k+{\beta }_k{p}_{k-1}$

12:  $t_k\Leftarrow G{p}_k$

13:  $k\Leftarrow k+1$

14: end while

15: return $m_k$ (Final model)

4. Results

We tested the algorithm developed on a single synthetic model for calibration purposes. First, a non-linear inversion was performed for the depth estimation of the basement. Finally, we present the results due to satellite gravity data for the depth to the basement estimation of the area of interest.

4.1. Synthetic Case: Bishop Model

In order to test the developed algorithm, we performed the inversion by depth for the Bishop model, a synthetic basement model created by upscaling and shifting the topographic data of an area of Bishop, CA, USA [36]. This model contains a variety of structures that resemble the complexity of real basement relief: two long faults, en echelon arrays of smaller faults, transfer zones between faults, and an unfaulted deep basin area. This features allows the model to be exploited by various authors to test new methods for depth to the basement estimation using magnetic data [36,37] and gravity data [38,39].

In the north-east region of the model, the shallower parts of the basement (≈100 m) are presented and, in the south-west region of the model, the deepest parts of the basement are located up to ≈9300 m, as seen in the 3D view of Figure 3 and the 2D floor view of Figure 4a, which acts as an unknown true depth model in the inversion. A total of $n_x\times n_y\times 2=61\times 61\times 2=3162$ prisms are considered in a region of interest $0\le x\le 380$ km for the east and $0\le y\le 402$ km for the north, given a spacing of $\Delta x=6.33$ km and $\Delta y=6.7$ km (width and length) for each prism. The density for the upper prism is ${\rho }_1=2100$ g/cm${}^3$, corresponding to sedimentary rocks above the basement of density $\Delta {\rho }_2=2750$ g/cm${}^3$ for the lower prisms. 3D gravity forward modeling was performed, measuring along all of the grid points on the surface, giving a total of $n_x\times n_y\times 2=61\times 61=1581$ gravimetric stations (length of the observed data vector $d^{\mathrm{obs}}$). The gravity data show anomalies corresponding to the high and lower depths of the Bishop depth to the basement model (Figure 4b).

For geophysical data inversion, a priori information was not considered; therefore, ${\mathbf{C}}_m=0$ and the parameter of regularization were chosen using the $L-$curve method for Tikhonov regularization. We performed 10 iterations of the CGLS method (Equation (6)) starting from a constant depth model ${\mathbf{m}}_0^z$ of 1000 m. It is worth mentioning that the starting model does not affect the final solution, but it can help for a fast convergence of the problem. The depth model after 10 iterations of the CGLS is shown in Figure 4c, and the final synthetic response due to this model is shown in Figure 4d.

The result obtained clearly exhibits a successful recovery of the structure of the basement with the correct depths and, henceforth, the gravimetric anomalies were fitted properly, with few discrepancies. These differences can be observed in Figure 4e,f for the model error and the gravity residual, respectively. For the depth model, errors up to ≈±0.7 km appear, mainly located at regions where there is an abrupt change in lateral variation in density (see Figure 3); however, in regions where the basement behaves predominantly horizontal, the errors are almost negligible (≈±0.01 km or ≈±10 m) as seen in the scale of Figure 4a,c. Although this error happens where lateral variation is presented, we infer that this is a product of the Tikhonov regularization constraint; the smoothness interposed does not recover parts where the total variation in the model is presented. We address [40] for an optional regularization where a focused model is needed. Consequently, this error in the model also produces a gravity residual in such lateral-variation regions of up to 2.15 mGals, which is not harmful for the final results.

A single 1D depth to the basement profile is presented, located at $y=201$ km for an x-variation as seen in Figure 5a. Again, the structure of the Bishop model is well-recovered in regions where there is not an abrupt change in density ($x=280$ km). With respect to the convergence of the algorithm, the plot of the objective function (Misfit) for each iteration shows a fast convergence at the earliest iterations (Figure 5b). The CGLS algorithm encounters a stagnation in the misfit after 10 iterations, which is a sign that the minimum has been reached.

4.2. Depth to the Basement Estimation for Real Gravity Data

The study area covers a minima and maxima latitude of 22.58 and 23.69 and minima and maxima longitude of −103.525 and −102.625, respectively; therefore, this area surrounds the aquifer polygon as seen in Figure 1 and Figure 6a.

We recollected gravity data from The Satellite Geodesy research group at the Cecil H. and Ida M. Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California San Diego [41]. After geophysical processing, the observed gravity data (Figure 6a) of the area of study present consistencies with the geology and the topography information of the region of interest (Figure 6b). The interpretation of these gravity anomalies is straightforward: the sedimentary rock deposits, where the Aguanaval River streams, match with the lowest gravity anomalies, where the central part of the aquifer lies; on the other hand, the highest gravity anomalies are presented at the boundaries of the aquifer because these regions act as a source for the recharge of the aquifer due to high elevation.

The inversion result after 10 iterations of the algorithm is shown in Figure 6c. It is clear that the data have been fit for all zones, except for some regions with errors up to ≈2 mGal (see scale of Figure 6c), mostly at the regions where the gravity maximum and minimum are presented. Such errors can be reduced if the weight on the regularization is lowered; however, models with a higher roughness could be presented. As expected, the depth model presents a similar shape to the observed gravity data. A depth to the basement between 1.8 and 5.5 km is estimated, where the maximum depth is located at the central part of the aquifer limits as seen in Figure 6d for a 2D view and Figure 6e for a 3D view of the basement relief. The model does not present zero depths, which means that the basement does not crop out to the surface. Such depths are consistent with the topography and the geology of the study zone. The convergence of this result is well-behaved and fast in the same way as the synthetic example; therefore, the cost function will not be presented.

One-directional profiles are picked for interpretation, centered at the domain of the aquifer, horizontally and vertically. The north–south profile, along a constant latitude of 103${}^{\circ }$W, suggests the presence of an inland basin filled with sedimentary rocks according to the geology of the study zone. In the same way, the west–east profile covers the aquifer perpendicularly, where such a basin is narrower. This suggests that the aquifer is elongated in the north–south direction and slightly inclined in the same direction as the aquifer’s polygon.

5. Conclusions

We developed an algorithm for depth to the basement estimation using the iterative CGLS method for the optimization, employing large sparse storage for the system of equations. The algorithm was calibrated using the Bishop synthetic model, successfully recovering the depth model and the observed gravity data, with neglectable errors due to Tikhonov regularization.

For satellite gravity data, the algorithm was applied for an unconfined aquifer located in the north-central region of Zacatecas, Mexico. The estimated relief has a depth between 1.8 and 5.5 km, where the deepest region is presented at the central part of the aquifer and shallower regions are approximately at the boundaries of the aquifer limits. This relief matches with the geology and the topography of the zone; however, the calibration remains as pending due to the lack of borehole data information (for an a priori model in the inversion) and the null acquaintance of other geophysical data available as an academic work or field report. Hence, we catalogued this depth to the basement model as an estimation, subject to modifications if a priori information could be obtained, but the relief will keep the same structure.

Finally, we conclude that this is an important step for a better understanding of the subsurface of the area of interest since the processes that occur within the aquifer need to be detailed, such as groundwater flow, for the sustainable management of this natural and vital resource.