Seismicity related to the hydraulic stimulation of GRT1, Rittershoffen, France

The Rittershoffen deep geothermal reservoir, in Northeastern France, is well characterized and has been extensively studied by a multidisciplinary approach. A hydraulic stimulation for the development of the geothermal reservoir was performed in June 2013. This injection of fluid led to seismic activity which was closely monitored by a dedicated set of seismic stations. The seismic sequence started during the injection but showed an unusual long quiet period of 4 d after shut-in before the occurrence of a second swarm of events. Here we take the opportunity of this well-monitored activity to gain insight into the geomechanical factors favouring the development of induced earthquakes. We apply a template matching approach and a relative relocation procedure to obtain a precise estimate of the geometries of the activated structures. Our approach shows that the induced events during the injection took place on two parallel planar structures. It shows that details of the seismicity generally obtained from borehole seismic network are achievable from surface network when an appropriate analysis is performed. The development of this induced seismicity is in good agreement with the known stress field and failure criterion proposed for the reservoir. In particular, the orientation of the activated structure, associated focal mechanisms and the overpressure needed to initiate the seismic activity are all in line with the geomechanical model of the area. The swarm of delayed events, 4 d after shut-in, can be well explained by considering an aseismic slip on the imaged fault and the related static stress transfer. We therefore suggest that the ability to monitor local slow aseismic movements at depth, in conjunction with precise tracking of the seismicity, is of primary importance to understand induced earthquake activity.


INTRODUCTION
Understanding how the injection of fluids in the crust is linked to the occurrence of earthquakes is of primary importance for the assessment of induced seismicity. Tracking and understanding this seismicity is important as it poses a risk to the nearby population and the infrastructures (Giardini 2009) and also provides information about the reservoir mechanics and the fluid pathway (Zang et al. 2014).
However, modeling the direct impact of fluid injection on induced events is not straightforward as it involves the coupling of various processes. Indeed the nucleation of earthquakes is affected both by the existing fractures and faults within the reservoir, the in situ stress field, the strength of the rock reservoir, the mechanical and chemical impact of the fluid flow, and the temperature variations, factors that all interact with one another (e.g. Blanpied et al. 1995). Most of these parameters are usually not known precisely nor monitored continuously at depth during the injection such that it is difficult to observe their influence on the development of the earthquake activity. The clear identification of the role of each of these parameters on the seismicity is also restricted by the limited resolution that can be achieved on earthquake detection and locations. Indeed, the structures activated by fluid injection are usually of very small scale compared to the uncertainties associated with the earthquake locations. It thus requires a dedicated seismological network and processing in order to record and locate precisely the associated seismic activity and highlight the geometry of the active structures.
Despite these difficulties some mechanisms have been tested and proved as having an influence on the development of the seismicity. The mechanical influence of the pore-pressure, P , on the triggering of earthquakes is generally evaluated from a Mohr-Coulomb failure criterion, where slip on an interface occurs when where σ n is the normal stress on the fault, τ the shear stress, µ i the internal friction coefficient of the interface and C 0 is the cohesion of the material. Note that for pre-existing weak interfaces like joints or faults, the friction coefficient and the cohesion are typically smaller than for intact materials (Scholz 2002). The experiment performed in Rangely oil field, Colorado, USA, showed that the onset of the Seismicity related to the hydraulic stimulation of GRT1 3 earthquake activity was linked to a threshold injection overpressure P that makes the left hand term in Eq. 1 positive (Raleigh et al. 1976). Cornet et al. (2007) also showed a nice example from Soultzsous-Forêts where seismicity onset and location are well predicted by the simple description of Eq. 1.
It shows that knowing the state of stress and the strength of the activated faults allows good prediction of the onset and the fate of the seismicity. However, other examples of induced seismicity exist where the earthquake activity presents a more complex dynamics that the one described by the Eq. 1 alone.
Indeed, in some cases, the diffusion of the pore pressure, elastic stress redistribution caused by the earthquakes or the aseismic movement within the reservoir have shown to influence the seismicity (Shapiro et al. 2002;Schoenball et al. 2012;Bourouis & Bernard 2007). It therefore appears that the induced earthquakes could be linked to various mechanisms.
Here we analyze the earthquake activity that develops during the hydraulic injection in the well GRT1 in Rittershoffen, France. This well is part of a deep geothermal project located at 6 km east of Soultz-sous-Forêts in Northern Alsace and developed by the ECOGI joint-venture. For this purpose, at the end of 2012, a first well (GRT1) was drilled to 2580 m depth through Triassic-sediments and into the crystalline basement. In order to enhance the reservoir permeability, a hydraulic stimulation was performed in the GRT1 well in June 2013. The hydraulic stimulation in GRT1 lasted 2 days (27 and 28 June 2013) and was recorded by a dedicated seismic network (Maurer et al. 2015;Baujard et al. 2017). The seismic activity related to the GRT1 hydraulic stimulation was processed in real-time and gave rise to a first seismicity catalog composed of a total of 212 events, from June 27 to July 4, 2013 (Maurer et al. 2015). The catalog reveals that the seismicity stopped shortly after injection but started again after 4 completely quiet days on July 2, in the form of an intense seismic swarm that lasted less than one day. The occurrence of this second swarm after a quiet period is relatively rare and poses the question of its mechanical origin. In order to understand how this second swarm developed several days after the injection was stopped, we apply a dedicated set of tools to recover and locate the earthquakes that occurred during this sequence as precisely as possible. We see that, given the current state of stress of the reservoir, the strength of the rock mass forming the reservoir and the orientation of the main activated structures, we can propose a logical geomechanical model for the development of the observed induced and triggered seismicity.

DATA
The seismic network around the Rittershoffen geothermal field at the time of the injection in June 2013 was composed of 18 surface stations in total (Maurer et al. 2015) (Figure 1). The recording sampling frequency of these stations is distributed as follows: 9 at 150 Hz, 4 at 100 Hz and 5 at 300 Hz. All stations were continuously recording the seismic signal. Most of the sites have 3 components sensors and, depending on the site, these are broadband or short period sensors. The material differences are the result of the various equipments installed in the area over time: the Soultz permanent network, the Rittershoffen permanent network and the Rittershoffen temporary network. These networks provide a dense station coverage with 12 stations located less than 5 km away from the injection point. The region at the south of the injection point, however, is not as well uniformly covered by seismic sensors because of permit issues that limited their installation.
We first process the seismic signals by performing a Short Term Average / Long Term Average (STA/LTA) approach (Earle & Shearer 1994) to detect possible seismic events over the whole time period. All signals are first filtered between 10 and 40 Hz, the STA is computed over 1s and the LTA over 5s of signal. We set a detection criterion when at least 3 stations reach the STA/LTA threshold fixed at 2. This leads to 946 detections. We then review manually all detections to discard false detections and only keep clear seismic events, resulting in a total of 686 earthquakes. We finally pick visible P-and S-wave arrival times for these events and obtain on average 8 P-wave picks per event. The location of these events, based on the picking arrival times, is then obtained by the software HYPOINVERSE (Klein 2002) using a 1D gradient velocity model of the area (Table 1). This velocity model is derived from a sonic log that was performed in GRT1 well (Maurer et al. 2015 Cornet et al., 2007 from imaging in these well are shown with thick black arrows. The name of the seismic stations are indicated above each instrument (triangles). The seismicity (orange and blue dots) defines a 1.5 km long cloud mostly located west of the injection well. We observe that events that occurred during the second swarm (blue dots) are mostly present in the northern end of this cloud.

TEMPLATE MATCHING DETECTION
The processing described above: (STA/LTA) detection combined with the location based on absolute arrival time picks, shows that almost all of the seismicity occurs less than 1 kilometer from the injection point ( Figure 1). We also observe that events of the second swarm show an offset compared to the events that occurred during the injection. However, the lack of resolution in the location of the detected earthquakes prevent any further detailed analysis at this stage. The low resolution is partly related to the difficulty of picking the arrival times of the noisy seismic signals and the supposedly small scale of the activated structures.
The seismicity we are investigating took place in a densely populated area (Cf. Figure 1) with numerous noise sources masking the earthquake signals . Not only can this cause imprecise picking arrival times but this could also lead to a reduced number of events in the final catalog. Indeed, a low signal-to-noise ratio at a sufficient number of stations, or events occurring close in time, will not result in a positive detection by the STA/LTA algorithm. It is also possible that some events that were detected by the STA/LTA procedure are actually not located confidently because of too few or inconsistent picks resulting in the removal of these events during latter processing stages. All these factors result in numerous missed events. In order to improve the seismic catalog and to recover events that may have been missed during the previous procedure, we employ a template matching approach. The template matching approach has proved to be successful in recovering previously missed events in various contexts e.g. (Peng & Zhao 2009;Helmstetter et al. 2015;Lengliné et al. 2016). It has also been successfully employed to track induced seismicity (Skoumal et al. 2014;Huang & Beroza 2015).
The template matching approach uses a known seismic signal in order to detect newer events similar to the tested one. This technique has the ability to detect events even when the signal to noise ratio at an individual station is lower than 1. It however requires the knowledge of some template events to compare the continuous seismic signal with. The choice of these events can affect the final set of newly detected events as we will detect mainly similar events to the tested ones. It implies that the newly detected events are not located too far from the template event, otherwise waveforms will not be similar and stacking the correlation signals at different stations after correcting from travel times between stations will not result in a coherent stack. One could consider all the detected events as pos-Seismicity related to the hydraulic stimulation of GRT1 7 sible template events. This approach will be quite unnecessarily time consuming as many detected events are similar such that it will akin to processing the same data several times thus obtaining the same results. Furthermore, because some identified events have a low signal to noise ratio, it is possible that at some stations the correlation detector will mainly correlate with the noise source and not the earthquake signal. In this case it will produce multiple false detections. To overcome these difficulties, we consider here a different approach. We first group all the 674 located events into clusters of similar waveforms. We compute the correlation matrix between all events around the P-wave on 128 samples long window (i.e. 1.28s) filtered between 8 and 25 Hz. The frequency range 8-25 Hz has been chosen based on the analysis of the spectrogram of the continuous data. This frequency range is best for filtering out most of the noise present in the data while keeping the earthquake signature and thus maximizing the signal to noise ratio. This range typically corresponds to the range where the corner frequency of the events is expected to fall within. All signals are first re-sampled at 100 Hz for all stations. If the mean normalized correlation coefficient at least at 3 stations is higher than 0.8, we associate the two events into the same cluster. We tested that reducing the correlation coefficient threshold for linking events in a same group do not lead to an improvement of the final number of detected events. We finally keep 13 clusters that all have at least 2 events. For each group we stack the waveforms of all events in the group in order to create a synthetic waveform representing the average waveform of the events of the group. Events are first aligned on the P-wave arrival before stacking is performed at all possible stations. We thus create synthetic template signals for the 13 groups of events. This approach is somewhat similar to the subspace detector, but here we consider only the first vector of the singular value decomposition that represent the stack of the events ( The template matching approach is performed by correlating the template signals with the continuous signal at 4 stations (E15, BETS, RITT and KUHL) representing 10 channels in total (for the first three stations, we consider the signals on the 3 components while at station KUHL only the vertical component is available) at each time step (0.01s). Template window signals are 2.56s long and start 0.5 s before the P-wave pick. Correlations are performed after filtering the template signal and the windowed continuous signal between 8-25 Hz. We use the overlap-add algorithm and the FFTW library for a fast implementation of the correlation computation (Harris & Paik 2006;Frigo & Johnson 1998) The selected stations and channels are the 4 closest stations from the injection point with the best signal to noise ratio and are distributed with a good azimuthal coverage relative to the injection point.
After the correlation coefficient is obtained at all individual channels for a given template event, these vectors are time-shifted according to the travel time differences between the stations. This step ensures that the resulting stack of the correlation coefficient vectors from these different stations is only coherent if the detected signal originates from a source proximal to the template event. In order not to be too restrictive using this procedure and possibly excluding new events not exactly co-located with the template event, we apply a maximum filtering over a duration of 0.1s before stacking the correlation signals at the different stations. The maximum filtering simply transforms the correlation signal by replacing every value by the maximum in its 0.1s range (Van Herk 1992; Gil & Werman 1993). It thus permits the detection of somewhat similar waveforms to the template but not located exactly in the same area. Stacking is performed over the 10 channels for each time step. We consider all parts of this averaged correlation signal where the correlation coefficient is higher than 0.4 as a possible detection.
We visually check that extracted waveforms do correspond to earthquake signals. Reducing this correlation coefficient threshold produce some new events that could not be ambiguously distinguished from noise. If multiple templates detect the same event we only consider one detection associated with the template for which the correlation coefficient is the highest over the considered time-window. We run this template matching approach for the seven days covering the injection period and the second swarm (from the 27 of June 2013 until 03 July 2013). This leaves us with 1395 events from the original set of 13 earthquake templates. We notably recover 93% of the original catalog's events. For all of these new detections we extract the seismic waveform on all stations around the P-wave arrival for 2.56 seconds. Our template detection approach confirms that only one earthquake occurred in the injection area during the nearly 4 days following the shut-in. Seismicity resumed on the 2nd of July 2013 with a strong burst with 225 events in this improved version of the catalog (See Figure 2).

RELATIVE RELOCATIONS
Because of the high similarity between the detected events ( Figure 3), it is possible to obtain a precise relative relocation of these events based on double-difference algorithm. We compute travel time delays between all the 1395 detected events. Delays are computed on 128 sample long windows filtered between 8 and 25 Hz and centered on the P-wave arrivals for vertical components, and on S-wave arrivals for horizontal components. We keep all travel time delays associated with a correlation coefficient higher than 0.6 for the relocation process. In total we observe a much higher signal to noise ratio and similarity of the signal on the horizontal components. This is attested by the number of delays finally retained: 7 255 735 from S waves (both E and N components) and 2 306 002 from P-waves.
The relocation is performed from the Software HYPODD (Waldhauser & Ellsworth 2000) where we set the initial position of all the events to be located on the open hole borehole trajectory at 2300m depth. We are able to relocate confidently 1393 out of the 1395 detected events.
The relocated dataset shows with finer details that the seismicity took place on two geometrically independent faults in the reservoir (Figure 4 and S1). We acknowledge that the overall geometries of the apparent structures we are resolving with the relocation can be slightly affected by the velocity model as this will modify the computed ray geometries (Kinnaert et al. 2016). There is a clear offset between these two quasi-parallel structures. The earthquakes taking place on the structure in the south are well fitted by a plane of azimuth N025 • E and dipping 74 • towards west ( Figure 5). The second structure cannot be as well fitted by a single plane. We find that the azimuth varies between N005 • E and N025 • E. The sparsity of the events on this second structure prevents us to obtain a reliable estimation of its dip (Figure 4). Events that occurred during the injection (induced events) are mostly located on the fault to the south and the events of the second swarm on the fault to the north. The open-hole section of the borehole extends between 1920 m to 2560 m but we observe that most earthquakes occur between 2200 and 2400 m depth. There exists a slight shift of nearly 20 m between the location of the first earthquakes and the borehole. Such an absolute shift is however not constrained from the relative relocation procedure and it is likely that the imaged structure actually intersects the borehole.
Based on borehole imaging of well GRT1 we constrain the first earthquake that occurred during the injection to coincide with the well trajectory at 2368 m depth where a major structure was observed (see section 7.1). The earthquake activity during the injection progressively migrates upward and to the south along the identified fault plane (See Figure 6). It thus mostly corresponds to an asymmetric migration relative to the injection point. At the end of the injection, the seismicity illuminates a planar 300 m long structure with a vertical extent of 200 m. The average along strike migration speed of the seismicity is of the order of 12 m.hr −1 . We observe that some events at the end of the injection, before shut-in, are located in the northern structure mainly activated 4 days later (see Figure 4 around latitude 48.8985 • N). Earthquakes that occurred during the second swarm appear to occur on a plane with quite a similar orientation as those activated during the injection. We also observe that there exists an approximately 100m shift between this plane and the one inferred from the location of earthquakes during injection (see Figures 4 and 6). It is possible that the two faults imaged by the seismicity are part of en-echelon system as commonly observed for faults in the Rhine graben (Schumacher 2002).
We also observe that earthquakes during both crises are all located around the same depth interval and that almost very few events are located below 2400 m. The second burst of events that occurred 4 days after shut-in started at 48.900 • N, at the middle of the visible structure where we observe a change of direction of the event locations. The earthquakes during this second burst appear to migrate quickly from the initiation point. This migration occurs in both directions from the first earthquake of this second burst. The depth range of the events in this second burst is similar to the events of the first swarm.

FOCAL MECHANISM
Obtaining a focal mechanism based on first motion polarities for the detected events is not an easy task as most of these events show a low signal to noise ratio. Furthermore, the station coverage is not optimal for deducing the focal mechanisms of the recorded earthquakes owing to the absence of stations in the south. However, it is still possible to estimate the mechanisms of a few large events with a higher signal amplitude. We focus on the two events that occurred during the injection with the largest signal to noise ratio at most stations and with the most clear first motions. At each station we visually determine the polarity of the signals. We compute the ray parameters for the two events based on their relocation. We observe that the inferred orientation of the fault plane from the spatial

RELATIVE MOMENT
We estimate the relative moments between all earthquakes of the sequence to gain insight into any possible moment variations during the investigated period. For all relocated events we employ a Singular Value Decomposition to compute earthquake's relative moment (Rubinstein & Ellsworth 2010).
The coefficients of the first basis vector for each event represent the amplitude of the coherent part between each tested waveform and can then be taken as representative of a relative moment. We employ this approach to compute the relative moments at 4 stations with a good azimuthal coverage around the injection well (E17, E15, BETS and RITT). The waveforms are not filtered and the singular value decomposition is performed on both horizontal components at each station. The final estimate of the relative moment is taken for each event as the median value obtained for all traces. All relative moments are normalized to the event with the lowest moment, here fixed arbitrarily to 1. We observe that the moment distribution for the events that occurred during the injection (first swarm) and during the second swarm are clearly distinct (Figure 7). Earthquakes that took place during the injection have lower moments. Their distribution can be described above the completeness moment by a high B value where B is defined as and N (M 0 ) is the number of earthquakes with a relative moment higher than M 0 and a is a prefactor.
We observe that the distribution of moments for the events that occurred during the second swarm is well compatible with a typical B = 2/3 value as observed for earthquakes worldwide (Scholz

Borehole Fault Zone and Seismicity Fault Plane
Well logging was performed in the well GRT1 both before and after the hydraulic injection. It reveals that a major structure associated with the largest temperature anomaly and shows that a significant permeability enhancement is located at 2368m depth. Acoustic imaging allows the determination of the orientation of this fault zone which strikes N175 • E, and dips 65 • W (Vidal et al. 2016). Geological logs show that this major fault zone intersects the borehole in the granitic basement close to the bottom of the upper altered zone of the granite. This major permeable structure (its estimated width is 24 cm) is then supposedly the main fluid pathway for the injected water. We thus hypothesize that the main structure oriented N25 • E and dipping 70 • W as evidenced by the earthquake locations could correspond to this main permeable structure imaged in the well bore. This is also in agreement with

Seismicity related to the hydraulic stimulation of GRT1 15
the first motion polarities. The slight variation of the azimuth of the inferred fault plane between seismicity and acoustic imaging could result from the uncertainty associated with the strike determination in the acoustic image but also from the different scale of the two types of measurement. Indeed, the orientation given by fitting the seismic cloud provides a large-scale global orientation while the borehole image only gives a local determination of the geometry of the fault structure. It confirms that the identified earthquakes do not occur as an isotropic 3D cloud but on a pre-identified planar fault.

Regional Stress Field and Induced Seismicity Onset
The state of stress in the investigated area has been extensively discussed in a number of studies mainly from measurements performed in the wells of the Soultz-sous-Forêts geothermal site (e.g. Klee & Rummel 1993). As this geothermal zone is located less than 7 km away from the current site and as no major active geological structure is located between the two sites, we will assume that they both are under the influence of the same regional stress field, i.e. as proposed in Cornet et al. (2007) P p = 0.9 + 0.0098z (6) where z is the depth (in meters) and stress are given in MPa. S H is the maximum horizontal stress, S h , the minimum horizontal stress and S v is the vertical stress. The direction of S H has been found to vary with depth in the area. As the analyzed earthquakes occurred over a limited depth range around 2300 m depth, we can consider that the orientation of S H is constant with an azimuth of N170 • E ±10 • (Cornet et al. 2007). The analysis of the orientation of the breakouts in the well GRT1 confirms that the direction of the maximum horizontal stress is nearly N-S and gives a direction of S H N180 • E (Vidal et al. 2016) The stress tensor of Eqs. 3-6 implies that it would require an overpressure in the well of 7.7 MPa for the fluid pressure to be higher than S h and thus create hydraulic fractures. The overpressure in the well did not exceed 3.1 MPa during the injection which ruled out such a scenario. Furthermore, the seismicity during this stimulation is not aligned in the direction of S H (perpendicular to S h ) which also confirms that the recorded events do not result from tensile hydraulic fractures. As we show below, the increase of overpressure first promotes the sliding on a pre-existing interface and at much lower stress that the one needed to propagate a mode I fracture.
The stress field of Eqs. 3-6 promotes the occurrence of strike-slip, normal faulting. For a typical friction coefficient µ = 0.85, the optimally oriented plane for rupture under a strike-slip environment has an azimuth of N25 • E. It therefore suggests that the imaged strike-slip fault corresponds to a nearly optimally oriented structure with a slip vector in agreement with the regional stress field. We can represent the Mohr diagram corresponding to the stress field at the depth of the fault plane where it intersects the borehole (2368m) (Figure 8). If the image log of the permeable structure at 2368 m actually corresponds to the same structure where seismicity is observed, as hypothesized, it then suggests that this structure was existing prior to the injection. We first consider the fault to be at the friction equilibrium defined by τ = µσ n , where µ = 0.85 is the typical friction coefficient at low normal stress (Byerlee 1978). We observe that even in the absence of perturbation caused by fluid injection, the fault is nearly critical. Only a very small amount of stress perturbation will promote slip on this interface. It suggests that the fault starts to slip at the very beginning of the injection, for a very small overpressure, as predicted by Byerlee's law. In the absence of seismicity it is suggested here that most of the slip occurred aseismically during this period.
We notice that the earthquake activity does not start immediately at the time of the injection. We Such a value of the internal friction coefficient close to 1.0 is typical for a wide variety of rocks (Carmichael 1982) . The orientation of the optimally planes for failure are 23.0 • from the direction of S H considering this internal friction coefficient. We compute that the fluid overpressure needs to increase to 2.65 MPa in order for the Mohr envelope to reach such a failure criterion (Figure 8).
The development of the induced seismicity during the injection very well obeys this failure criterion.

Seismicity Migration During the Injection
We propose two main hypotheses to explain the migration of the seismic activity that occurs during the injection. In a first model, we propose that the water pressure resulting from the injection is propagat- ing in an open fault. It leads to a stable crack propagation that is related to the volume of the injected water. If we suppose that this injected volume, V is expanding over a disk of radius r and thickness h, then V (t) = 2πhr 2 (t) where h is the thickness of the fault where the fluid is propagating and is supposed constant. The injected volume can be computed from the different flow rates Q(t) used in the various steps of the injection procedure, Combining Eqs. 7 and 8 leads to the expression of the radius of expanding injected volume Applying Eq. 9 to the actual injection history leads to results presented in Figure 9. The best fit of the model to the seismicity data is obtained when setting the thickness of the fault (the boundary of the injected volume) to 12 mm. We show that during the first part of the injection, when we observe the migration of the seismicity, the radius r defining the limit of the propagating crack in our model is expanding almost linearly with time.
Our second hypothesis to explain the migration of the seismicity is to consider that the pressure pulse is diffusing in a permeable medium. In this case, and hypothesizing that the applied perturbation to the medium can be approximated as a step function, Shapiro et al. (1997) showed that the pressure pulse is propagating away from the injection point as r(t) = √ 4πDt, where D is the hydraulic diffusivity of the surrounding medium. Applying this model to the seismicity recorded during the injection, we found that the best hydraulic diffusivity that fits the data is 7 10 −2 m 2 s −1 . This value is of the same order of magnitude to the value deduced for the nearby Soultz reservoir during the 1993 injection (5 10 −2 m 2 s −1 ) (Shapiro et al. 2002). We observe that the first hypothesis, despite its simplicity, explains the data at least as well as the second one. It suggests that the development of the seismicity on the fault plane can be well explained knowing the failure criterion of the fault, the actual state of stress and the injection history.

The triggering of the second swarm -Coulomb stress transfer from continuous slow aseismic slip
It is difficult to explain the occurrence of the earthquakes of the second swarm only from the diffusion of the injected fluid. Indeed, in this case we would rather expect a continuous seismic activity Seismicity related to the hydraulic stimulation of GRT1 19 over the entire period extending from the shut-in to the time of this second swarm. Such a prolonged activity after the end of injection is very often encountered in geothermal exploitation (Zang et al. 2014). Furthermore, as the over-pressure in the well reaches almost 0 MPa at the end of the 28th June it is difficult to propagate such a pressure front with no over-pressure at the source. This scenario is confirmed by the analysis of the moment distribution which points to a very fast decay of the moment distribution of earthquakes during the injection and to a much lower decay compatible with normal tectonic earthquakes for the events of the second swarm. It is tempting to link this change of the B value to the fluid overpressure as observed in various reservoirs (e.g. Bachmann et al. 2012). It then indicates that the events of the second swarm are certainly related to a low fluid pressure. These observations point out that another process is responsible of the observed activity. We test the hypothesis that a large scale aseismic slip on the fault plane imaged by the seismicity could have been responsible for the triggering of the activity during this second swarm. Aseismic slip in geothermal reservoirs have been reported from various observational evidences: shift of the borehole, velocity variations and repeating earthquakes (Cornet et al. 1997;Calò et al. 2011;Bourouis & Bernard 2007). The effect of fluid injection on natural interfaces has been found to promote aseismic sliding (Guglielmi et al. 2015;Wei et al. 2015;De Barros et al. 2016). Such a slow slip within the reservoir would be coherent with low amplitude events observed during the injection compared to the post-injection events (e.g. Lengliné et al. 2014)).
We suppose that the whole structure imaged by the seismicity during the injection defines a single plane that slips slowly during the injection. It thus defines a plane of surface A=200 m × 300 m, and we suppose that the movement of this aseismic slip is in agreement with the one reported by the focal mechanism of earthquake taking place on this plane (pure left lateral strike slip faulting).
From the dimension of this plane, an empirical relation for earthquakes gives a typical displacement of the order of 1 cm (Wells & Coppersmith 1994). We assume that such an empirical relation also applies to slow slip although the stress drop involved in such events might be lower than those of earthquakes (e.g Brodsky & Mori 2007). We note that it is difficult to constrain the actual size of the hypothesized aseismic slip with no broadband sensors at depth and that the approximate dimension of the plane given by the extension of the seismicity gives only a lower bound estimation. We thus look at the impact of uniform 1 cm of pure left-lateral slip on the defined fault in the variation of Coulomb stress changes on the nearby fault hosting the activity of the second swarm. As this second fault has a variable strike, we compute the Coulomb stress changes for 2 possible orientations of the receiver fault that spans the possible azimuth range (N005 • E-N025 • E) of the fault as observed from the seismicity.
In each case we consider that the receiving fault is vertical but we check that assuming a dip of 74 • , similar to the slipping fault, has no significant impact on the calculated stress changes. All stress changes are computed at a depth of 2300m, corresponding to the depth of the earthquakes of the second swarm. We observe that the second swarm is mostly located in a region of stress increase. This stress increase is the most important when computed on structures with an almost N-S orientation.
Such an orientation is similar to the strike of the northern branch of the fault that ruptured during the second swarm.
Although modest, these stress changes on the second fault can be responsible for the increase in seismicity observed 4 days later. Similar amplitude of stress increases are also found in the south of the slipping plane, but no earthquakes were recorded in this area. It is possible that no fault exists at that location or that faults are badly oriented for rupture such that the stress changes do not lead to an increase in seismicity. Several explanations can be invoked to explain the delay of the seismicity of the second swarm. First, the aseismic plane has been progressively slipping since the injection started and continued to slip after shut-in such that it reaches a 1 cm cumulative slip only after several days and then triggers the earthquake activity. However this explanation is not very plausible because it will be difficult to promote slip on this fault plane after shut-in as the fault pressure goes quickly to zero and thus moves the interface away from failure. A second hypothesis is that all the aseismic slip actually took place during the injection. We note that the stress increase on the second fault is only modest ( Figure 10). If this second fault is initially in the same state of stress as the one related to the injection, then it shows that the stress change of ∼ 2 bars is not large enough to promote instantaneous failure on this plane. However, we can consider that this stress step actually causes a clock advance of the nucleating asperities on the fault plane that is slowly loaded (Gomberg et al. 1998;Perfettini et al. 2003). This stress step makes the frictional unstable patches accelerate instantaneously, such that they get closer to failure but still need more time before they actually reach instability. Then as the first event on this second plane occurs we can suppose that most of this part of the fault is now close to failure such the extra stress change caused by the elastic stress redistribution of the event itself triggers this short lived burst of activity during the second swarm.

IMPLICATION FOR POST-INJECTION SEISMICITY
Post-injection seismicity occurs frequently in many cases of induced seismicity (e.g. Häring et al. 2008). In most of these instances, the seismic activity following the shut-in is continuous and does not stop and resume at latter times as observed during the GRT1 stimulation. It is then difficult to discriminate if the post-injection seismicity results from the only effect of the pore-pressure diffusion or if aseismic movements within the reservoir could also participate to the triggering of earthquakes. The unique GRT1 induced seismicity sequence suggests that the post-injection events are well explained as the result of such slow slip events. It is therefore appealing to propose that aseismic movements are ubiquitous features of stimulated reservoir and can be responsible for the observed post-injection seismicity.

CONCLUSION
Using dedicated tools (template matching, relative relocation), we were able to obtain a precise image of the fault structures activated during the injection. The induced seismicity related to fluid injection appears to agree with a geomechanical model but it requires detailed knowledge of the investigated area: the existing fracture and fault network, the regional stress field amplitudes and orientations, the failure criterion for the existing fractures and faults, and the pore pressure at depth during the injection.
It also requires a detailed analysis of the seismicity with precise location. Such information allows an estimation of when and where the seismicity would appear in the geothermal reservoir. As aseismic movements appear to have a significant role in the deformation of the reservoir and can redistribute the stress locally it is also important to monitor the displacement associated with such movements at depth.

22
O. Lengliné, M. Boubacar and J. Schmittbuhl Unfortunately, geodetic measurements at the surface might not be sensitive enough to resolve the small displacement on these structures. Some downhole broadband instruments or dedicated geodetic monitoring might be needed in the future to capture these slow movements.

ACKNOWLEDGMENTS
We deeply thank F. Cornet Dolatshahi for grammatical assistance. We thank two anonymous reviewers for their comments. We also thank ECOGI for funding part of the monitoring network and for sharing some geoscientific data.
This work was conducted in the framework of the Labex G-EAU-THERMIE-PROFONDE, which is co-funded by the French government under the program Investissements dAvenir, a cross-border collaboration between industrial companies, ECOGI, and ES-Gothermie, academics (KIT, EOST) and the EEIG Heat Mining of Soultz-sous-Forłts. We also thank the Geophysical Instrument Pool Potsdam from the GFZ German Research Centre for Geosciences (Germany) for providing a set of temporary seismological units.