# Modelling earthquake location errors at a reservoir scale: a case study in the Upper Rhine Graben

X. Kinnaert, E. Gaucher, U. Achauer, T. Kohl

► To cite this version:

X. Kinnaert, E. Gaucher, U. Achauer, T. Kohl. Modelling earthquake location errors at a reservoir scale: a case study in the Upper Rhine Graben. Geophysical Journal International, 2016, 206 (2), pp.861-879. 10.1093/gji/ggw184. hal-03514609

HAL Id: hal-03514609

https://cnrs.hal.science/hal-03514609

Submitted on 6 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modelling earthquake location errors at a reservoir scale: a case study in the Upper Rhine Graben

X. Kinnaert,<fnref n="1" />,<fnref n="2" /> E. Gaucher,<fnref n="1" /> U. Achauer<fnref n="2" /> and T. Kohl<fnref n="1" />

Accepted 2016 May 10. Received 2016 February 3; in original form 2015 August 27

## SUMMARY

Earthquake absolute location errors which can be encountered in an underground reservoir are investigated. In such an exploitation context, earthquake hypocentre errors can have an impact on the field development and economic consequences. The approach using the state-of-the-art techniques covers both the location uncertainty and the location inaccuracy—or bias—problematics. It consists, first, in creating a 3-D synthetic seismic cloud of events in the reservoir and calculating the seismic traveltimes to a monitoring network assuming certain propagation conditions. In a second phase, the earthquakes are relocated with assumptions different from the initial conditions. Finally, the initial and relocated hypocentres are compared. As a result, location errors driven by the seismic onset time picking uncertainties and inaccuracies are quantified in 3-D. Effects induced by erroneous assumptions associated with the velocity model are also modelled. In particular, 1-D velocity model uncertainties, a local 3-D perturbation of the velocity and a 3-D geostructural model are considered. The present approach is applied to the site of Rittershoffen (Alsace, France), which is one of the deep geothermal fields existing in the Upper Rhine Graben. This example allows setting realistic scenarios based on the knowledge of the site. In that case, the zone of interest, monitored by an existing seismic network, ranges between 1 and 5 km depth in a radius of 2 km around a geothermal well. Well log data provided a reference 1-D velocity model used for the synthetic earthquake relocation. The 3-D analysis highlights the role played by the seismic network coverage and the velocity model in the amplitude and orientation of the location uncertainties and inaccuracies at subsurface levels. The location errors are neither isotropic nor aleatoric in the zone of interest. This suggests that although location inaccuracies may be smaller than location uncertainties, both quantities can have a cumulative effect. Besides, small velocity uncertainties applied to the whole 1-D profile can lead to large increase of the location uncertainties. However, local variations of the velocity field around the well may have negligible effects that would make such a feature undetectable with an absolute location method. Although the reference 1-D velocity model was built from well log data, the results show that it is not a good representative of a more realistic 3-D model including a fault and its associated block shift. The amplitude and distribution of the induced location inaccuracies are such that the positioning and the orientation of features delineated by seismicity are distorted and may be difficult to correctly interpret.

**Key words:** Probability distributions; Earthquake source observations; Seismic monitoring and test-ban treaty verification; Seismicity and tectonics; Fractures and faults; Europe.

## 1 INTRODUCTION

Earthquake hypocentres constitute a unique source of information for understanding the physical processes at the origin of earthquakes, describing the subsurface and quantifying earthquake seis-

mic hazard. It is the primary attribute of an earthquake without which other characteristics such as occurrence time, seismic moment, magnitude or focal mechanism cannot be determined. The merging of these individual properties into catalogues allows investigating the space-time interactions between earthquakes or the

<footnote n="1">¹Karlsruhe Institute of Technology, Institute of Applied Geosciences, Division of Geothermal Research, Kaiserstrasse 12, Geb 50.40, D-76131 Karlsruhe, Germany. E-mail: xavier.kinnaert@kit.edu</footnote>

<footnote n="2">²EOST-IPGS UMR7516, 5 rue René Descartes, F-67000 Strasbourg, France</footnote>
dynamic and spatial scale of the seismic ruptures. Tectonic interpretation or fault identification using the spatial distribution of earthquakes is also the purpose of many studies. However, earthquake location errors exist and need to be properly quantified for reliable result interpretation.

The earthquake location error can be described as the combination of two quantities: the location inaccuracy and the location imprecision. The latter will be confused with the *a posteriori* location uncertainties that are directly linked to the uncertainties which are taken into account during the location process (e.g. Tarantola 2005). Typically, *a priori* picking uncertainties of seismic body waves are integrated in the inverse problem and lead to *a posteriori* location uncertainties which are therefore part of the location result. On the contrary, the location inaccuracy is defined as the wrong positioning of the earthquake hypocentre due to all effects that have been ignored in the inverse problem, either for practical reasons or because of a lack of knowledge. These simplifications of the reality introduce systematic errors or bias in the computation of the earthquake location. In most cases, the use of a velocity model not representative of the effective seismic propagation medium would lead to earthquake location inaccuracies (Pavlis 1986; Bardainne & Gaucher 2010; Husen et al. 2013).

The earthquake hypocentres together with their associated errors determine the scale at which the subsurface can be investigated. The use of earthquake hypocentres as direct markers of geological structures illustrates this. For example, if uncertainties of several tens of kilometres are suitable to describe subduction zones at 200 km depth (e.g. Pegler & Das 1998), they are inadequate for local scale studies covering areas smaller than $10 \times 10 \text{ km}^2$. Interpretation of the interaction between ongoing activity and geological structures in mines may require earthquake location errors of the order of a few metres (e.g. Kwiatek et al. 2010; Plenkers et al. 2010). To assess the efficiency of successive hydraulic fracturing stages for enhanced oil and gas recovery, microearthquake location errors smaller than 25 m are often necessary (e.g. Bardainne & Gaucher 2010). Megies & Wassermann (2014) show that, in a geothermal field, depth uncertainties larger than 500 m may be problematic for unequivocal interpretation of the structures and physical processes at the origin of the earthquakes induced at ca. 5 km depth. These examples emphasize that, generally, more than one order of magnitude between earthquake location errors and description of earthquake-based spatial features is necessary for valuable use of the information.

In this work, we focus on uncertainties and inaccuracies of absolute location of earthquakes in underground reservoirs. Therefore, relative location algorithms of earthquakes based on double-differences (e.g. Waldhauser & Ellsworth 2000) or master-slave events (e.g. Fréchet 1985; Poupinet et al. 1985) are not considered, although they often constitute a further step to better constrain hypocentres of clustered earthquakes. Location techniques based on waveform stacking and migration from dense array recordings (Kao & Shan 2004; Gharti et al. 2010; Drew et al. 2013; Grigoli et al. 2013) are not examined either, despite used routinely in hydraulic fracturing monitoring for oil and gas fields. The modelling of the hypocentre errors associated with these location methods would differ from the one we apply here and remains beyond the scope of this work.

It is common practice to deploy seismic networks over underground reservoirs in applications such as geothermal energy production, oil and gas production, underground storage or mining. In these cases, the monitored volume is typically of the order of $10 \times 10 \times 10 \text{ km}^3$. Within these industrial contexts, earthquake location errors may have a major impact on the field development and

economic consequences. In particular, experience in the geothermal domain shows that induced seismicity always contributed to the identification of faults in the reservoir (Sausse et al. 2010; Kraft & Deichmann 2014; Edwards et al. 2015; Frietsch et al. 2015). In many cases, this led to optimize the positioning of several wells, and sometimes to evaluate the field economic performances (e.g. Held et al. 2014).

The understanding of the physical processes at the origin of the induced seismicity in geothermal contexts greatly increased with the observation and analysis of this phenomena (Cornet et al. 2007; Lengliné et al. 2014; Zang et al. 2014; Gaucher et al. 2015b), as well as the development of forecasting and mitigation approaches (GEISER 2013; Gaucher et al. 2015a). In the domain of enhanced oil and gas recovery similar use of recorded seismicity is often made. However, to our knowledge, only few studies provide, model or discuss earthquake location errors (e.g. Pavlis 1986; Lomax et al. 2009; Husen & Hardebeck 2010), especially at such a spatial scale and with respect to inaccuracies (or bias). Earthquake locations are usually provided with uncertainties but without inaccuracies, which, by definition, must be evaluated separately from the standard location algorithm. Yet, it remains of importance to quantify them because they may be much larger than the location uncertainties. In other words, the location uncertainty domain may not contain the true earthquake location and is therefore not sufficient to correctly interpret the earthquake locations. Since inaccuracy implies that the real earthquake location is known, only synthetic models may be considered, unless a few controlled seismic sources have known locations (Bardainne & Gaucher 2010; Husen et al. 2013). Consequently, dedicated procedures, which are relatively time consuming and which require *a priori* knowledge to deliver beneficial results, must be developed. Our approach is described in the next section. This lack of *a priori* modelling of the earthquake location uncertainties and inaccuracies in reservoir contexts motivated this work and its application to the Rittershoffen (France) geothermal field. This site is located in the Upper Rhine Graben (URG) and surrounded by several other enhanced geothermal systems (EGS) like Soultz-sous-Forêts (France), Landau, Insheim and Bruchsal (all in Germany). Hence, although the quantitative results will be related to the Rittershoffen site, the existing or future geothermal fields in the URG could benefit from the general conclusions.

The paper is organized as follows: in Section 2, we first recall the general problem of absolute earthquake location in 3-D and then present the multi-step approach based on the state-of-the-art techniques that we applied to model the location errors. The case study is described in Section 3 which gives the basic information necessary to apply the methodology (geothermal field context, seismic monitoring configuration and reference seismic velocity model). In Section 4, realistic scenarios for seismic wave picking and velocity model errors are applied to model the associated hypocentre location errors. Each case is presented and discussed before giving concluding remarks and perspectives to the work.

## 2 METHODOLOGY

### 2.1 Earthquake absolute location

The earthquake absolute location is a classical nonlinear problem in geophysics, whose objective is to minimize discrepancies between seismological observations—presently only seismic wave onset times—and the solution of the forward problem, in a given physical system, for a given space–time source. The minimization
process is referred to as the inverse problem, whereas the modelling of the observations is referred to as the forward problem (Tarantola 2005). Accordingly, the earthquake hypocentre and origin time depend, among others, on the seismic onset times and associated uncertainties (i.e. the observations), on the velocity model (i.e. the physical system), and on the location algorithm applied (i.e. the minimization process).

In practice, the earthquake location numerical code NonLinLoc (NLL) developed by Lomax et al. (2000) was used in the process of modelling the location errors. This code, which is widely used in the seismological community, has several advantages. First, it keeps the non-linearity of the location inverse problem. Second, it is able to locate earthquakes in 3-D velocity models, like we intend to do; and third, it is appropriate to a reservoir scale. Thanks to a Bayesian formalism, the location probability density of an earthquake hypocentre can be estimated everywhere in the location zone expected to include the event. The absolute location method proposed by Wittlinger et al. (1993) and implemented in NLL is applied. It combines the traveltime computation algorithm of Podvin & Lecomte (1991) to solve the forward problem, with the least-square formalism of Tarantola & Valette (1982) to define the misfit function and to compute the location probability density function (PDF). Accordingly, the PDF σ at location X is defined by

$$ \sigma(X) = K \cdot \exp \left\{ -\frac{1}{2} \left[ (\tilde{T}_{\text{Obs}} - \tilde{T}_{\text{Calc}}(X))^T \cdot C^{-1} \cdot (\tilde{T}_{\text{Obs}} - \tilde{T}_{\text{Calc}}(X)) \right] \right\} \quad (1) $$

with K a normalizing constant, $\tilde{T}_{\text{Obs}}$ the vector of the observed arrival times minus their weighted mean, $\tilde{T}_{\text{Calc}}$ the vector of the computed arrival times minus their weighted mean, and $C^{-1}$ the weight matrix. The t superscript indicates the transpose operation.

The $\tilde{T}_{\text{Calc}}$ vector is obtained by the Podvin & Lecomte (1991) numerical code which solves the Eikonal equation using the Huygen's principle and a finite-difference approximation. This numerical code computes the traveltimes between a seismic sensor and any source located on a Cartesian grid, for the first arrivals of the body-waves. It works in any type of seismic velocity model, in particular in strongly heterogeneous 3-D models, without losing its robustness. Therefore, it is well adapted to the reservoir scale and to velocity models we wish to consider. We do not include in the weight matrix $C^{-1}$ any factor due to incorrect theoretical computation of the traveltimes by the Podvin & Lecomte (1991) code. This is justified by the fact that only this traveltime computation code will be applied throughout the study (see Section 2.2). However, the $C^{-1}$ diagonal matrix contains the inverse of the uncertainties (standard deviations) of the observed arrivals times $T_{\text{Obs}}$ which are independent from each other (see Section 4.1 for details).

The earthquake hypocentre is located where the PDF is the highest, in other words where the misfit function (right term in brackets in eq. 1) is the smallest. From the determination of the PDF in the location zone, the true a posteriori earthquake location uncertainty is available in the 3-D space. Assuming a Gaussian distribution of the PDF around the best location, the location uncertainty can be described by the so-called confidence ellipsoid which includes all grid points of confidence level larger than 68.3 per cent. The analysis of the confidence ellipsoid using principal components gives the expected location (gravity centre of the grid-point distribution) and the three orientations and lengths of the uncertainty orthogonal axes (Lomax et al. 2000). In the following, the location uncertainty will be quantified by the half-length of the largest confidence ellipsoid

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic of the multistep approach applied to model earthquake location errors.</figcaption></figure>

axis and, therefore, the location should be understood within ± the uncertainty length.

Because $\tilde{T}_{\text{Calc}}$ is calculated once for a given velocity model, with this formalism, the earthquake location uncertainties do not contain any factor associated with the errors in the velocity model. This explains why uncertainties and inaccuracies have to be discriminated although they both contribute to the earthquake location error. We will see in Section 4.2.1 that the formalism may be adapted under specific assumptions to integrate velocity model uncertainties into location uncertainties.

To guarantee an earthquake hypocentre associated with the best absolute fit and a correct spatial estimate of the a posteriori location uncertainties, eq. (1) is calculated for all grid-points in the location zone (i.e. grid-search option). Once the hypocentre has been obtained, the earthquake origin time can be computed since it is equal to the weighted mean of the observed arrival times minus the computed ones (Moser et al. 1992).

For several reasons, commonly used linearized iterative methods such as Hypo71 (Lee & Lahr 1975), Hypoinverse (Klein 1978) or Hypoellipse (Lahr 1999) have not been retained to perform the analysis. Although they allow fast earthquake location, the linearization of the location inverse problem around an estimated initial solution was considered as a possibly too strong assumption with regards to the seismic network coverage. Moreover, we did not want to be limited to 1-D velocity-depth models for the different configurations we intend to investigate. However, few characteristics of these linearized iterative methods will be reproduced to examine their effects on the location errors. In particular, the use of observed seismic arrival times rounded at 10 ms and the use of 1-D velocity-depth models will be shown for Rittershoffen.

Interested readers can find more details about the absolute location of earthquakes in seismological books (e.g. Lee & Stewart 1981, chap. 6, Lay & Wallace 1995, chap. 6), or in Lomax et al. (2009), Husen & Hardebeck (2010) and Pavlis (1986) for problematics closer to ours.

## 2.2 Synthetic modelling

In order to quantify earthquake location errors (inaccuracies and uncertainties), a multi-step approach is applied (Fig. 1). In short, it consists in (1) defining synthetic earthquake hypocentres and computing, in a given velocity model, the associated seismic arrival times at the seismic stations of a network (the synthetic modelling step), (2) relocating the earthquake hypocentres using the simulated times but after perturbing the initial conditions (the relocation step), and (3) comparing the relocated hypocentres with the initial ones.
In the synthetic modelling step, the sources are positioned in the subsurface, where they are expected to occur within the reservoir. Then, the arrival times of the *P* and *S* waves corresponding to the synthetic sources and observed on the seismic network are computed in a given velocity model. This model is supposed to represent the reality of the seismic wave propagation in the Earth. Therefore, the associated traveltimes will be taken as the observed arrival times at the seismic network and used to relocate the synthetic earthquakes. As presented in detail in Section 4.2, several models will be simulated. Each model may represent a reality of the subsurface associated with an *a priori* information. For the same reasons as described previously, the numerical code of Podvin & Lecomte (1991) is applied at this step to compute the seismic traveltimes.

The next step of the methodology consists in relocating the earthquakes from the synthetic arrival times, which were previously generated, by applying the NLL standard location code described in the previous subsection. Therefore, it represents a real processing flow, especially as the initial hypotheses will be changed to mimic the intrinsic lack of knowledge of the reality. One unique velocity model will be used throughout the study to relocate the synthetic earthquakes. This model is a 1-D depth-velocity model consisting of a series of horizontal layers of constant *P*- and *S*-wave velocities. Such a flat-layer model geometry is applied because it is used in a vast majority of real cases even at a local scale (e.g. Dorbath *et al.* 2009; Bönnemann *et al.* 2010; Maurer *et al.* 2015). This is sometimes justified by the lack of information available to better characterize the subsurface velocities, but it is also a convenient way to quickly obtain earthquake locations because these model geometries are easy to handle in linearized location numerical codes, especially in real-time processing software. This is precisely the effects of such simplifications that we want to quantify, especially when *a priori* information gives evidences that it may not apply.

Once the synthetic earthquake relocation is obtained, it is compared to the initial location to quantify the location error. In the error, we distinguish the location inaccuracy from the location uncertainty. The former is defined as the spatial distance separating the initial source hypocentre with the relocated one. The latter corresponds to the spatial domain delimited by the event relocation probability at the 68.3 per cent confidence level. This level is chosen to allow comparison with the widely used standard deviation output by many location numerical codes.

# 3 CASE STUDY: RITTERSHOFFEN GEOTHERMAL FIELD

## 3.1 Geological settings

The about 300 km long URG is part of the European Cenozoic rift system, which extends over more than 1000 km in Central Europe (Ziegler 1992; Schumacher 2002). The deep geothermal field of Rittershoffen (Alsace, France) is located at the western margin of the NE-SW-striking central segment of the URG. The geological settings in Rittershoffen are rather comparable with those described for the geothermal field of Soultz-sous-Forêts (Dezayes *et al.* 2011), which is located about 6 km to the NW (Fig. 2). Polyphase reactivations of a complex pattern of Late Variscan and Permo-Carboniferous crustal discontinuities (Schumacher 2002; Ziegler *et al.* 2006) and intense Cenozoic faulting during evolution of the URG resulted in complex subsurface structures. Hence, several graben and horst structures, which dissect the approximately 2000-m thick Mesozoic and Cenozoic succession, are delimited

by synthetic and antithetic faults, predominantly N-S to NNE-SSW striking (Fig. 3). The granitic basement is highly fractured and faulted and, at its uppermost part below the boundary to the overlying sediments, it shows distinct indications of alteration and probably represents a former erosional discordance (Genter 1989). Subsurface fluid flow along permeable fault zones under an extensional or transtensional stress state probably explains the occurrence of large-scale convection cells and positive heat anomalies in the whole area (Pribnow & Schellschmidt 2000; Baillieux *et al.* 2013). Such anomalies contribute to the deep exploitation of this renewable energy resource in several active fields such as Landau, Insheim, Bruchsal (all in Germany) or Soultz-sous-Forêts (France) and furthermore highlight the geothermal potential in this part of the URG (Meixner *et al.* 2016).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** General view of the Upper Rhine Graben and of the deep geothermal fields (red circles) close to Rittershoffen (black circle).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Vertical cross-section of the geological structure between Soultz-sous-Forêts and Rittershoffen geothermal fields (adapted from GeORG Project Team 2015).</figcaption></figure>
## 3.2 Field development and induced seismicity

Together with the acquisition of new geophysical data, the long hydrocarbon production history and the proximity of the Soultz-sous-Forêts EGS field contributed to the detailed characterization of the Rittershoffen subsurface and the field development plan. The site benefits from one of the largest thermal gradient observed in the URG and temperatures of approximately 165 °C were measured at 1800 m depth (Baujard et al. 2014). Unlike the Soultz-sous-Forêts geothermal field, but like the Landau and the Insheim fields, the Triassic sandstone and the underlying Palaeozoic granite constitute the exploited reservoir formations at Rittershoffen. To increase the chance to access permeable zones, the normal fault delimiting the two blocks below Rittershoffen was targeted by the geothermal doublet (Fig. 3). This normal fault, which is approximately oriented N–S and dipping 60°W, accounts for approximately 200 m vertical shift. Both wells of the doublet are drilled down to ca. 2400 m mean-sea-level, that is 200 m below the Buntsandstein–granite interface.

To enhance the connectivity between GRT1 and the reservoir, several stimulations were carried out in 2013 (Baujard et al. 2014). All stimulations induced seismicity which was recorded by a surface seismic network. Several hundreds of events were detected and located, all with local magnitude smaller than $M_L = 1.6$ (Maurer et al. 2015). Most of the seismicity is located in the southwest of GRT1 and centred close to the bottom of the hole at ca. 2500 m depth. The seismic cloud is roughly oriented N–S to NNE–SSW and is approximately 2 km long, 1 km wide and 2 km high (Maurer et al. 2015).

Accordingly, to model the earthquake location errors, we will generally distribute the synthetic earthquake hypocentres in a cube of approximately 4 km side centred on the GRT1 open-hole mid-depth. The sources will be regularly spaced either in the 3-D volume, or on several planes included in the volume (see grey boxes and segments in Fig. 4). The detailed synthetic source positions will be presented in each of the examined scenarios.

## 3.3 Seismic monitoring network

The paper focuses on the period covering the chemical and mechanical stimulations of the well GRT1 which was carried out in June 2013 (Baujard et al. 2014). The seismic network taken as reference for the location error analysis is made of 12 permanent and 5 temporary stations actively monitoring at that time (Maurer et al. 2015). Among the 12 permanent stations, 7 have three-component short-period seismometers (L-4C-3-D) and the remaining have a vertical short-period seismometer only (L-4C). The signal is sampled either at 100 Hz or at 150 Hz. All temporary stations have three-component short-period seismometers (L-4C-3-D) whose signal is sampled at 300 Hz. As shown in Fig. 4, this network only covers the northern part of GRT1, and the associated effect on the earthquake location errors will be determined. Permitting issues prevented from deploying seismic stations in the southern part of GRT1 before the stimulation. Among the 17 active stations, the furthest two (GUNS and LAMP) are discarded from the analysis. According to N. Cuenot (private communication, 2015), the signal to noise ratio of the recorded induced seismicity at these stations was too low to enable body-wave picking and to include them in the location process. So, it was decided to keep the same assumption for the synthetic test. We further assumed that the synthetic P-wave arrival times could be used for all remaining 15 stations but that the S-wave arrivals could be used only for the 12 three-component stations. Hence, the P- and S-wave arrival times are always used in the (re)location process.

This configuration did not change during the stimulation of GRT1 and is used over the whole study. However, more seismic stations were deployed afterwards and the effect of adding one surface station (E3316) in the southern part of the network (Fig. 4) is discussed in the Section 4.3.

## 3.4 Seismic velocity model

At Rittershoffen, a 1-D velocity-depth profile was created from several measurements taken in the well GRT1. This model, which was used to process the induced seismicity (Maurer et al. 2015), is considered as the 1-D reference model in the study and will always be used for relocating the synthetic sources. The P-wave interval velocities were computed from a zero offset VSP for the principal stratigraphic layers, which were identified from the geological log of the well (Aichholzer et al. 2015). Then, using a compression-and shear-velocity log, the Vp/Vs ratio was averaged within each layer from a three-component sonic log performed in GRT-1 and applied to obtain the S-wave velocity profile. Fig. 5 shows the P- and S-wave velocity profiles at Rittershoffen.

Within this 13-layers model, the P-wave velocities range between $1320 \text{ m s}^{-1}$ at surface and $5815 \text{ m s}^{-1}$ in the granitic formation, and the S-wave velocities between $620 \text{ m s}^{-1}$ at surface and $3275 \text{ m s}^{-1}$ in the granitic formation. The Vp/Vs ratio varies between layers, from 1.68 to 2.12, the highest values being observed for the tertiary formations. Three embedded low-velocity layers exist for the depth intervals 1025–1300 m, 1630–2000 m and 2100–2200 m. Two large velocity contrasts are also observed at the top of the Lias layer (1365 m) and at the top of the granite (2200 m).

In the NLL location numerical code, this 1-D velocity model is discretized on a 10-m mesh size in the east, north and depth directions.

# 4 RESULTS AND DISCUSSIONS

In this section, we present several scenarios applied to investigate the location uncertainties and inaccuracies, which could be expected at Rittershoffen, based on the reference seismic network. Several effects are examined: first the seismic wave picking uncertainties and inaccuracies, and second the velocity model uncertainties and inaccuracies, for a few realistic scenarios. For all scenarios, the reference 1-D velocity model is used to relocate the synthetic earthquakes.

For clarity purpose, we used the GRT1 wellhead (1 010 653.18 m east, 2 447 831.75 m north, Lambert II extended) as the geographical origin for the latitude and longitude presentation of the results. Unless specified, the depth is a true vertical depth from mean sea level given in meter.

## 4.1 Location errors driven by the seismic onset time

### 4.1.1 Picking uncertainties

Earthquake hypocentre errors may be due to a variety of effects. First, the effect of the arrival time uncertainties on the amplitude and geometry of the hypocentre uncertainties, in the reference velocity model, using the reference seismic network, has to be quantified. Accordingly, the 1-D-velocity model is used both for the synthetic modelling step and the relocation step. The only difference between both steps is the introduction of uncertainties in the synthetic arrival times. Consequently, no location inaccuracy is expected, which means that the hypocentre is positioned at the initial source location.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Map of the seismic network deployed at Rittershoffen. The Soultz-sous-Forêts permanent stations (blue squares), the Rittershoffen permanent stations (red squares) and the temporary stations (green triangles) are shown as well as GRT1 wellhead (red crossed circle). Station E3316 (reversed green triangle) was unavailable during the period of interest and is used only in Section 4.3. The largest rectangle delimits the velocity model zone used in this study. The smallest rectangle delimits the area in which the seismic sources are simulated; the two segments are the projections of the vertical sections of simulated sources. The projection of the zone with perturbed velocity (see Section 4.2.1) is shown as a grey circle. All coordinates are in Lambert II extended system.</figcaption></figure>

The *P*- and *S*-wave arrival times were given uncertainty values from manual picking of several seismograms observed by the network during the stimulation of GRT1 (N. Cuenot, private communication, 2015). For the eight stations the closest to GRT1 (BETS, RITT, KUHL, E3301, E3308, E3313, E3315 and E3317, see Fig. 4), *P*- and *S*-picking uncertainties were set to ±20 ms, and for the seven remaining stations (FORA, KEFF, OBER, OPS, SCHW, STUN and SURB, see Fig. 4), they were set to ±50 ms. So, no picking uncertainty difference was noted between the *P*- and the *S*-wave arrivals. These picking uncertainties constitute the matrix *C* of eq. (1) and will apply in the rest of the study. The synthetic sources were placed every 200 m in the 3-D cube centred at the bottom of the GRT1 well.

Fig. 6 shows a 3-D view of the relocation uncertainties at the 68.3 per cent and 99.7 per cent confidence levels, for a subset of 27 sources around GRT1 well separated by 1200 m. These levels are often used and correspond to one and three standard deviations respectively. As observed, most of the uncertainties at these confidence levels look ellipsoidal and their shape can be simply described

by three orthogonal directions and three lengths (e.g. Lomax *et al.* 2009). This also means that the uncertainties can be assumed Gaussian distributed around the highest probability location. To check whether the confidence ellipsoid quantifies correctly the uncertainties, we compared the hypocentre relocation with the expected location (i.e. the gravity centre of the distribution of all grid points included in the 68.3 per cent confidence level, see Section 2.1) for all synthetic sources, and we assumed that similar values were enough to validate the hypothesis. For 86 per cent of the cases, the median and the 75 per cent percentile of the location discrepancies are below the 10-m mesh size thereby negligible. However, median discrepancies up to 30 m (and 40 m for the third quartile) exist for sources located close to the large velocity contrasts at 1400 m and 2200 m depth, as could be expected (e.g. Lomax *et al.* 2009). Increase of the distortion of the uncertainties from a Gaussian distribution is also observed when moving southwards from GRT1, because the seismic network does not cover this zone. Hence, although we will assume in the following that the confidence ellipsoid of the location uncertainty is reliable in most of the volume below Rittershoffen, it
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Profiles of the reference *P*-wave (blue curve) and *S*-wave (black curve) velocities at Rittershoffen overlaying the 150 random profiles for the *P*-wave (red curves) and the *S*-wave (green curve) velocities (see Section 4.2.1 for details). The bottom depth of the well GRT1 (horizontal dashed line) and the upper depth of the injection zone (horizontal dotted-dashed line) are also shown.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** 3-D view of the relocation uncertainties at 68.3 per cent (dark blue) and 99.7 per cent (light blue) confidence levels of 27 sources distributed around GRT1 well (red curve). Part of the seismic stations is also displayed (grey cones).</figcaption></figure>

is important to keep in mind that this is not correct close to the high velocity contrasts at 1400 and 2200 m depth (Fig. 5) and southwards from GRT-1.

Fig. 7 shows a horizontal section, at the bottom of the well GRT1 (2414 m), of the length of the largest axis of the 68.3 per cent confidence ellipsoid. In the rest of the paper, this quantity will be associated with the location uncertainty. As observed, the uncertainty increases towards south and east of the GRT1 well. Such a horizontal variability is kept over depth (Fig. 8 bottom).

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Horizontal section of the location uncertainty at 2414 m. The stations above the location zone are displayed (white triangles) as well as the GRT1 well trajectory (black curve).</figcaption></figure>

Boxplots of the uncertainty variation as a function of latitude, longitude and depth (Fig. 8) emphasize these spatial variations for the set of sources in the 3-D volume. From a median value of 140 m north of GRT1, the uncertainty increases up to 190 m at 2 km south of GRT1. Uncertainty also slightly increases towards east: from 140 to 150 m. However, as shown by the results as a function of depth, the uncertainty variations depend mainly on the earthquake depth. Starting from about 50 m 2 km above GRT1 bottom hole, the location uncertainties increase to 150 m at the top of the granite (at approximately 2500 m, close to GRT1 bottom hole) and then remain almost constant. At the depth of the stimulated zone (2089 m), the uncertainty represents about 6 per cent of the distance to the surface. The orientation of the uncertainty is also anisotropic (but not shown). It roughly points towards the KUHL station, which is located about 4 km NNW from GRT1 well-head (see Fig. 4), with an inclination between 40° and 45°. Consequently, for a given confidence level, the rough tendency for the hypocentre location is to either become deeper away from GRT1 or shallower closer to GRT1.

All seismic stations used for the analysis are located north of GRT1 or at similar longitude, with a higher station density in the western part of the investigated zone. Such a seismic coverage explains very well the spatial distribution of the uncertainties both in amplitude and direction. The uncertainty variation with depth is also consistent with the lack of down-hole station and with the strong velocity contrasts: the more south, and/or the more east and/or the deeper the earthquake sources are from the GRT1 well-head, the larger are the uncertainties.

To conclude, this analysis shows that even in the correct velocity model, with realistic picking uncertainties and with optimistic use of all 15 stations of the network, earthquake location uncertainties of minimum 150 m are to be expected in the granite (from 2400 m) and of approximately 50 m at 400 m depth.

### 4.1.2 Picking precision

As mentioned previously, many linearized iterative algorithms work with initial seismic arrival times of 10 ms precision. We wish to investigate here whether or not what can be viewed as a random
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Boxplots of the location uncertainties as a function of latitude (top left), longitude (top right) and depth (bottom left). The uncertainties median (points), 25th and 75th percentiles (crosses), and range (dotted segment) are shown. The 10-m relocation-mesh size is displayed (grey dashed-line) as reference.</figcaption></figure>

rounding of the arrival times at the stations plays a role in the hypocentre inaccuracy, at the spatial scale of the Rittershoffen field (i.e. few kilometres on surface and a depth around 2.5 km).

To examine this effect, synthetic sources were placed every 200 m in the 3-D cube centred at the mid-depth of the GRT1 open-hole section. Then, between the synthetic modelling step and the relocation step, nothing was changed but the *P* and *S* traveltimes which were rounded to the closest 10 ms.

The relocation of the synthetic sources using NLL leads to inaccuracies, which may exceed the 10-m relocation-mesh size. The larger inaccuracies are observed to the south and to the east of GRT1, as well as in depth and, therefore, have a distribution comparable to the one of the relocation uncertainties.

Fig. 9 displays the horizontal, depth and total location inaccuracies as a function of the latitude. As observed, no relevant depth variation is observed since the median, 25th and 75th percentiles of the inaccuracies remain all within the relocation-mesh size (10 m). The median of the horizontal inaccuracy slightly decreases from south to north but is still smaller than 10 m. In fact, however, values larger than the mesh size are observed only when total inaccuracy is considered. Then, from approximately 12 m 2 km north of GRT1 the median inaccuracies increase to 15 m 2 km south of GRT1. The

75th percentile also increases from 15 m to 20 m. These largest values remain, however, relatively small; but, considering the range, distances between synthetic and relocated hypocentres may reach 45 m. Fig. 9 also displays the variation of the location inaccuracies as a function of depth, independently from the latitude and longitude. Inaccuracies increase with depth and the median becomes larger than the mesh size below the first strong velocity contrast, at about 1400 m depth. The median inaccuracy can reach 18 m, 2 km below the GRT1 injection zone and the 75th percentile, 20 m.

When locating using only the *P*-wave arrival time, all other things remaining similar, the inaccuracies were multiplied by a factor of 2 to 3, the spatial distribution being unchanged.

To conclude, for the case of Rittershoffen, rounding the *P*- and *S*-wave arrival times at 10 ms can lead to inaccuracies of the order of 10–20 m, with a maximum of 45 m. From 1400 m depth, the largest bias is approximately 40 m. The detailed study highlighted that despite relatively small, these inaccuracies are not isotropic, neither in depth nor in latitude. Their spatial distribution, which depends on the network coverage and on the velocity model, is similar to the distribution of the location uncertainties, which remain similar to those presented in the previous subsection. As a consequence, the resulting location inaccuracies should be added to the location
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Boxplots of the depth (top-left), the horizontal (top-right) and the total (bottom left) location inaccuracies of the synthetic earthquakes as a function of the latitude, and location inaccuracies as a function of depth (bottom right). On this last subplot, the 1-D velocity profile is shown in grey. The relocations were computed by NLL using the *P*- and *S*-wave arrivals rounded at the closest 10 ms. For details on the symbols, see Fig. 8.</figcaption></figure>

uncertainties rather than randomly included inside the location uncer-
tainties (although the latter are larger). There is a cumulative
effect in the location error due to the 10-ms picking accuracy.

**4.2 Location errors driven by the velocity model**

Velocity model errors lead to earthquake location inaccuracies,
which may be so large that the true earthquake location may even not
be included in the relocation uncertainties (Bardainne & Gaucher
2010; Gesret et al. 2015). Such an effect may be dramatic for inter-
pretation when location uncertainty is mistaken for location error.
In this subsection, we model the effects induced by three differ-
ent types of velocity uncertainties and inaccuracies which could
be encountered at Rittershoffen. First, we introduce uncertainties
in the original 1-D velocity profile; second, the velocity model is
perturbed in 3-D around the stimulated zone; and third, a fault and
its associated velocity shift are introduced leading to another 3-D
model. One after each other, these models will represent the real
Earth in which the synthetic traveltimes will be calculated. How-

ever, in each case, the relocation will be computed in the initial 1-D
velocity model still considered as the reference (Fig. 5).

4.2.1 1-D velocity model uncertainty

The 1-D velocity model proposed for Rittershoffen (see Section
3.4) is certainly not perfect beyond the hypothesis of the intrin-
sic lateral invariability. It has been created by merging information
from the GRT1 Z-VSP for the *P*-wave velocity, and the GRT1 sonic
log, for the Vp/Vs ratio, which are not error-free measurements and
which do not characterize the ground properties at similar spatial
scale. In particular, the sonic log provides velocities which can dif-
fer from the Z-VSP and which deviate notably from their moving
average over depth. To model the uncertainty in the 1-D veloc-
ity profile, 150 Vp and Vs profiles were generated by randomly
drawing values of a Gaussian distribution with 5 per cent standard
deviation, centred on the original velocities (Fig. 5). This 5 per cent
level is chosen because almost all velocity values measured by the
sonic log fall within the interval defined by plus or minus three
times this standard deviation. The random variations were taken
<figure><img src="image_12.png" /><figcaption>**Figure 10.** 3-D view of the relocation uncertainties at 68.3 per cent confidence levels for the 81 sources distributed every 300 m in the E–W vertical plane (light grey) and in the N–S vertical plane (dark grey) around GRT1 well (red curve). Part of the seismic stations is also displayed (grey cones) and the initial 1-D velocity model is shown as an N–S vertical section.</figcaption></figure>

independent per seismic phase and per layer (whose depth remained unchanged), but only those preserving the increase or the decrease of velocities between two layers were retained. The depths of the layers were kept unchanged because they are well constrained from the chronostatigraphic log of the GRT1 borehole (Aichholzer et al. 2015).

Assuming that the numerous velocity profiles correctly sample the density distribution of the real Rittershoffen one, they can be used in a probabilistic framework to propagate the velocity model uncertainty to the earthquake location uncertainty (Gesret et al. 2015). According to the Bayesian formalism developed by Gesret et al. (2015), the probability density distribution of the earthquake location is equal to the sum of the probability density distributions computed within each velocity profile. In such a formalism, the propagation of the picking uncertainties (as per eq. 1) is also kept; hence, the final hypocentre location and its uncertainty depend on both the velocity and the picking uncertainties. This approach, in which the non-linearity of the location problem is kept, provides a final relocation uncertainty domain always including the real location.

In practice, the synthetic modelling step is carried out in each of the 150 perturbed velocity profiles (for the *P* and *S* waves) and the probability density function of the earthquake location computed within the initial 1-D velocity profile using the grid-search algorithm of NLL. The maximum of the PDF resulting from their summation gives the hypocentre relocation and the 68.3 per cent confidence level, their uncertainty. Two vertical planes of synthetic sources have been created, one in the N–S direction and one in the E–W direction, both crossing the bottom of GRT1 well. The N–S direction corresponds roughly to the direction of the initially located earthquakes in the area (Maurer et al. 2015) and to the direction of the main horizontal stress (Cornet et al. 2007). Each plane contains 9 × 9 sources spaced by 300 m. During the relocation step, the reference seismic network and the *P*- and *S*-wave arrival times with their reference uncertainties are used.

Fig. 10 gives a 3-D view of the relocation uncertainties at 68.3 per cent confidence level for each of the 81 sources of the E–W and N–S vertical planes. The associated volumes correspond roughly to ellipsoids in the granitic formation, deeper than 2200 m; however, this is less valid for the sedimentary cover in which sev-

eral velocity contrasts exist and are potentially emphasized during the randomization process. Such an observation is consistent with the one made considering no velocity model variation (see Section 4.1) but the effect is strengthened. Nevertheless, to simplify the description of the location uncertainties due to velocity model uncertainties, we will suppose that the confidence ellipsoids are good enough representatives. The introduction of 5 per cent uncertainty in the velocity model leads to larger variations of the earthquake hypocentre uncertainties compared to the fixed-model case. Location uncertainties up to 650 m are observed, which is more than 2.5 times larger (Fig. 11) than the ones in Section 4.1.1. As previously observed, the deeper, the more east and the more south the events are of the well GRT1, the larger are the uncertainties.

The perturbations of the initial 1-D-velocity model also introduce inaccuracies of the hypocentre relocations. As a matter of fact, despite the initial model is randomly perturbed according to a Gaussian distribution, the non-linearity of the location problem with regards to the velocity models leads to hypocentre relocations different from the original locations. In the sediments, these inaccuracies (median value for the 2 × 81 sources) are of the order of 10 m to 20 m but can reach 60 m in the granitic formation, and a maximum inaccuracy of 155 m is observed. In all of these cases, however, the inaccuracies are lower than the associated uncertainties. No systematic difference between the initial and the relocated hypocentres is observed along depth. However, it exists along the longitude and latitude and is, once again, spatially distributed like the uncertainties along these directions.

To conclude, Gaussian perturbation of the Rittershoffen initial velocity profile by 5 per cent can lead to more than 250 per cent increase of the location uncertainty and displacement of the original hypocentre. In the granitic formation, uncertainties in the range of 200 m to 650 m may be expected, which is not negligible with regards to the depth of the investigated zone: 2200 m to 3300 m.

### 4.2.2 Local 3-D perturbation of the initial velocity model

In this second case, variation of the initial 1-D velocity model is applied to the area around the injection zone of the GRT1 well. This zone may be seen as a volume in which the seismic velocity properties may change due to fluid injection. It is known that seismic velocities may differ between dry and fluid saturated rock samples, but also by varying the effective stress of rock samples (Spencer & Nur 1976; Lockner et al. 1977; Stanchits et al. 2003). In the Soultz-sous-Forêts enhanced geothermal system, Calò et al. (2011) highlighted *P*-wave seismic velocity decreases of up to 10 per cent during massive water injection, at 5 km depth, into the granite. This example illustrates our motivation to perform the present test and to quantify the effects of such localized variations on the seismic event absolute locations.

Therefore, over the 650 m depth interval of the GRT1 fluid injection and within a radius of 325 m around the well, the *P*- and *S*-wave velocities were decreased by a factor of 10 per cent, whatever the geological formation was (granite, Buntsandstein and Muschelkalk). The velocity perturbed zone corresponds to a low-velocity cylinder embedded in the original 1-D model and the velocity model representative of the real Earth becomes 3-D. Synthetic sources are located on two vertical planes, which are crossing GRT1 at the average injection depth, one striking N–S and one striking E–W. The sources are distributed on a regular 40-m mesh including the perturbed zone and are relocated on a 20-m mesh.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Vertical section at 2 m north (left) and 215 m west (right) of the largest location uncertainty interpolated on a 10-m mesh. The GRT1 well trajectory is displayed (black curve).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Projection of the initial locations and of the relocations of the synthetic earthquake hypocentres for the vertical E–W plane (left) and for the vertical N–S plane (right). The relocated hypocentres (blue points) should be shifted along the red line to come back to their initial location. The shadowed area represents the perturbed velocity zone.</figcaption></figure>

Fig. 12 shows the displacements of the earthquake locations in-
duced by the use of the 1-D velocity model instead of the perturbed
3-D model. Both vertical planes are presented. As observed, the
largest location inaccuracies occur in the depth range correspond-
ing to the deeper half of the perturbed velocity zone, below the top of
the granite. Moving away from the perturbed zone, the inaccuracies
decrease. Events in the western part are located at the west of their
initial position whereas events in the eastern part are located more
to the east. Similar relative behaviour exists for the events on the
N–S vertical plane. Hence, earthquake relocations have a tendency
to move away from their original locations relatively to the centre
of the perturbed zone. We also observe higher inaccuracies in the
eastern and southern parts than in the western and northern parts.
All these effects can be well explained by the seismic coverage and

the seismic ray bending, as for the earthquakes inaccuracies above
the perturbed zone, which result from the propagation of refracted
seismic waves in the inaccurate velocity model.

The median of the horizontal location inaccuracy can reach al-
most 40 m for the E–W plane (Fig. 13), and is larger than 20 m (the
relocation-mesh size) in the deeper half of the perturbed velocity
zone. A maximum of 55 m horizontal shift is observed. For the
N–S plane, the median is smaller but the largest uncertainty can
reach 85 m. In most cases, the vertical inaccuracy is smaller than
the horizontal one, whatever the vertical plane is, and in 50 per cent
of the cases they are not significant taking into account the mesh
size. Locally, however, inaccuracies up to 75 m exist.

No significant differences are observed between the length
and the orientation of the relocation uncertainties in this
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Boxplots of the relocation inaccuracies for the seismic sources located on the E–W vertical plane (left column) and on the N–S vertical plane (right column). The horizontal (top row), vertical (middle row) and total (bottom row) relocation inaccuracies as a function of the event depth are presented. The vertical grey lines at ±20 m indicate the mesh size of the relocation model, the 1-D velocity profile is also displayed as a grey curve. For details on the symbols, see Fig. 8.</figcaption></figure>
perturbed velocity test compared to the unperturbed velocity model (Section 4.1).

This test highlights that a decrease of 10 per cent both in Vp and Vs, in a localized zone around the well GRT1, is difficult to detect from the absolute location of earthquakes. We stress here that we are only considering absolute locations and that further modelling should be carried out to quantify the effects associated with relative location techniques (outside the scope of this work). The resulting inaccuracies are one order of magnitude smaller than the location uncertainties and have no systematic trends which would notably change the shape of the vertical planes. Besides, the average *P*- and *S*-wave time residuals at the seismic stations after relocation in the perturbed model are all smaller than $2 \pm 4$ ms, for the events, on both planes. In other words, the time variations induced by the 10 per cent velocity decrease inclusion would be very difficult if not impossible to detect using absolute location techniques because they are smaller than the picking uncertainties. Furthermore, because of their similar amplitudes, the inaccuracies induced by this local velocity decrease or by the use of 10 ms arrival-time precisions may be mistaken (see Section 4.1.2). Finally, since the average residuals are smaller than their associated uncertainties, it does not make any difference by applying the corresponding station corrections to improve the event relocations.

Testing a Vp and Vs decrease of 20 per cent instead of 10 per cent led to location inaccuracies distributed similarly but with larger amplitudes. Despite larger, the inaccuracies remain smaller than the location uncertainties and the time residual still smaller than the picking uncertainties. The effect of fluid injection on the decrease of the *P*- and *S*-wave velocities is generally not the same and may be associated with a reduction of the Vp/Vs ratio. Accordingly, we also tested 10 per cent and 5 per cent decreases of the *P*- and *S*-wave velocities respectively. This corresponds to approximately 5 per cent decrease of the Vp/Vs ratio. In such a case, the location inaccuracies are very small and the median over depth remains almost smaller or of the order of the mesh. Consequently, the perturbed velocity zone associated with each of these configurations would also be undetected from absolute location of earthquakes.

### 4.2.3 Geo-structural effect: 3-D velocity model with a fault

Heat exploitation at Rittershoffen was motivated by the presence of an unusual high-temperature anomaly and a fault which could drive geothermal fluids (see Section 3.1, Fig. 3). For these reasons, the well GRT1 was drilled to cross the major fault below surface. This fault separates, from the granitic formation up to the Tertiary formations, two blocks of the graben which are vertically shifted by 200 m. Consequently, although assuming a flat velocity model built from data acquired in the well GRT1 is locally representative, this is certainly not the case away from the well.

To investigate the effect of this structural setting on the earthquake locations, a 3-D velocity model was built for Rittershoffen. An N–S fault dipping 60° W and crossing the well at 2200 m was included. According to the existing seismic profiles, this is an acceptable planar representation of the fault, which surface is obviously more complex. The initial 1-D-velocity profile was kept for the western block of the fault but was shifted 200 m upwards for the eastern block. So, the eastern block is on average faster than the western block.

Four planar surfaces with synthetic events were created. Two of them are vertical, one striking N–S the other one striking E–W, another one is horizontal, and the three of them cross GRT1 at

2089 m depth (injection range mid-depth). The last planar surface containing synthetic sources is the fault plane. For all planar surfaces, the synthetic sources are positioned on a 50-m squared mesh of 2400 m side length and centred on the injection point. As usual, to quantify the absolute location errors, the synthetic modelling was performed in the 3-D velocity model and the relocation computed in the reference 1-D velocity model.

On Fig. 14, the four initial planar surfaces are displayed with the four associated surfaces obtained after hypocentre relocation. One can first note that, since the location inaccuracies are varying in direction and in amplitude, the initial planar surfaces become curved surfaces. The sources initially located on the fault (Fig. 14, top) are systematically shifted towards east, slightly down and slightly southwards for the sources located north of GRT1 well and northwards for those located south of GRT1 well. This observation is consistent with the results obtained by Pavlis (1986) for a comparable simulation of a two-block velocity model in the Morgan Hill area (California). The median horizontal shift dominated by the eastern shift varies between 300 and 350 m. The median depth shift varies from 55 m for the deepest sources to 200 m for the shallowest. Despite these spatial inaccuracies, the average azimuth and dip of the relocated source surface are similar to the initial planar surface. However, a more complex transformation is observed for the three other surfaces of synthetic sources which are not parallel to the fault. In these cases, the initial planes are deformed in a continuous manner which depends on the initial source locations, but the fault plane and the associated block shift are neither delineated nor directly visible in the results. For the horizontal plane at 2089 m (Fig. 14, top), with increasing Easting coordinate of the synthetic sources the eastern shift increases, and the initial southward shift becomes a northward shift. The eastern shift is however dominating and the median horizontal location inaccuracies vary from 250 to 350 m from the western to the eastern side of the plane (Fig. 15). Still from west to east, the depth difference between the relocated and the initial sources decreases from 230 to –30 m (Fig. 15), which means that, although the sources are relocated deeper than expected on the western side of the plane, they are shallower on the eastern side. The sources of the E–W vertical plane (Fig. 14, bottom) are relocated eastwards from their original position with increasing shifts from the upper side of the original plane to the lower one. An original southward shift from the western side of the E–W plane becomes a northward shift at the eastern side of the plane. This change always occurs to the east of the modelled fault but for smaller easting coordinates as depth increases. The median of the horizontal inaccuracies ranges between 250 and 320 m. With regards to the depth, the relocated events are deeper and deeper than expected as a function of depth for the upper-western part of the plane but are shallower for the lower-eastern part (located to the east of the fault). Median depth inaccuracies range from –90 to 240 m. At last, the N–S vertical plane (Fig. 14, bottom) moves eastwards by a larger offset at depth. A southward-down shift from the upper-northern corner of the N–S plane becomes a northward-up shift at the lower-southern corner of the plane. Therefore, the median of the horizontal inaccuracies increases regularly with depth from 200 m in the sediments to 400 m in the granite and the vertical inaccuracies from –50 m in the bedrock to 230 m in the sediments (Fig. 15).

To summarize the inaccuracy results, not considering the fault and the associated block shift leads to systematically locate the earthquakes eastwards from expected when the 1-D velocity model built from GRT1 well data is used. The offsets are in average close to 350 m. Moreover, the sources at the NW-upper corner of the
<figure><img src="image_16.png" /><figcaption>**Figure 14.** 3-D view of the relocation of the synthetic events initially located on the horizontal plane and along the fault (top) and on the N-S and E-W vertical plane (bottom). The initial locations were on the grey planes whereas the relocations are on the coloured planes. The colour scale is associated with the largest uncertainty value which varies between 50 and 250 m. The location of the fault is shown by the brown plane (bottom).</figcaption></figure>

investigated zone are relocated deeper and more southwards than originally but this tendency is reversed when considering sources located at the east of the fault. However, the fault position is not associated with a sharp variation of the location inaccuracies. In general, the inaccuracies are enhanced with the depth. Since the inaccuracies are variable in space, sources initially distributed on planar surfaces are relocated on curved surfaces. This means that planar features delineated by earthquakes are deformed. In other words, the azimuth and dip directions of such features are not kept constant (except for sources distributed on the fault).

Fig. 14 also displays the location uncertainties associated with the relocated sources (colour scale of the planes). They range between 50 and 250 m and are bigger but still of the order of magnitude of the uncertainties observed in the original 1-D velocity model, both in amplitude and in spatial distribution (see Section 4.1).

To conclude, this test clearly shows the very strong location bias induced by neglecting the fault and the associated shifted blocks at Rittershoffen for events initially located within a radius of at least 1200 m around the injection point. The location inaccuracies, of the order of 200–400 m, are larger than the location uncertainties which vary between 50 and 250 m. Consequently, they dominate and cannot be considered as included in the location uncertainties. As shown, inaccuracy values and directions vary in space in a continuous manner although eastward shift is dominating. This prevents from applying a systematic and constant correction factor to the relocated hypocentres, although a westward shift of minimum 200 m could be used. This also makes difficult the delineation of the fault separating the two blocks. Moreover, expected delineation of seismic events in any direction or on any plane would be deformed in most cases and, as a consequence, spatial
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Boxplots of the horizontal (left column) and vertical (right column) inaccuracies for the horizontal plane along east (top row) and along north (middle row) and for the N–S vertical plane along depth (bottom row). For details on the symbols, see Fig. 8.</figcaption></figure>

interpretation should be revised. This anisotropic bias is controlled by the fault geometry and the block shift as well as the coverage of the network.

This test shows that the reference 1-D velocity model built from GRT1 well data is not representative of the 3-D model with the

fault and it quantifies the induced errors. It does not seek for the best 1-D velocity model which would be representative of the 3-D model with the fault; or, in other words, for the 1-D velocity model minimizing the location errors. This work would imply a travel-time tomography which is a domain we do not consider here.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Horizontal section of the location uncertainty at 2414 m using the 16-station network. The stations above the location zone are displayed (white triangles) as well as the GRT1 well trajectory (black curve). Same display as Fig. 7.</figcaption></figure>

## 4.3 Addition of one surface station

The previous discussions apply to the seismic network which was monitoring during the seismogenic phase of the GRT1 hydraulic stimulation. As shown, the network, which is mainly deployed northwest from the well, has a strong impact on the analysed earthquake hypocentre errors. Here, we shortly present and discuss the principal effects of adding one surface station.

To keep a realistic approach, the chosen supplementary station, E3316 (Fig. 4), is taken from the temporary network installed after the stimulation of the GRT1 well and before the drilling of the GRT2 well. It is located in the forest, approximately 3 km SSE from the GRT1 well, on the eastern block delimited by the fault described in Section 4.2.3. Therefore, the seismic network coverage of the zone of interest is improved. All origins of location error investigated in the previous two subsections were also applied to this extended network and the processing flow remained the same. We focus, however, on the main results and their differences with those of the 15-station network.

At station E3316, the *P*- and *S*-picking uncertainties were set to ±20 ms, like the other temporary stations. Following the approach described in Section 4.1.1, Fig. 16 shows, as presented in Fig. 7, a horizontal section of the location uncertainty, at 2414 m, driven by the picking uncertainties in the reference 1-D velocity model. The section illustrates the general variation of the spatial distribution of the uncertainties which becomes more homogeneous horizontally; the depth distribution being still controlled by the layered velocity model. Boxplot analysis of the uncertainties calculated in the volume delimited by ±2 km around the well highlights a reduction of the location uncertainties, between the 16-station and the 15-station networks, which reaches approximately 20 m in the granite. These results show the interest of adding one surface station to diminish the location uncertainties and to get more uniform uncertainty distributions, at least horizontally. The other tested configurations, which are considering different velocity models between the synthetic modelling phase and the relocation phase, also agree with this conclusion.

The addition of the station in the test investigating the effect of the 1-D velocity model uncertainty (see Section 4.2.1) reduces the cor-

responding location inaccuracies by maximum 50 m, in the granitic formation. However, an opposite effect is observed when considering the 3-D velocity model with the fault (see Section 4.2.3). In that case, an increase of the location inaccuracies is observed. Fig. 17 presents 3-D views of the relocation of the synthetic events initially located on the four tested source planes, with similar point of views and colour scale than in Fig. 14. The median horizontal shift, which was ranging between 200 and 350 m with the 15-station network, is now ranging between 250 and 425 m. Besides, the vertical median inaccuracies vary between -30 m and 275 m. So, the inaccuracies have increased and the planes are still distorted, thus preventing from direct interpretation of seismicity alignments in terms of location, azimuth and dip of major structures. This apparent contradictory behaviour between the impact of the 1-D velocity model uncertainty and that of the model with the fault can be explained. In the latter case, the addition of the station E3316 located in the eastern block delimited by the fault is strongly in contradiction with the assumed 1-D velocity model used to relocate the synthetic earthquakes. As a matter of fact the seismic ray paths between one event and the station E3316 will be very different in the 1-D model and in the 3-D model with the fault. Therefore, this station alone brings a lot of inconsistency in the (least-square) inverse problem since the associated seismic arrival times are outliers with regards to most of the other observed times. A more detailed study, outside the scope of this work, should be conducted to see if this effect may be counter balanced by adding more and more stations to the southern and eastern parts of the monitored area.

## 5 CONCLUSIONS AND OUTLOOK

In this paper, the impact of several factors on the absolute location of earthquakes in a reservoir was investigated. The methodology, which used state-of-the-art techniques, consisted in relocating synthetic hypocentres under hypotheses different from the data modelling step, and comparing the initial and relocated hypocentres. Hence, the effects of the *P*- and *S*-wave onset time uncertainties and inaccuracies were examined as well as the effects of the velocity model uncertainties and inaccuracies. In particular, we looked at the location errors driven by using a 1-D velocity model instead of several 1-D or 3-D velocity models. The analysis was applied to the Rittershoffen geothermal field, where seismicity was induced between 1000 and 5000 m, under the seismic monitoring conditions existing during the chemical and mechanical stimulations of the well GRT1. The 3-D analysis of the location uncertainties and inaccuracies covers a zone of a couple of kilometres around the open-hole section of the well.

At Rittershoffen, in an ideal case where, first, a 1-D velocity model is well representative of the propagation medium, and second, the *P*- and *S*-wave pickings and their uncertainties are representative of the recorded induced seismicity, the earthquakes location uncertainties will range between 50 m at 400 m depth and 150 m from 2400 m depth (in the granitic formation). The coverage of the seismic network, which is deployed in the north of the GRT1 well, controls the spatial distribution of the uncertainties in amplitude and direction. This layout leads to uncertainties roughly pointing towards the KUHL station, which is located about 4 km NNW from GRT1 well-head. Although a picking precision of 10 ms leads to inaccuracies smaller than the uncertainties, their similar spatial distribution have a cumulative effect and the location inaccuracy cannot be totally considered has contained in the uncertainty. Such an effect would be typical of the use of a location algorithm requiring 10 ms picking precision as input. The location errors driven by the
<figure><img src="image_19.png" /><figcaption>**Figure 17.** 3-D view of the relocation, by the 16-station network, of the synthetic events initially located on the horizontal plane and along the fault (top) and on the N–S and E–W vertical plane (bottom). The initial locations were on the grey planes whereas the relocations are on the coloured planes. The colour scale is associated with the largest uncertainty value which varies between 50 and 250 m. The location of the fault is shown by the brown plane (bottom). Same display as Fig. 14.</figcaption></figure>

hypotheses associated with the velocity model may be very large. As shown, only 5 per cent uncertainty in the reference 1-D velocity profile of Rittershoffen can multiply by a factor of 2.5 the uncertainties of the ideal case, thus leading to uncertainties up to 650 m in the granite at 3200 m. On the contrary, decreasing by 10 per cent the *P*- and *S*-velocities in the neighbourhood of the GRT1 open-hole section has negligible effect and makes such a feature undetectable with an absolute location method. Not considering the Rittershoffen fault and the associated block shift in the velocity model induces very strong location biases, larger than the uncertainties. Although the reference 1-D velocity model is based on well data and centred on the zone of interest, it is not a good representative of the 3-D model with the fault and does not minimize the location errors. While the location uncertainties range between 50 and 250 m, the expected sources are shifted by 200 to 400 m, mainly eastward. The

location inaccuracies, however, vary continuously in space, thus making difficult any reliable interpretation of directions or surfaces delineated by the located seismicity.

As emphasized, the location errors driven by the velocity model may be dominating at the considered reservoir scale. This confirms the common sense which recommends using any *a priori* information to better constrain the initial velocity model or to define its uncertainties, which can be later integrated into the location uncertainty. The 3-D analysis of the results highlighted that the location errors are neither aleatoric nor isotropic but clearly driven by the seismic network coverage and the velocity model. This suggests that although location inaccuracies may be smaller than location uncertainties, both quantities may have a cumulative effect.

The test consisting of adding one surface station to the network in an area so far uncovered by the seismic network showed that the
location uncertainties are reduced, as may be expected. However, it also highlighted that the hypocentre inaccuracies can, on the contrary, increase due to the addition of inconsistency between the real Earth and its model. This underlines once more the distinction which must be made between imprecision and inaccuracy.

In this study, we did not perform an exhaustive analysis of all parameters which could influence the earthquake location at Rittershoffen. We focused on those which are usually not taken into account during standard location processing because they require synthetic modelling; in particular the discrepancies between the velocity model and the real seismic propagation medium. Moreover, we selected scenarios which looked the most relevant with regards to the *a priori* information available for Rittershoffen. To allow clear identification of each effect, the associated hypotheses for relocating the earthquakes were changed one by one. However, among numerous other options, it could be reasonable, for example, to jointly consider uncertainties of the 1-D velocity profile and the inclusion of the fault. When studying the location errors driven by erroneous assumptions on the velocity models, we did not seek for a better velocity model which would minimize the location errors. This work, which would imply travel-time tomography, could constitute a future analysis.

For Rittershoffen, the quantitative results of this study can constitute an *a priori* knowledge useful for interpretation or processing of seismological data. Nearby geothermal fields such as Landau, Insheim, Bruchsal (all in Germany) or Soultz-sous-Forêts (France) have geological settings very similar to Rittershoffen geothermal reservoir. To some extent the present results can be used in these contexts, in particular at the Landau and the Insheim fields where seismicity was induced and for which the Triassic sandstone and the Palaeozoic granite also constitute the reservoir. At Soultz-sous-Forêts, the geothermal reservoirs were developed deeper into the granite (below 3500 m); however, major faults also delimit lifted blocks, which lead to velocity contrasts in the 1400 m-thickness sedimentary cover. Several geothermal fields in the URG are currently under development or explored and can benefit from this study. Since the described methodology is independent from the induced seismicity recorded at the site, it can also help in quantifying the location capabilities of a given network at a given site, even prior to the network deployment.

The earthquake hypocentre constitutes the primary attribute of a seismic event. Error of this attribute on the determination of secondary ones should be investigated in future works. The impact on the focal mechanisms, which are used to better describe the reservoir structure, is of special interest. Finally, this work focused on quantifying errors of earthquake absolute locations. The extensive use of relative location methods to obtain earthquake hypocentres indicates that adapting and applying the developed methodology to this processing frame is necessary.

## ACKNOWLEDGEMENTS

This work was conducted in the framework of the excellence laboratory 'Labex G-EAU-THERMIE PROFONDE' (University of Strasbourg). It was funded by the French National Research Agency, as part of the French 'Investments for the future' program, by the Energy Baden-Württemberg AG (EnBW), and by the French-German University (DFH-UFA). We wish to thank the ECOGI joint venture and the Électricité de Strasbourg – Géothermie company (ESG) for sharing data. We are grateful to V. Maurer, A. Genter and C. Baujard (ESG) for numerous discussions about this work, to N. Cuenot (EEIG 'Heat Mining') for sharing his experience on the

Rittershoffen raw seismic data and to J. Schmittbuhl (EOST, Strasbourg University/CNRS) for a fruitful review of an early version of the manuscript. The authors would like to thank F. Grigoli and T. Plenefisch for their valuable comments and suggestions for the improvement of the manuscript.

## REFERENCES

Aichholzer, C., Duringer, P., Orciani, S. & Genter, A., 2015. New stratigraphic interpretation of the twenty-eight-year old GPK-1 geothermal well of Soultz-sous-Forêts (Upper Rhine Graben, France), in *4th European Geothermal Workshop*, Strasbourg, France.

Baillieux, P., Schill, E., Edel, J.-B. & Mauri, G., 2013. Localization of temperature anomalies in the Upper Rhine Graben: insights from geophysics and neotectonic activity, *Int. Geol. Rev.*, **55**(14), 1744–1762.

Bardainne, T. & Gaucher, E., 2010. Constrained tomography of realistic velocity models in microseismic monitoring using calibration shots, *Geophys. Prospect.*, **58**(5), 739–753.

Baujard, C., Genter, A., Maurer, V., Dalmais, E., Graff, J.-J. & Schmittbuhl, J., 2014. The ECOGI EGS project in Rittershoffen, France, in *GRC 38th Annual Meeting*, Portland, OR, USA.

Bönnemann, C. et al., 2010. Das seismische Ereignis bei Landau vom 15. August 2009, in *Abschlussbericht der Expertengruppe, Seismisches Risiko bei hydrothermaler Geothermie*, pp. 55, Hannover.

Calò, M., Dorbath, C., Cornet, F.H. & Cuenot, N., 2011. Large-scale aseismic motion identified through 4-D P-wave tomography, *Geophys. J. Int.*, **186**(3), 1295–1314.

Cornet, F.H., Bérard, T. & Bourouis, S., 2007. How close to failure is a granite rock mass at a 5 km depth?, *Int. J. Rock Mech. Min. Sci.*, **44**(1), 47–66.

Dezayes, C., Beccaletto, L., Oliviero, G., Baillieux, P., Capar, L. & Schill, E., 2011. 3-D vizualization of a fractured geothermal field: the example of the EGS Soultz site (Northern Upper Rhine Graben), in *36th Workshop on Geothermal Reservoir Engineering*, Stanford University, Stanford, California.

Dorbath, L., Cuenot, N., Genter, A. & Frogneux, M., 2009. Seismic response of the fractured and faulted granite of Soultz-sous-Forêts (France) to 5 km deep massive water injections, *Geophys. J. Int.*, **177**(2), 653–675.

Drew, J., White, R.S., Tilmann, F. & Tarasewicz, J., 2013. Coalescence microseismic mapping, *Geophys. J. Int.*, **195**(3), 1773–1785.

Edwards, B., Kraft, T., Cauzzi, C., Kästli, P. & Wiemer, S., 2015. Seismic monitoring and analysis of deep geothermal projects in St Gallen and Basel, Switzerland, *Geophys. J. Int.*, **201**(2), 1022–1039.

Fréchet, J., 1985. Sismogénèse et doublets sismiques, *Thèse d'État*, Université J. Fourier, Grenoble.

Frietsch, M., Groos, J. & Ritter, J.R., 2015. Detection and Delineation of a Fracture Zone with Observation of Seismic Shear Wave Anisotropy in the Upper Rhine Graben, SW Germany, *Pure Appl. Geophys.*, **172**(2), 267–282.

Gaucher, E., Schoenball, M., Heidbach, O., Zang, A., Fokker, P.A., van Wees, J.-D. & Kohl, T., 2015a. Induced seismicity in geothermal reservoirs: a review of forecasting approaches, *Renew. Sustainable Energy Rev.*, **52**, 1473–1490.

Gaucher, E., Schoenball, M., Heidbach, O., Zang, A., Fokker, P.A., van Wees, J.-D. & Kohl, T., 2015b. Induced seismicity in geothermal reservoirs: physical processes and key parameters, in *World Geothermal Congress*, Melbourne, Australia.

GEISER, 2013. 'Geothermal Engineering Integrating Mitigation of Induced Seismicity in Reservoirs, Final Reports'. Available at: http://www.geiser-fp7.fr/ReferenceDocuments/Pages/ReferenceDocuments.aspx, last accessed 26 August 2015.

Genter, A., 1989. Géothermie Roches Chaudes Sèches: le granite de Soultz-sous-Forêts (Bas-Rhin, France), Université d'Orléans/BRGM.

GEORG Project Team, 2015. 'INTERREG IV Upper Rhine - Geopotentials of the deep Upper Rhine Graben'. Available at: http://maps.geopotenziale.eu/?app=georg&lang=en, last accessed 26 August 2015.
Gesret, A., Desassis, N., Noble, M., Romary, T. & Maisons, C., 2015. Propagation of the velocity model uncertainties to the seismic event location, *Geophys. J. Int.*, **200**(1), 52–66.

Gharti, H.N., Oye, V., Roth, M. & Kuhn, D., 2010. Automated microearthquake location using envelope stacking and robust global optimization, *Geophysics*, **75**(4), MA27–MA46.

Grigoli, F., Cesca, S., Vassallo, M. & Dahm, T., 2013. Automated seismic event location by travel-time stacking: an application to mining induced seismicity, *Seismol. Res. Lett.*, **84**(4), 666–677.

Held, S., Genter, A., Kohl, T., Kölbel, T., Sausse, J. & Schoenball, M., 2014. Economic evaluation of geothermal reservoir performance through modeling the complexity of the operating EGS in Soultz-sous-Forêts, *Geothermics*, **51**, 270–280.

Husen, S. & Hardebeck, J.L., 2010. *Earthquake Location Accuracy*, Community Online Resource for Statistical Seismicity Analysis, doi:10.5078/corssa-55815573.

Husen, S., Kissling, E. & Deschwanden, A., 2013. Induced seismicity during the construction of the Gotthard Base Tunnel, Switzerland: hypocenter locations and source dimensions, *J. Seismol.*, **17**(1), 63–81.

Kao, H. & Shan, S.-J., 2004. The Source-Scanning Algorithm: mapping the distribution of seismic sources in time and space, *Geophys. J. Int.*, **157**(2), 589–594.

Klein, F.W., 1978. Hypocenter location program HYPOINVERSE. Part 1. User's guide to versions 1,2,3 and 4. Open-File Rep. 78-694, U.S. Geol. Surv., pp. 113.

Kraft, T. & Deichmann, N., 2014. High-precision relocation and focal mechanism of the injection-induced seismicity at the Basel EGS, *Geothermics*, **52**, 59–73.

Kwiatek, G., Plenkers, K., Nakatani, M., Yabe, Y. & Dresen, G. & JAGUARS-Group, 2010. Frequency-magnitude characteristics down to magnitude –4.4 for induced seismicity recorded at Mponeng Gold Mine, South Africa, *Bull. seism. Soc. Am.*, **100**(3), 1165–1173.

Lahr, J.C., 1999. HYPOELLIPSE: a computer program for determining local earthquake hypocentral parameters, magnitude, and first-motion pattern, Open-File Rep. 99-23, Version 1.1, U.S. Geol. Surv., pp. 119.

Lay, C. & Wallace, T.C., 1995. *Modern Global Seismology*, Academic Press.

Lee, W.H.K. & Lahr, J.C., 1975. HYPO71 (revised): a computer program for determining hypocenter, magnitude, and first motion patterns of local earthquakes, Open-File Rep. 75-311, U.S. Geol. Surv., pp. 116.

Lee, W.H.K. & Stewart, S.W., 1981. *Principles and Applications of Microearthquake Networks, Advances in Geophysics*, Supplement 2, Academic Press.

Lengliné, O., Lamourette, L., Vivin, L., Cuenot, N. & Schmittbuhl, J., 2014. Fluid-induced earthquakes with variable stress drop, *J. geophys. Res.*, **119**, 8900–8913.

Lockner, D.A., Walsh, J.B. & Byerlee, J.D., 1977. Changes in seismic velocity and attenuation during deformation of granite, *J. geophys. Res.*, **82**(33), 5374–5378.

Lomax, A., Michelini, A. & Curtis, A., 2009. Earthquake location, direct, global-search methods, in *Encyclopedia of Complexity and Systems Science*, pp. 2449–2473, ed. Meyers, R.A., Springer.

Lomax, A., Virieux, J., Volant, P. & Berge-Thierry, C., 2000. 5 Probabilistic earthquake location in 3D and layered models, in *Advances in Seismic Event Location*, pp. 101–134, eds Thurber, C. & Rabinowitz, N., Springer.

Maurer, V., Cuenot, N., Gaucher, E., Grunberg, M., Vergne, J., Wodling, H., Lehujeur, M. & Schmittbuhl, J., 2015. Seismic monitoring of the Ritter-shoffen EGS project (Alsace, France), in *World Geothermal Congress*, Melbourne, Australia.

Megies, T. & Wassermann, J., 2014. Microseismicity observed at a non-pressure-stimulated geothermal power plant, *Geothermics*, **52**(0), 36–49.

Meixner, J., Schill, E., Grimmer, J.C., Gaucher, E., Kohl, T. & Klingler, P., 2016. Structural control of geothermal reservoirs in extensional tectonic settings: an example from the Upper Rhine Graben, *J. Struct. Geol.*, **82**, 1–15.

Moser, T.J., van Eck, T. & Nolet, G., 1992. Hypocenter determination in strongly heterogeneous Earth models using the shortest path method, *J. geophys. Res.*, **97**(B5), 6563–6572.

Pavlis, G.L., 1986. Appraising earthquake hypocenter location errors: a complete, practical approach for single-event locations, *Bull. seism. Soc. Am.*, **76**(6), 1699–1717.

Pegler, G. & Das, S., 1998. An enhanced image of the Pamir–Hindu Kush seismic zone from relocated earthquake hypocentres, *Geophys. J. Int.*, **134**(2), 573–595.

Plenkers, K., Kwiatek, G., Nakatani, M., Dresen, G. & Group, T.J., 2010. Observation of Seismic Events with Frequencies $f > 25$ kHz at Mponeng Deep Gold Mine, South Africa, *Seismol. Res. Lett.*, **81**(3), 467–479.

Podvin, P. & Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, *Geophys. J. Int.*, **105**(1), 271–284.

Poupinet, G., Fréchet, J., Ellsworth, W.L., Frémont, M.J. & Glangeaud, F., 1985. Doublet analysis: improved accuracy for earthquake prediction studies, *Earthq. Predict. Res.*, **1**, 147–159.

Pribnow, D. & Schellschmidt, R., 2000. Thermal tracking of upper crustal fluid flow in the Rhine graben, *Geophys. Res. Lett.*, **27**(13), 1957–1960.

Sausse, J., Dezayes, C., Dorbath, L., Genter, A. & Place, J., 2010. 3D model of fracture zones at Soultz-sous-Forêts based on geological data, image logs, induced microseismicity and vertical seismic profiles, *C. R. Geosci.*, **342**(7–8), 531–545.

Schumacher, M.E., 2002. Upper Rhine Graben: role of preexisting structures during rift evolution, *Tectonics*, **21**(1), 6-1-6-17.

Spencer, J.W. & Nur, A.M., 1976. The effects of pressure, temperature, and pore water on velocities in westerly granite, *J. geophys. Res.*, **81**(5), 899–904.

Stanchits, S.A., Lockner, D.A. & Ponomarev, A.V., 2003. Anisotropic changes in P-wave velocity and attenuation during deformation and fluid infiltration of granite, *Bull. seism. Soc. Am.*, **93**(4), 1803–1822.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, SIAM.

Tarantola, A. & Valette, B., 1982. Inverse problems = quest for information, *J. Geophys.*, **50**, 159–170.

Waldhauser, F. & Ellsworth, W.L., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward fault, California, *Bull. seism. Soc. Am.*, **90**(6), 1353–1368.

Wittlinger, G., Herquel, G. & Nakache, T., 1993. Earthquake location in strongly heterogeneous media, *Geophys. J. Int.*, **115**(3), 759–777.

Zang, A., Oye, V., Jousset, P., Deichmann, N., Gritto, R., McGarr, A., Majer, E. & Bruhn, D., 2014. Analysis of induced seismicity in geothermal reservoirs – an overview, *Geothermics*, **52**, 6–21.

Ziegler, P.A., 1992. European Cenozoic rift system, *Tectonophysics*, **208**(1–3), 91–111.

Ziegler, P.A., Schumacher, M.E., Dèzes, P., Van Wees, J.-D. & Cloetingh, S., 2006. Post-Variscan evolution of the lithosphere in the area of the European Cenozoic Rift System, *Geological Society, London, Memoirs*, **32**(1), 97–112.