3D Convolution Conjugate Gradient Inversion of Potential Fields in Acoculco Geothermal Prospect, Mexico

Potential field data have long been used in geophysical exploration for archeological, mineral, and reservoir targets. For all these targets, the increased search of highly detailed three-dimensional subsurface volumes has also promoted the recollection of high-density contrast data sets. While there are several approaches to handle these large-scale inverse problems, most of them rely on either the extensive use of high-performance computing architectures or data-model compression strategies that may sacrifice some level of model resolution. We posit that the superposition and convolutional properties of the potential fields can be easily used to compress the information needed for data inversion and also to reduce significantly redundant mathematical computations. For this, we developed a convolution-based conjugate gradient 3D inversion algorithm for the most common types of potential field data. We demonstrate the performance of the algorithm using a resolution test and a synthetic experiment. We then apply our algorithm to gravity and magnetic data for a geothermal prospect in the Acoculco caldera in Mexico. The resulting three-dimensional model meaningfully determined the distribution of the existent volcanic infill in the caldera as well as the interrelation of various intrusions in the basement of the area. We propose that these intrusive bodies play an important role either as a low-permeability host of the heated fluid or as the heat source for the potential development of an enhanced geothermal system.


INTRODUCTION
Potential field data such as gravity and magnetics are among the first geophysical data used in mineral and hydrocarbon exploration. Their continued use has resulted in a historical improvement on data acquisition and interpretation methodologies as well as in the development of surveying instruments. One example is the development of modern airborne gravity gradiometers (Zhdanov et al., 2004;Nabighian et al., 2005;Dransfield and Zeng, 2009;Jekeli, 2006) and the increased use of unmanned aerial vehicles for aeromagnetic surveys (e.g., Aleshin et al., 2020;Jiang et al., 2020;Parshin et al., 2020;Walter et al., 2020). These technological developments have yielded the assimilation of the larger potential field data sets needed to achieve higher subsurface detail for reservoir, mineral, and archaeological studies.
Most algorithms for the numerical computation of the potential fields due to the arbitrarily shaped volumes are based on analytical responses of defined volumes such as rectangular prisms (Sorokin, 1951). Some examples of these algorithms are the 2D inversion linear programming developments by Mottl and Mottlová (1972), Cuer and Bayer (1980), and Safon et al. (1997). A complete review of existing algorithms for the computation of gravity and gravity gradient effects due to some geometric bodies such as rectangular prisms and polyhedrae can be found in Li and Chouteau (1998).
In order to carry out a large-scale data inversion, we need a fine three-dimensional discretization of the subsoil with a large number of parameters; thus, we need an efficient algorithm. Whereas exact analytical expressions allow the computation of arbitrarily shaped models, the reiterated computation of trigonometric and logarithm functions bears some computational cost. To face this challenge, several approaches have been developed. They include the extensive use of computational resources using parallel computing (e.g., Moorkamp et al., 2010) as well as various data-model compression strategies such as fast Fourier transform (e.g., Pilkington, 1997;Shin et al., 2006;Caratori Tontini et al., 2009), wavelet transform for sensitivity matrix compression (e.g., Li and Oldenburg, 2003;Davis and Li, 2011;Martin et al., 2013;Sun et al., 2018), or mesh refinement (e.g., Ascher and Haber, 2001).
While we concede that all the computational approaches described above have a large history of success when applied to the inversion of potential field data for large data sets, we consider that the especial properties of potential field data can be used further to reduce the industriousness of potential field data inversion problems. We note that potential fields follow the same principles: they are conservative (i.e., result on harmonic fields that can be described by a scalar field) and also depend on the relative position between the source and the measurement point. Both features are fundamental for Fourier transform processing and inversion approaches that have long been in place (e.g., Parker and Huestis, 1974). They have also proven advantageous for multiple tools for data processing and analysis (cf. Blakely, 1996).
The convolution property of the potential fields have been used by Caratori Tontini et al. (2009) to get a simple expression in the Fourier domain and perform a 3D forward modeling for both gravity and magnetic anomalies of a given distribution of density or magnetization contrasts and provide a faster tool for modeling anomalies; similarly, Chen and Liu (2019) express the gravity field like a convolution integral and introduce an optimized algorithm using the FFT to compute the gravity response along a plane; nevertheless, none of these algorithms is applied to the inversion of gravity and magnetic data.
In this work, we propose that given the regular accommodation of both large spatial data grids and discretized three-dimensional volumes, we can take the advantage of the superposition principle inherent to potential fields as well as the convolution-based property of their associated integral equations to establish a general framework for an exact sensitivity matrix compression useful for an efficient 3D inversion of potential field data such as gravity, magnetics, and gravity gradient data. We first show the theoretical foundations and the resolution power through conventional singular value decomposition. We then propose a computational framework for a convolution-based conjugate gradient 3D inversion algorithm for potential field data and prove the algorithm for various potential field data combined in a test example. We apply our algorithm to gravity and magnetic data from the Acoculco geothermal zone in Mexico.

Computation of Potential Fields for 3D Volumes of Rectangular Prisms
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 dU(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 where c 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 (dg) is FIGURE 1 | Illustration of the contribution to the potential U(r) associated with the three-dimensional distribution of monopolar (ρ(r′)) or dipolar (M(r′)) sources.
conservative, irrotational, and according to the Helmholtz theorem, it can be described by the gradient of the differential gravity potential (e.g., Blakely, 1996), For a finite volume V′, the total potential can be described in the form of a volume integral (e.g., Zhdanov et al., 2004) 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: where κ 4π × 10 −7 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 for the gravity potential and for the magnetic potential, where ·p denotes 3D convolution (p) 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: where dxdydz. Using a discrete 2D convolution (*) for a Δx, Δy equally spaced grid: With this notation, we may compute g z as: where W gz k zW D k zz , i.e., (following Banerjee and Das Gupta, 1977) FIGURE 2 | Illustration of convolution filters (W k ) for the g z , T xx , and T xy elements when acting on two shifted rectangular prisms located on the k − th layer of the model. Note that the anomalies result on a simple displacement in space.
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 where W T k ∇∇W T k (α, β, c), i.e., We may also use the same function (W T k ) for the magnetic field as: 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 twodimensional domain for each layer of cells ( Figure 2), which can also be easily stored to avoid repeated computations when applying iterative procedures for inversion.

Convolution-Based Conjugate Gradient Inversion
In order to solve the inverse problem, we use a quadratic norm to define the following objective function: where d obs are the observed data, A is the sensitivity matrix, m accommodates all the model parameters (in our case either Initialization compute filters W k , and propose initial m m 0 , FIGURE 3 | Three-dimensional configuration of the test density model. The model is composed of two anomalous rectangular prisms embedded in a homogenous and finite three-dimensional volume. Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 759824 density or magnetization values per cell), C dd is the covariance matrix of the observed data d obs (assumed diagonal), C 00 is the covariance matrix of the a priori model m 0 , 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.  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 759824 5 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 optimalm estimator given by: 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 Am) and 2) when updating the search direction r b − A T C −1 dd A. They are computed as: and We can see the full adaption of the algorithm in Table 1. M-preconditioned convolution-based. conjugate gradient algorithm to solve

Test Model
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 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  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 759824 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.

SVD Resolution Analysis
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. • For both cases, the decrease of the singular values becomes significant at discrete numbers. In this example, roughly every 441 element, which is the number of cells for each individual depth slice, indicating that resolution at depth decays faster when using a map of any potential field data.
• The magnitude of the sharp decreases at every 441 element is at least one order of magnitude: by example, the S 1 singular value for the sensitivity matrix of g z and GGT together is 11.50 times greater than that of S 442 and 140.32 times greater than that of S 923 ; as comparison, for the first singular value of g z alone, (S gz1 ) is 136.19 times greater than S gz442 and 5,571.42 times greater than S gz923 . This reflects a rapid loss of resolution when descending through each layer of the model.

Inversion Experiment
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.0 g/cm 3 , a standard deviation of the model σ 00 0.01 g/cm 3 , 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 3 , 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 Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 759824 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 3 , 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.

Geological Framework
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-SEtrending 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 Ma), and Pliocene ( ∼ 3.9, −, 3 Ma) 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.

Gravity Data Inversion
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,   (Figures 8, 10). Note that the distinction between the proposed intrusives is mainly based on their differences in magnetization contrast.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 759824 9 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 3 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.

Magnetic Data Inversion
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 Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 759824 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.

INTEGRATED INTERPRETATION
The resulting 3D distributions of density contrast ( Figure 8) and magnetization contrast ( Figure 10) are somewhat limited in resolution. This is a natural consequence of the divergent decay of both potential fields reflected not only at depth but also laterally due to the separation of the original aeromagnetic flight lines and the sparsity of the actual gravity stations. Nevertheless, there are several geological features that can be confidently inferred from the interpreted models. From the point of view of the density contrast model of Figure 8, the clearest feature to identify corresponds to the distribution of the volcanosedimentary infill in the area, which is delimited by the shallow largest negative density contrast. Differently, the magnetization contrast distribution depicts various located units of both magnetic and non-magnetic materials, most of them inferred at basement depths that may be originated by various intrusive events.
In general, we may identify at least four units, which clearly match the reported lithological groups in their stratigraphic and structural disposition (see, for e.g., López-Hernández et al., 2009;Avellán et al., 2020). They include (see Figure 11): Unit I. A low-density contrast zone with various magnetization contrast values. These combined values imply unconsolidated materials with a range of mineralogical compositions. These materials correlate with the alluvial and volcaniclastic deposits with a low level of compaction that came to occupy the large caldera structure. Unit II. A zone with high-density contrast and high-positive magnetization contrast. This unit may correlate to various intrusions dominated by ferromagnetic compositions as reported in the late events of the volcanic sequences of the area (cf. Sosa-Ceballos et al., 2018;Avellán et al., 2020). Unit III. A zone with high-density contrast and high-negative magnetization contrast. This zone may correspond to nonmagnetic intrusions. According to the apparent intrusion sequence, these intrusions precede the magnetic intrusions of Unit II, since it is intersected by them at various locations. Unit IV. A high-density contrast and null magnetization contrast zone. This unit matches the position of the cretaceous limestone that conforms to the regional basement in the area. This basement is largely intruded by units II and III.
The structure of the caldera is noticeable in Figure 8A by the lower density contrast of the deposits that came to fill the formed topographic depression. Differently, the basement shows an evident structural control, having dominant nonmagnetic intrusives (Unit III) along a preferable trend NW-SE which runs along the east of the mapped volume. One of these intrusives seems to correspond to the hornblende granite reached by the exploratory well EAC-1 (López-Hernández et al., 2009). The intrusives of the more ferromagnetic composition occur in various events aligned in a preferential SW-NE direction ( Figure 10) and are likely to be associated with the latest events of magmatism in the area that resulted in several basalt-dominant deposits of more recent age (Avellán et al., 2020). Both alignments are clearly associated with more regional lineaments identified in the larger extension of the TMVB (see, for e.g., Avellán et al., 2019).
From a geothermal potential point of view, it seems clear that the area around the exploratory well EAC-1 is largely extended in between Cretaceous limestone and the various intrusives which are naturally highly impermeable rocks. The evidence of high temperature in the well, however, indicates the influence of a nearby heat source. Considering the various potential ages of the intrusive events and the closeness of the intrusive of Unit II to EAC-1, it seems very likely that this material or their younger magmatic feed may be responsible of both fracture motivation and heat transfer as suggested by Sosa-Ceballos et al. (2018). In this scenario, the mapped flanks of this Unit II may be a suitable place to explore temperature potential to evaluate the feasibility of developing an enhanced geothermal system. It is noted that this high temperature evidence only occurs at depth in the well surroundings and does not seem to be correlated to the shallower basin of the caldera filled with volcanic deposits that may conform to upper aquifers for the zone.

CONCLUSION
We have developed a convolution-based conjugate gradient algorithm for the inversion of potential field data to produce three-dimensional volumes of density or magnetization contrasts. The algorithm computes exact convolutional filters that permit the storage of a highly compressed but exact sensitivity matrix. By using regular data grids and models, the convolution prevents not only the repeated computation of costly mathematical instructions but also the use of interpolated values from the filters.
Using singular value decomposition and a synthetic test, we show that the proposed methodology bears neither loss of data information nor model resolution achieving highly detailed inversion models in small computational infrastructures.
We prove the algorithms applicability on the Acoculco geothermal area in Mexico, where we successfully inverted land gravity and aeromagnetic data. From the combined density and magnetization contrasts values and the lithological information provided by an exploratory well, we could distinguish a group of intrusive bodies at depth as potential low-permeability geothermal reservoirs and their interaction with younger intrusive bodies as potential heat sources. We found no apparent connection of the deep basement with the volcanosedimentary cover and thus no direct connection with local aquifers in the area. The combined models yielded a coherent interpretation of a complicated volcanic caldera and helped to elucidate their implications for the development of an actual enhanced geothermal system.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subjected to the following licenses/restrictions: private data for copyright. Requests to access these data sets should be directed to lgallard@cicese.mx.

AUTHOR CONTRIBUTIONS
JC contributed to investigation, conceptualization, methodology, visualization, and writing-original draft. LG contributed to conceptualization, methodology, supervision, resources, and writing-review and editing.