Gravimetric and morpho‑structural analyses in the superhot geothermal system Los Humeros: an example from central Mexico

The influence of deep and regional geological structures is becoming increasingly important in superhot geothermal systems due to their proximity to the transition between brittleness and ductility. In the Los Humeros geothermal field in Mexico, where subsurface fluids reach temperatures of over 350 °C, the surface structures resulting from the collapse of calderas have so far only been interpreted at the local scale. The aim of this work is to place the recent tectonic and volcano‑tectonic geo‑ morphologic evolution and structures in the Los Humeros volcanic area in a regional context. NE‑ and NW‑striking dominant structures resulting from a morpho‑structural analysis on a regional scale are confirmed by negative and positive anomalies, respectively, after Butterworth filtering of gravity field data with different wave‑ lengths over a local area of about 1000 km 2 . By analyzing the slip and dilation trends of the observed directions, we show the relevance of the regional context for res‑ ervoir exploration. The magnitudes of the principal stresses we estimate indicate a trans‑tensional fault regime, a combination of strike‑slip and normal faulting. The structures derived from the gravity and morpho‑structural analyses, which are parallel to the maximum horizontal stress, have the highest potential for tensile and shear fail‑ ure. Therefore, the corresponding negative gravity anomalies could be related to frac‑ ture porosity. Consequently, we hypothesize that these structures near the transition between brittleness and ductility control fluid flow in the Los Humeros geothermal field.


Introduction
In geothermal systems with temperatures above 350 °C, referred to as superhot geothermal systems (SHGS), higher productivity can be achieved with fewer wells than conventional ones (Cladouhos et al. 2018).SHGS have been successfully drilled at depths near or within the brittle-ductile transition zone (BDT) in continental and oceanic crust (Kruszewski and Wittig 2018;Watanabe et al. 2019;Reinsch et al. 2017).To understand SHGS, the base of the permeable crust, which is controlled by fractures and in which geothermal fluids flow, albeit transiently, i.e., the depth of this regional rheological boundary, must be considered.This depends on various mechanical and physical parameters, such as tectonic context, pressure, temperature, strain rate, lithology and fluid pressure (e.g., Kohlstedt et al. 1995).
The Reykjanes geothermal area in Iceland, for example, illustrates how the interaction of mechanical and physical parameters promotes brittle deformation: aseismic zones, such as those in the Reykjanes field, have been associated with low resistivity zones determined by inversion of magnetotelluric (MT) data (Saemundsson and Einarsson 2014;Guðnason et al. 2015Guðnason et al. , 2016)).Similar correlations have also been observed in other regions, e.g., in the Taupo Volcanic Zone in New Zealand and in the Villarrica low temperature field in southern Chile (Bertrand et al. 2012;Bryan et al. 1999;Pavez et al. 2020).
The Los Humeros geothermal field (LHGF), located in the eastern part of the Trans-Mexican Volcanic Belt (TMVB), does not have a clearly defined BDT.The northern part of the current production area has the highest temperatures of up to 380 °C at a depth of 2 km and is, therefore, a SHGS that has not yet been fully exploited.Inversion of MT data (Benediktsdóttir et al. 2019) shows anomalies with low resistivity trending N-S to NNW-SSE to a depth of 8 km.Seismic activity indicates brittle deformation to a depth of at least 3 km, as well as a heterogeneous stress state in different sectors of the field and at different depths (Lorenzo-Pulido 2008;Rodríguez et al. 2012;Lermo et al. 2016;Toledo et al. 2019Toledo et al. , 2020;;Jousset et al. 2019).Assuming a continental crust mechanically controlled by the rheology of quartz and a geothermal gradient of 100 °C/km, Calcagno et al. (2019) interpreted a BDT at about 4-4.5 km below the surface.
Given the relevance of deep structures for SHGS, in this study we complement the existing data from Los Humeros with gravity measurements, which can support integrated subsurface models (Calcagno et al. 2008;Guillen et al. 2008), and also provide valuable information on deep structures.Sequential Butterworth filters applied to the Bouguer anomalies reveal structures at different depth levels (Abdelfettah et al. 2014;Baillieux et al. 2014).Due to the spatial correlation of thermal and gravity anomalies (Baillieux et al. 2013), the latter has also been used to investigate fracture porosity in geothermal reservoirs (Baillieux et al. 2014;Altwegg et al. 2015).As such, gravity data physically characterize geologic lineaments.In the LHGF, the joint inversion shows a striking link between low density, low magnetization, and fault zones possibly related to hydrothermal alteration (Carrillo et al. 2022).
In addition, regionally distributed tectonic lineaments have proven to reflect the accumulation over time of deformation extending into the lower crust (Regenauer-Lieb and Petit 1997).Our morpho-structural approach involves a combination of tectonic geomorphology and structural geology, with the analysis of remote sensing imagery, digital terrain topography, multi-scale maps, local field mapping, and geological information to identify features of tectonic significance.This approach has been applied in active mountain chains (Piccardi et al. 1997(Piccardi et al. , 1999;;Boccaletti et al. 2002Boccaletti et al. , 2001;;Delle Donne et al. 2007;Nirta et al. 2021) and in volcanic areas (Boccaletti et al. 1998;Giaquinta et al. 1999;Pierotti et al. 2017;Piccardi et al. 2017) including the Los Humeros Volcanic Complex (LHVC) (Norini et al. 2019).
To investigate the relevance of the regional context to the reservoir exploration, the slip and dilation tendencies of the structures derived from gravity and morphostructural results are analyzed based on a generic fault model (Morris et al. 1996;Ferrill et al. 1999).Thus, we show which faults are favorably oriented to transmit geothermal fluids (Siler et al. 2016).
So far, the structures resulting from the caldera collapse in the LHGF have only been interpreted on a local scale (Norini et al. 2015;Carrasco-Núñez et al. 2017a).The present work aims to consider the recent tectonic and volcano-tectonic geomorphology of the Los Humeros volcanic area at a regional scale, characterize deep structures by gravity anomalies, and infer their geothermal relevance by analyzing their reactivation potential.In this way, we also contribute to defining of the local and regional BDT depth.

The Trans-Mexican Volcanic Belt
The TMVB is an almost 1000 km-long Neogene continental volcanic arc resulting from the subduction of the Cocos and Rivera plates beneath the North American plate.The evolution of the subduction system, as well as the nature of the crust beneath the belt, mainly control its large spatial and temporal variation, as well as its magma composition, encompassing several ancient fault systems that were partially reactivated at different times (Ferrari et al. 2012).
The main fault systems within the TMBV are described by Ferrari et al. (2012) (Fig. 1).The age of deformation suggests a coherent pattern of southward migration of tectonic activity.However, the fault system between Querétaro and Toluca in the center of the TMVB does not clearly show such as pattern of tectonic activity during the Pleistocene-Holocene. East of this fault system, the Neogene deformation is less intense.
The western part of the TMVB is characterized by the intersection of the E-Wtrending Chapala, the N-S-trending Colima, and the NW-SE-trending Tépic-Zocoalco grabens, southwest of Lake Chapala.The central part is affected by W-E-trending Quaternary normal faults with known historical seismicity, whereas no major active faults have been documented in the eastern part of the TMVB (Suter 1991).
In central Mexico, results from Suter (1991) indicate a normal faulting stress regime with a minor left-lateral strike-slip component.The direction of the

The Los Humeros Volcanic Complex
The LHVC is a Pleistocene basaltic andesite-rhyolite explosive caldera volcano formed by two explosive eruptions encompassed between 460 and 70 ka (Ferriz and Mahood 1984;Carrasco-Nunez et al. 2012), followed by two minor post-caldera eruptive events, dated at 44.8 ± 1.7 and 2.8 ± 0.03 ka, respectively, the latter affecting the southern border of the caldera rim (Carrasco-Núñez et al. 2017b;Lucci et al. 2020).The LHVC holds an active geothermal system, the LHGF, located in the eastern sector of the TMVB (Fig. 1), northeast of the Pico de Orizaba-Cofre de Perote range (Carrasco-Núñez et al. 2017a).The location of the caldera is controlled by the intersection of pre-existing fault systems, NNW-and NE-oriented, respectively, determining an asymmetric caldera, with the deepest part settled in its south-western sector (Maestrelli et al. 2021, with references therein).Within the Los Humeros caldera, the volcano-tectonic structural system is controlled by the smaller Los Potreros caldera, where faults deriving from the caldera collapse are oriented from NNW-SSE to E-W (Norini et al. 2015).A schematic geological map of the LHVC is presented in Fig. 2.
The caldera structure is well captured in the 3-D integrated geological model of Los Humeros proposed by Calcagno et al. (2022) (Fig. 3).This integrated scale of the model is based on the extent of the geophysical surveys performed during the GEMex project.
Holocene volcanism records the recurrent injection of compositionally distinct magma batches uprising from different crustal depths.This indicates the transition from the caldera stage magmatic system dominated by a single, large, and shallow magmatic body to a polybaric post-caldera stage plumbing system made up of a lower crust mafic reservoir feeding smaller magma batches that are vertically distributed (Carrasco-Núñez et al. 2022).The Los Humeros area is interpreted to be affected by the resurgence of the caldera floor (Norini et al. 2015;Toledo et al. 2020), induced by an inferred magmatic intrusion, representing the heat source of the geothermal system, exploited in the central and northern caldera areas (Arellano et al. 2003;Gutierrez-Negrin and Izquierdo-Montalvo 2010;Lelli et al. 2021).
Heterogeneous stress state in different sectors of the LHGF and at maximum depths of 3 km have been inferred (Lorenzo-Pulido 2008;Rodríguez et al. 2012;Lermo et al. 2016;Toledo et al. 2019;Jousset et al. 2019).However, the age of the caldera formation, the post-caldera plumbing system, and the hypothesized magmatic resurgence allow us to assume that the local variation of stress induced by the caldera collapse (Holohan et al. 2013;Corbi et al. 2015;Oliva et al. 2022) has stopped.In this study, we consider that the Los Humeros area is affected by regional tectonic stress as the primary driver, mainly controlling the present geothermal fluids pathways.
The exhumed Miocene geothermal system of Las Minas (shown in Fig. 1 with a green circle) is located 20 km to the northeast of Los Humeros.It is considered as an analog of the LHGF since it exposes rocks belonging to the Los Humeros pre-volcanic basement

Morpho-structural analysis
With the aim of improving our understanding of the main trends in recent tectonic structures and volcanic-tectonic geomorphology, a morpho-structural analysis was carried out at the LHVC at local and regional scales.The analysis was based on remote sensing data such as topographic maps, rectified satellite images from Landsat, Google Earth, and Bing, digital elevation models, and geological maps at different scales.
This analysis was completed by fieldwork to fix pin-points for the remote sensing analysis.The field survey was performed mainly within the Pleistocene caldera area, the main target for geothermal exploration, as well as along other structures outside the caldera.Faults and fractures, as well as evidence of recent faulting affecting the volcanic deposits, were collected to verify some of the inferred structures and to learn about the local geomorphologic response of recent faulting in this volcanic area.Alignments of volcanic structures, especially within monogenic volcanic fields, were also mapped while paying attention to their morpho-structural characteristics.The criteria used to distinguish tectonic faults, although may be considered subjective, are based on tectonic-geomorphological considerations.For instance, a fault exhibiting a high scarp and longer and more accentuated lateral continuity is considered a major fault with respect to one that features subdued scarps and is laterally discontinuous (Piccardi et al. 1999).
As morpho-structures indicative of recent tectonics, we considered those geomorphological lineaments in laterally continuous scars whose features suggest that they could be related to tectonics affecting younger geological or geomorphological elements, such as lava deposits, volcanic features (such as crater rims, alignment of monogenic cones or volcanic centers) or dated paleosurfaces.The simultaneous presence of several pieces of geomorphologic, volcanic-tectonic, or field mapping evidence confirm their tectonic origin.Examples of the methodology described can be found in the Additional file 1.

Gravity data
In the following, we describe the procedures of data acquisition and processing, and the forward modeling of the gravity response of the 3-D geological model (Fig. 3).Section "Filtering analysis" describes the concept of Butterworth filtering.It should be noted that the comparison with the results of the morpho-structural analysis is based on the Butterworth-filtered data.

Survey and processing
Gravity data in the LHGF was acquired during two field surveys (November-December 2017 and April 2018) in an area of 24 km × 23.5 km during the GEMEX project (www.gemex-h2020.eu; Fig. 2) using a CG-5 Autograv with an instrumental accuracy of 0.001 mGal.The positions were measured by differential GPS (DGPS) using two Trimble 5700 receivers and two Trimble antennas with an instrumental vertical accuracy of 5 mm.From the 344 gravity measurements, 263 were acquired along ten E-W profiles of 5.5 km length with typical inter-station and inter-profile distances of 200 and 500 m, respectively.The survey was completed by a NE-SW oriented, 31-km-long profile, including 81 gravity measurements with an inter-station distance of about 375 m.Three measurements were rejected due to unreliable DGPS records.
Standard deviations of the raw data measurements range between 0.008 and 0.119 mGal and between 0.0012 and 0.3893 m for the gravity (Fig. 4a) and the vertical component of the DGPS data (Fig. 4b), respectively.The application of thresholds of 0.05 mGal and 0.2 m results in a final number of 318 stations for further processing, i.e., 92.4% of the acquired data.
Instrumental drift, earth tide, latitude, and elevation corrections were carried out using GravProcess Matlab codes (Cattin et al. 2015).To estimate the instrumental drift, repeated gravity measurements were carried out at the beginning and end of each day of measurements at a base station (red star in Fig. 2).The instrumental drift correction was estimated by a least-squares fit of the weighted time series at this station with a first-order polynomial function.The earth tide correction was calculated for each gravity measurement according to Agnew (2007Agnew ( , 2012)).Free-air and simple Bouguer corrections were made using the DGPS-processed data.
A complete Bouguer anomaly (Kane 1962;Nagy 1966) was obtained using a Bouguer density of 2670 kg m −3 .The terrain correction was carried out using the Oasis Montaj ™ software (Geosoft) up to a maximum distance of 166.7 km from each gravity station (LaFehr 1991).Since the gravitational effect of a prism decreases with the square of the distance, digital elevation models (DEM) were used with the following resolutions, which decreases with increasing distance from the gravity stations: Note that for the 14 north-easternmost stations the 1 m DEM was not available and thus the 15 m resolution DEM was the highest resolution used in the terrain correction.
Given the standard deviation threshold of the GPS measurements, the uncertainty introduced by the terrain correction is estimated to be below 0.06 mGal.The uncertainty in elevation introduced by the DEM is approximately 0.6 mGal, summing up with the standard deviation threshold of the gravity measurements (0.05 mGal) to a mean uncertainty of about 0.7 mGal.
The complete Bouguer anomaly in the LHGF (Fig. 5) was complemented with a regional complete Bouguer anomaly with a reference density of 2670 kg m −3 including 13,737 measurements (courtesy of Comisión Nacional de Hidrocarburos).The regional trend represents the increase in crustal thickness from the Gulf of Mexico towards the center of the continent (Urrutia-Fucugauchi and Flores-Ruiz 1996).From this comprehensive data set, the residual anomaly was determined by subtracting the regional trend (Fig. 5b) in the area of the geological model.The latter was determined using a Gaussian filter with a cutoff wavelength of 20.5 km.The respective radially averaged power spectrum of the complete Bouguer anomaly is included in Additional file 1: Fig. S6.
Due to the uniform geology below 0 m.a.s.l., where only basement occurs, we limit the forward modeling to this depth.The geological model was discretized into 3-D voxels with a minimum cell size of 210 m in the X and Y directions, resulting in a horizontal cell distribution of nx = 133 and ny = 104.In the vertical dimension, a variable cell size between 131 and 151 m with a geometric factor of 1.1 was used, resulting in a number of vertical cells of nz = 32.Finally, each voxel of the geological model is associated with a stratigraphic unit and its respective density value.

Filtering analysis
To place the filter analyses in a regional context, a subset of 1083 regional gravity field stations was added to the local survey (see outer black rectangle in Fig. 5).The spatial data of the full Bouguer anomaly were transformed to the wavenumber domain by a fast Fourier transform and Butterworth filters (Butterworth 1930) were defined and applied using Oasis Montaj ™ software (Geosoft).In this context, the spatial data of the study area are extended to ensure continuity and periodicity.It is then interpolated to make it periodic along both coordinates.For reference and analysis purposes, the radial average spectrum is calculated, and high-pass and band-pass filters are defined and applied.The full wavelength ranges of the high-pass and band-pass Butterworth filters used are listed in Table 2.As a first approximation, it is assumed that a longer wavelength represents deeper structures.

Slip and dilation tendency analysis
Generic slip and dilation tendencies (Morris et al. 1996) were calculated and compared to the observed orientations of the residual gravity anomalies and the structures identified during the morpho-tectonic analysis.
The regional stress in Los Humeros and surrounding areas is defined by a NE-trending maximum horizontal stress with an azimuth, θ S Hmax , of 42°E ± 8° (Heidbach et al. 2016) and a dominant NNW-trending crustal extension direction (García-Palomo et al. 2018).This stress field interacts with the crustal uplift indicated by the high heat flow in the TMVB (Prol-Ledesma et al. 2018;Prol-Ledesma and Morán-Zenteno 2019;Espinoza-Ojeda et al. 2023).Seismic tomography (Toledo et al. 2020) suggests a caldera resurgence in the Los Humeros area, which implies that uplift is dominant with respect to crustal extension today.
The stresses are therefore calculated in two ways: (i) inside the caldera, where the geothermal fluid pressure is active (pore fluid factor = 0.9 ; cfr.Liotta and Ranalli 1999), and (ii) outside the caldera, where the fluid pressure is assumed to be hydrostatic ( = 0.4 ).
In this study, we estimate the horizontal effective principal stresses based on linear elasticity assumptions (Thiercelin and Plumb 1994) where ν is the rock's static Poisson's ratio and E is the rock's static Young's modulus, both taken from Weydt et al. (2021).Total vertical principal stress, S v , and pore pressure, P p , gradients to calculate the effective vertical stress (σ v ) were obtained from Kruszewski et al. (2022).
The minimum, ε hmin , and maximum tectonic strains, ε Hmax , were determined on the basis of a 1-D geomechanical model using a strain-corrected method (Thiercelin and Plumb 1994;Blanton and Olson 1999) and neglecting thermal effects.The geomechanical model developed was based on the results of the H 64 well (red circle in Fig. 2), including the lithologic column and elastic rock properties, and was calibrated using the stress values obtained from a leakage test in the same well (Kruszewski et al. 2022).The resulting tectonic strain amounted to ε hmin of −4.0 ± 9.0 • 10 −5 and ε Hmax of 127 ± 34.1 • 10 −5 .
We estimate the effective principal stresses to the top of the BDT as this is the base of crustal permeability controlled by fractures and where the flow of geothermal fluids is channeled.However, this is an approximation for a gradual transition that likely extends over several kilometers.By equating the critical stress for frictional failure and the creep stress in the ductile state, the BDT depth can be estimated as follows (Ranalli 1995(Ranalli , 1997;;Liotta and Ranalli 1999): where α is a parameter that depends on the tectonic regime and the frictional properties of faults (here,α = 0.68 , assuming a coefficient of friction of 0.6 and normal faulting), ρ is the density of overlying rocks, g the acceleration of gravity, z the depth, λ the pore fluid factor (a ratio of pore fluid pressure to the lithostatic stress), ϵ the ductile strain rate, R the gas constant, T the absolute temperature, and A, n, and E the creep parameters of the rock (see Table 3).
Frictional parameters were chosen assuming a normal faulting regime (Sibson 1974) and a coefficient of friction of 0.6 (Byerlee 1978).Creep parameters were chosen assuming a quartz-dominated crust (Calcagno et al. 2019).Different pore fluid factors were considered in the analysis, from 0.4 (low pore fluid pressure) to 0.9 (high pore fluid pressure).In the ductile regime, we considered tectonic strain rates between 10 -16 s −1 and 10 -14 s −1 as references.
Regarding the temperature gradient considerations, more than 3800 new geothermal gradients were compiled by Espinoza-Ojeda et al. ( 2023) from deep boreholes.In the Los Humeros area, an average geothermal gradient of about 120 °C km −1 was estimated (data can be downloaded from: https:// github.com/ omesp inoza ojeda/ Mexico_ HF_ datab ase.git), and in the whole TMVB province, a geothermal gradient above 80 °C km −1 is expected in most areas (Prol-Ledesma and Morán-Zenteno 2019).The BDT depth was then calculated taking into account two geothermal gradients, one representing the geothermal area and the other the regional environment.
As slip and dilation tendencies strongly depend on the azimuth of the fault plane, the analysis was performed for a set of generic faults (Meixner et al. 2016) with fault plane azimuths from 0 to 180° and mean dip angles of 70° (Additional file 1: Table S1) under the estimated stress field.

Results
In the following, the lineaments obtained from the morpho-structural analysis is presented to show the predominant tectonic structures and their orientation.In Section "From the complete Bouguer to the residual anomaly", we present the residual anomaly after subtraction of the regional trend.We compare this residual anomaly with the gravity forward model of the integrated 3-D geologic model (Section "Forward modeling results compared to the residual anomalies") to demonstrate significant agreement within the caldera.Finally, in Section "High-and band-pass Butterworth filter results", we present sequential Butterworth filtering to decipher structures outside the caldera from the gravity field data.The results from the gravity field data point in sub-parallel directions to the morpho-tectonic results, and are thus used for the petrophysical characterization of the same.

Table 3 Frictional and creep parameters used in the BDT calculation
Assuming a continental crust that is mechanically governed by quartz's rheological behavior, parameters for wet quartzcontrolled rheology were chosen (Ranalli 1995(Ranalli , 1997) )

Morpho-structural analysis results
The regional structures of the LHVC, which were derived from a morpho-structural analysis of the eastern sector of the TMBV, are shown in Fig. 6.A detailed description of the recognized morpho-structures can be found in Additional file 1: Figs.S1-S5.The following structural trends were identified: 1. NE-SW oriented structures visible in the Las Minas valley and Los Humeros (e.g., the Perote fault).In the volcanic complex, NE-SW oriented major faults such as the Hidalgo Fault controls the southeastern side of the Los Humeros caldera (Fig. 6).
Other similarly oriented structures are present inside and outside the volcano.2. NW-SE oriented structures, visible in the southwest of the Cofre de Perote volcano, delimit the extent of the Perote fault.In the volcanic complex, the NW-SE oriented structures are the Maxtaloya Fault, one of the main volcanic-tectonic structures controlling the LHGF and the Arenas Fault.A relevant buried volcanic-tectonic structure northeast of the Los Potreros caldera also has this orientation.3. The most significant of the N-S oriented smaller structures are located in the western part of the Los Humeros caldera in the area of Oyameles and in the eastern part of the Los Potreros caldera.In addition, N-S oriented structures appear to control the emplacement of the N-S elongated lava dome of C. Arenas and the Xalapazco crater south of Los Humeros.4. East-west oriented structures: These are manifested at the surface by minor faulting.
This trend appears to have some control over the northern margin of the Los Humeros volcanic edifice.Compared to the NW-SE, the NE-SW oriented structures may be older, as the former appear to control the regional tectonic evolution of the area.The Perote and Las Minas faults, as well as the subsided basin in the southwest of the Las Minas valley are examples.The N-S oriented structures generally have a limited extent and abut the major tectonic faults.They are therefore younger and determined by the latter.This applies, for example, to the lava dome of C. Arenas, the Xalapazco crater and the western structures of the Los Humeros caldera, which abut the NW-SE oriented faults.E-W oriented minor faults could be controlled by subsidence; however, subsidence due to soil compaction is not the subject of this paper as it requires further investigation for proper evaluation.Faulting in the central area of the LHVC is more likely to reflect local subsidence due to magmatic fluctuations or other effects related to the geothermal system.The fault kinematic indicators show that the Maxtaloya/Humeros fault and other major faults are mostly normal with an average dip of 70°.An overview of the kinematic indicators observed on the main fault structures in the Los Potreros caldera is provided in Additional file 1: Fig. S7 and Table S1.

From the complete Bouguer to the residual anomaly
The regional trend of increasing gravity values from SW to NE in the complete Bouguer anomaly (Fig. 7a) is eliminated in the residual anomaly of the area (Fig. 7b).In the latter, lowest residual gravity values are observed in the North-Eastern part of the LHGF and the Los Humeros caldera (Fig. 7b).Medium residual anomalies are observed to the SW of the caldera, trending NW-SE and SW-NE.

Forward modeling results compared to the residual anomalies
The forward model of gravity (Fig. 8) roughly reproduces the low anomaly in the central part of the Los Humeros caldera, but there is no further correlation between the calculated and observed gravity.This is especially true for the southwest of the forward modeled area.In contrast to a relatively high anomaly in the forward modeling, the residual anomaly map (Fig. 7b) shows a low anomaly that crosses the caldera and extends in the NE-SW direction.To the southwest of the central region of the caldera, the residual anomaly shows intermediate anomaly values in the NW-SE direction, in contrast to the low anomaly observed in forward modeling (Fig. 8).Less pronounced alignments are observed in the NE-SW direction (Fig. 8).
The observed discrepancies are further investigated by sequential Butterworth filtering to understand the influence of potentially deep structures on the gravity data.

High-and band-pass Butterworth filter results
Considering the accuracy of our gravity measurements of 0.7 mGal, the 3-km high-pass filter was not included in the analysis.High-pass and band-pass filters with increasing wavelengths from 6 to 20 km and 6 to 30 km, respectively, applied to the complete Bouguer anomaly are shown in Fig. 9a-c and d, e.To investigate the potentially deep structures, we show a wavelength band of 10-30 km (Fig. 9f ).
Randomly distributed low-and high-gravity anomalies are observed in the 6-km high-pass filter (Fig. 9a).Although the dynamics of the 6-km filters extend over several mGal and the size of the anomalies is confirmed in the range of the data collected in this study (inner black rectangle in Fig. 5), it is not interpreted further.At a corner wavelength of 10 km, the previously randomly distributed anomalies begin to align (Fig. 9b).A group of negative anomalies surrounded by positive anomalies forms in the center of the Los Humeros caldera.In the southwest and northeast of the study area, the negative anomalies are oriented in a NE-SW direction and the positive anomalies are preferentially oriented in a NW-SE direction at a higher corner wavelength of 20 km (Fig. 9c).The anomalies observed at this filter level seem to correspond to the anomalies observed after subtracting the regional trend determined with the Gaussian filter (Fig. 7b).This indicates that the regional structures are represented by wavelengths ).Red and black lines represent major faults and the boundaries of the caldera, respectively (Calcagno et al. 2022) longer than 20 km.The annular alignment appears to open up at this corner wavelength and turn into NE-SW and NW-SE alignments.Thus, a large low-gravity anomaly runs NE-SW across the entire caldera and beyond (Fig. 9c).However, it is intersected to the south by a NW-SE-striking high-gravity anomaly that coincides with major NW-SE-striking surface faults associated with the LHGF.
Band-pass filtering confirms the NE-SW directed positive and negative anomalies as well as the transformation of the ring-shaped structure into NW-SE-trending high anomalies and NE-SW low anomalies.Note that the NE-SW and NW-SE alignments persist through the different wavelength ranges and are assumed to represent the dominant structural setting.

Slip and dilation tendency analysis results
In the following, we first present the assumptions we have made for a generic estimate of the slip and dilation tendency analysis.It should be noted that neither the gravity nor the morpho-structure analyses provide information on the depth of the corresponding geological or tectonic structures.However, it must be assumed that they occur above or at most in the BDT.The rheological profiles of the TMVB and LHGF areas for estimating the depth of the BDT are shown in Fig. 10.Assuming a high pore fluid pressure in the LHGF, the depth of the BDT can be estimated to be between 2.5 Fig. 9 Residual anomalies obtained from high-pass and band-pass filtering with corner wavelengths of 6, 10, 20, and 6-10, 6-30, and 10-30 km, respectively.Red and black lines represent major faults and the boundaries of the caldera, respectively (Calcagno et al. 2022).Black crosses represent the gravity measurement locations and 3.3 km as a first approximation (Fig. 10a); in comparison, the depth of the BDT in the TMVB would be between 2.9 and 4.6 km for low and high pore fluid pressure, respectively (Fig. 10b).(Toledo et al. 2020), and b density variation with depth from the regional density model in Los Humeros (Carrillo et al. 2022) To further substantiate our estimate, we refer to the results of seismic events recorded over a 12-month period in Los Humeros (Toledo et al. 2020).The peak of the event distribution is at a depth of about 3 km below ground level (Fig. 11a) and then decreases to zero at a depth of about 5-6 km.Similarly, the density values from the 3-D inversion (Carrillo et al. 2022) show a decrease from a depth of around 3 km (Fig. 11b).Lithologically, this decrease occurs within a basement that can be assumed to be homogeneous.The density decrease can thus be explained by decreasing porosity as an alternative to lithological changes.Experimental results support this assumption (Violay et al. 2017).
These observations, together with the calculations of the rheological profiles (Fig. 10), indicate that the BDT zone in the LHGF begins at a depth of about 3 km below ground level.The stress estimation results show a trans-tensional faulting regime from strikeslip to normal at 3.2 km depth and at a P p of 23.3 MPa (Additional file 1: Fig. S8).It implies a regional trans-tensional deformation to be partitioned in the main fault trends affecting the areas.The magnitudes of the total and effective principal stresses are shown in Table 4.
Since the transition depth of 3.2 km is estimated as the upper limit of the BDT and deeper levels indicate a normal faulting regime, the generic slip and dilation tendency  analysis were performed for such a regime.Inclinations of 70° are taken into account, as the kinematic indicators show this average inclination.The maximum of dilation tendency at θ SHmax is accompanied by maxima of the slip tendency at 22° and 62° (Fig. 12).Note that increasing dip leads to a minimum of slip tendency around θ SHmax .The minimum slip and dilation tendencies both occur perpendicular to θ SHmax .This indicates that there is no preferential reactivation on NW-SE oriented structures, i.e., such as observed in combination with gravity highs.In contrast, structures dominated by gravity minima and NE-SW oriented may be reactivated for both opening and shearing failure if the dip tends to angles of about 70°.

Discussion
The eastern sector of the TMVB is characterized by two coeval main systems of faults, respectively, NNW-and NE-oriented (Norini et al. 2015;García-Palomo et al. 2018;Olvera-García et al. 2020a, 2020b;Gómez-Alvarez et al. 2021;Mennella et al. 2022).These two orientations were recognized in the morpho-tectonic lineaments in the study area.The kinematic meaning is suggested by the analysis of fault surfaces carried out in the Las Minas area, considered the analogue, exhumed geothermal system of Los Humeros.Here, Olvera-García et al. (2020a) recognized that the NW-oriented structures are dominantly oblique to normal, with the latter often being the last kinematic event overprinting the previous kinematic indicators, whereas the NE-oriented faults resulted in a dominantly normal movement.This different evolution has been explained considering two kinematic events: during the first, the NW-oriented crustal stretching is the prominent force, thus the NW-oriented structures play the role of transfer faults (oblique kinematics), acting with the NE-oriented normal faults; during the second event, uplift is dominant, and preexisting structures (i.e., both NW-and NE-oriented faults) are forced to play the role of normal faults.
These oriented NE-and NW-oriented structures were also observed both in the residual anomaly map (Fig. 7b) and in the residual anomalies obtained by filtering analysis (Fig. 9). Figure 9a-c, shows how the major NE-SW as well as the NW-SE-oriented structures control the major residual anomalies, and while other minor structures (and caldera structures as well) play a subdued role in respect to those.The same relationships are highlighted in Fig. 9d, e.However, a quantitative interpretation of the anomalies is not presented in this study.
The calculated gravity carried out by forward modeling in the Los Humeros area (Fig. 8) shows a sub-circular gravity low anomaly, while the observed gravity (Fig. 7b) shows a lineal trend of a low gravity in the NE-SW direction, where the geothermal field is located.To the SW, the low gravity is interrupted by a NW-SE medium-high anomaly: an abrupt change in density between the Maxtaloya and Arenas faults (Fig. 6).This change is not only observed in density, but also in magnetization and resistivity (Corbo-Camargo et al. 2020) and limits the most productive zone of the geothermal field.
The discrepancies between the forward-modeled gravity data using an existing geological model and the observed gravity data, even after applying Butterworth filters, underscore the complexity of the subsurface structure.The geological model may not be able to represent the complex subsurface structures and variations in gravity anomalies.
Our estimation of the stress state revealed a trans-tensional fault regime at depth in the basement, ranging from strike-slip to normal faulting.According to Kruszewski et al. (2022), a normal fault regime based on focal mechanisms predominates at depth.This suggests that deformation in the basement is more influenced by the regional stress state and less by caldera collapse, which is consistent with the resurgence hypothesized by Norini et al. (2015) and Toledo et al. (2019).These faulting regimes are in accordance with the dominant kinematic movements of the NNW-and NE-striking fault systems proposed by Olvera-García et al. (2020a).Furthermore, the latter is in agreement with the predominant NE-and NW-oriented structures obtained by our gravimetric and morpho-structural analyses.
Our analysis of slip and dilatation tendencies shows that NE-SW oriented structures have the highest reactivation potential.This orientation coincides with one of the main orientations identified in the morpho-structural (Fig. 6) and gravity analyses (Figs. 7 and  9), and with the orientation of some of the main faults present in the geological model (Fig. 3).The relationship between high reactivation potential, fracture porosity and lowgravity anomalies has already been demonstrated in various geothermal environments (Altwegg et al. 2015;Guglielmetti et al. 2013;Baillieux et al. 2014).
In contrast, other considerations can be made for the NW-SE trending faults.The latter run almost parallel to the main tectonic stretching direction and are therefore hardly exposed to dilation.This is obviously a consequence of the boundary conditions imposed by the calculation, which consider the fault orientation as the decisive factor for a favorable dilatancy.However, the parallelism between the tectonic stretching direction and the fault trend supports the interpretation that the NW-SE trending faults, which were active at the same time as the normal faults, played the role of transfer faults within the same extensional regime, as proposed by Olvera-García et al. (2020a), who studied the Los Humeros analog.
The contribution of the transfer fault zones, i.e., with strike-slip to trans-tensional kinematics, in a geothermal context, has been investigated in different areas (e.g., Alçiçek et al. 2013;Liotta and Brogi 2020;Liotta et al. 2021).The key factor is linked to the orientation of the intermediate stress axis, i.e., from vertical to oblique, determining the kinematics of each single transfer fault.As described by Sibson (2000), the intermediate stress axis mimes the orientation of the structural channels that drive fluid flow.The more vertical it is, the more efficiently the deep fluids are distributed laterally into the dilated normal faults.
Finally, it is important to mention that our analysis was based on a generic fault model since the maximum depth of each fault in the LHGF is unknown.For this reason, we estimated the BDT depth as the top of the transition zone from brittle to ductile rheological behavior, and not all faults necessarily reach that depth.They could go even deeper, as the BDT is not a sharp boundary but a transition zone, as shown by the fact that earthquakes can occur to some extent at deeper structural levels.
In this framework, the inversion of gravity field data shows that density starts to decrease at a depth of about 3 km (Fig. 11b, Carrillo et al. 2022) and that local seismic activity peaks at this depth before decreasing (Fig. 11a, Toledo et al. 2020), as expected when the BDT is reached.

Conclusions
The morpho-structural approach and the gravimetric analysis show that the structural setting of the LHVC is mainly dominated by NE-SW and NW-SE oriented structures, which also favored the development of the Los Humeros caldera and controlled the geometry.The alignments are not only observed inside the volcanic complex, but also in the Las Minas Valley, the Perote fault in the Cofre de Perote, and to the SW of Los Humeros.These structures appear to be independent of the presence of the volcano and calderas and are therefore considered to be of tectonic origin.Furthermore, these oriented tectonic structures appear to have direct control over the emplacement of the Los Humeros volcano, as also shown by regional geology and residual anomalies.
Because the caldera collapse is a Pleistocene event and today resurgence, i.e., uplift, is considered effective, our slip and dilation tendency analysis was focused on the regional tectonic structures.The combination of high reactivation potential and low gravity along regionally significant structures that are NE-SW oriented in the study area indicates fracture porosity as the cause of the low gravity.
observed here is coherent with the dextral displacement of the Rhyolitic domes depicted in Fig. S3a.Men give scale.Figure S6.Radially average power spectrum of the Complete Bouguer Anomaly (Figure 5a of the main manuscript).A cutoff wavelength of 20.5 km was selected to separate the lower and higher frequencies, related to regional and local anomalies, respectively.The cutoff wavelength was chosen where the slope of the power spectrum changes, indicating a change in the depth of the anomalous source.Figure S7.Location of the kinematic indicators collected on the main structures within Los Potreros caldera (Table 1) displayed on a 1 m high-resolution DEM (kindly provided by CeMie-Geo project).See the kinematic indicators associated to these structures in Table S1. Figure S8.Average total principal stresses and pore pressure with depth in Los Humeros area.A strike-slip to normal faulting regime transition is observed at 3.3 km depth at a P p of 23.83 MPa.Table S1.Summarized kinematic indicators collected on dominant normal faulting main structures within the Los Potreros caldera.Average fault dip: 70°.

Fig. 2
Fig. 2 Schematic geological map of the LHVC compared to the Digital Elevation Model (from Instituto Nacional de Estadística y Geografía, INEGI, Mexico) after Norini et al. (2015).The location of the LHVC within the TMVB is shown in Fig. 1

Fig. 3 a
Fig. 3 a Integrated 3-D geological model of Los Humeros (Calcagno et al. 2022) with four geological groups: basement (green), pre-caldera (blue), rocks from the caldera (purple), post-caldera rocks (brown).b Visualization of the faults: regional faults (purple) and local faults (red) Fig. 4 Standard deviation of the a gravity and b vertical component of the differential GPS measurements.Standard deviation windows of 0.05 mGal and 0.2 m, respectively

Fig. 5 a
Fig. 5 a Interpolated regional complete Bouguer anomaly using a density of 2670 kg m −3 provided by the Comisión Nacional de Hidrocarburos (http:// www.gob.mx/ cnh).b Regional trend of the complete Bouguer anomaly using a Gaussian filter.Inner/outer black rectangles delineate the areas of the processed gravity data in the LHGF and Butterworth filtering analysis, respectively.The white rectangle shows the area of the 3-D geological model

Fig. 6
Fig.6Regional structures of the LHVC inferred from a morpho-structural analysis of the eastern sector of the TMBV.The black dashed rectangle delineates the area of the Butterworth filtering analyses (Fig.5).The green dashed rectangle delineates the area of the morpho-structural analysis shown in detail in Additional file 1: Figs.1-5

Fig. 7 a
Fig. 7 a Complete Bouguer anomaly of the modeled area.b Residual anomalies of the LHGF and surroundings.Black crosses indicate gravity stations.Red and black lines represent major faults and the boundaries of the caldera, respectively(Calcagno et al. 2022)

Fig. 8
Fig. 8 Gravity response obtained from forward modeling of the 3-D geological model of Los Humeros (limited to 0 m.a.s.l.).Red and black lines represent major faults and the boundaries of the caldera, respectively(Calcagno et al. 2022)

Fig. 10
Fig. 10 Rheological profiles for the LHGF (a), and the TMVB (b).Red and blue lines delineate the geothermal fluid pressure with pore fluid factor λ inside and outside of the caldera, respectively.The latter reveals hydrostatic condition.Black lines show assumed tectonic strain rates in the ductile regime

Fig. 11 a
Fig. 11 a Local earthquake distribution in Los Humeros from 1-D inversion(Toledo et al. 2020), and b density variation with depth from the regional density model in Los Humeros(Carrillo et al. 2022)

Fig. 12
Fig. 12 Slip and dilation tendencies of 70° and 90° dipping structures as a function of the fault plane azimuth under a normal faulting regime at 3.2 km depth and a maximum horizontal stress of N42°E

Table 1
(Weydt et al. 2021ty values, ρ bulk, and standard deviation of analogue samples of the geological groups of LHGF and Las Minas(Weydt et al. 2021) used in the forward modeling

Table 2
Wavelength ranges of the high-and band-pass filters applied to the complete Bouguer anomaly

Table 4
Total and effective principal stresses are estimated for the depth of 3.2 km, where the transition of strike-slip to normal faulting regime occurs