Local earthquake tomography at the Los Humeros geothermal field, Mexico

About: The article was published on 2020-03-09 and is currently open access. It has received 6 citations till now. The article focuses on the topics: Geothermal gradient.

Summary

1. 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.
• 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.
• The authors image the 3-D Vp and Vp/Vs structures, along with the seismicity distribution at Los Humeros geothermal field.
• In section 5, the authors compute and average the 3-D tomography of several parametrized models using the minimum 1-D velocity model as initial input.

2. Geologic and Tectonic Setting

• A second set of faults parting from the main system runs N-S, NE-SW, and E-W. Both sets disappear at the surface when approaching Los Potreros caldera rim .
• Here, the authors briefly describe the lithologies found in the four major groups.
• The lower portion of the basement, also called the Teziultlán Massif, is mainly composed of old intrusive igneous and metamorphic rocks (Paleozoic granites, greenschists, and granodiorites) (Quezadas-Flores, 1961; Viniegra, 1965; Yáñez & García, 1982).
• This unit is characterized primarily by eruptive events which resulted in the formation of Los Humeros and Los Potreros calderas (Carrasco-Núñez & Branney, 2005; Carrasco-Núñez et al., 2012; Ferriz & Mahood, 2009; Norini et al., 2019).

3.1. Seismic Network

• From September 2017 to September 2018, the authors 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 (Toledo et al., 2019).
• The array consisted of two complementary subnetworks each configured to characterize shallow and deeper structures using different seismic processing techniques .
• 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; TOLEDO ET AL.
• This subnetwork was primarily designed for local microseismicity retrieval (Gaucher et al., 2019), 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).

3.2. Local Earthquake Detection

• The authors focused the event detection mainly on Los Potreros caldera (Gaucher et al., 2019) using Python tools based on the ObsPy library (Beyreuther et al., 2010).
• The authors 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).
• 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.
• The authors 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).
• After picking P and S phases, these earthquakes were located using an oct-tree search (Lomax et al., 2000, 2009) in a homogeneous 3-D volume TOLEDO ET AL.

4. The 1-D Velocity Model

• The authors 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).
• 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 following layers are then defined roughly every 0.5 km at shallow depths and progressively increase to 1 and 2 km for deeper levels.
• Finally, the authors manually select the model with the lowest misfit value that best coincides with the main trend of most frequent velocity values as the minimum 1-D velocity model.
• 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.

5. The 3-D Seismic Tomography

• After selecting the reference 1-D velocity model and hypocenter locations, the authors 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.

5.1. Model Parametrization

• An appropriate model parametrization is suggested by Evans et al. (1994) and Husen et al. (2000, 2003) as that with the finest possible node spacing, which allows inversions without strong derivative weighted sum (DWS) heterogeneities.
• 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.
• Iterations are performed through finer grids using a coarse model output as input for a new inversion with a finer grid.
• Taking into account the raypath configuration, the authors chose an initial lateral parametrization of 1 × 1 km2 and 0.5 km internode spacing with depth within Los Humeros caldera region .
• DWS values are, as expected, larger within the regions above the seismic clusters (higher ray density) and extend to 0 to−1 km depth.

5.3. Model Quality and Uncertainty

• An adequate quality assessment of the solution is typically carried out to validate inversion performance.
• Figure 12 shows the recovered velocity anomalies for Vp and Vp/Vs across several depth slices.
• These areas coincide with those of lower spread and high RDE values in Figures D1 and C1, respectively.
• Overall, the uneven seismicity distribution (raypath configuration) marks a very limited area (shallow north central portion of Los Potreros caldera) of good recovery.
• 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.

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.
• In a similar manner, Cluster C2 reflects the position of Los Humeros fault further south.
• Given their vicinity to injection wells , most of them could probably be induced/triggered events.
• 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.

6.3.1. Main Results

• The authors results provide new detailed insights into the 3-D P velocity structure of Los Humeros geothermal field.
• Nevertheless, the authors show absolute velocity values with contour lines in the cross sections to ease interpretation.
• If the authors observe the Vp variations in cross sections , the high-velocity anomaly is located mostly toward the east of Los Humeros normal fault.
• This feature barely reaches the limits of the imaging capabilities of their data set and must therefore be interpreted with caution.

6.3.2. Discussion

• Some of the velocity anomalies at −2.6 km depth accurately follow the surface geology .
• To determine possible unit boundaries, the authors 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.
• 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.

6.4.1. Main Results

• Average Vp/Vs values varied between 1.50 and 1.77 throughout the studied region , with standard deviations in the order of ±0.02 .
• 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.
• Further in depth, the high Vp/Vs anomalies could hint at regions with increased liquid content (Gassmann, 1951).
• 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.

7. Conclusions

• A new seismological analysis using a dense temporary seismic network was undertaken at Los Humeros geothermal field.
• 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.
• 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 low Vp/Vs in combination with low Vp values could indicate gas bearing regions (Gassmann, 1951; Husen et al.

Figures

Figure 4. Results of the 1-D inversions using VELEST. The 35 best (a) P wave and (b) S wave final velocity models (gray lines). Heat maps for the same (c) P wave and (d) S wave set of final models. The yellow lines indicate all initial velocity models used. The selected minimum 1-D models are indicated in black lines in panels (a) and (b), and in red lines in panels (c) and (d).

Figure 10. Averaged DWS distribution at different depth levels. Panels (a) and (c) show two depth slices for the Vp model DWS distribution at −2.6 and −1.10 km depths. Panels (b) and (d) show the depth slices for the Vp/Vsmodel DWS distribution at −2.6 and −1.10 km depths. Darker shading indicates regions of higher ray density.

Figure 14. Vpmodel variations with respect to the minimum 1-D velocity model at different depth levels. Panels (a)–(d) show depth slices for the resulting Vp model variations at −2.60, −2.10, −1.60, and −1.10 km depths, respectively. Green circles mark the location of earthquakes ±150m away from slice. Dashed red lines indicate the boundary at which spread values are less than or equal to 1.5. Gray areas mark the regions where the DWS is less than or equal to 5.

Figure 5.Minimum 1-D model showing (a) the selected Vp and Vsmodels along with two available models (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 8. Cross Sections (a) A1-A1′ and (b) B1-B1′ for the Vpmodel DWS distribution (Figure 7). Darker shading indicates regions of higher ray density. Black crosses indicate the node positions.

Figure 9. Nodes of the 228 inversion grids used to estimate the average velocities.

Figure 13. Cross Sections A2-A2′ , B2-B2′ , C2-C2′ , and D2-D2′ (Figure 12) for the retrieved Vp and Vp/Vsmodel variations. Panels (a), (c), (e), and (g) show the recovered Vp anomalies, and panels (b), (d), (f), and (h) show the recovered Vp/Vs anomalies of the checkerboard test. Gray areas mark the regions where the DWS is less than or equal to 5.

Figure 1. (a) Surface geology, (b) main structures, and well locations at LHVC (modified from Carrasco-Núñez, Hernández, et al., 2017; Norini et al., 2015). (c) Locations of the Trans-Mexican Volcanic Belt (TMVB) and LHVC (red triangle).

Figure 17. (a–d) Cross sections for the Vp/Vsmodel (Figure 16). Green circles mark the locations of earthquakes ±200m away from the slice. Dashed red lines indicate the boundary at which spread values are less than or equal to 1.5. Gray areas mark the regions where the DWS is less than or equal to 5. Approximate locations of main structures are indicated in black. Vertical green lines indicate the positions of neighboring injection wells.

Figure 12. Checkerboard recovery. Panels (a), (c), and (e) show three depth slices for the recovered Vp anomalies at −2.6, −2.10, and −1.10 km depths. Panels (b), (d), and (f) show the recovered Vp/Vs anomalies at −2.6, −2.10, and −1.10 km depths. The red and blue squares mark the positions of the synthetic high- and low-velocity anomalies. Gray areas mark the regions where the DWS is less than or equal to 5.

Figure 6. Raypath distribution after the 1-D inversion: (a) map view and (b) N-S and (c) E-W projections. Seismic stations are represented as purple triangles, and local events as green circles. The projections of three injection wells are marked as red lines in the cross sections. Dark solid lines in the map view indicate structures inferred at the surface, and the gray lines correspond to topographic contours.

Figure 2. Topographic map and temporary seismic network at Los Humeros geothermal field. Blue and red triangles mark the positions of three component short-period (Mark L-4C-3D) and three-component broadband (Trillium Compact 120s) sensors, respectively. The reference station for the 1-D inversions (also a three-component broadband Trillium Compact 120s sensor) is marked as a red circle. Several identified and inferred structures are delineated in black (modified from Carrasco-Núñez, Hernández, et al., 2017; Norini et al., 2015).

Figure 16. Vp /Vs structure at different depth levels. Panels (a)–(d) show depth slices for the resulting Vp/Vsmodel at −2.60, −2.10, −1.60, and −1.10 km depths, respectively. Green circles mark the locations of earthquakes ±150m away from the slice. Dashed red lines indicate the boundary at which spread values are less than or equal to 1.5. Gray areas mark the regions where the DWS is less than or equal to 5.

Figure 11. RMS misfit variation for the 228 inverted models.

Figure 7.Model parametrization and DWS distribution at different depth levels for an initial (unrotated) inversion grid. Panels (a) and (c) show two depth slices for the Vpmodel DWS distribution at −2.6 and −1.10 km depths. Panels (b) and (d) show the depth slices for the Vp/Vsmodel DWS distribution at −2.6 and −1.10 km depths. Darker shading indicates regions of higher ray density. Gray crosses indicate node positions.

Figure 3. Distribution of the detected local earthquakes after a nonlinear localization in a homogeneous 3-D volume with a P wave velocity of 3.5 km/s and a Vp/Vs ratio of 1.73. Triangles mark the station positions, and dark solid lines indicate structures inferred at the surface. Red stars mark the positions of three injection wells. C1, C2, and C3 indicate the positions of three main seismic clusters. Depths are defined relative to sea level.

Figure 15. (a–d) Cross sections for the Vpmodel variations (Figure 14). Green circles mark the locations of earthquakes ±200m away from the slice. Dashed red lines indicate the boundary at which spread values are less than or equal to 1.5. Gray areas mark the regions where the DWS is less than or equal to 5. Dashed gray lines indicate different absolute velocity levels, and solid gray lines mark approximate unit boundaries. Approximate locations of main structures are indicated in black. Vertical green lines indicate the positions of neighboring injection wells.

Full text

Local Earthquake Tomography at Los Humeros
Geothermal Field (Mexico)
T. Toledo
1,2
, E. Gaucher
3
, P. Jousset
1
, A. Jentsch
1
, C. Haberland
1
, H. Maurer
4
,
C. Krawczyk
1,2
, M. Calò
5
, and A. Figueroa
6
1
GFZ German Research Center for Geosciences, Potsdam, Germany,
2
Faculty VI – Planning Building Environment:
Institute of Applied Geosciences, TU Berlin, Berlin, Germany ,
3
Institute for Applied Geosciences: Geothermal Energy
and Reservoir Technology, Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany,
4
Institute of Geophysics:
Applied and Environmental Geophysics, ETH Zurich, Zurich, Switzerland,
5
Instituto de Geofísica, UNAM, Mexico City,
Mexico,
6
Instituto de investigaciones en Ciencias de la Tierra, Universidad Michoacana de San Nicolás de Hidalgo,
Morelia, Mexico
Abstract A passive seismic experiment using 25 broadband and 20 short-period stations was
conducted between September 2017 and September 2018 at Los Humeros geothermal field, an important
natural laboratory for superhot geothermal systems in Mexico. From the recorded local seismicity, we
derive a minimum 1-D velocity model and obtain 3-D Vp and Vp/Vs structures of Los Humeros. We
improved the classical local earthquake tomography by using a postprocessing statistical approach. Several
inversions were computed and averaged to reduce artifacts introduced by the model parametrization
and to increase the resolution of the investigated region. Finally, the resulting Vp and Vp/Vs structures
and associated seismicity were integrated with newly acquired geophysical and petrophysical data for
comprehensive interpretation. The recorded seismicity is mainly grouped in three clusters, two of which
seem directly related to exploitation activities. By combining new laboratory measurements and existing
well data with our Vp model, we estimate possible geological unit boundaries. One large intrusion-like
body in the Vp model, together with neighboring high Vp/Vs anomalies, hints at a region of active
resurgence or uplift due to the intrusion of new magma at the northern portion of the geothermal field. We
interpret high Vp/Vs features as fluid bearing regions potentially favorable for further geothermal
exploitation. Deep reaching permeable faults cutting the reservoir unit could explain fluid flow from a
deeper local heat source in the area.
1. 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 elec-
tric 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 num-
ber of studies on the geochemical (e.g., Martinez & Alibert, 1994), geological (e.g., Carrasco-Núñez,
Hernández, et al., 2017; 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., 2001, 2016, 2008; Norini 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
RESEARCH ARTICLE
10.1029/2020JB020390
Key Points:
• High-quality earthquake data were
collected to image the Vp and
Vp/Vs models for the first time
at Los Humeros geothermal field
(Mexico)
• Inversions were performed by
extending the classical earthquake
tomography using a postprocessing
statistical approach
• Geological unit boundaries and
fluid and gas bearing zones were
interpreted considering new
geological, geophysical, and
petrophysical data
Correspondence to:
T. Toledo,
taniat@gfz-potsdam.de
Citation:
Toledo, T., Gaucher, E., Jousset, P.,
Jentsch, A., Haberland, C.,
Maurer, H., et al. (2020). Local
earthquake tomography at Los
Humeros geothermal field (Mexico).
Journal of Geophysical Research:
Solid Earth, 125, e2020JB020390.
https://doi.org/10.1029/2020JB020390
Received 13 JUN 2020
Accepted 26 OCT 2020
Accepted article online 30 NOV 2020
The copyright line for this article was
changed on 18 DEC 2020 after original
online publication.
©2020. The Authors.
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits
use, distribution and reproduction in
any medium, provided the original
work is properly cited.
TOLEDO ET AL. 1of29

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 ade-
quate 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 differ-
ent 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 struc-
tural 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 aver-
age 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 geologi-
cal 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).
2. Geologic and Tectonic Setting
LHVC is a Quaternary geological complex constituted by two nested calderas: the older (ca. 460 kyr) outer
18–20 km wide Los Humeros caldera and the younger (70 kyr) subordinate 5–8 km wide Los Potreros caldera
(Carrasco-Núñez, Hernández, et al., 2017; Carrasco-Núñez, López-Martínez, et al., 2017; Carrasco-Núñez
et al., 2018; Calcagno et al., 2018), where most of the injection and geothermal production activities take
place (Figure 1). An extensive fault network crosses the main production zone of the geothermal field
and is responsible for secondary permeability in the reservoir. Several faults (e.g., Los Humeros fault and
the Loma Blanca fault) favor fluid flow and present strong hydrothermal alteration at the surface (Norini
et al., 2015, 2019). The main fault system runs around 8 km in a NNW-SSE direction and includes the
Maztaloya fault and Los Humeros fault. A second set of faults parting from the main system runs N-S,
NE-SW, and E-W. Both sets disappear at the surface when approaching Los Potreros caldera rim (Figure 1).
From a geological perspective, Los Humeros geothermal field can be divided into four distinct groups:
(1) regional metasedimentary basement, (2) precaldera, (3) caldera, and (4) postcaldera volcanic phases,
which can be subdivided into nine local lithostratigraphic units (Carrasco-Núñez, López-Martínez, et al.,
2017; Calcagno et al., 2018). Here, we briefly describe the lithologies found in the four major groups. The
lower portion of the basement, also called the Teziultlán Massif, is mainly composed of old intrusive igneous
and metamorphic rocks (Paleozoic granites, greenschists, and granodiorites) (Quezadas-Flores, 1961;
Viniegra, 1965; Yáñez & García, 1982). These rocks are covered by an up to 3 km thick Mesozoic sedimen-
tary basement mostly constituted of limestones, with some silts and shales. The basement is overlain by
the precaldera group (10.5–0.155 Ma) mainly composed of andesitic, dacitic, and to a minor extent, basaltic
lavas also known as Teziutlán andesites. The Teziutlán volcanic unit hosts the active geothermal reservoir
and has a thickness larger than 1,500 m in some of the geothermal wells within LHVC (Arellano et al., 2003;
Carrasco-Núñez, Hernández, et al., 2017; Carrasco-Núñez, López-Martínez, et al., 2017; Cedillo-Rodríguez,
1997, 1999; Ferriz & Mahood, 2009; Gutierrez-Negrin & Izquierdo-Montalvo, 2010; Lorenzo-Pulido, 2008;
Norini et al., 2019; Yáñez & García, 1982). The basalts and andesites are sealed above by low-permeability
Quaternary ignimbrites of variable thickness belonging to the caldera stage (Arellano et al., 2003;
Cedillo-Rodríguez, 1997, 1999; Gutierrez-Negrin & Izquierdo-Montalvo, 2010; Lorenzo-Pulido, 2008;
TOLEDO ET AL. 2of29

Journal of Geophysical Research: Solid Earth 10.1029/2020JB020390
Figure 1. (a) Surface geology, (b) main structures, and well locations at LHVC (modified from Carrasco-Núñez, Hernández, et al., 2017; Norini et al., 2015).
(c) Locations of the Trans-Mexican Volcanic Belt (TMVB) and LHVC (red triangle).
Norini et al., 2019). This unit is characterized primarily by eruptive events which resulted in the formation
of Los Humeros and Los Potreros calderas (Carrasco-Núñez & Branney, 2005; Carrasco-Núñez et al., 2012;
Ferriz & Mahood, 2009; Norini et al., 2019). The postcaldera stage (0.05 – < 0.003 Ma) was influenced by
different intracaldera eruptive phases (effusive and explosive). Rhyodacitic, andesitic, and basaltic lavas and
pyroclastic material (Carrasco-Núñez et al., 2018) were produced by various monogenetic volcanic centers
which are scattered between Los Potreros and Los Humeros caldera rims (Norini et al., 2015, 2019). During
that time, another significant eruption took place which resulted in the 1.7 km oval shaped Xalapazco crater
in the south of the complex (Carrasco-Núñez et al., 2018).
3. Seismic Monitoring and Data Processing
3.1. Seismic Network
From September 2017 to September 2018, we deployed and maintained a temporary seismic network com-
prising 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
(Toledo et al., 2019). 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;
TOLEDO ET AL. 3of29

Journal of Geophysical Research: Solid Earth 10.1029/2020JB020390
Figure 2. Topographic map and temporary seismic network at Los Humeros geothermal field. Blue and red triangles
mark the positions of three component short-period (Mark L-4C-3D) and three-component broadband (Trillium
Compact 120s) sensors, respectively. The reference station for the 1-D inversions (also a three-component broadband
Trillium Compact 120s sensor) is marked as a red circle. Several identified and inferred structures are delineated in
black (modified from Carrasco-Núñez, Hernández, et al., 2017; Norini et al., 2015).
Lermo et al., 2001, 2016, 2008; Urban & Lermo, 2013) and where most of the producing and injecting wells
are located. This subnetwork was primarily designed for local microseismicity retrieval (Gaucher et al.,
2019), 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.
3.2. Local Earthquake Detection
We focused the event detection mainly on Los Potreros caldera (Gaucher et al., 2019) 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, 2009) in a homogeneous 3-D volume
TOLEDO ET AL. 4of29

Journal of Geophysical Research: Solid Earth 10.1029/2020JB020390
Figure 3. Distribution of the detected local earthquakes after a nonlinear localization in a homogeneous 3-D volume with a P wave velocity of 3.5 km/s and a
Vp/Vs ratio of 1.73. Triangles mark the station positions, and dark solid lines indicate structures inferred at the surface. Red stars mark the positions of
three injection wells. C1, C2, a nd C3 indicate the positions of three main seismic clusters. Depths are defined relative to sea level.
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.
4. 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 com-
putes 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 k m for deeper levels. The depth intervals
were chosen taking into account well data interpretation (Norini et al., 2015, 2019) and exhaustive testing.
TOLEDO ET AL. 5of29

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "Local earthquake tomography at los humeros geothermal field (mexico)" ?

The authors improved the classical local earthquake tomography by using a postprocessing statistical approach. The authors interpret high Vp/Vs features as fluid bearing regions potentially favorable for further geothermal exploitation.