Local earthquake tomography at the Los Humeros geothermal field, Mexico

The GEMex * project is a recently finalized European-Mexican collaboration that aimed to improve the understanding of two geothermal fields: Acoculco and Los Humeros Volcanic Complex . These sites are located in the Trans-Mexican Volcanic Belt, a region that hosts numerous active volcanoes and is favorable for geothermal exploitation. Currently, the Los Humeros Volcanic Complex is one of Mexico’s main geothermal systems with an installed capacity of ~95MW. Many studies have been performed at this site since the 70s highlighting several features and characteristics of the shallow subsurface. However a thorough knowledge of structures and behavior of the system at greater depths is still quite sparse. Hence one main objective of the GEMex project was to conduct several geological, geochemical, and geophysical experiments to investigate deeper structures for future development of local and regional geothermal resources.


Introduction
Los Humeros Volcanic Complex (LHVC) is a superhot geothermal system located at the eastern edge of the Trans-Mexican Volcanic Belt (TMVB), a volcanically active region favorable for geothermal energy exploitation. It is one of the oldest producing fields in the region, with more than 60 wells drilled up to ∼3 km deep since the early 1980s (Arellano et al., 2003;Cedillo-Rodríguez, 1999;Gutierrez-Negrin & Izquierdo-Montalvo, 2010;Rocha-López et al., 2010). Currently, it has an installed capacity of ∼95 MW electric power and is administered by the Comisión Federal de Electricidad (CFE) (Romo-Jones et al., 2018). Temperatures as high as 400 • C have been measured in several producing wells at ∼2.5 km depth. However, geothermal fluids at these temperatures are presently not being exploited. Despite the large number of studies on the geochemical (e.g., Martinez & Alibert, 1994), geological (e.g., Carrasco-Núñez, López-Martínez, et al., 2017), structural (e.g., Norini et al., 2015), and geothermal (e.g., Gutierrez-Negrin & Izquierdo-Montalvo, 2010) properties of the reservoir, a solid understanding of the conditions and underground structures at depth is still rather sparse. Only a few deep probing geophysical studies (resistivity 2-D profiles and seismic surveys) in recent years have provided notions of the local stress field and structures of the geothermal field (Arzate et al., 2018;Gutierrez-Negrin & Quijano-Leon, 2004;Lermo et al., 2001Lermo et al., , 2016Lermo et al., , 2008Norini et al., 2019;Urban & Lermo, 2013).
One objective of this study is to investigate the deeper structures of the geothermal system and to locate and better understand the deep superhot fluids for their exploitation. Passive seismic methods are to this purpose widely exploited in geothermal prospecting (e.g., Calò & Dorbath, 2013;Jousset et al., 2011;Muksin et al., 2013). Seismic properties such as the compressional P (Vp) and shear S (Vs) wave velocities, and Journal of Geophysical Research: Solid Earth 10.1029/2020JB020390 the Vp/Vs ratio structures have proven reliable tools to describe lithologies and possible variations due to changes in fluid composition, rock porosity, and temperature (e.g., Gritto & Jarpe, 2014;Husen et al., 2004;Ito et al., 1979;Mavko & Mukerji, 1995). These are key features in geothermal exploration and monitoring.
One conventional approach to obtain the seismic properties of a target area is the 3-D tomographic inversion of P and S wave arrival times from local earthquakes, as observed in records of seismometers deployed in the area of interest. The 3-D velocity structure is typically obtained through a joint inversion of hypocenter locations and velocity structures using an a priori parameterized 3-D grid model of the subsurface. Classical tomographic results are strongly influenced by the inversion grid or node spacing choice, and hence, its adequate selection is fundamental to retrieve the main features of the subsurface for reliable interpretation. A too fine model could, for example, lead to poor-resolution values and/or artifacts such as grid oscillations, whereas a too coarse model (especially a coarse fixed grid) could overlook smaller underground features. In addition, significant smearing can be introduced when the chosen grid does not follow the orientation of the main anomalies. In this work, we extend the conventional tomographic method of a single fixed model grid by using a postprocessing statistical approach. We compute and average several inversions using different model parametrizations to achieve higher spatial accuracy, reduce the effects of poor parametrization selection, and overall increase model resolution.
In this study, we image the 3-D Vp and Vp/Vs structures, along with the seismicity distribution at Los Humeros geothermal field. In the first part of this study, we compile information on the geological and structural setting of Los Humeros area. Later, we describe the passive seismic experiment and the data processing workflow followed to detect and locate the local seismic events. We use the retrieved earthquake catalog to derive a new so-called minimum 1-D velocity model in section 4. In section 5, we compute and average the 3-D tomography of several parametrized models using the minimum 1-D velocity model as initial input. Finally, section 6 proposes a first interpretation of the obtained results in relation to existing geological information and newly acquired petrophysical, geochemical, and geophysical data (Bär & Weydt, 2019;Benediktsdóttir et al., 2019;Jentsch et al., 2020;Lucci et al., 2020;Urbani et al., 2020).

Seismic Network
From September 2017 to September 2018, we deployed and maintained a temporary seismic network comprising 25 three-component broadband (Trillium Compact 120s) and 20 three-component short-period (Mark L-4C-3D) sensors recording continuous seismic data at sampling rates of 200 and 100 Hz, respectively . The array consisted of two complementary subnetworks each configured to characterize shallow and deeper structures using different seismic processing techniques ( Figure 2). A denser (∼1.6-2 km interstation distance) pseudorhomboidal array was located mainly in the inner Los Potreros caldera where previous studies have identified the occurrence of local seismicity (Gutierrez-Negrin & Quijano-Leon, 2004;  Lermo et al., 2001Lermo et al., , 2016Lermo et al., , 2008Urban & Lermo, 2013) and where most of the producing and injecting wells are located. This subnetwork was primarily designed for local microseismicity retrieval , local earthquake tomography, beamforming of ambient noise (Löer et al., 2020), time-reverse imaging (Werner & Saenger, 2018), and autocorrelation techniques (Verdel et al., 2019). The second much sparser network (∼5-10 km interstation distance) was placed around the outer Los Humeros caldera and was mainly intended for imaging deeper large-scale structures with techniques such as ambient noise tomography (Granados et al., 2020;Martins et al., 2020), among others.

Local Earthquake Detection
We focused the event detection mainly on Los Potreros caldera  using Python tools based on the ObsPy library (Beyreuther et al., 2010). We calibrated a recursive short-time-average through long-time-average (STA/LTA) detection algorithm (Trnkoczy, 2012;Withers et al., 1998) on several days of the recently acquired seismic data set (2017-2018) and on a set of local seismic events recorded between 2005 and 2006 by the permanent network operated by the CFE (Lermo et al., 2008). We exhaustively tested the detection performance through several days of the recent seismic database using a wide range of parameter combinations. The optimum parameters selected were a combination of bandpass filter between 10 and 30 Hz, STA and LTA windows of 0.2 and 2 s, respectively, and on and off trigger thresholds of the computed STA/LTA function at 3.5 and 1.0, respectively. To account for the P and S wave arrivals, the STA/LTA function was computed from a single amplitude trace determined by the square root of the sum of the three single component squared traces for each station. Finally, a detection was declared as such when the triggering window of at least five stations from the dense subnetwork coincided in time (Trnkoczy, 2012;Withers et al., 1998). We reviewed each triggered detection and manually picked P and S wave arrivals of local events and their associated empirical uncertainty range using the Python Obspyck tool (Megies, 2016).
From a total of 1,586 detections, 488 were identified as local events. After picking P and S phases, these earthquakes were located using an oct-tree search (Lomax et al., 2000(Lomax et al., , 2009) in a homogeneous 3-D volume with a P wave velocity of 3.5 km/s and a Vp/Vs ratio of 1.73 ( Figure 3). Later, we reselected the seismic events with a greatest angle without observation (GAP) of less than 180 • and at least three P and three S wave arrivals (333 events in total) for the calculation of a minimum 1-D velocity model and their relocation. The recorded seismicity is mostly located below the dense array within Los Potreros caldera and mainly grouped into three distinctive clusters, marked as C1, C2, and C3 in Figure 3.

The 1-D Velocity Model
We use the retrieved travel time data from the filtered catalog (333 events with 2,146 P wave and 2,146 S wave picks) as input for a joint inversion to determine the so-called minimum 1-D Vp and Vs models and the hypocenter relocations using the code Velest (Kissling et al., 1994). The code Velest iteratively computes forward modeled data (predicted travel times), using a ray tracer in an initial model (1-D velocity model, hypocenter locations, and station corrections), compares the synthetic data to the observed data set, and updates the model such that the RMS (root-mean-square) misfit between the two is minimized (Tarantola, 2005). This procedure uses regularization parameters (damping factors) to tackle instabilities due to data uncertainties (Levenberg, 1944;Marquardt, 1963) and continues until a maximum number of iterations is reached.
The estimation of a minimum 1-D model consists of a trial and error process in which typically a broad range of plausible initial models is tested to ensure covering as many potential solutions as possible and select the best fitting model. This procedure is necessary because the inversions are based on linearization and thus strongly depend on the initial model. In this work, we performed the inversion of 10,648 initial models with varying P wave velocities at the surface, vertical velocity gradients, and Vp/Vs ratios (thus also varying Vs models) over five iterations. The software Velest allows tracing rays to the true station elevations. However this option poses the limitation of locating all stations within the first layer. With this in mind, we set the uppermost layer thickness to more than 1 km, which corresponds to the approximate elevation difference between the highest and lowest recording stations. The following layers are then defined roughly every 0.5 km at shallow depths and progressively increase to 1 and 2 km for deeper levels. The depth intervals were chosen taking into account well data interpretation (Norini et al., 2015(Norini et al., , 2019 and exhaustive testing. We define depths relative to sea level throughout this manuscript. The uncertainty of P and S arrivals was defined by the weighting factors 0, 1, 2, and 3 corresponding to the estimated picking uncertainties of up to 0.03 s, 0.08 s, 0.12 s, and higher, respectively. Finally, we selected Station DB13 as the reference station given its location at approximately the center of the dense array and geothermal field and its high number of recorded P and S arrivals (red circle in Figure 2). Figures 4a and 4b depict all the initial (yellow lines) and the 35 resulting velocity models with lowest associated RMS misfit values (gray lines) for Vp and Vs, respectively. The misfit values for this set of best models range from 95.0 to 96.1 ms. Note the good agreement for several layers of the resulting velocity models, particularly between −1.0 and 1.5 km depths where most of the seismicity is located. This coherence becomes less obvious at shallower and greater depths, where only a few hypocenters are located. Few models show values with large deviations from the more obvious trend. We create a heat map with this set of final models in Figures (Lermo et al., 2008;Löer et al., 2020), (b) the resulting Vp/Vs ratio, and (c) the earthquake distribution over depth after the 1-D inversion. Solid lines indicate the depth intervals with best sensitivity for each model. Figure 5 depicts the resulting P and S wave velocity models, the associated Vp/Vs ratio, and the final event distribution for the selected minimum 1-D velocity model. The model is best resolved between −2 and 2 km depths approximately, which is consistent with the hypocenter distribution shown in Figure 5c. The seismic events are restricted up to approximately 4 km depth from the surface, with a maximum number of events between −0.5 and 0.2 km depths. The seismicity presents a systematic shift toward ∼0.5-0.8 km greater depths and some improvements in clustering after the 1-D inversion. Two velocity models proposed by Lermo et al. (2008) (P wave) and Löer et al. (2020) (S wave) are marked in green and magenta, respectively, in Figure 5a. The Vp model proposed by Lermo et al. (2008) was derived using a seismic reflection profile and reaches an approximate depth of −0.5 km, below which a default value of 5.18 km/s is assigned. The Vs model derived by Löer et al. (2020) was obtained using three-component ambient noise beamforming and is most sensitive in the interval between −0.5 and 10 km depths. Notice the good correlation between the derived minimum 1-D Vp model and the model obtained by Lermo et al. (2008), especially between −2.0 and −0.5 km. Similarly, there is a good agreement between the derived minimum 1-D Vs model and the model obtained by Löer et al. (2020) between −1.0 and 2.2 km depths.
The station delays associated with the selected minimum 1-D velocity model are also consistent with the local geology and are further described in Appendix A.

The 3-D Seismic Tomography
After selecting the reference 1-D velocity model and hypocenter locations, we used the 3-D travel time inversion code SIMUL2000 (Eberhart-Phillips, 1990;Eberhart-Phillips & Michael, 1998;Evans et al., 1994;Thurber, 1983) to estimate the 3-D velocity structure of the geothermal field. Forward calculations are computed using a pseudo bending method (Uhrhammer, 1987), and inversions are performed using an iterative damped least squares scheme. The software SIMUL2000 allows for the simultaneous inversion of Vp and Vp/Vs ratio instead of Vs to account for the generally lower resolution of Vs models due to larger uncertainties of S wave arrival determination, most of them being hampered by the coda of the P wave (Thurber, 1993;Thurber & Eberhart-Phillips, 1999). Inversions in this section are computed using the minimum 1-D Vp model and a homogeneous Vp/Vs of 1.71 obtained from a Wadati diagram analysis of the stacked events.

Model Parametrization
An appropriate model parametrization is suggested by Evans et al. (1994) and Husen et al. (2000Husen et al. ( , 2003 as that with the finest possible node spacing, which allows inversions without strong derivative weighted sum (DWS) heterogeneities. The DWS is a measure for ray density which takes into account the number of crossing rays, ray-node separation, and raypath length in the vicinity of each node (Evans et al., 1994;Husen et al., 2000). Kissling et al. (2001) advise taking into account both a priori information about the underground structure and resolution capabilities of a data set when selecting appropriate model parametrization. Accounting for known subsurface features could lead to a too fine model node spacing which in turn could result in lower-resolution values and sparse imaging. On the other hand, a coarse model parametrization, although yielding higher-resolution values due to increased ray density, could potentially overlook smaller features in the subsurface. To avoid this effect, some authors (e.g., Abers & Roecker, 1991;Bijwaard et al., 1998) use uneven node spacing in their inversions. This technique could, however, complicate interpretations given some velocity changes may be inadvertently interpreted as underground features. One methodology often used is a graded inversion approach (Evans et al., 1994;Husen et al., 2003). Iterations are performed through finer grids using a coarse model output as input for a new inversion with a finer grid.
Some novel postprocessing techniques include the averaging and weighted averaging of several inversion results using different initial model parametrizations (e.g., Haberland et al., 2009). This technique helps enhancing velocity anomalies, reduces possible smearing effects and model noise due to the parametrization choice, overcomes the fixed coarse parametrization limitation of the inversion code, and improves the model resolution. In this study, similar to Calò and Dorbath (2013), we compute several inversions using different model parametrizations which we later average on a finer grid. Figure 6 shows the event to station raypath configuration in the initial minimum 1-D velocity model. Raypaths are unevenly distributed and mostly located within Los Potreros caldera, reaching ∼1.5 km depth for the most part. Taking into account the raypath configuration, we chose an initial lateral parametrization of 1 × 1 km 2 (Figure 7) and 0.5 km internode spacing with depth within Los Humeros caldera region (Figure 8). We then progressively increased the node spacing in regions outside Los Humeros caldera and below the seismicity. Figures 7 and 8 show the DWS distribution after a 3-D inversion using the chosen grid. DWS values are, as expected, larger within the regions above the seismic clusters (higher ray density) and extend to 0 to −1 km depth. After the inversion, the seismicity presents a systematic shift toward ∼0.5-0.8 km shallower depths which could be attributed to initially setting the station delays of the minimum 1-D model to 0. To avoid decreases in retrieved velocity amplitudes, we fixed the hypocenter locations during the first iteration and updated the new relocations after each velocity inversion (Husen et al., 2003).
We constructed several inversion grids starting off with the previous grid and then rotated it in 15 • steps. A similar procedure was carried out after displacing the inversion grid center 0.2 and 0.5 km toward the north, south, east, west, northwest, northeast, southeast, and southwest. Figure 9 depicts the nodes of the 228 inversion grids used. Although Calò (2009) proposes the use of fewer models to reach a decent benefit of the postprocessing technique, we opted for a larger number of inversions to yield a more statistically stable final model.
Note that both Los Potreros and Los Humeros calderas are densely covered by nodes. To construct a new averaging grid, we interpolated all outputs post inversion onto a regular finer grid with 0.1 km spacing using the same interpolation scheme as applied in SIMUL2000 (Thurber, 1983). Then we averaged each point of the new grid. Figure 10 shows several depth slices for the average DWS using the new grid. The covered gray areas (regions with any ray density) become larger than when using a single initial inversion grid, which could be attributed to both the averaging of DWS values, but also to taking into account model parametrizations that favor different ray orientations.

Regularization
We performed a trade-off test (Appendix B) to select adequate damping values for a single inversion grid (Eberhart-Phillips, 1990). Given the node spacing remained the same for all 228 inversion grids, we kept these damping factors fixed. Figure 11 shows the misfit value changes with each progressive iteration for the 228 inversions. The misfits start at 0.2 s for all models and progressively reduce to a mean value of ∼0.102 s, which falls within the range of the picking uncertainties, thus confirming the successful inversion of the models used with the chosen regularization parameters.

Model Quality and Uncertainty
An adequate quality assessment of the solution is typically carried out to validate inversion performance. The main objective is the identification of poorly resolved areas and unrealistic model perturbations or artifacts that may have been introduced and affect interpretation. Several parameters and procedures that analyze ray distribution and density include the evaluation of the hit count, the diagonal element of the model resolution matrix (MRM), the spread function (Michelini & McEvilly, 1991;Toomey & Foulger, 1989), the smearing information of each node provided by the MRM, and tests with synthetic data such as checkerboard (Zelt, 1998) and recovery tests. In this study we calculate and display the results of a checkerboard test. Averaged spread values and resolution diagonal elements (RDE) are shown in Appendices C and D, respectively.
We carried out a checkerboard test to determine what kind of anomalies can be retrieved with the seismic network and catalog used. We perturbed the minimum 1-D Vp model with alternating ±12% anomalies and the starting Vp/Vs model (constant value of 1.71) with ±10% perturbations. These anomalies are in agreement with the range of velocities obtained after the 3-D inversion. Positive and negative anomalies are indicated in blue and red colors, respectively (Figures 12and13). They comprise four nodes of the single grid (Figure 7) in the horizontal directions (∼1.5 km × 1.5 km) and two nodes in the vertical direction (∼0.7-1 km). These perturbations have the approximate size of some of the main anomalies obtained in the real data inversion. We computed synthetic traveltimes using the retrieved seismic catalog and added Gaussian noise with ±0.065 s standard deviation, which corresponds to the standard deviation of the uncertainty distribution of the manual picks. Similar to the real data inversion, we determined the damping parameters by performing a trade-off test. Then, we carried out inversions for the 228 grids of Figure 9 and averaged the results.  Figure 12 shows the recovered velocity anomalies for Vp and Vp/Vs across several depth slices. Regions with higher ray density, particularly those closer to the identified seismic clusters, appear to be best recovered. These areas coincide with those of lower spread and high RDE values in Figures D1 and C1, respectively. Depths between ∼−2.10 to −1.6 km seem best resolved (Figures 12a-12d), with errors of approximately ±2-3% in the best cases for both Vp and Vp/Vs. At shallower and deeper levels (Figures 12e and 12f) the polarity of the anomalies are recovered only toward the center and smeared toward the edges of Los Potreros caldera. Velocity uncertainties in these regions vary around ±5-8%. The checkerboard is also well reproduced in some regions of the vertical sections ( Figure 13). The recovery looks best at the center of the Cross Sections A2-A2 ′ , B2-B2 ′ , and C2-C2 ′ to a depth of −1 km, after which anomalies become smeared once more. The polarities of some velocity anomalies, however, are reproduced to 0 to −0.2 km depth toward the center of Section B2-B2 ′ . Overall, the uneven seismicity distribution (raypath configuration) marks a very limited area (shallow north central portion of Los Potreros caldera) of good recovery. This region coincides with the area where spread values fall below a value of 1.5.
It is worth noting that a single inversion using the same grid configuration as the one used to produce the synthetic model perturbations is able to retrieve these anomalies accurately. However, when the parametrization is considerably different (e.g., a shifted or rotated grid), the reconstruction is worsened and significant smearing is introduced, distorting the anomalies. In reality it is complicated to know beforehand the location and direction of the perturbations. By averaging several inversions using different model configurations, the dependence of the results on a single model parametrization is reduced, therefore significantly enhancing the accuracy of the final model.

Results and Discussion
The retrieved Vp (Figures 14 and 15) and Vp/Vs (Figures 16 and 17) models together with the seismicity distribution provide new insights on the geometry of the different geological units of LHVC and the relation between some of the main structures and the geothermal system. In this section, we describe and discuss the results of the 1-D and 3-D velocity inversions in relation to existing geological and newly acquired petrophysical, geochemical, and geophysical data at Los Humeros geothermal field.

The 1-D Velocity Model
The selected minimum 1-D velocity model ( Figure 5) shows not only a good agreement with alternative studies on the region (Lermo et al., 2008;Löer et al., 2020) but also one possible geological boundary observed in retrieved well data (Norini et al., 2019). A prominent discontinuity is seen at ∼−1.2 km depth in Figure 5 for both Vp and Vs models. This discontinuity could potentially mark an average transition between the precaldera group and the sedimentary basement, as it is also seen in well data at ∼−1.5 to −0.2 km depth (Norini et al., 2019). Vp/Vs ratio values are fairly constant but rather low at different depth levels, with the exception between −0.5 and 0 km depths, where most of the seismicity is concentrated (hence the larger average value of 1.71). These unusually low values could be the consequence of lateral averaging and the irregular distribution of sources and receivers. For this reason, we consider that the minimum 1-D velocity model should not be overinterpreted.

Seismicity Distribution 6.2.1. Main Results
The final earthquake catalog was obtained by averaging the coordinates and origin times of the output catalogs of each inversion. The standard deviation of each component was used to quantify the location uncertainties which were on average 131 m, 127 m, 214 m, and 0.027 s for x, y, z, and origin time, respectively. These values represent in some way the intrinsic trade-off between the errors of the model and hypocenter estimations.

Discussion
Similar to the catalog obtained after the initial relocation, the seismicity is grouped in three clusters (Figure 3). The northernmost cluster (C1 in Figure 3) has already been evidenced in Lermo et al. (2008) and is situated close to the main production area, where two out of three neighboring injection wells are located (red stars in Figure 3). The southwestern cluster (C2 in Figure 3) is located to the west of Los Humeros fault, close to a third injection well. Finally, a deeper cluster (C3 in Figure 3) is located toward the east, between Las Papas and Las Viboras faults. The remaining events are scattered within the geothermal field, with some following major structures such as the Maztaloya fault.
Cluster C1 (Figures 15a and 17a) has a narrow subvertical distribution relating to Los Humeros Fault Zone at the surface. We define Los Humeros Fault Zone as the combination of very closely spaced (∼100-250 m) N-S fault strands (Loma Blanca fault, Los Humeros fault, and Los Conejos fault). In a similar manner, Cluster C2 (Figures 15c and 17c) reflects the position of Los Humeros fault further south. Given their vicinity to injection wells (see green vertical lines in Figures 15 and 17), most of them could probably be induced/triggered events. The third cluster (C3) is located at a deeper level toward the east (Figures 15b, 15c, 17b, and 17c). At the surface, this region coincides with the area between the E-W trending Las Papas and Las Viboras faults but does not seem directly associated with any geothermal wells. However, the proximity of this cluster to C2 could potentially indicate a deeper fluid pathway toward the east. These two faults show no hydrothermal alteration along their strike at the surface (Norini et al., 2019). However, the presence of this cluster could hint to an increase of permeability of these faults at greater depths.

Vp Structure 6.3.1. Main Results
Our results provide new detailed insights into the 3-D P velocity structure of Los Humeros geothermal field. We present our results in a series of horizontal depth slices ( Figure 14) followed by E-W vertical sections ( Figure 15) across Los Potreros caldera. Average P wave velocities range from ∼2-4.1 km/s, with standard deviations in the order of ±0.06 km/s ( Figure E1). We show velocity perturbations relative to the minimum 1-D Vp model, as it is easier to observe the physical properties of rocks (e.g., presence of fluids and temperature) in this form. Nevertheless, we show absolute velocity values with contour lines in the cross sections to ease interpretation.
Close to the surface (Figure 14a), a large low-velocity anomaly (∼−8 to −12%) marked as a (Figure 14a) is located at a highly faulted region toward the center of Los Potreros caldera. This anomaly is surrounded by high-velocity anomalies (∼+6 to +10%) to the northeast, south, and west (marked as b). Further in depth (Figure 14b), a clear division between high (∼+10%) and low (∼−5%) velocity anomalies (c and d) is separated by the main Los Humeros fault. This velocity contrast remains at depth, although the low-velocity anomaly to the east appears attenuated (∼−3%) at deeper levels. One large high (∼+10 to +13%) velocity anomaly (h) appears at −1.60 km depth (Figure 14c) almost following Los Humeros and Maztaloya faults, which extends in a much narrower corridor at −1.10 km depth (j in Figure 14d). If we observe the Vp variations in cross sections (Figures 15a-15c), the high-velocity anomaly is located mostly toward the east of Los Humeros normal fault. At larger depths, a second minor high (∼+3%) velocity anomaly (mainly visible in Sections B-B' and C-C') appears further to the east where Cluster C3 is located. This feature barely reaches the limits of the imaging capabilities of our data set and must therefore be interpreted with caution.
A smaller low (∼−5 to −3%) velocity anomaly appears west of the northern portion of Los Humeros Fault and is traceable at depth (f , g, and i in Figure 14). This anomaly is also seen in the cross sections  ( Figures 15a-15c), mostly west of the buried La Antigua fault. The low-velocity region visible to the east in Cross Sections B-B' and C-C' is associated to d in Figure 14b.

Discussion
Some of the velocity anomalies at −2.6 km depth (Figure 14a) accurately follow the surface geology (see Figure 2). The low-velocity anomaly marked as a is located where undefined pyroclastics belonging to the postcaldera stage are deposited (Figure 1). In a similar manner, the high-velocity anomalies marked as b are related to regions with rhyodacitic, andesitic, and basaltic volcanic rocks.

Note.
Modified from Bär and Weydt (2019) 10.1029/2020JB020390 The high-velocity anomalies observed at depth correspond to intrusive-like bodies in the absolute velocity contour lines. Similar structures have been indicated by the combined interpretation of structural field analysis, and analog modeling at Los Potreros caldera and interpreted by Norini et al. (2019) and Urbani et al. (2020) as discontinuous resurgence associated with the intrusion of multiple magma bodies rather than a single magma chamber.
The shallow low-velocity anomalies in Cross Sections B-B' and C-C' indicate the increased thickness of the postcaldera unit in these areas and reveal the variable deposition of volcanic materials during the complex volcanic history of the caldera.
To determine possible unit boundaries, we marked the positions of several neighboring interpreted wells (Carrasco-Núñez, López-Martínez, et al., 2017) in the cross sections shown in Figure 15. We compared the depth ranges of the units seen in the interpreted wells with ultrasonic pulse velocity measurements of collected core and outcrop rock samples (Table 1) and our retrieved velocities. Then, we marked approximate unit boundaries with solid gray lines in Figure 15. We deduced the postcaldera stage (pyroclastics) with Vp around ∼2.2-2.4 km/s at shallow depths. This velocity range is well within the range of laboratory measurements of collected rock samples (2.0-3.7 km/s). The caldera stage (mainly ignimbrites) lower boundary was interpreted at around −2.0 km depth with average Vp of 2.8 km/s. Laboratory measurements range between ∼1.8-3.5 km/s for rocks found in this unit. Below, we interpreted the precaldera unit (andesites) with Vp up to ∼3.8-3.9 km/s. After that the marbles and limestones belonging to the basement start in our sections with Vp ≥ 3.9 km/s. The basement boundary is shallower close to Los Humeros fault zone.

Vp /Vs Structure 6.4.1. Main Results
Average Vp/Vs values varied between 1.50 and 1.77 throughout the studied region (Figures 16 and 17), with standard deviations in the order of ±0.02 ( Figure E1).
At shallow depth (Figure 16a), a prominent low Vp/Vs anomaly (≤1.60) is located at the northern portion of Los Humeros fault system with the lowest values concentrated between Los Humeros, La Cuesta, and Cueva Ahumada faults (k). This anomaly extends at depth toward the southeast (l and n in Figure 16). This behavior is also seen in the cross sections of Figure 17, where in many cases the anomaly extends at depth mostly toward the east of Los Humeros fault.
On the other hand, two higher Vp/Vs anomalies (≥1.71) are observed at the northern portion of Los Potreros caldera (o in Figure 16) at −1.10 km depth. These anomalies are particularly evident in Figures 17a and 17b. In Figure 17a the high Vp/Vs anomalies are divided by the Los Humeros fault and the anomaly toward the east is higher.

Discussion
Although Vp/Vs values appear low in Los Humeros (minimum ∼1.5), they are not uncommon in volcanic and geothermal regions (e.g., Husen et al., 2004, minimum ∼1.57;Muksin et al., 2013, minimum ∼1.47). The shallow low Vp/Vs anomaly coincides in shape and position with a conductive body (≤10 Ωm) imaged by a new magnetotelluric (MT) survey at Los Humeros geothermal field (Benediktsdóttir et al., 2019), hinting at the location of the cap rock composed of ignimbrites. If we consider porous media (and low Vp), this region of low Vp/Vs values can be inferred as a gas bearing chamber (Gassmann, 1951). Such interpretation is supported by the modeling of low Vp and low Vp/Vs anomalies in porous volcanic rocks in Husen et al. (2004). This hypothesis is confirmed by a new survey of CO 2 emissions at the surface (Jentsch et al., 2020), where higher flux regions coincide with k in Figure 16a.
Further in depth, the high Vp/Vs anomalies could hint at regions with increased liquid content (Gassmann, 1951). The anomaly to the east of Los Humeros fault zone (Figure 17a) coincides with a region close to the bottom of a neighboring injection well. West of Los Humeros fault zone, the second high Vp/Vs anomaly coincides with generally lower Vp values (3.2-3.4 km/s), which could be an indication of rocks, namely, the andesites, influenced by the presence of liquid. These areas could potentially be considered for further exploration and exploitation of the geothermal field.
A local heat source could be assumed as located at greater depths transporting heat along permeable faults especially in the region close to Los Humeros Fault zone. Such a hypothesis would be in correlation with the analog and structural work by Urbani et al. (2020)

Conclusions
A new seismological analysis using a dense temporary seismic network was undertaken at Los Humeros geothermal field. We collected high-quality earthquake data to image the Vp and Vp/Vs models for the first time in this region. These models were obtained by extending the classical local earthquake tomography using a postprocessing statistical approach. Several models were inverted and averaged to reduce the potential bias introduced by the choice of model parametrization and enhance the final spatial resolution. The results were then carefully integrated with new geophysical, geological, and petrophysical data for interpretation.
The statistical approach reduces the potential smearing resulting from selecting model parametrizations that do not align with anomaly location and orientations, which are in many cases unknown prior to computing a tomography. The consideration of different initial grids allows for a much finer solution and helps overcome the code restriction of using a fixed coarse grid. It also allows assessing and reducing the error bars of the final solution.
From this analysis, we identified three seismogenic areas within Los Potreros caldera, one of which does not appear in direct relation to any geothermal wells. A deep earthquake cluster was located between Las Papas and Las Viboras faults. Although these faults show no hydrothermal alteration at the surface (Norini et al., 2019), the presence of this seismic cluster suggests that these faults may increase permeability at depth. The vicinity of this new cluster to another one close to an injection well could potentially highlight a deeper fluid pathway toward the east.
The main geological boundaries found from well and new petrophysical data were also found in our Vp model. The presence of two instrusive bodies supports the idea of resurgence at Los Potreros caldera (Norini et al., 2019;Urbani et al., 2020). Urbani et al. (2020) suggest that intrusions in the region are the result of the inflation of the magma chamber at depth and may represent locations of local heat source(s).
The Vp/Vs model also supports the resurgence or uplift due to the intrusion of new magma at Los Potreros caldera. High Vp/Vs ratio anomalies are located to each side of Los Humeros Fault zone, where the hypocenters in between have a subvertical configuration. This could hint at a deeper heat source transporting hot fluids upward along permeable faults. A new petrological study ) also suggests such a system, with several (ephemeral) magma pockets in the crust being fed by multiple magma transport and storage layers. The high Vp/Vs values in this region could potentially indicate higher fluid content. Therefore, this area could be further studied.
Above this anomaly, a low Vp/Vs region coincides with the conductive clay cap seen in a new MT study (Benediktsdóttir et al., 2019). The low Vp/Vs in combination with low Vp values could indicate gas bearing regions (Gassmann, 1951;Husen et al., 2004). This hypothesis is also supported by a new CO 2 emissions survey at the surface (Jentsch et al., 2020). The shallower portions with the lowest Vp/Vs value coincide with regions of higher CO 2 fluxes.
Further steps to be considered for better understanding of the geothermal system include an attenuation tomography and the imaging of deeper structures with techniques such as ambient noise tomography. In addition, a more quantitative approach such as a cluster analysis of different physical properties will be performed to improve the accuracy of the interpretation and to build a conceptual model.

Appendix A: Station Corrections Associated With the 1-D Velocity Model
Station corrections are defined as scalar terms accounting for near-surface velocity variations below each seismic station. In other words, they are potential indicators of surface geology and/or site conditions. Figure A1 shows the (a) P wave and (b) S wave station corrections associated with the minimum 1-D velocity model. There are slightly higher P delays at stations located toward the southeast of Los Potreros caldera. This is a region characterized by several fault outcrops and undefined pyroclastic deposits. The delays decrease toward the northwestern edge of Los Humeros caldera. This area is characterized by basalts and andesites at the surface (Figure 1). Station delays for stations further away from Los Humeros caldera show higher values than those within the dense array, which could be associated with larger picking uncertainties. S delays ( Figure A1b) are also relatively balanced within Los Potreros caldera. Table A1 shows the retrieved station correction values. Locations, elevations, and sensor types are available as the seismic network associated metadata in Toledo et al. (2019).

Appendix B: Trade-Off Test Sample for a Single Model Parametrization
Damping parameters for an inversion using a single model parametrization are chosen such that the data variance is minimized at a moderate model variance. We first determine the damping factor for Vp by testing several values while fixing the damping factor for Vp/Vs. In a similar manner we select a damping factor for Vp/Vs by testing a range of values in combination with the selected Vp damping factor. Through this approach, the damping parameters chosen for this experiment were 7 and 10 for Vp and Vp/Vs models, respectively ( Figure B1). Given the node spacing did not vary when inverting for the different inversion grids, damping factors remained the same throughout all inversions performed.