Let dm be a physical property of a particle or elemental volume of matter (e.g., mass or electric charge) located at the r′ position in space; a set of these particles will interact with a certain force depending on the associated properties as gravitational or electric force. Historically, the mathematical description of these forces was given independently in what it is described as some fundamental laws of physics, for instance the law of universal gravitation of Newton. For this law, the differential gravity potential d U r is inversely proportional to the distance between a source point at r ′ = x ′ , y ′ , z ′ and the measurement point at r = x , y , z and can be expressed as d U r = − γ g d m r ′ r , (1) where γ _g is the gravitational constant, r = | r − r ′ | = x − x ′ 2 + y − y ′ 2 + z − z ′ 2 is the source–receiver distance, and dm is the fraction of mass as shown in Figure 1. The differential contribution on the gravity vector d g is conservative, irrotational, and according to the Helmholtz theorem, it can be described by the gradient of the differential gravity potential (e.g., Blakely, 1996), d g r = − γ g d m r ′ ∇ 1 r . (2) For a finite volume V′, the total potential can be described in the form of a volume integral (e.g., Zhdanov et al., 2004) U r = − γ g ∭ V ′ 1 | r − r ′ | ρ r ′ d v ′ . (3) A similar treatment can be performed for dipolar sources such as the magnetic field B . A stationary magnetic field B can be described by a scalar magnetic potential V r .

For a finite volume V′, the total magnetic potential V r can be computed using the magnetization vector M r ′ by: V r = − κ ∭ V ′ M r ′ ⋅ ∇ 1 | r − r | ′ d v ′ , (4) where κ = 4π × 10⁻⁷ H/m is the magnetic permeability of free space. Similar equations are also applicable to gradients of monopole fields including the gravity tensor.

It is important to note that all these fields depend on the relative position (e.g., Equations 3 and 4); these equations depend on a common factor that can be defined as a new function W 0 r = 1 | r | . This change of variable allows us to express Equations 3 and 4 in terms of three-dimensional convolution integrals, resulting in U r = − γ g ρ r ∗ W 0 r , (5) for the gravity potential and V r = − κ M r ⋅ ∗ ∇ W 0 r , (6) for the magnetic potential, where ⋅∗ denotes 3D convolution ∗ for every inner product ⋅ .

We can see that Equations 5 and 6 have the form of a convolution; this convolution can be solved for a finite volume using the appropriate W filter, computed by the field response of an individual prism for each depth layer.

In general, for a discretized volume divided in various layers of homogenous rectangular prisms with x − and y − regular dimensions (Figure 1), Equation 5 yields: U x l , y m , z n = γ g ∑ i ∑ j ∑ k ρ i j k W k D x l − x i , y m − y j , z n − z k , (7) where W k D α , β , γ = ∫ α α + Δ x ∫ β β + Δ y ∫ γ γ + Δ z k 1 x 2 + y 2 + z 2 d x d y d z . Using a discrete 2D convolution * for a Δx, Δy equally spaced grid: U = γ g ∑ k ρ k ∗ W k D . (8) With this notation, we may compute g _z as: g z = − γ g ∑ k ρ k ∗ W k g z , (9) where W k g z = ∂ W k D ∂ z , i.e., (following Banerjee and Das Gupta, 1977) W k g z α , β , γ = α ⁡ ln β + r + β ⁡ ln α + r − γ ⁡ tan − 1 α β γ ρ | α α + Δ x | β β + Δ y | γ γ + Δ z k . (10) The gradient of the gravity field (gravity gradient tensor or GGT) contains the information of the vertical and horizontal gradients as well as the gravity field curvature (Jekeli, 2006); this means that it defines better lateral contrasts and discriminate depths and improves structural or geometrical indicators of the field (Butler, 1995). Similarly, we may compute the full GGT T as T = − γ g ∑ k ρ k * W k T , (11) where W k T = ∇ ∇ W k T α , β , γ , i.e., W k T α , β , γ = tan − 1 γ β r α − ln r + γ − ln r + β − ln r + γ tan − 1 α γ r β − ln r + α − ln r + β − ln r + α tan − 1 α β r γ | α α + Δ x | β β + Δ y | γ γ + Δ z k . (12) We may also use the same function W k T for the magnetic field as: B = κ ∑ k M k ⋅ ∗ W k T . (13) The discrete convolution implies that if we have a fixed z position and a constant cell width, the derivative is the convolution of a characteristic 2D filter, which does not need to be evaluated individually. This allows for an efficient computation of sensitivities for thousands of observed data and horizontal cells. In fact, due to its symmetry, we only need to compute a single quadrant of the filter to cover the complete two-dimensional domain for each layer of cells (Figure 2), which can also be easily stored to avoid repeated computations when applying iterative procedures for inversion.

In order to solve the inverse problem, we use a quadratic norm to define the following objective function: ϕ = ‖ d o b s − A m ‖ C d d − 1 2 + ‖ m − m 0 ‖ C 00 − 1 2 + α p 2 ‖ D m ‖ 2 , (14) where d _obs are the observed data, A is the sensitivity matrix, m accommodates all the model parameters (in our case either density or magnetization values per cell), C _dd is the covariance matrix of the observed data d _obs (assumed diagonal), C ₀₀ is the covariance matrix of the a priori model m ₀, and D is a discrete derivative operator that depends on the α _p penalty term and gears the search toward smoothed property distributions as a regularizing constraint. Am is the 3D model response that includes all the convolution coefficients indicated by Equations 9, 11, and 13.

We use a least-squares minimization, i.e., deriving (14) with respect to m and equating to zero for each component (Tarantola and Valette, 1982). We then obtain a linear system of equations for the optimal m ̂ estimator given by: A T C d d − 1 A + C 00 − 1 + α p 2 D T D m ̂ = A T C d d − 1 d o b s + C 00 − 1 m 0 , (15) where T means transposed.

This system of linear equations can be solved using various linear algebra strategies. Direct solutions, however, may be prohibitive for large-scale three-dimensional problems, whereas the use of iterative schemes have to either face the repeated computation of the sensitivity matrix A or the storage of their compressed versions, usually at the cost of losing some resolution power.

To solve (15), we adapt the preconditioned conjugate gradient method (CG) described by Nocedal and Wright (2006) and take full advantage of the convolution property of the gravity and magnetic fields for an ensemble of equally sized cells. With this strategy, we avoid storing the sensitivity matrix or its Hessian and also the repetition of costly mathematical computations. As usual, preconditioners are recommended to reduce the number of iterations needed to solve the inversion problem.

To incorporate the convolution approach described in section 2.1, we modify the algorithm and illustrate the basic change using g _z as an example. The two main changes to incorporate the discrete convolution in the CG search are given as follows: 1) when computing the model response g = A m and 2) when updating the search direction r = b − A T C d d − 1 A . They are computed as: g l m = ∑ k ∑ j ∑ i m 0 i j k W k x l − x i , y m − y j g z , (16) and r α , β , γ = b α , β , γ − ∑ l ∑ m g l m σ l m 2 ⋅ W γ α − l , β − m g z . (17)

We can see the full adaption of the algorithm in Table 1.

M-preconditioned convolution-based.

conjugate gradient algorithm to solve

A T C d d − 1 A m ̂ = A T C d d − 1 d .

We performed a rather conventional singular value decomposition (SVD) resolution analysis using the complete gravity tensor information to characterize exclusively the data sensitivity matrix. In this theory, a general rectangular m × n matrix A is factorized into two orthogonal vector basis U _m×n and V _n×n in the form A = U S V T , (18) where S _m×n is a partially diagonal matrix. With the aim of comparing the lateral and depth resolution of g _z and the six elements of GGT, the matrix A is analyzed for two cases: a) using g _z data along with the six elements of the gravity gradient tensor (GGT) and b) using only g _z .

Our synthetic model consists of a 2000, ×, 2000, ×, 1000 m volume discretized with 21 × 21 × 20 equally spaced rectangular prisms. Two density heterogeneities are formed with 3 × 3 × 3 elemental prims as shown in Figure 3. For the sole purpose of this sensitivity analysis, the A matrix is reassembled column by column using W filters as needed for each one of the 8,820 individual prisms that comprise the model.

In this section, we only discuss the resolution analysis results in terms of the singular values S, leaving the singular vector (U and V) inspection in the appendix for the interested reader. As in conventional plots, the singular values are ordered from the highest to the lowest (Lanczos, 1996) and plotted in a log–log scale. Figure 4 shows the singular values normalized by the highest value for both g _z plus GGT data and g _z data alone cases.

For this plot, we may observe the following:

• Given both sets of singular values are normalized, it is clear that singular values when using only g _z decay faster than when adding GGT data. This accounts for the effective information supplied by the various forms of potential field data.

We use the synthetic model shown in Figure 3 to test the 3D inversion of synthetic gravimetric and magnetic data through our convolution-based CG method. Using this test model, we created synthetic data for both the vertical component of gravity g _z and the GGT and added normal random noise of 5% of the maximum value of the corresponding data type. This resulted in σ _dd = 0.01 mGal for g _z and σ _dd = 0.1 mGal/cm for all the components of the GGT data. Additionally, we start the inversion program with a density contrast of m ₀ = 0.0 g/cm ³, a standard deviation of the model σ ₀₀ = 0.01 g/cm ³, and a regularization parameter of α _p = 10.

The inversion process was stopped after 100 iterations, and the resulting model is shown in Figure 5. We note that when inverting only g _z data, the recovered density contrast values range between − 0.204 and 0.605 g/cm ³, and the location of the largest positive contrast corresponds exactly to that of the heterogeneities of the original test model. The negative values accommodate around the recovered positive density contrast and may reflect the smoothed response when facing the original sharp contrast existent in the test model. As expected from the resolution analysis, the models bear a large smearing effect at depth when data resolution decreases.

We also note that when inverting GGT and g _z data together, the recovered density contrast values range between − 0.243 and 0.817 g/cm ³, thus achieving a closer match to the values of the original test model heterogeneities. The distribution of these positive heterogeneities also depicts better the original boundaries of the test blocks and shows a reduced relative smearing at depth. In all these scenarios, the algorithm proves successful at retaining the expected resolution power in iterative schemes.

The Acoculco caldera (AC) belongs to the Trans-Mexican Volcanic Belt (TMVB), the largest Neogene volcanic arc in North America with a length of almost 1000 km between 18°30′ and 21°30′ N in central Mexico, where volcanic activity is reported to have started about 16 Ma ago (Ferrari et al., 2012; Calcagno et al., 2018). The AC is located approximately 35 km southeast of the city of Pachuca, in the Mexican state of Hidalgo, in the eastern part of the TMVB. This caldera locates between the coordinates UTM 14 570 000 & 610,000 E and 2,190,000 & 2,220,000 N (López-Hernández et al., 2009); it has a semicircle shape and covers an area of ca. 40 × 30 km (Figure 6).

The AC is interesting for geothermal potential because of the presence of an extensive surface hydrothermal alteration, cold acid springs, and gas dischargers. Since 1981, the area has been analyzed by the Comision Federal de Electricidad (CFE) of Mexico that has drilled an exploratory well near Los Azufres with a depth of 2000 m in 1994, where stabilized temperatures rise above 300°C. The main objective of the well was to determine whether high temperatures and permeabilities existed at depth (López-Hernández et al., 2009). A second exploratory well (EAC-2) was drilled in 2008 with a depth of 1900 m and a maximum temperature of 264°C. None of the wells produced fluids; so, the zone is currently considered as a prospect to develop an enhanced geothermal system (EGS).

The caldera complex sits at the intersection of two regional fault systems with NE–SW and NW–SE orientations. A NE–SW alignment of volcanic cones and medium-sized composite volcanoes can be observed. These volcanoes could be related to the NE-striking Apan–Piedras Encimadas Lineament (López-Hernández et al., 2009; Calcagno et al., 2018). The NW–SE-trending fault system is represented by subtle morphological lineaments between the Pachuca and Apan regions, NW of the Acoculco zone (López-Hernández et al., 2009). The AC was built atop Cretaceous limestones, the Zacatán basaltic plateau of unknown age, early Miocene domes ( ∼ 12.7 − 10.98 M a ), and Pliocene ( ∼ 3.9 , − , 3 M a ) lava domes (Avellán et al., 2020).

The oldest outcropping rocks are located in the eastern area and correspond to Cretaceous sedimentary rocks from the Sierra Madre Oriental. These rocks cannot be found in other zones inside the region; nevertheless, various exploratory wells in the AC show their existence at depth. These rocks are affected for various calco-alkaline events, which occurred in most of the area, resulting in several sequences of volcaniclastic deposits, lava flows, and intrusives. The existing igneous materials vary both in their composition (from basaltic to rhyolitic) and their structural arrangement, making the area a difficult target of exploration (López-Hernández et al., 2009).

According to López-Hernández et al. (2009), the stratigraphic sequence of the Tulancingo–Acoculco complex can be described from the exploratory well EAC-1 with the following lithology (from bottom to top): 340 m of an intrusive body who is responsible for metamorphism of the boxing rock, 870 m of an intensely metamorphosed sedimentary sequence (skarn), and 790 m of a volcanic sequence related with the activity of the complex.

The analyzed area inside the AC is located between the coordinates UTM 14 582 000 & 602,000 E and 2,198,000 & 2,210,000 N covering an area of 20 × 12 km; the gravity data were collected in stations located through existing roads and highways in the area, to make a network of stations distributed in the form of polygons with a mean diameter of 5 km and measurement points separated every 250 m (Figure 7). The measurements were made using Worden Master (Texas Instruments) gravity meters with an accuracy of tenths of μGal (López-Hernández et al., 2009).

The Bouguer anomaly data were re-sampled on a grid spaced 200 m and processed to obtain a regular grid for analysis. A linear regional trend using the values at the ends of the map was removed to capture the gravimetric effect of the overlying material to basement. The gravity values were shifted to produce a residual negative gravity anomaly suitable for inversion (Figure 7). We see that the minimum gravity value locates near the center of the studied map and forms a semicircular structure located north of the borehole EAC-1. This negative gravity anomaly extends, with lower intensities, in the preferential regional directions (NW–SE and NE–SW). The maximum gravity values are located to the NE and NW of the studied area denoting a decreased vulcanosedimentary thickness.

The Acoculco density model was composed of a 20 × 12 km volume using 101 × 61 × 100 equally sized (200m × 200m × 30m) rectangular prisms. We assigned a standard deviation for the observed data of σ _dd = 0.01 mGal.

Our program starts with an initial model m _ini = 0.0 g/cm ³ and was stopped when the rms reached a value less than 2% with respect to the maximum gravity value; the average time of the iterations was of 458.35 s compared with the 1,250.1 s of the computation time of the matrix A. The model response is displayed in Figure 7.

The estimated density contrast model is shown in Figure 8. For the deepest model slice (2,800 m), it becomes clear that the high-density contrast material surrounds (E, NE, NW, and SW) the less dense material to the north of borehole EAC-1, which goes in agreement with the overall structural caldera. There are, however, clear preferential directions of the less dense materials (NW–SE, NE–SW, E–W, and N–S), which agrees with the known regional tectonic scenario. These features also imprint in the shallower depth slice shown in Figure 8B, which indicates the action of the normal smearing of the model at depth. This depth smearing is confirmed by the selected vertical sections included in Figure 8.

The shallowest model slice (Figure 8A) is almost completely dominated by low-density contrast materials that more likely belong to the upper vulcanosedimentary sequence that fill in the structural crater left by the various calderic episodes.

Figures 8D,E are vertical sections, which mainly denote the varied thicknesses of the vulcanosedimentary cover.

As in the case of the gravimetric data, the area of the AC for the magnetic data is between the coordinates UTM 14 582 000 & 602,000 E and 2,198,000 & 2,210,000 N. The aeromagnetic total magnetic field (TMI) data were provided by the Servicio Geologico Mexicano for the Acoculco area. The distance among flight lines was 1 km with direction N–S; the E–W control lines are 5 km separated, and the flight height was 300 m above ground level. The data were processed with MagMap tools in Oasis Montaj^® to produce the 200 m-spaced reduced-to-the-pole (RTP) anomaly map used for inversion (Figure 9).

To perform the RTP data inversion from Acoculco and maintain a magnetization contrast model consistent to the previous density model for Acoculco, we used a magnetization contrast model with exactly the same dimensions and spacing to that used for the density contrast model. We then assumed a standard deviation to the observed RTP data of 0.1 nT and started the process with an initial model of m _ini = 0.0 A/m. The program stopped at the 100th iteration taking an average time per iteration of 450.98 s. The model reached an acceptable normalized data misfit (rms = 3.34), and their data response is shown in Figure 9B.

Figure 10 shows various slices of the resulting magnetization contrast model. At the deepest horizontal slice, at 2,800 m depth below ground, we can observe several major magnetization contrast zones with trends in different directions. The largest zone with minimum magnetization contrast intensities orients NW–SE. This zone is flanked east by a narrower zone with large magnetization contrast. The borehole EAC-1 locates exactly between two of these positive magnetization contrast zones.

At 2000 m depth (Figure 10B), the largest negatively magnetized heterogeneity trending NE–SW and the positive magnetization contrast volume shown at the largest depth remain but are narrower; in the north-east corner, a positive magnetization contrast region with trend NW–SE is found. All these anomalous volumes are subdivided in smaller regions depicting an increased geological heterogeneity. At 1,500 m depth (Figure 10A), we observe several isolated volumes with local positive or negative magnetization contrast. The borehole EAC-1 is flanked NE–SW by two local negative magnetization contrast regions.

The EW section (Figure 10D) evidences the continuity at depth of the positive (west) and negative (east) magnetization contrast heterogeneities. These regions are overlain by several local magnetization contrast areas. Differently, the S–N section (Figure 10E) shows two deeper high-magnetization contrast volumes flanking the low-magnetization contrast heterogeneity around the borehole EAC-1.