TY - JOUR
AU - Saenger, E H
AB - SUMMARY The estimation of the source–location accuracy of microseismic events in reservoirs is of significant importance. Time-reverse imaging (TRI) provides a highly accurate localization scheme to locate events by time-reversing the recorded full wavefield and back propagating it through a velocity model. So far, the influence of the station geometry and the velocity model on the source–location accuracy is not well known. Therefore, sensitivity maps are developed using the geothermal site of Los Humeros in Mexico to evaluate the spatial variability of the source–location accuracy. Sensitivity maps are created with an assumed gradient velocity model with a constant vp–vs ratio and with a realistic velocity model for the region of Los Humeros. The positions of 27 stations deployed in Los Humeros from September 2017 to September 2018 are used as surface receivers. An automatic localization scheme is proposed that does not rely on any a priori information about the sources and thus negates any user bias in the source locations. The sensitivity maps are created by simulating numerous uniformly distributed sources simultaneously and locating these sources using TRI. The found source locations are compared to the initial source locations to estimate the achieved accuracy. The resulting sensitivity maps show that the station geometry introduces complex patterns in the spatial variation of accuracy. Furthermore, the influence of the station geometry on the source–location accuracy is larger than the influence of the velocity model. Finally, a microearthquake recorded at the geothermal site of Los Humeros is located to demonstrate the usability of the derived sensitivity maps. This study stresses the importance of optimizing station networks to enhance the accuracy when locating seismic events using TRI. Earthquake location, Numerical modeling, Earthquake source observations, Induced seismicity, Wave propagation 1 INTRODUCTION A crucial part of any reservoir exploration and monitoring scheme is to locate and characterize microseismic events using passive-seismic data recordings, since the spatial and temporal distribution of these events provides insight into geomechanical deformation, subsurface stress changes and fluid mitigation (Shapiro et al. 2002; Duncan & Eisner 2010; Maxwell et al. 2010). Passive-seismic station networks are often optimized for tools such as ray-based methods to locate seismic events, tomographic methods to obtain the velocity structure and noise-based methods to estimate the anisotropy. Recently, migration-based localization methods such as time-reverse imaging (TRI; Fink et al. 1999; Gajewski & Tessmer 2005; Saenger 2011) have been applied to these data sets as well since they promise high-accuracy localizations without the time-intensive identification of phases in the seismic traces. However, recent studies indicated that the station geometry significantly influences the achievable source–location accuracy with TRI (Werner & Saenger 2018; Franczyk 2019). Therefore, a quick test is proposed to evaluate the spatial variability of source–location accuracy by generating sensitivity maps for a specific network and velocity model. Microseismicity is usually characterized by numerous events that occur close in space and time and often occur in regions with elevated ambient seismic noise levels. Therefore, obtaining accurate hypocentral parameters for the majority of events is challenging. Currently, ray-based techniques (Aki & Richards 1980) are the most common methods to obtain reliable and accurate event locations. However, these methods rely on the identification of events and their onset times in the seismic traces and are thus hampered by low signal-to-noise ratios and a high frequency of events occurring in a small area and in a short period of time. In contrast, TRI uses the whole waveform and seismic wave-propagation solvers to locate the events. Thereby, TRI relies on the time invariancy of the wave equation and the source–receiver reciprocity. The time-reversed wavefield is back propagated through an adequate elastic velocity model until it ideally focuses on the original source locations. Therefore, TRI is exceptionally well-suited for events with low signal-to-noise ratios or events that occur close in time and space. The scientific interest for TRI as an alternative localization method has increased significantly in the past decade. Since computational power increased and high-resolution 3-D velocity models became available, TRI has been applied to a number of case studies (Horstmann et al. 2015; Kim & Lees 2015; Kocur et al. 2016; Ge & Han 2018). Very recently, microseismic events have been located using TRI (Li et al. 2019; Yang & Zhu 2019). However, a common challenge is the necessity of a well-constraint 3-D elastic velocity model. Köhn et al. (2016) showed that improving the velocity model also increases the source–location accuracy. The next step in accurately imaging the subsurface is to include viscoelasticity in the velocity structure and the TRI method (Zhu 2015; Oren & Shragge 2019). Another common challenge is the inadequacy of the station network. Attempts at constraining an optimized network revealed that the parameters of the station geometry interact non-linearly (Werner & Saenger 2018; Franczyk 2019), which makes theoretical predictions of the performance of a network for TRI very challenging. The accuracy of a located seismic source represents the spatial and temporal deviation from the true source location. In synthetic studies, all source characteristics are well-constraint. Therefore, synthetic studies enable the systematic investigation of the source–location accuracy in dependence of other factors, such as the station network and velocity model. Sensitivity maps are defined as 3-D volumes containing the source–location accuracy. Sensitivity maps are created by locating numerous synthetic sources simultaneously and evaluating the obtained deviation from the initial sources. Since it is not feasible to identify numerous source localizations by visually inspecting the TRI images, a workflow is developed for the automatic localization of sources. With slight modifications, this workflow could also be applied to locate microearthquakes. In the second part of this study, the workflow is expanded to include numerous sources and develop guidelines for the positioning of the sources. Finally, sensitivity maps are generated for the geothermal site of Los Humeros, Mexico, which was investigated in the scope of the GeMex project (Jolie et al. 2018). The complexity of the sensitivity maps underlines the importance of evaluating the source–location accuracy before locating microearthquakes in the reservoir. Furthermore, the significant influence of the station network on the spatial variation of the source–location accuracy is shown. To demonstrate the usability of the derived sensitivity maps, a microearthquake recorded in Los Humeros is located. 2 TEST SITE: LOS HUMEROS GEOTHERMAL FIELD IN MEXICO Los Humeros is a superhot geothermal system situated in the eastern part of the Trans-Mexican Volcanic Belt. The field is operated by the Comisión Federal de Electridad (CFE). A detailed description of the multicaldera geological complex can be found in Calcagno et al. (2018). Previous studies investigated the seismicity in the area (Lermo et al. 2007, 2016; Urban & Lermo 2013). However, only a limited number of seismic stations were available at that time. To investigate the seismicity in more detail, 45 seismic stations were deployed from September 2017 to September 2018 in Los Humeros as part of the GeMex project (Jolie et al. 2018). The station network is designed as a combination of a dense network composed of 27 stations with a mean interstation distance of |$1.5 \, \textrm{km}$| and a sparse network composed of 23 stations with a mean interstation distance of |$8.5 \, \textrm{km}$| (Toledo et al. 2019). In this study, the locations of the 27 stations of the dense network are used (Fig. 1). This dense network is composed of 17 broad-band stations and 10 short-period stations. These stations are in the centre of the field close to the geothermal wells. All receivers are placed at the free surface of the model since altitude differences are minor compared to the grid spacing and can be neglected. The maximum altitude difference between stations is |$340 \, \textrm{m}$|⁠, which is smaller than the dominant S-wave wavelength (about |$460 \, \textrm{m}$|⁠) in the models. All simulations in this study are performed using a model size of |$20 \, \textrm{km}$| × |$20 \, \textrm{km}$|⁠. The model depth varies depending on the depth of the synthetic sources. The origin point is set to be at the centre of the model (Fig. 1) and corresponds to |$19^\circ \, 40.2^\prime \, \textrm{N}$| and |$97^\circ \, 27^\prime \, \textrm{E}$|⁠. The grid spacing is chosen to be |$20 \, \textrm{m}$|⁠. Figure 1. Open in new tabDownload slide Seismic station network deployed in Los Humeros between September 2017 and September 2018 plotted in the model domain as used in this study. Blue triangles depict the 10 short-period stations and red triangles depict the 17 broad-band stations. The four sources used for demonstrating the method are marked as black circles. Figure 1. Open in new tabDownload slide Seismic station network deployed in Los Humeros between September 2017 and September 2018 plotted in the model domain as used in this study. Blue triangles depict the 10 short-period stations and red triangles depict the 17 broad-band stations. The four sources used for demonstrating the method are marked as black circles. For the creation of sensitivity maps it is crucial to use the same velocity model that will be used to locate genuine microearthquakes since the velocity model itself is expected to influence the source–location accuracy. At the time of this study, the most sophisticated velocity model is a 1-D velocity model (Löer et al. 2020). It is a combination of an S-wave velocity model derived from the analyses of seismic noise and a P- and S-wave model from traveltime tomography (Toledo 2019). The S-wave velocity model is available down to a depth of |$15 \, \textrm{km}$|⁠, while the P-wave velocity model is available to a depth of |$10 \, \textrm{km}$|⁠. However, below a depth of |$6 \, \textrm{km}$| the P-wave velocity is constant. Therefore, in this study, the vp–vsratio is assumed to be constant at depths ranging from |$6$| to |$15 \, \textrm{km}$|⁠. The value of the vp–vs ratio is taken as the average of the vp–vs ratio in the depth above |$6 \, \textrm{km}$|⁠. The P-wave velocity in the depth from 6 to |$15 \, \textrm{km}$| is calculated from the S-wave velocity and the constant vp–vs ratio. Furthermore, the model is interpolated to the used grid spacing and some of the strong velocity contrasts are thereby removed (Fig. 2). Figure 2. Open in new tabDownload slide P-wave velocity vp, S-wave velocity vs and vp–vs ratio for an assumed gradient velocity model with a constant vp–vs ratio and a velocity model based on Löer et al. (2020). Figure 2. Open in new tabDownload slide P-wave velocity vp, S-wave velocity vs and vp–vs ratio for an assumed gradient velocity model with a constant vp–vs ratio and a velocity model based on Löer et al. (2020). 3 AUTOMATIC SOURCE LOCALIZATION USING TRI An automatic and user-independent localization scheme is essential for the efficient and reproducible localization of numerous seismic events. Such a scheme may consist of three parts (cf. Fig. 3): a) The time-reverse simulation of a recorded seismic wavefield and application of imaging conditions to obtain 3-D image volumes. b) The optimization of these images to enhance the contrast between spots of focused energy and the background. c) The identification of focused energy spots and, in the synthetic case, comparison to the original sources. These three parts are described in detail in the following sections. Figure 3. Open in new tabDownload slide Outline of steps used in this study to locate synthetic sources and derive the sensitivity map. Colours represent the corresponding sections. Figure 3. Open in new tabDownload slide Outline of steps used in this study to locate synthetic sources and derive the sensitivity map. Colours represent the corresponding sections. 3.1 Time-reverse imaging A rotated staggered-grid finite-difference scheme (Saenger et al. 2000) is used to propagate the seismic wavefield in the forward and the time-reversed simulations. A finite-difference operator of second order is used in time as well as in space. The sides and bottom of each model are absorbing boundaries (Clayton & Engquist 1980) and the top surface of the model is implemented as a free surface using two layers of vacuum. These vacuum layers mimic the impedance of the solid earth–atmosphere boundary. A typical simulation with a model depth of 6 km consists of 300 million grid points and a simulation with 4000 time steps takes about 6.5 hr on three nodes on a mid-size cluster computer. For all sources, a Ricker wavelet with a central frequency of |$f_\mathrm{ c} = 5 \, \textrm{Hz}$| is used. Sources are implemented as moment tensor sources. Two source types are used: an explosion source (only Mxx = Myy = Mzz ≠ 0) and a strike-slip source (only Mxy ≠ 0). The time increment is set to |$0.001 \, \textrm{s}$| to ensure stability. For demonstration purposes, four synthetic explosion sources are placed in the model (see Fig. 1). Two of the sources are placed at a depth of |$1.8 \, \textrm{km}$| and the other two sources are placed at a depth of |$4.2 \, \textrm{km}$|⁠. Their horizontal positioning is staggered to minimize the interference between sources. The locations of the 27 stations deployed in the geothermal field of Los Humeros are used as receiver positions. A homogeneous velocity model with a P-wave velocity of |$v_\textrm{p} = 4000 \, \rm{m\,s}^{-1}$| is used in this demonstration case. The S-wave velocity is estimated as |$v_\textrm{s} ={v_\textrm{p}}/{\sqrt{3}} = 2309 \, \rm{m\,s}^{-1}$|⁠. The density is derived from the P-wave velocity using an empirical equation from Brocher (2005) to be ρ = 2393 kg m−3. The four sources were excited simultaneously. At the surface of the model, at the receiver positions, the full wavefield is recorded. Very similar example traces can be seen in Werner & Saenger (2018) in fig. 12 and in Saenger (2011) in fig. 11. The complete traces recorded at the surface of the model during the forward simulation are reversed in time, normalized and reinserted at the same positions they were recorded at. Receivers become sources and the time-reversed wavefield is back propagated through the same velocity model used for the forward simulation. Since synthetic data is used, no altering of the traces besides normalization is necessary. Normalizing the traces ensures equal energy inputs at all receivers and therefore enhances the localization results (Kim & Lees 2015). During the time-reversed simulation, imaging conditions are computed to quantify the convergence of the wavefield. When using synthetic explosion sources, the P-wave energy density imaging condition |$\mathbf {I}_\mathrm{ p}$| is used (Saenger 2011): $$\begin{eqnarray*} \mathbf {I}_\mathrm{ p}(\vec{x}) = \max _{t \in [0,T]} (\lambda + 2 \mu ) \left[\nabla \cdot \vec{u}(\vec{x},t) \right]^2. \end{eqnarray*}$$(1) Here, λ and μ represent the Lamé parameters and |$\nabla \cdot \vec{u}(\vec{x},t)$| describes the divergence of the displacement field. When using strike-slip sources, the total energy density imaging condition |$\mathbf {I}_\mathrm{ e}$| is calculated as the sum of the stresses σij times strains εij at every point (Saenger 2011): $$\begin{eqnarray*} \mathbf {I}_\mathrm{ e}(\vec{x}) = \max _{t \in [0,T]} \sum _i \sum _j \left[\sigma _{ij}(\vec{x},t) \varepsilon _{ij}(\vec{x},t) \right] . \end{eqnarray*}$$(2) Furthermore, the time steps |$t_{\mathrm{ max}} (\vec{x})$|⁠, when the maximum amplitudes in |$\mathbf {I}_\mathrm{ e}$| or |$\mathbf {I}_\mathrm{ p}$| occurs at each grid point, are observed. While clusters of high amplitudes in the imaging conditions itself represent focusing spots and may reveal the location of seismic sources, tmax at the same point in the time volume reveals the time the wavefield converged at this point and can later be used to determine the temporal error of the localization by comparing it to the excitation time of the initial sources. In |$\mathbf {I}_\mathrm{ p}$| as obtained during the time-reversed simulation (Fig. 4), the position of the four initial explosion sources are marked as red circles. In addition to the amplitude of |$\mathbf {I}_\mathrm{ p}$| in the 2-D depth slices, the maximum amplitude in each of the three dimensions is plotted in Figs 4(d)–(f). The highest amplitudes are observed close to the surface. However, at a depth of |$1.8 \, \textrm{km}$| a small increase is visible in the maximum amplitude against depth (Fig. 4d). This depth relates to the depth of the shallower two sources. There is no peak visible at the second source depth. The maximum amplitude in east and north direction seems to be dominated by the high amplitudes at the surface and cannot be related to the source positions but rather to the station positions. Figure 4. Open in new tabDownload slide P-wave energy density Ip after back propagating the time-reversed wavefield (eq. 1). The four explosion sources are marked with red circles in (b) and (c) and red dashed lines in the profiles of maximum amplitude (d)–(f). The high amplitudes at the surface visible in (a) are dominating the depth profile in (d). This is also the case for the two horizontal profiles (e) and (f). The peaks are induced by the stations at the surface. In (b) and (c), the four sources can be visually identified. Figure 4. Open in new tabDownload slide P-wave energy density Ip after back propagating the time-reversed wavefield (eq. 1). The four explosion sources are marked with red circles in (b) and (c) and red dashed lines in the profiles of maximum amplitude (d)–(f). The high amplitudes at the surface visible in (a) are dominating the depth profile in (d). This is also the case for the two horizontal profiles (e) and (f). The peaks are induced by the stations at the surface. In (b) and (c), the four sources can be visually identified. 3.2 Image processing The four sources are visually distinguishable in Figs 4(b) and (c). However, high amplitudes at the surface close to the stations hinder an automatic localization. Due to geometrical spreading, the wavefield is attenuated in depth and the amplitudes close to the stations relatively higher. Additionally, these high amplitudes are caused by the generation of surface waves. Additionally, energy can accumulate in parts of the model solely due to the velocity structure. Witten & Artman (2011) introduced the correction of these artefacts by dividing the TRI image by an image obtained during the propagation of a noise wavefield from the same surface receivers through the same velocity model. The noise wavefield highlights areas where energy accumulates without converging on source locations. Therefore, the image obtained with a noise wavefield is called illumination map in this study. To have a noise wavefield with the same amplitude and frequency content as the used data, the non-time-reversed wavefield is used to create the illumination map. Here, non-time-reversed traces mean the traces as recorded by the stations without time-reversing them. The same set-up used for the time-reversed simulation is used again with non-time-reversed traces and a significantly longer runtime. The runtime was chosen to ensure that the whole model is illuminated. The results from the time-reversed simulation are then divided by these illumination maps to dampen areas of high amplitudes. In |$\mathbf {I}_\mathrm{ p}$| after division by the illumination map (Fig. 5), the source locations appear more distinct while the high amplitudes at the surface have vanished. The maximum amplitude in east and north directions (Figs 5d–f) shows peaks corresponding to the positions of the four sources. After division by the illumination map, the sources can be identified clearly in the 2-D depth slices of the imaging conditions as well as the 1-D plots of the maximum amplitude. However, it is not possible to pick localizations from the 1-D plots without any a priori information about the sources. From the 1-D plots, it would also be possible that there are eight sources, four in each depth. Figure 5. Open in new tabDownload slide Ip after division by the illumination map. The four explosion sources are marked with red circles in (b) and (c) and red dashed lines in the profiles of maximum amplitude (d)–(f). The two source depths are visible as peaks in the depth profile in (d). The horizontal coordinates of the sources are visible in the horizontal profiles (e) and (f). The sources dominate compared to Fig. 4 where the surface stations dominate. The vertical dashed blue line in (d)–(f) marks the threshold Th (eq. 3) used for eliminating background noise from the imaging conditions. In (b) and (c), the four sources can be visually identified. Figure 5. Open in new tabDownload slide Ip after division by the illumination map. The four explosion sources are marked with red circles in (b) and (c) and red dashed lines in the profiles of maximum amplitude (d)–(f). The two source depths are visible as peaks in the depth profile in (d). The horizontal coordinates of the sources are visible in the horizontal profiles (e) and (f). The sources dominate compared to Fig. 4 where the surface stations dominate. The vertical dashed blue line in (d)–(f) marks the threshold Th (eq. 3) used for eliminating background noise from the imaging conditions. In (b) and (c), the four sources can be visually identified. After the division by the illumination map, some background noise remains visible in the images. Therefore, the images need to be cleaned to be able to automatically identify the source locations. All points with an amplitude beneath a threshold are muted from the 3-D volume. The threshold Th is found by calculating the mean of the maximum amplitude for all three directions and then taking the mean value from those: $$\begin{eqnarray*} T_{\mathrm{ h},\mathrm{ k}} = \frac{1}{3} \sum _{i = 1}^3 \left( \frac{1}{n_i} \sum _{j = 1}^{n_i} \mathbf {I}_{\mathrm{ k}}^{\mathrm{ max}}\right) . \end{eqnarray*}$$(3) The index i represents the three coordinate directions and ni is the number of sample points in each of the three directions. The index k represents the type of imaging condition. Muting everything below the threshold Th will remove the lower amplitudes in the background but retain the peaks corresponding to the source locations. After eliminating the background noise, the only non-zero values are close to the original source positions in this four-source case. Fig. 6 shows the four recovered explosion sources in more detail after eliminating the background noise. The grey-scale background shows the normalized amplitude of |$\mathbf {I}_\mathrm{ p}$| with everything below the threshold set to zero amplitude. Since explosion sources were simulated, the expected shape of the located sources would be spherical. However, the localizations have different shapes and sizes which may be a result of the station distribution since the velocity model is homogeneous and all sources were equal in strength. Additionally, source number 2 (Fig. 6b) is hardly visible in the grey-scale image. Figure 6. Open in new tabDownload slide Close-up view of depth slices of Ip of four explosion sources after the division by the illumination map and elimination of background noise. (a)–(d) show sources one to four as labelled in Fig. 1. The original source positions are marked as red circles and are in the depth of the shown slices. The localizations found in the vicinity of the sources are marked as blue circles with larger circles representing localizations closer to the original source location. The localization chosen from the algorithm to fit to the source is marked in orange and the depth of this localization is given. The S-wave wavelength is plotted for scale. The details of the orange localizations can be found in Table 1. Figure 6. Open in new tabDownload slide Close-up view of depth slices of Ip of four explosion sources after the division by the illumination map and elimination of background noise. (a)–(d) show sources one to four as labelled in Fig. 1. The original source positions are marked as red circles and are in the depth of the shown slices. The localizations found in the vicinity of the sources are marked as blue circles with larger circles representing localizations closer to the original source location. The localization chosen from the algorithm to fit to the source is marked in orange and the depth of this localization is given. The S-wave wavelength is plotted for scale. The details of the orange localizations can be found in Table 1. 3.3 Obtaining source localizations The automatic identification of localizations is done using the normalized 3-D imaging condition after the division by the illumination map and after the elimination of the background noise (cf. Fig. 6). To identify all localizations, the point with the absolute maximum amplitude in the whole 3-D imaging condition is found and a sphere with a radius of one S-wave wavelength around this point is muted. Afterwards, the point with the new absolute maximum amplitude is found and the process repeated until the whole model has zero amplitudes. This iterative approach enables to find localizations independent of the number of initial sources and thus does not imply any a priori knowledge. During the iterative process, the focus size is determined as the number of points inside the sphere with non-zero amplitude. The focus size is representative of the rate of convergence of the wavefield and may be used as an estimate for the uncertainty of the localization (Li et al. 2019). The maximum focus size in the test case with the four sources was about |$2 \times 10^4 \, \textrm{grid~points}$|⁠. Assuming spherical localizations, this would relate to a sphere with a radius of about |$340 \, \textrm{m}$| (Table 1), which is significantly smaller than the S-wave wavelength of |$460 \, \textrm{m}$| for this case with the homogeneous velocity model. This confirms that using the S-wave wavelength results in spheres large enough to cover the localizations. However, some localizations may be elongated in one direction and therefore deviate from a sphere. These localizations may be selected multiple times by the algorithm, which can be seen in Fig. 6. The blue circles represent all initially found localizations after applying some general filters to the localizations. The larger circles represent localizations closer in depth to the depth slice shown in the grey-scale background image. Table 1. Characteristics of the source localizations found for the four explosion sources using a homogeneous velocity model. . Initial source . Localization . Amplitude . Focus . Temporal . Spatial . . (x, y, z) (km) . (x, y, z) (km) . (-) . radius (m) . error (ms) . error (m) . source 1 −3.6, 3.6, 1.8 −3.58, 3.58, 1.74 0.99 296.1 9 66.3 source 2 3.6, −3.6, 1.8 3.78, −3.92, 2.08 0.88 254.7 102 461.7 source 3 3.6, 3.6, 4.2 3.64, 3.62, 4.24 0.88 335.3 16 60 source 4 −3.6, −3.6, 4.2 −3.5, −3.56, 4.06 1 347.1 34 176.6 . Initial source . Localization . Amplitude . Focus . Temporal . Spatial . . (x, y, z) (km) . (x, y, z) (km) . (-) . radius (m) . error (ms) . error (m) . source 1 −3.6, 3.6, 1.8 −3.58, 3.58, 1.74 0.99 296.1 9 66.3 source 2 3.6, −3.6, 1.8 3.78, −3.92, 2.08 0.88 254.7 102 461.7 source 3 3.6, 3.6, 4.2 3.64, 3.62, 4.24 0.88 335.3 16 60 source 4 −3.6, −3.6, 4.2 −3.5, −3.56, 4.06 1 347.1 34 176.6 Open in new tab Table 1. Characteristics of the source localizations found for the four explosion sources using a homogeneous velocity model. . Initial source . Localization . Amplitude . Focus . Temporal . Spatial . . (x, y, z) (km) . (x, y, z) (km) . (-) . radius (m) . error (ms) . error (m) . source 1 −3.6, 3.6, 1.8 −3.58, 3.58, 1.74 0.99 296.1 9 66.3 source 2 3.6, −3.6, 1.8 3.78, −3.92, 2.08 0.88 254.7 102 461.7 source 3 3.6, 3.6, 4.2 3.64, 3.62, 4.24 0.88 335.3 16 60 source 4 −3.6, −3.6, 4.2 −3.5, −3.56, 4.06 1 347.1 34 176.6 . Initial source . Localization . Amplitude . Focus . Temporal . Spatial . . (x, y, z) (km) . (x, y, z) (km) . (-) . radius (m) . error (ms) . error (m) . source 1 −3.6, 3.6, 1.8 −3.58, 3.58, 1.74 0.99 296.1 9 66.3 source 2 3.6, −3.6, 1.8 3.78, −3.92, 2.08 0.88 254.7 102 461.7 source 3 3.6, 3.6, 4.2 3.64, 3.62, 4.24 0.88 335.3 16 60 source 4 −3.6, −3.6, 4.2 −3.5, −3.56, 4.06 1 347.1 34 176.6 Open in new tab The found localizations are filtered using a maximum allowed spatial error and a maximum allowed temporal error based on the initial sources. The maximum allowed spatial error is defined as half the distance to the next source when using multiple sources. In this test case with only four sources, the maximum allowed spatial error is set to about four S-wave wavelengths. The maximum allowed temporal error is set to one period |$T = {1}/{f_\mathrm{ c}} = {1}/{5 \, \textrm{Hz}} = 0.2 \, \textrm{s}$|⁠. When locating events in the field data, this filtering of localizations may be performed manually using the 3-D images of the imaging condition. From the blue localizations in Fig. 6, the one with the largest imaging-condition amplitude in the vicinity of each source is chosen to be the localization for that specific source. If there are no localizations in the vicinity of a source, this source could not be located. This removes multiple localizations of the same source and removes artefacts remaining from the elimination of the background. Choosing the localization with the largest amplitude results in the algorithm choosing the centre of the localization. Since the intersource distance is set during the forward simulation, the range in which to look for localizations is adapted accordingly. For genuine microearthquakes, this would need to be performed more carefully. In our test case with four explosion sources, there were initially 38 localizations found. After eliminating those with a large temporal and spatial error, 28 localizations remained. All of the four sources could be paired with one of the localizations. The four localizations that were associated to the sources were the four with the highest overall amplitude of all 38 localizations, the smallest temporal error of all localizations and the largest focus radius. In Table 1, the characteristics of the four source localizations are given. As already seen, source number 2 has the largest temporal error and the largest spatial error of all the sources. However, the focus radius is comparably small. A quite small focus radius, close to half the S-wave wavelength therefore may indicate that not enough information about the wavefield is present at this position. Therefore, the focus radius may not be used directly as the uncertainty of the localization. 4 LOCATING NUMEROUS SYNTHETIC SOURCES SIMULTANEOUSLY To obtain the source–location accuracy in the whole model, numerous synthetic sources need to be distributed in the model space. This significantly reduces the simulation time compared to testing each source location in an individual simulation. However, the source–location accuracy must not be influenced by the source distribution and the number of sources. In this section, different source distributions are investigated using the same homogeneous velocity as before. Using a homogeneous velocity model ensures that the described changes in the source–location accuracy are due to the source distribution and the receiver distribution. As receivers, the locations of the 27 stations deployed in Los Humeros are used. All tests were performed once with explosion sources and once with strike-slip sources. However, results are only shown for the explosion sources since the strike-slip sources resulted in the same guidelines for the source distribution. 4.1 Optimal intersource distance All sources are placed with equal intersource distances in north and east directions to allow for a fast set-up of simulations. The optimal intersource distance was found to be |$2.4 \, \textrm{km}$|⁠, which relates to slightly more than five S-wave wavelengths (⁠|$5 \times 460 \, \textrm{m} = 2.3 \, \textrm{km}$|⁠) and three P-wave wavelengths (⁠|$3 \times 800 \, \textrm{m} = 2.4 \, \textrm{km}$|⁠). This intersource distance is also slightly larger than the maximum interstation distance of the network which is about |$2.1 \, \textrm{km}$|⁠. The focus radius, spatial error and temporal error for the found localizations of 36 uniformly spaced sources with an intersource distance of |$2.4 \, \textrm{km}$| can be seen in Figs 7(a)–(c). The homogeneous velocity model and all of the short-period and broad-band stations were used. The localizations are plotted at the position they were located at. The original source positions are marked as grey circles. As expected, sources beneath the network can be located while the sources outside the network are located with a larger error or not at all. The focus radius of the explosion sources is about 115 to |$320 \, \textrm{m}$| large, which relates to a fourth to a half of the S-wave wavelength. In Fig. 7(d), the results are sorted into three categories: (1) High-accuracy source locations have a spatial error smaller than the focus radius and therefore these localizations cover the initial source locations. Furthermore, the spatial error is smaller than half of the S-wave wavelength, smaller than |$230 \, \textrm{m}$|⁠, and the temporal error is smaller than half of the period, smaller than |$0.1 \, \textrm{s}$|⁠. (2) Medium-accuracy localizations are all localizations that do not fit the criteria of the high-accuracy localizations. (3) The categories marked in red are sources that cannot be located with the used set-up. Figure 7. Open in new tabDownload slide Focus radius, spatial error and temporal error for sources in a depth of |$1.8 \, \textrm{km}$| when using a homogeneous velocity model and all of the short-period and broad-band stations. Initial source locations are shown as open circles and the localizations are filled circles at the position of their localization. Different intersource distances and number of sources are compared. The right column shows the achieved accuracy for each of the initial sources. Localizations with a smaller spatial error than focus radius and with a spatial error smaller than |$230 \, \textrm{m}$| and a temporal error smaller than |$0.1 \, \textrm{s}$| are categorized as high accuracy. Sources that could not be located are marked in red. See Fig. 1 for the station geometry. Figure 7. Open in new tabDownload slide Focus radius, spatial error and temporal error for sources in a depth of |$1.8 \, \textrm{km}$| when using a homogeneous velocity model and all of the short-period and broad-band stations. Initial source locations are shown as open circles and the localizations are filled circles at the position of their localization. Different intersource distances and number of sources are compared. The right column shows the achieved accuracy for each of the initial sources. Localizations with a smaller spatial error than focus radius and with a spatial error smaller than |$230 \, \textrm{m}$| and a temporal error smaller than |$0.1 \, \textrm{s}$| are categorized as high accuracy. Sources that could not be located are marked in red. See Fig. 1 for the station geometry. Decreasing the intersource distance (Figs 7e–l) results in some sources not being located, although the general area where sources are located stays the same. Intersource distances of |$1.8$| and |$2.4 \, \textrm{km}$| seem to give similar results. However, some sources beneath the network can no longer be located. Intersource distances of |$1.2 \, \textrm{km}$| do not yield satisfying results. The source locations seem to interfere significantly with each other hampering the location of some sources. The observed decrease in location accuracy with a decrease in the intersource distance may be due to the increased number of stations. To keep the area covered with sources similar, the number of sources was increased when decreasing the intersource distance. In Figs 7(m)–(p), the localization results of 64 sources with an intersource distance of |$2.4 \, \textrm{km}$| is shown. These result look more similar to Figs 7(a) and (e), where the intersource distance is the same, than to Figs 7(e)–(h), where the number of sources is the same. Therefore, the influence of the intersource distance seems to be more significant than the influence of the total number of sources. However, Figs 7(m)–(p) also have some sources that cannot be located directly beneath the network. This may be due to absorbing boundaries that are not capable of absorbing the complete wavefield and therefore introduce higher amplitudes at the sides of the model. The sources in the centre experience smaller amplitudes and may not be located because their amplitude is too small in comparison to the threshold used for eliminating the background noise from the images. To investigate whether sources need to be placed |$2.4 \, \textrm{km}$| apart due to the wavelength or the interstation distance, additional simulations were performed only with the 17 broad-band stations (Fig. 8). The largest interstation distance for the broad-band stations is about |$3.9 \, \textrm{km}$|⁠. Figs 8(a)–(d) show the results of 36 sources with an intersource distance of |$2.4 \, \textrm{km}$|⁠. There are some sources that are not located although they are underneath the network and were expected to be located. Figs 8(e)–(h) show the results of locating 25 sources with an intersource distance of |$4 \, \textrm{km}$|⁠, which is larger than the maximum interstation distance. Here, all sources underneath the network can be located. However, the sources are quite far apart. Therefore, we conclude that the previously found optimal intersource distance of |$2.4 \, \textrm{km}$| for the network with all of the short-period and broad-band stations may be due to the maximum interstation distance of the network. Figure 8. Open in new tabDownload slide Using only the 17 broad-band stations, the focus radius, spatial error and temporal error for sources in a depth of |$1.8 \, \textrm{km}$| when using a homogeneous velocity model are shown. Initial source locations are shown as open circles and the localizations are filled circles at the position of their localization. Different intersource distances and number of sources are compared. The right-hand column shows the achieved accuracy for each of the initial sources. Localizations with a smaller spatial error than focus radius and with a spatial error smaller than |$230 \, \textrm{m}$| and a temporal error smaller than |$0.1 \, \textrm{s}$| are categorized as high accuracy. Sources that could not be located are marked in red. See Fig. 1 for the station geometry. Figure 8. Open in new tabDownload slide Using only the 17 broad-band stations, the focus radius, spatial error and temporal error for sources in a depth of |$1.8 \, \textrm{km}$| when using a homogeneous velocity model are shown. Initial source locations are shown as open circles and the localizations are filled circles at the position of their localization. Different intersource distances and number of sources are compared. The right-hand column shows the achieved accuracy for each of the initial sources. Localizations with a smaller spatial error than focus radius and with a spatial error smaller than |$230 \, \textrm{m}$| and a temporal error smaller than |$0.1 \, \textrm{s}$| are categorized as high accuracy. Sources that could not be located are marked in red. See Fig. 1 for the station geometry. 4.2 Investigate each source depth individually Previous simulations were performed with numerous sources in one source depth only. However, the four-source test case showed that sources can be located simultaneously in multiple depths. Figs 9(a)–(d) show the result of 36 explosion sources in a depth of |$1.8 \, \textrm{km}$| and Figs 9(e)–(h) show 36 explosion sources in a depth of |$4.2 \, \textrm{km}$|⁠. Those sources were located using two simulations. One simulations with sources only in the shallower depth and one simulation with sources only in the larger depth. Figure 9. Open in new tabDownload slide Focus radius, spatial error and temporal error for sources in a depth of |$1.8\,$| and |$4.2 \, \textrm{km}$| are shown when locating sources in one depth per simulation (a–h) or in two depths simultaneously (i–p). A homogeneous velocity model and all of the short-period and broad-band stations are used. Initial source locations are shown as open circles and the localizations are filled circles at the position of their localization. Different intersource distances and number of sources are compared. The right-hand column shows the achieved accuracy for each of the initial sources. Localizations with a smaller spatial error than focus radius and with a spatial error smaller than |$230 \, \textrm{m}$| and a temporal error smaller than |$0.1 \, \textrm{s}$| are categorized as high accuracy. Sources that could not be located are marked in red. See Fig. 1 for the station geometry. Figure 9. Open in new tabDownload slide Focus radius, spatial error and temporal error for sources in a depth of |$1.8\,$| and |$4.2 \, \textrm{km}$| are shown when locating sources in one depth per simulation (a–h) or in two depths simultaneously (i–p). A homogeneous velocity model and all of the short-period and broad-band stations are used. Initial source locations are shown as open circles and the localizations are filled circles at the position of their localization. Different intersource distances and number of sources are compared. The right-hand column shows the achieved accuracy for each of the initial sources. Localizations with a smaller spatial error than focus radius and with a spatial error smaller than |$230 \, \textrm{m}$| and a temporal error smaller than |$0.1 \, \textrm{s}$| are categorized as high accuracy. Sources that could not be located are marked in red. See Fig. 1 for the station geometry. Figs 9(i)–(p) show the same sources but all of the 72 sources located in one single simulation. Most of the deeper sources cannot be located and some of the shallower sources can also not be located. The focus radius of the deeper sources is small. This indicates that not enough information of the wavefield reaches the larger depth. The reason could be that there are too many sources that overlap in the recorded traces. Alternatively, the very regular distribution of sources could cause an ambiguity in the source location that results in deeper sources being located at the shallower source depth. Additionally, the overall amplitude decreases with depth due to geometrical spreading and the amplitude of the deeper sources may fall below the threshold. Simulating each depth individually has the additional advantage that the depth can be sampled with a higher resolution than the interstation distance generally allows. Therefore, it is recommended to create sensitivity maps for sources in each desired depth individually and then combine them to a 3-D image afterwards. An alternative method to circumvent the need for multiple simulations is to excite the sources in each depth at different times. This staggering in time may separate the recorded source signals and allow the location of sources in each depth at different simulation times permitting to use different thresholds at different times. However, tests showed that the shallow sources will be located more easily since the amplitude generally decreases with depth. Therefore, staggering the sources in time potentially still underestimates the source–location accuracy in depth. 5 SENSITIVITY MAPS FOR LOS HUMEROS Sensitivity maps are created for the geothermal field of Los Humeros to estimate the source–location accuracy in the reservoir. This can be used to identify if TRI is a suitable method for locating events close to known faults or geothermal wells. As receiver positions, the locations of the 27 stations deployed in Los Humeros are used. The velocity model based on Löer et al. (2020) is used as shown in Fig. 2. Additionally, a gradient velocity model is assumed as a first approximation (Fig. 2) that is more realistic than a homogeneous velocity model. 36 sources were spaced |$2.4 \, \textrm{km}$| apart and in depths of 1 to |$12 \, \textrm{km}$|⁠. Explosion sources are used with both velocity models and strike-slip sources are used with the realistic velocity model. The accuracies of the source localizations were sorted into the same three categories as in the previous section. The simulations for the individual depth slices were put together and interpolated in-between the source locations to obtain the 3-D sensitivity map. The sensitivity map for the gradient velocity model and explosion sources is shown in Fig. 10. High-accuracy localizations are found beneath the network and in depths of up to |$8 \, \textrm{km}$|⁠. The medium accuracy localizations are found in a much wider area covering almost the whole space beneath the network. With increasing depth, the area is decreasing. However, there are small patches where the accuracy is decreased. Since a gradient velocity model is used here, these complex pattern of accuracies can be assumed to be caused by the station network. Figure 10. Open in new tabDownload slide Sensitivity map for the Los Humeros geothermal field using the gradient velocity model and all of the short-period and broad-band stations. There are 36 explosion sources in each of the nine depths: 1, 2, 3, 4, 5, 6, 7, 8, 10 and |$12 \, \textrm{km}$|⁠. The colour scheme shows the categorization of the localization results into high accuracy (spatial error is smaller than focus radius, spatial error is smaller than |$230 \, \textrm{m}$| and temporal error is smaller than |$0.1 \, \textrm{s}$|⁠), medium accuracy and sources that could not be located. Figure 10. Open in new tabDownload slide Sensitivity map for the Los Humeros geothermal field using the gradient velocity model and all of the short-period and broad-band stations. There are 36 explosion sources in each of the nine depths: 1, 2, 3, 4, 5, 6, 7, 8, 10 and |$12 \, \textrm{km}$|⁠. The colour scheme shows the categorization of the localization results into high accuracy (spatial error is smaller than focus radius, spatial error is smaller than |$230 \, \textrm{m}$| and temporal error is smaller than |$0.1 \, \textrm{s}$|⁠), medium accuracy and sources that could not be located. The sensitivity map for the realistic velocity model and explosion sources can be seen in Fig. 11. The general distribution of high and medium accuracies is smaller than the sensitivity map with the gradient velocity model. The area of high accuracies at the shallowest investigated depth seems larger than for the gradient velocity model since the velocities in the gradient model are larger in the shallowest part of the model. Furthermore, at a depth of about |$6 \, \textrm{km}$|⁠, the accuracies seem to decrease quite rapidly. Figure 11. Open in new tabDownload slide Sensitivity map for the Los Humeros geothermal field using the velocity model based on Löer et al. (2020) and all of the short-period and broad-band stations. There are 36 explosion sources in each of the eleven depths: 1, 2, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 10 and |$12 \, \textrm{km}$|⁠. The colour scheme shows the categorization of the localization results into high accuracy (spatial error is smaller than focus radius, spatial error is smaller than |$230 \, \textrm{m}$| and temporal error is smaller than |$0.1 \, \textrm{s}$|⁠), medium accuracy and sources that could not be located. Figure 11. Open in new tabDownload slide Sensitivity map for the Los Humeros geothermal field using the velocity model based on Löer et al. (2020) and all of the short-period and broad-band stations. There are 36 explosion sources in each of the eleven depths: 1, 2, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 10 and |$12 \, \textrm{km}$|⁠. The colour scheme shows the categorization of the localization results into high accuracy (spatial error is smaller than focus radius, spatial error is smaller than |$230 \, \textrm{m}$| and temporal error is smaller than |$0.1 \, \textrm{s}$|⁠), medium accuracy and sources that could not be located. The sensitivity map for strike-slip sources with the realistic velocity model is shown in Fig. 12. The distribution of accuracies appears to be different than for explosion sources. This is expected since the two source types have very different radiation patterns. Sources can only be located to a depth of about |$5 \, \textrm{km}$|⁠. If a detailed interpretation is done in context with, for example, well locations, multiple source types should be considered. However, since seismic events are rarely of one pure source types, the number of sensitivity maps could be endless when considering all source types. Therefore, we propose to use sensitivity maps with explosion sources and keep in mind that other source types may be located in parts where explosion sources cannot and vice versa. Figure 12. Open in new tabDownload slide Sensitivity map for the Los Humeros geothermal field using the velocity model based on Löer et al. (2020) and all of the short-period and broad-band stations. There are 36 strike-slip sources in each of the eleven depths: 1, 2, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 10 and |$12 \, \textrm{km}$|⁠. The colour scheme shows the categorization of the localization results into high accuracy (spatial error is smaller than focus radius, spatial error is smaller than |$230 \, \textrm{m}$| and temporal error is smaller than |$0.1 \, \textrm{s}$|⁠), medium accuracy and sources that could not be located. Figure 12. Open in new tabDownload slide Sensitivity map for the Los Humeros geothermal field using the velocity model based on Löer et al. (2020) and all of the short-period and broad-band stations. There are 36 strike-slip sources in each of the eleven depths: 1, 2, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 10 and |$12 \, \textrm{km}$|⁠. The colour scheme shows the categorization of the localization results into high accuracy (spatial error is smaller than focus radius, spatial error is smaller than |$230 \, \textrm{m}$| and temporal error is smaller than |$0.1 \, \textrm{s}$|⁠), medium accuracy and sources that could not be located. The number of sources located in each depth for the three sensitivity maps (Fig. 13a) is compared with the vp–vs ratio of the two velocity models. From a depth of 4 to |$5 \, \textrm{km}$|⁠, the vp–vs ratio significantly increases for the realistic velocity model. At this depth, the number of sources located with the realistic velocity model decreases. This effect cannot be seen when viewing only the number of high-accuracy localizations (Fig. 13b) and the mean spatial error and mean focus radius (Figs 13c and d). However, the mean temporal error (Fig. 13e) decreases slightly in this depth for the realistic velocity model compared to the gradient velocity model. The focus radius of the strike-slip sources is generally smaller than that of the explosion sources (Fig. 13d) and the number of high-accuracy localizations decreases rapidly at a depth of |$5 \, \textrm{km}$|⁠. Otherwise, no significant difference between the two source types can be identified. Figure 13. Open in new tabDownload slide Depth profiles of (a) the number of successfully located sources, (b) the number of localizations with a high accuracy, (c) the mean spatial error, (d) the mean focus radius and (e) the mean temporal error. The vp–vs ratio of the gradient velocity model and the velocity model based on Löer et al. (2020) is plotted for reference. Strike-slip sources are shown in red and explosion sources in black. Results produced with the gradient velocity model are shown with a dashed line and results with the velocity model based on Löer et al. (2020) are shown as a solid line. Figure 13. Open in new tabDownload slide Depth profiles of (a) the number of successfully located sources, (b) the number of localizations with a high accuracy, (c) the mean spatial error, (d) the mean focus radius and (e) the mean temporal error. The vp–vs ratio of the gradient velocity model and the velocity model based on Löer et al. (2020) is plotted for reference. Strike-slip sources are shown in red and explosion sources in black. Results produced with the gradient velocity model are shown with a dashed line and results with the velocity model based on Löer et al. (2020) are shown as a solid line. 6 LOCATION OF A MICROEARTHQUAKE USING THE DERIVED SENSITIVITY MAP The derived sensitivity maps may be used to expand the workflow described for locating synthetic sources to locate events in the field data set recorded in Los Humeros (Toledo et al. 2019). Specifically, instead of using the synthetic waveforms as an input, the recorded waveforms are resampled to fit the time increment, filtered between 2 and 12 Hz, normalized and time-reversed. Filtering was done using a fourth-order Butterworth filter twice to retain all phase information. The selected frequency range is chosen to keep the finite-difference simulation accurate with respect to numerical dispersion. Higher frequencies are possible if the grid spacing is reduced and thus results in increased computation times. By limiting the frequencies to the lower range, the lower-frequency part of the microearthquake is located. A prominent low-frequency content is not uncommon for shallow induced earthquakes (Cesca et al. 2011). Future studies will exploit the use of multiple frequency bands. The pre-processed waveforms were inserted into the model and the workflow (Fig. 3) was followed. After dividing the TRI images by the illumination map and before applying thresholding, one of the derived sensitivity maps (Figs 11 and 12) is used as a mask. The areas that are red and yellow in the sensitivity map are muted in the TRI images. The area with medium accuracy (yellow) is muted here as well due to the high level of scattering but may be included in other cases. The result is a list of points representing clusters of high amplitudes. These need to be evaluated manually. To demonstrate this workflow, a microearthquake recorded on 2018 April 22 at around 21:03:50 (Fig. 14a) is located using the steps described above. At the time of the event, 22 of the 27 available stations were actively recording (Fig. 14b). Since the geology in Los Humeros is very heterogeneous, the 1-D velocity model based on Löer et al. (2020) is too simplified to represent the subsurface correctly. In combination with the sparse station network, this leads to an abundance of artificial focusing spots when trying to locate events. Adding the sensitivity map to the workflow, mutes artificial focusing spots in areas of the model where no localization is possible. The source locations are, therefore, more prominently visible and can be picked more easily from the TRI images. Since mostly shear dislocations are expected in this context, the sensitivity map derived with synthetic strike-slip sources (Fig. 12) is used. This sensitivity map was produced with the total energy density imaging condition Ie but may be used for all imaging conditions since all of them image the same earthquake. The importance is to use a sensitivity map produced with a source type that is expected to occur in the reservoir to correctly account for the influence of the radiation pattern. Figure 14. Open in new tabDownload slide Localization of a microearthquake that occurred on 2018 April 22. (a) Example time-reversed seismic traces recorded at stations DB02 and DB25 after pre-processing. The traces recorded at all available stations are used as input for the localization. Arrows point to the corresponding stations of the shown example traces. (b) Localizations in Ie, Ip and Is plotted together with used station network. (c) Cleaned and normalized TRI image of Ie with the determined source location. Figure 14. Open in new tabDownload slide Localization of a microearthquake that occurred on 2018 April 22. (a) Example time-reversed seismic traces recorded at stations DB02 and DB25 after pre-processing. The traces recorded at all available stations are used as input for the localization. Arrows point to the corresponding stations of the shown example traces. (b) Localizations in Ie, Ip and Is plotted together with used station network. (c) Cleaned and normalized TRI image of Ie with the determined source location. To locate the microearthquake three imaging conditions are used: the total energy density imaging condition Ie (cf. eq. 2), the P-wave energy density imaging condition Ip (cf. eq. 1) and the S-wave energy density imaging condition Is which is sensitivity to the shear wavefield (Saenger 2011): $$\begin{eqnarray*} \mathbf {I}_\mathrm{ s}(\vec{x}) = \max _{t \in [0,T]} \mu \left[\nabla \times \vec{u}(\vec{x},t)\right]^2 . \end{eqnarray*}$$(4) Here, |$\nabla \times \vec{u}(\vec{x},t)$| denotes the curl of the displacement field, the displacement of the shear wavefield. The automatic search for clusters of high amplitudes found 166 localizations in the total energy density imaging condition Ie. However, manual inspection of these localizations revealed that most of them are small in size and amplitude. Four localizations with a normalized amplitude greater than 0.9 are identified in Ie (see red circles in Fig. 14b). Similarly, two localizations are identified in Is (see orange squares Fig. 14). The localizations with a smaller amplitude that are excluded here, are considered to be artificial or secondary scatterers and may be suppressed further with an updated velocity model. The larger amplitude localizations are all in an area with a total extent of one to one and a half dominant S-wave wavelengths in depth and northern direction and an extend of about eight dominant S-wave wavelengths in eastern direction. When the fourth localization that is farthest from the others in Ie is excluded, the extent in eastern direction reduces to about four and a half dominant S-wave wavelengths. One localization in Ie and one localization in Is are only about half the dominant S-wave wavelength apart (Fig. 14b). Localizations found in Ip are farther away from Ie and Is and are considered to be artificial. The large extent of the localizations suggests that there may be lateral inhomogeneities in the seismic velocity that are not resolved with the 1-D velocity model. If this particular event would have been in the part of the model muted by the sensitivity map, no localizations would be expected to stand out as larger in size and amplitude than the bulk of localizations. The source location of this microearthquake is determined to be in the centre between the two closest localizations found in Ie and Is. The source was located at |$\mathrm{ east }= -1.41 \, \textrm{km}$|⁠, |$\mathrm{ north} = 1.85 \, \textrm{km}$| and |$\mathrm{ depth} = 1.43 \, \textrm{km}$|⁠. This location is in good agreement with preliminary results derived from ray-based methods (E. Gaucher, 2019, personal communication). Further studies will work on enhancing this result and understanding the source process. 7 DISCUSSION The sensitivity maps created here for the geothermal field of Los Humeros enable the estimation of the source–location accuracy in the investigated model domain. As expected, sources can generally be located beneath the network in shallower depths and in larger depths in a smaller area. However, the pattern of source locations with high- and medium-accuracy localizations are rather complex. This highlights the importance of such an accuracy study to properly understand the influence of the station geometry and the velocity model on the accuracy. With this kind of accuracy study, the localization of (micro)earthquakes may be significantly enhanced. The location accuracy with existing networks may be estimated. Furthermore, the developed workflow can be tested for different station geometries to optimize existing or planned networks. The automatic source localization used in this study was developed for synthetic point sources but can be adapted to true events. If the true events can no longer be represented as point sources, the TRI images itself may be used to characterize the source. No a priori information about the sources is used during the localization. Therefore, the method excludes user bias and may be used to locate any number of randomly distributed point sources. However, sources too close together or close to the model boundaries may hinder individual localizations. Typical spatial errors obtained in this study were smaller than one S-wave wavelength (⁠|$460 \, \textrm{m}$|⁠) and typical temporal errors were smaller than one period (⁠|$0.2 \, \textrm{s}$|⁠) of the source signal. The sources categorized as high-accuracy sources have even smaller spatial and temporal errors. This accentuates the future importance of TRI as a high-accuracy localization method. For generating sensitivity maps, it proved useful to locate sources only in one depth per simulation to accurately locate the individual sources. However, the test case with four sources demonstrated that in general sources in different depths can be located, albeit with a weaker peak amplitude. Therefore, it is expected that true events can be located in multiple depths as long as there are not too many sources at the same time. Furthermore, in this study, all sources were generated with the source signal independent of the source depth. In reality, source strength may vary and therefore the accuracy may be increased. In most of the simulations presented in this study, an explosion source was used. However, an anisotropic source type, such as strike-slip source, was investigated as well and showed in general a similar pattern in the source location accuracy. The details varied substantially between sensitivity maps for strike-slip sources and explosion sources. This is the result of the different radiation patterns of the two significantly different source types. However, although the source types are very different, the general structure of the sensitivity maps was quite similar. If the aim of the sensitivity maps is a general idea about the area in which sources may be located, it is sufficient using an isotropic source type. For a detailed interpretation of the accuracy along a certain fault or in the vicinity of a well, multiple source types should be investigated. However, this could require a significant amount of testing to incorporate most of the source types occurring in such a reservoir setting. Sensitivity maps were created with a synthetic gradient velocity model and with a velocity model based on Löer et al. (2020) for Los Humeros. The gradient velocity model revealed a maximum depth for localizations of about |$8 \, \textrm{km}$| while the velocity model based on Löer et al. (2020) only enabled to locate sources to a depth of about |$6 \, \textrm{km}$|⁠. Strike-slip sources could only be located accurately to a depth of |$5 \, \textrm{km}$| with the this velocity model. Since the velocity models are only 1-D the influence of the velocity model on the source–location accuracy is only in vertical direction. It is expected that a similar influence can be seen for 3-D variations in the velocity model. The influence of the velocity model may be due to the change in wavelength. This can be seen in very shallow depths where the velocity is higher for the gradient velocity model than for the velocity model based on Löer et al. (2020) and more sources could be located with the latter. When both velocity models are very similar, the number of located sources is very similar. However, it has to be noted that in this study it is assumed that the correct velocity model is known. The same velocity model was used for the forward and the time-reversed simulation. An erroneous velocity model may further increase uncertainties and reduce the overall accuracy, as can be seen in the location of a microearthquake. Nevertheless, using an erroneous velocity model that hinders localizations in a part of the model may result in sources not being located that otherwise could be located. Therefore, it is important to be aware of the influence the velocity structure itself has on the accuracy. The sensitivity maps show that the influence of the station geometry is larger than the influence of the velocity model. Although the velocity models are only 1-D, the accuracy shows complex 3-D variations. The conclusion is that the interstation distance and the general positioning of the stations influence the accuracy throughout the model. This effect was observed by Franczyk (2019) and validated here using numerical experiments. In comparison to Werner & Saenger (2018), it is shown that different parameters of the station geometry seem to interact non-linearly to produce the complex sensitivity maps, which is also mentioned by Franczyk (2019). All simulations in this study were performed noise-free and therefore the sensitivity maps need to be understood as a best-case scenario of the source–location accuracy. As seen in other studies (Gajewski & Tessmer 2005; Witten & Artman 2011; Werner & Saenger 2018) an increased noise level decreases the location accuracy but does not hinder the localization. However, when sensitivity maps are used for the design or evaluation of station networks, sensitivity maps should be created including the estimated noise for the study region to obtain a sensitivity map that shows in which regions source localizations is very likely to be accurate. Noise-free sensitivity maps show in which regions no source localizations can be expected. Most studies using time-reverse imaging focus on enhancing the localization capabilities of the method by, for example, introducing new imaging conditions or by optimizing the localization scheme in general (Anderson et al. 2009; Witten & Artman 2011; Horstmann et al. 2015; Kim & Lees 2015; Franczyk et al. 2017). However, this study emphasizes that an optimization of the station network may significantly enhance the source–location accuracy. 8 CONCLUSIONS In this study, a workflow was developed to evaluate the spatial variation of the source–location accuracy that can be achieved with TRI. Numerous synthetic sources were located simultaneously while using the station distribution and the velocity model of the geothermal field of Los Humeros. Using numerous sources simultaneously enables the evaluation of the accuracy in a larger part of the model while keeping computational costs low. The sensitivity maps for the geothermal field of Los Humeros highlight the complexity of the spatial variation in the source–location accuracy and the necessity to assess the accuracy for each project individually. The usability of the derived sensitivity maps was demonstrated by locating a microearthquake. To locate numerous sources simultaneously, the identification of source locations in the TRI images was performed automatically and without a priori information about the sources. It was found that the distribution of the synthetic sources needs to follow some guidelines to ensure that each of the sources can be evaluated individually: (1) Sources can be located if they are at least as far apart as the largest interstation distance. (2) Sources should not be closer than four S-wave wavelengths to the model boundaries. (3) Sources should be located in each depth separately and then combined for a 3-D sensitivity map. Applying the developed work flow to the station network and velocity model for the geothermal field in Los Humeros reveals two major influences on the location accuracy: The station network seems to control the horizontal extent of the area in which sources can be located. The 1-D velocity model seems to dictate the maximum depth in which sources can be located. In the geothermal field of Los Humeros sources can be located underneath the network and to a depth of about 6 km with a high accuracy. However, there appear to be areas in the model where a localization is challenging. If a microearthquake occurs in one of these zones, it may not be detectable. In general, the influence of the station geometry on the source–location accuracy is larger than the influence of the velocity model. This means that the most significant increase in accuracy for TRI may be achieved by optimizing the station network. ACKNOWLEDGEMENTS The authors would like to acknowledge the Comisión Federal de Electricidad (CFE) for kindly providing support and advice and for granting access to their geothermal fields for the deployment of the seismic stations. We also acknowledge our Mexican and European colleagues for their help and collaboration in deployment, maintenance and retrieval of the stations. This project received funding from the European Union’s Horizon 2020 - Research and Innovation Framework Programme under grant agreement no. 727550 (http://www.gemex-h2020.eu, last access: 2019 October 16). The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu, last access: 2019 October 16) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). Author Contributions: CF was the lead author and wrote the paper, performed the numerical experiments and analysed the results. EHS provided the finite-difference code, gave advice and support on the scope of the project, provided background knowledge and helped with editing the paper. REFERENCES Aki K. , Richards P., 1980 . Quantitative Seismology: Theory and Methods , WH Freeman & Co . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Anderson B.E. , Guyer R.A., Ulrich T.J., Bas P.-Y.L., Larmat C., Griffa M., Johnson P.A., 2009 . Energy current imaging method for time reversal in elastic media , Appl. Phys. Lett. , 95 ( 021907 ), 1 – 3 ., doi.org/10.1063/1.3180811 . 10.1063/1.3180811 Google Scholar OpenURL Placeholder Text WorldCat Crossref Brocher T.M. , 2005 . Empirical relations between elastic wavespeeds and density in the Earth’s crust , Bull. seism. Soc. Am. , 95 ( 6 ), 2081 – 2092 ., doi.org/10.1785/0120050077 . 10.1785/0120050077 Google Scholar OpenURL Placeholder Text WorldCat Crossref Calcagno P. , Evanno G., Trumpy E., Gutiérrez-negrín L.C., Macías J.L., Carrasco-núñez G., Liotta D., 2018 . Preliminary 3-D geological models of Los Humeros and Acoculco geothermal fields ( Mexico ) – H2020 GEMex Project , Adv. Geosci. , 45 , 321 – 333 ., doi.org/10.5194/adgeo-45-321-2018 . 10.5194/adgeo-45-321-2018 Google Scholar OpenURL Placeholder Text WorldCat Crossref Cesca S. , Dahm T., Juretzek C., Kühn D., 2011 . Rupture process of the 2001 May 7 Mw 4.3 Ekofisk induced earthquake , Geophys. J. Int. , 187 ( 1 ), 407 – 413 ., doi.org/10.1111/j.1365-246X.2011.05151.x . 10.1111/j.1365-246X.2011.05151.x Google Scholar OpenURL Placeholder Text WorldCat Crossref Clayton R.W. , Engquist B., 1980 . Absorbing boundary conditions for wave-equation migration , Geophysics , 45 ( 5 ), 895 – 904 ., doi.org/10.1190/1.1441094 . 10.1190/1.1441094 Google Scholar OpenURL Placeholder Text WorldCat Crossref Duncan P.M. , Eisner L., 2010 . Reservoir characterization using surface microseismic monitoring , Geophysics . 75 , 75A139 – 75A146 . Google Scholar OpenURL Placeholder Text WorldCat Fink M. , Cassereau D., Derode A., Prada C., Roux P., Tanter M., Thomas J.-L., Wu F., 1999 . Time-reversed acoustics , Rep. Prog. Phys. , 281 ( 5 ), 91 – 97 . Google Scholar OpenURL Placeholder Text WorldCat Franczyk A. , 2019 . Using the Morris sensitivity analysis method to assess the importance of input variables on time - reversal imaging of seismic sources , Acta Geophys. 67 , 1525 – 1533 ., doi.org/10.1007/s11600-019-00356-5 . 10.1007/s11600-019-00356-5 Google Scholar OpenURL Placeholder Text WorldCat Crossref Franczyk A. , Leśniak A., Gwiżdż D., 2017 . Time reversal seismic source imaging using peak average power ratio (PAPR) parameter , Acta Geophys. , 65 ( 2 ), 299 – 308 ., doi.org/10.1007/s11600-017-0022-0 . 10.1007/s11600-017-0022-0 Google Scholar OpenURL Placeholder Text WorldCat Crossref Gajewski D. , Tessmer E., 2005 . Reverse modelling for seismic event characterization , Geophys. J. Int. , 163 ( 1 ), 276 – 284 ., doi.org/10.1111/j.1365-246X.2005.02732.x . 10.1111/j.1365-246X.2005.02732.x Google Scholar OpenURL Placeholder Text WorldCat Crossref Ge Q. , Han L., 2018 . Time reverse locating of noise source based on isotime point imaging , J. Appl. Geophys. , 159 , 38 – 51 ., doi.org/10.1016/j.jappgeo.2018.07.015 . 10.1016/j.jappgeo.2018.07.015 Google Scholar OpenURL Placeholder Text WorldCat Crossref Horstmann T. , Harrington R.M., Cochran E.S., 2015 . Using a modified time-reverse imaging technique to locate low-frequency earthquakes on the San Andreas Fault near Cholame, California , Geophys. J. Int. , 203 ( 2 ), 1207 – 1226 ., doi.org/10.1093/gji/ggv337 . 10.1093/gji/ggv337 Google Scholar OpenURL Placeholder Text WorldCat Crossref Jolie E. et al. ., 2018 . GEMex—a Mexican-European research cooperation on development of superhot and engineered geothermal systems , in Proc. 43rd Work. Geotherm. Reserv. Eng. , Stanford, California, pp. 1 – 10 . Google Scholar OpenURL Placeholder Text WorldCat Kim K. , Lees J.M., 2015 . Imaging volcanic infrasound sources using time reversal mirror algorithm , Geophys. J. Int. , 202 ( 3 ), 1663 – 1676 ., doi.org/10.1093/gji/ggv237 . 10.1093/gji/ggv237 Google Scholar OpenURL Placeholder Text WorldCat Crossref Kocur G.K. , Saenger E.H., Grosse C.U., Vogel T., 2016 . Time reverse modeling of acoustic emissions in a reinforced concrete beam , Ultrasonics , 65 , 96 – 104 . Google Scholar OpenURL Placeholder Text WorldCat Köhn D. , De Nil D., al Hagrey S.A., Rabbel W., 2016 . A combination of waveform inversion and reverse-time modelling for microseismic event characterization in complex salt structures , Environ. Earth Sci. , 75 ( 18 ), doi:10.1007/s12665-016-6032-4. Google Scholar OpenURL Placeholder Text WorldCat Lermo J. , Yanet A., Quintanar L., Lorenzo C., 2007 . Estudio sismológico del campo geotérmico de Los Humeros, Puebla, México. Parte I: Sismicidad, mecanismos de fuente y distribución de esfuerzos , Geotermia , 21 ( 1 ), 25 – 41 . Google Scholar OpenURL Placeholder Text WorldCat Lermo J. , Lorenzo C., Antayhua Y., Ramos E., Jiménez N., 2016 . Sísmica pasiva en el campo geotérmico de los Humeros, Puebla-México y su relación con los pozos inyectores , in XVIII Congreso Peruano de Geología, Sociedad Geológica del Perú, Lima, Perú, p. 5 . OpenURL Placeholder Text WorldCat Li M. , Li H., Tao G., Ali M., Guo Y., 2019 . Microseismic event location using multi-scale time reversed imaging , J. Pet. Sci. Eng. , 174 , 144 – 160 . Google Scholar OpenURL Placeholder Text WorldCat Löer K. , Toledo T., Norini G., Zhang X., Curtis A., Saenger E.H., 2020 . Imaging the deep structures of the Los Humeros geothermal field, Mexico, using three-component ambient noise beamforming , Prep. Google Scholar OpenURL Placeholder Text WorldCat Maxwell S.C. , Rutledge J., Jones R., Fehler M., 2010 . Petroleum reservoir characterization using downhole microseismic monitoring , Geophysics , 75 ( 5 ), doi:10.1190/1.3477966. Google Scholar OpenURL Placeholder Text WorldCat Oren C. , Shragge J., 2019 . 3D anisotropic elastic time-reverse imaging of surface-recorded microseismic data , in SEG Int. Expo. 89th Annu. Meet. , pp. 3076 – 3080 . Google Scholar OpenURL Placeholder Text WorldCat Saenger E.H. , 2011 . Time reverse characterization of sources in heterogeneous media , NDT E Int. , 44 ( 8 ), 751 – 759 . Google Scholar OpenURL Placeholder Text WorldCat Saenger E.H. , Gold N., Shapiro S.A., 2000 . Modeling the propagation of elastic waves using a modified finite-difference grid , Wave Motion , 31 ( 1 ), 77 – 92 . Google Scholar OpenURL Placeholder Text WorldCat Shapiro S.A. , Rothert E., Rath V., Rindschwentner J., 2002 . Characterization of fluid transport properties of reservoirs using induced microseismicity , Geophysics , 67 ( 1 ), 212 – 220 . Google Scholar OpenURL Placeholder Text WorldCat Toledo T. , 2019 . Local Earthquake tomography of the Los Humeros geothermal field , Prep. Google Scholar OpenURL Placeholder Text WorldCat Toledo T. et al., 2019 . Dataset of the 6G seismic network at Los Humeros, 2017-2018 , GFZ Data Services, doi:10.14470/1T7562235078 OpenURL Placeholder Text WorldCat Urban E. , Lermo J.F., 2013 . Local seismicity in the exploitation of Los Humeros Geothermal Field, Mexico , in Proc. Thirty-Eighth Workshop on Geothermal Reservoir Engineering , Stanford, California , pp. 1 – 11 . Google Scholar OpenURL Placeholder Text WorldCat Werner C. , Saenger E.H., 2018 . Obtaining reliable source locations with time reverse imaging: limits to array design, velocity models and signal-to-noise ratios , Solid Earth , 9 ( 6 ), 1487 – 1505 . Google Scholar OpenURL Placeholder Text WorldCat Witten B. , Artman B., 2011 . Signal-to-noise estimates of time-reverse images , Geophysics , 76 ( 2 ), MA1 – MA10 . Google Scholar OpenURL Placeholder Text WorldCat Yang J. , Zhu H., 2019 . Locating and monitoring microseismicity, hydraulic fracture and earthquake rupture using elastic time-reversal imaging Jidong Yang and Hejun Zhu , Geophys. J. Int. , 216 , 726 – 744 . Google Scholar OpenURL Placeholder Text WorldCat Zhu T. , 2015 . Viscoelastic time-reversal imaging , Geophysics , 80 ( 2 ), A45 – A50 . Google Scholar OpenURL Placeholder Text WorldCat © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Sensitivity maps for time-reverse imaging: an accuracy study for the Los Humeros Geothermal Field (Mexico)
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa160
DA - 2020-07-01
UR - https://www.deepdyve.com/lp/oxford-university-press/sensitivity-maps-for-time-reverse-imaging-an-accuracy-study-for-the-osV2FOW2FR
SP - 231
EP - 246
VL - 222
IS - 1
DP - DeepDyve
ER -