Seismic asperity size evolution during ﬂuid injection: case study of the 1993 Soultz-sous-Forˆets injection

SUMMARY The injection of ﬂuid in the upper crust, notably for the development or exploitation of geothermal reservoirs, is often associated with the onset of induced seismicity. Although this process has been largely studied, it is not clear how the injected ﬂuid inﬂuences the rupture size of the induced events. Here we re-investigate the induced earthquakes that occurred during an injection at Soultz-sous-Forˆets, France in 1993 and studied the link between the injected ﬂuid and the source properties of the numerous induced earthquakes. We take advantage that deep borehole accelerometers were running in the vicinity of the injection site. We estimate the moment and radius of all recorded events based on a spectral analysis and classify them into 663 repeating sequences. We show that the events globally obey the typical scaling law between radius and moment. However, at the scale of the asperity, ﬂuctuations of the moment are important while the radii remain similar suggesting a variable stress drop or a mechanism that prevents the growth of the rupture. This is conﬁrmed by linking the event source size to the geomechanical history of the reservoir. In areas where aseismic slip on pre-existing faults has been evidenced, we observed only small rupture sizes whereas in part of the reservoir where seismicity is related to the creation of new fractures, a wider distribution and larger rupture sizes are promoted. Implications for detecting the transition between events related to pre-existing faults and the onset of fresh fractures are discussed.


I N T RO D U C T I O N
Induced seismicity has significantly increased during the recent years notably because of fluid injection at depth in particular related to hydrofracking and wastewater disposal (Ellsworth 2013). This has led to a significant increase concerning the risk posed by large induced or triggered earthquakes (Dahm et al. 2013;Grigoli et al. 2017;Candela et al. 2018). There is for example a seismic hazard related to the exploitation of deep geothermal system (Giardini (2009), as after the recent M5.4 Pohang earthquake (Grigoli et al. 2018;Kim et al. 2018). What are the controlling factors that favour the development of large, potentially destructive induced or triggered earthquakes, is still an open question that is very important for many industrial operations involving injection of fluids into subsurface (e.g. conventional and non-conventional oil production, geothermal energy exploitation, waste disposal, etc., McGarr 2014; Guglielmi et al. 2015;Elst et al. 2016).
While large induced events have to be avoid in order to guarantee safe operations of reservoir exploitation, small induced events still offer the opportunity to study the mechanical response of these fractured reservoirs submitted to intense perturbations. For instance, deep EGS geothermal reservoirs present a specific context for fault behaviour analysis since they are at high temperature and relatively at low normal stress (reservoirs are typically shallower than 5 km depth at 150-200 • C, Cornet et al. 2007;Huenges & Ledru 2011;Olasolo et al. 2016;Vallier et al. 2019). Indeed, this context favours the transition from frictionally stable or conditionally stable regime with velocity strengthening behaviour to zones frictionally unstable with velocity weakening regimes that are responsible for seismic ruptures (Scholz 1998). It thus allows the investigation of the link between earthquakes and slowly deforming interfaces that are typically studied usually at much greater depth at the bottom of the seismogenic zone where the brittle-creep transition occurs (Beroza & Ide 2011). Furthermore, because the imposed loading on the reservoir can be controlled from the wellhead injection rate or pressure and because a dense instrumentation can be deployed at proximity of the reservoir, it provides both information on the forcing term and high resolution observations on induced effects, both being favourable for a better knowledge of the fracture and fault behaviour. Fluid-injection operations at many industrial sites can be seen as large scale experiments (between the lab-scale and the scale of a developed fault system) that can help to decipher the Seismic asperity size evolution during fluid injection 969 links between fluid flow, slow aseismic movement and earthquake nucleation (Bernard 2001).
For example, investigating the induced earthquakes during an almost 1-yr-long circulation experiment in the Soultz-sous-Forêts geothermal site in 2010, Lengliné et al. (2014) evidenced that rapid and large changes of the seismic slip was observed on some repeating seismic asperities while their size remains almost constant. The related variable stress drop was attributed to the role of fluid changing the effective normal stress and the frictional stability condition over the interface where the seismic rupture is taking place. Other instance of variable stress drop events have also been reported (Lin et al. 2016;Staszek et al. 2017). In each case, the link with a potential role of fluid cannot be unambiguously drawn. This is because hydraulic data might not be present and because of the limited number of earthquakes analysed in each case.
Because most of the induced events related to injection in a geothermal reservoir have a small magnitude, capturing these events and achieving a detailed analysis of this seismicity, requires an appropriate instrumentation. Notably, borehole instruments, when available, provide a wealth of data compare to seismic surface network. This allows a better description of the ongoing deformation process that happens during the injection. Here we revisit the 1993 Soultz-sous-Forêts stimulation of the deep GPK1 well. This injection has been already extensively studied. It produced an abundant seismic activity (e.g. . A particular characteristic of this injection experiment is that repeated borehole imagery, before and after the injection, identified planar structures intersecting the borehole that slip aseismically . This was also confirmed from the analysis of the repeating earthquakes that took place on one of this structure whose cumulative slip matches the one observed from borehole measurement (Bourouis & Bernard 2007). A dedicated network of borehole sensors was operating at the time of the injection providing an unique opportunity to see the details of the seismicity and to investigate the role of fluid on the source properties of the induced events (Baria et al. 1996).
Here we study the sources of the microseismic events that were recorded during the 1993 episode in order to detect any possible links between the source parameters of the seismic events and the injection history. In particular, we investigate if the events possibly associated with aseismic slips within the reservoir have a distinct signature from the events that result from the creation of new fractures that appeared at later times and in a different part of the reservoir. For this, we first computed source properties for all events recorded, using a typical Brune's spectral model (Brune 1970). We then isolated repeating event sequences on individual asperities and obtained refined source parameters for these events. We show that the scaling of these events obeys globally the typical scaling law between radius, r and moment, M 0 . Interestingly, when events are analysed at the scale of the asperity, this scaling law is not fulfilled and the moment M 0 is evolving independently of the radius r over a range that explains the typical broad dispersion of such (r, M 0 ) scaling law (Kanamori & Anderson 1975). It suggests that an external factor is influencing the size of the seismic rupture at the asperity scale. We also observe that events size is on average changing over the course of the injection and for different parts of the reservoir. We interpret this control of the event size as a result of the fluid pressure on the rupture mode.

S E I S M I C I T Y
We analyse the induced seismicity associated with the hydraulic stimulation of Soultz-sous-Forêts, France that occurred between September and October 1993 in two phases. The hydraulic stimulation was performed through the injection well GPK1. The injection process lasted from 01 September 1993 to 16 October 1993 with an interruption between September 17th and October 11th. The injection flow rate was increased by steps of 6 l s -1 every 2 days, reaching 36 l s -1 at the end of September episode and about 50 l s -1 at the end of October episode (Fig. 1a). In total, approximately 44 000 m 3 of fluid was injected at a depth between 2850 and 3400 m (3550 m in October). We show the daily number of events together with the injection flow rate in Fig. 1(b), where the link between seismicity and the injection process is clearly evidenced.
About 15 000 events were recorded by 3 accelerometers and 1 hydrophone installed in 4 deep wells (at depths between 1400 and 2000 m), in close proximity to the injection well, over the 2 months of the experiment (Figs 1c and d). While the gap of seismicity between September 17th and October 11th is due to the injection shut in, several gaps of seismicity are observed on September 7th, 11th and 12th and are instead due to a loss of data. Each accelerometer is a 4-component instrument, that has a flat response up to 1 kHz and is oriented with one vertical component and the other three are 109.5 • away from each other and from the vertical direction. Conversion from 4 components to 3 components was performed during post-processing by Gaucher (1998). In this study, we use the converted 3 components for analysis. The data acquisition was performed in triggered mode and all seismic events were recorded on 2-s-long triggered windows sampled with a frequency of 5 kHz.
Event locations were firstly determined by the Camborne School of Mines Associates (Jones et al. 1995) from P and S waves arrivals and assuming an homogeneous velocity structure (V p = 5.85 km s -1 , V s = 3.34 km s -1 ). In a second step, relative location of events was performed based on travel time differences by Bourouis (2004). These latter locations are used in this study. As the stimulation starts, the seismicity is generated all around the injection source, at a depth of about 2950 m. The seismicity cloud grows gradually during the stimulation, with a distribution of events along subvertical planar structures Moriya et al. 2003;Evans et al. 2005). A few days after the stimulation started, the seismicity is observed to migrate principally to shallower depths (Cornet & Morin 1997). The final extent of the seismicity cloud reaches a distance of 800 m from the injection well at the end of the stimulation.

Events characterization
We first identify the P-wave arrival times for all events and at all stations through a STA/LTA approach (Earle & Shearer 1994) applied to all signals filtered in a common frequency band between 80 and 270 Hz (STA = 0.02 s, LTA = 0.3 s and threshold = 3.5). We then select a 0.25-s-long signal window surrounding the P-arrivals (0.05 s before and 0.20 s after) and computed the acceleration amplitude spectra. The spectral energy of the events ranges between 80 and 500 Hz for the stations 4550 (see Fig. 2) and 4616 and hydrophone EPS1, and up to 270 Hz for the furthest station 4601. At higher frequencies, the spectral amplitudes decay rapidly due to the attenuation effect.
We model the acceleration amplitude spectra of the events with the Brune's model (Brune 1970) to estimate absolute values of the seismic moments M 0 and the corner frequencies f c of each event (see Fig. S1). We use data from the seismic stations 4550 and 4616, the closest to the seismicity cloud and apply the analysis to the 80-500 Hz frequency band to hinder the noise contamination at   Table S1 and Fig. 1). On the right the associated acceleration amplitude spectra using the same colour code for the time occurrence of the event and an arbitrary absolute scale (A.U.). The strong similarity of the spectra shows the similarity of the event waveform. The dashed black line corresponds to the averaged background noise spectrum computed on the 0.25-s-long windows preceding the events. lower frequencies and the attenuation effects affecting the higher frequencies. We fit the acceleration spectra with the function S, where f, is the frequency, o is the high frequency level of acceleration amplitude spectrum. The term H n (f), in eq. (1), represents a noise term estimated from the analysis of signal windows preceding the events. The values of f c and o that best fit the amplitude spectra are found through a grid search algorithm. We iterate for the whole dataset of seismic events and provide for each of them the estimation of the absolute values of o and f c and the associated uncertainties. These uncertainties are approximated from the width of the objective function, which, under the hypothesis of a normal distribution, Seismic asperity size evolution during fluid injection where ρ is the density of the propagation medium, V p is the Pwave propagation velocity. We here use ρ = 2350 kg m -3 and V p = 5850 m s -1 as taken from Gaucher (1998). The source-receiver distances, r, are estimated from Bourouis (2004) and R p is the radiation pattern, we consider for this last parameter, a mean value for P waves (R p = 0.52).
Assuming a circular rupture model, we then infer the source dimensions from the corner frequency estimations through the following relationship (Madariaga 1976): where V r is the rupture velocity, V s is the shear wave velocity, V s = 3340 m s -1 taking V p /V s = 1.75 from Gaucher (1998) and c is a constant (c = 0.32, Madariaga 1976).

R E P E AT I N G E A RT H Q UA K E S I D E N T I F I C AT I O N
In this section, we identify families of seismic events that are generated by repeated ruptures on the same asperity. These repeating events are not only characterized by highly correlated seismic waveforms but also by sharing a common source location and source mechanism. Forming families of repeating events allows us to use relative methods leading to more accurate estimates of the difference of source parameters. All changes that are resolved at the scale of the repeating earthquake sequence, are uniquely associated to temporal variation of the event sources. We identify those sequences of repeating events in a two step procedure. We first gather events on the basis of high waveform cross-correlation (i.e. multiplets). We then select events in each family of similar waveform events, those that are colocated and reject events that are not. For this second step, we (i) estimate the source dimension through a spectral analysis; (ii) assess the inter-event distance from t S − t P time delays and (iii) verify that the event rupture areas of the sequence do overlap by checking that the source dimensions are larger than the inter-event distances.

Cross-correlation analysis and clustering
We build a cross-correlation matrix where elements are given by the maximum of the cross-correlation function (CC) computed between two 0.55 s-long-windows (0.05 s preceding the STA/LTA picking time and 0.5 s after, or 2750 points long) of waveforms filtered between 80 and 270 Hz. This frequency band is chosen to maximize the signal-to-noise ratio of the recorded waveforms. The use of 0.55 s-long windows enables to encompass the P and S waves in Reccurence times (days) Figure 5. Probability density functions (pdf) of differential times for all events in the catalogue (green circles), for all events in a repeating sequence (purple triangles). The orange crosses represent the pdf of differential time within repeating sequences. The colour lines represent the best-fitting gamma distribution for each population. the seismic waveforms. Since the analysed signals are acceleration records, the resulting correlations emphasize the highest frequencies of the signal, compared to the use of correlations of velocity measurements. We then cluster events according to the correlation coefficient computed from the vertical components for the station 4550 that is the closest to the injection source. The clustering threshold for the normalized correlation coefficient is set to 0.75. We test several correlation coefficient thresholds and keep 0.75 as a good compromise between a too weak selection of events that actually merges events of different sequences and a too tight selection that separates events that are part of a same sequence. We chose an open and non exclusive algorithm for building clusters: in a set of events, only one connection (event pair with CC > 0.75) is enough to cluster with other events. In this manner 6970 events were classified into 1941 sequences of similar events.

Event size estimation
We apply a spectral ratio method to compare corner frequencies and seismic moments between pairs of events of each sequence . The method provides relative information which has a higher accuracy than the comparison of absolute estimates using a Brune's spectral fitting of all individual event waveforms. The spectral ratio for a pair of events is estimated through the following expression: where X 1 and X 2 are the Fourier transform of the signals of two events, estimated on 1024 sample-long window taken around the P-wave arrival and smoothed with a Tukey window. The * symbol represents the complex conjugate. In the case of two events with f c1 f c2 f c , the difference in corner frequencies is deduced from the slope a of log G(f) through the relationship (Got & Frechet 1993;Lengliné et al. 2014 For each pair of repeating earthquakes, we therefore compute log G(f) from the spectral ratio averaged over the 3 components (see an example in Fig S2). We then estimate the corner frequency differences and the associated uncertainties by a linear fitting of the slope of log G(f) in the domain [f c1 , f c2 ]. Using eq. (3), we then infer the variation in source dimensions r = R 1 − R 2 according to: Similarly, the ratio of seismic moments between two events in a pair is inferred from the spectral ratio: in the low frequency band when noise is negligible. At frequencies lower than f c , the spectral ratio is expected to be flat and its level is then used to estimate the relative moment when neglecting the second term in eq. (7). However, we observe that the noise contribution is relatively important at low frequencies (typically below 80 Hz), hiding the signal information and preventing us from using the simplification of the second term in eq. (7). We thus use the whole expression in eq. (7) and compute the spectral ratios in a common frequency range, 80-200 Hz, for all the repeating earthquakes. In this frequency range, the signal is indeed well above the noise level. We then infer from this ratio and the estimates of f c1 and f c2 , the moment ratios of each pair of events. We apply the spectral ratio technique to all pairs of events belonging to a common sequence. This produces a large amount of relative measurements ( r and moment ratio) since radius and moment differences are computed for each possible pair of events in a sequence. To take advantage of this large set of data, we develop a minimization approach. Estimating the source properties (radius or moment) m from the relative measurements is done by resolving the linear problem d = Gm, where d contains the relative measurements and G is a matrix composed of all combinations of event pair differences. This minimization approach (in the least square sense) is performed by using the following general solution for inverse problems (Tarantola 2005): where [m prior , C M ] is the a priori information (source properties and covariance matrix) obtained from the Brune's spectral fitting, [d, C D ], the relative measurements and associated uncertainties from the spectral ratio analysis. The relative measurements contribute to reduce the uncertainties on the estimation of the source properties (see Fig. S3).

Interevents distances
We estimate distances between earthquakes in each sequence of similar events through a high resolution cross-correlation analysis to get a fine measurement of S−P traveltimes (Poupinet et al. 1984;Bouchon et al. 2011). We assume that for each event pair of the sequence, their interdistance is much smaller than the distance to the receiver. We slide a 128 sample-long window on the pair of signals along the entire signal duration, with steps of 64 samples. At each step, we compute a cross-correlation function between both signals and identify the time lag associated to the maximum of the function. We then determine T P and T S as the variation in P and S arrivals between the pair of events. This variation is zero when the two investigated events have an exact common origin. T S−P = T S − T P is related to the distance d between both events through the following expression: assuming a common ray path for both P and S waves and where k V = V P −Vs V P V S and θ is the angle between the vector connecting the sources of the two events and the direction to the receiver.
We compute T S−P for each pair of events in each sequence and by averaging the results on the three components and iterating for the three stations 4550, 4616 and 4601. For each event pair, we thus obtain three measurements, one corresponding to each station. Finally, we select the highest distance between the three distance estimates in order to approach the case of cos θ = 1.

Overlapping criterion
Finally, in each earthquake sequence, we compare the sum of the radii and the inter-event distance for each pair of events. An overlapping criterion is introduced and defined as an inter-event distance smaller than the sum of event pair radii. We only keep events in each sequence when their source overlaps with at least two other events of the sequence. This last criterion is useful to ensure that the soformed repeating sequence represents a compact set of individual ruptures without mixing closely located but distinct repeating sequences. In this manner, we classified 4252 seismic events into 663 sequences comprising 3 or more repeating earthquakes (doublets are excluded in this selection).

Scaling of repeating events
We investigate here the moment-radius (R, M 0 ) scaling of all identified repeating events (see Fig. 3). We first observe in this figure that repeating events have a radius R, distributed between 2 and 13 m and a moment M 0 in the range [3 × 10 7 − 3 × 10 10 , N m −1 ]. We however keep in mind that these absolute values depend on arbitrary parameters such as the rupture velocity but also on hypothesis such as the circular shape of the rupture zone (Aki 1967;Kanamori & Anderson 1975). The scaling between both parameters is expected to be of the form M 0 ∝R 3 as observed from numerous earthquakes worldwide (e.g. Kanamori & Anderson 1975). Indeed, following Eshelby (1957) and Madariaga (1976), when assuming a circular crack model for the earthquake rupture, the seismic moment M 0 , scales as: where σ is the average stress drop. We see from Fig. 3 that globally the scaling (eq. 10) is fulfilled and most events fall in a range of average stress drops between 0.2 and 20 MPa. Again the absolute values of the stress drops are dependent on arbitrary choices like the rupture model and its parameters. We however observe that a dispersion of the stress drops exists. Interestingly, when we investigate at the scale of a single repeating earthquake sequence, we obtain a very different result. We show in Fig. 3 several examples of repeating earthquake sequences using coloured circles (one colour per sequence). For each sequence, we observe a small variation of R for a large evolution of M 0 , very differently from the classical scaling: M 0 ∝R 3 . It has to be noted that the variation in moment within a repeating sequence is not a monotonic time evolution but is randomly distributed in time with no clear history of the moment evolution (see Fig. 4). When we are analysing a single repeating sequence, we are dealing with relative measurements of the moment and radius variations that are well constrained and sensitive. It thus suggests that most of the dispersion around the theoretical moment-radius scaling is taking place at the scale of the repeating earthquake sequence itself.
To complement the description of the repeating sequences, Fig. 5 shows the distribution of recurrence times between consecutive events for all events of the catalogue and only events that belong to a sequence. Both populations are fitted by a gamma-distribution: with a slow decay [i.e. a large γexponent: γ ≈ 0.7 (Corral 2004;Davidsen et al. 2007)] and a short time cutoff (B ≈ 1 hr). We also computed the same, probability density function for recurrence times of events within repeating sequence but treating in this case, each sequence independently and not mixing times of the different sequences. We now recover a best-fitting gamma-distribution with a higher decay exponent (i.e. a smaller γ exponent: γ ≈ 0.13) and a larger time cutoff (B ≈ 8 d, comparable to the duration of the stimulation). This suggests that the events in sequence have a different temporal organization compared to the overall behaviour of the seismicity. Events within sequences are shown to have a strong and long-time correlation but the initiation of the sequence is similar to the whole seismicity with almost no memory, close to a Poissonian process.

Scaling within repeating sequence
In order to investigate more systematically the moment-radius relationship that we observe within a repeating sequence, we normalize seismic moments and source radii within each sequence by their respective average values in the sequence. We plot the relationship between these two parameters in Fig. 6. We observe that the normalized moments within a sequence span two orders of magnitude while at the same time there is no significant variations in the source radius. This was already observed during the 2010 water circulation test in Soultz-sous-Forêts ) and for repeating micro-earthquakes sequences in Taiwan (Lin et al. 2016). We do however observe a linear trend between these two parameters in Fig. 6, that is a non-vertical fit with a large slope (close to 20). It suggests that there might exist at the scale of the asperity, a different link between these two parameters. The dispersion that we observe in Fig. 6 is then explained by this non-classical behaviour at the scale of a repeating asperity. It possibly implies a mechanism that controls either the radius of the rupture event (quenched disorder on the fault plane, stress distribution), or/and the proportion of seismic slip taking place on a given asperity.

Two populations of asperities
In order to look at the possible mechanisms that control the variations of the source parameters, we investigated how they do change over the course of the injection. Here we concentrate our analysis on the first part of the stimulation in September 1993 (from September 1st to September 20th). We report in Fig. 7 the inferred radii of the repeating earthquakes as a function of time. We observe a spreading of the radii between 2 and 13 m. This spreading is actually increasing few days after the beginning of the injection (on September 6th, coincidentally with an increase of the flow rate) and then remains stationary during the rest of the stimulation. The average radius is also showing two regimes. We observed that during the first few days, most of the earthquakes have a small dimension (typically the inferred radius is rather constant around 4 m). In a second period, the average radius of the repeating earthquakes evolves to 6 m and shows a plateau for the rest of the stimulation. We identified two representative periods of the same duration (2 d) during these two regimes: one around the 5th of September and one around the 14th of September. In Fig. 7, we detailed two properties of the repeating events during these two periods: the distribution of their radius and the distribution of their magnitude. During the first period, the radius distribution is peaked around 4 m which is significantly different from the second period when it appears larger and broader. Interestingly we observe that the magnitude distributions are similar as a Gutenberg-Richter power law distribution with the same b-values [b ≈ 1.5 evaluated using the maximum-likelihood method of Ogata & Katsura (1993)] but a different prefactor (i.e. a-value). We conclude that there exists two populations of events that can be hardly distinguished by a b-value monitoring but better evidenced by a source radius characterization.
The transition between both periods dominated each by a specific population of repeating events (small and peaked for the first period and larger and broadly distributed for the second), is actually related to a specific transition in the flow rate history: the increase from 6 to 12 l s -1 on September 6th, the 5th day of the stimulation, as shown in Fig. 8. Indeed this transition marked by a red line in the figure, corresponds to the increase of the average source radius from 4 to 6 m (Fig. 8c) and can be readily observed in the waveforms of repeating sequence (Fig. 9).
In Fig. 8(b), we compare this time evolution to the depth evolution of all the events. We see from this figure the transition from 12 to 18 l s -1 on Sept 8th where the seismicity starts to rise from the injection depth (2900 m) to shallower zones (2000 m deep, cf. yellow vertical line in Fig. 8). The transition has been interpreted as the onset and growth of a quasi-vertical fresh hydraulic fracture Cornet 2016). Here we show that the population of repeating events existing during the hydraulic fracture growth was existing before this 12-18 l s -1 increase, that is to say at the 6-12 l s -1 transition highlighting that the process of nucleating a new fresh fracture has started earlier. As shown in Fig. 8(d), the average radius change is not related to a change in the proportion of repeating events in the whole population of seismic events. The Gutenberg-Richter distribution is characterized by two parameters for magnitudes M larger than the completeness magnitude: log 10 N(M) = a − bm. In Figs 8(e) and (f), is shown the time evolution of a and b using several sliding window varying between 300 and 600 events. We clearly see a sharp increase of the a-value at the 6-12 l s -1 transition and fluctuations of the b-value around an average value of b ≈ 1.5. The b-value for the sequence events is typically slightly larger than for the whole seismicity. It confirms the result of Fig. 7 that each population of sequence events share a similar b-value but a different a-value. Computation of the seismogenic index, as defined by Shapiro (2018) reveals an increase at the time of the flow rate increase. However this increase is not easily notable because of the lack of data just after this interval that produces an apparent decrease of the seismogenic index. We emphasize that our analysis of the evolution of the a-parameter corresponds to incremental Gutenberg-Righter distributions and not to the cumulative population of events as done for the seismogenic index.

Evolution of source parameters with time and depth
We see from Fig. 10 that at the beginning of the injection, most of the seismicity is clustered at a depth around 2900 m. This depth corresponds to the upper limit of the open-hole section (2850-3400 m) where injection occurred. We observe that as the injection progresses, most of the earthquakes take place around a restricted depth range (2850-3150 m). They present a small source dimension as confirmed by the radius distribution that is peaked at 4 m independently of time and significantly different from distributions in other (a) . From top to bottom: (a) Acceleration spectra computed on a window encompassing the P-wave arrival for events of two repeating sequences (2 shades of blue) at station 4550. For each sequence we observe that the corner frequency is very similar and is around 150 Hz for the events of the first group (light blue) and 340 Hz for the events of the second group (dark blue). All events in the second group are associated with a low radius and they all occurred in the first 3 d of the injection. All the events of the first group took place after the injection step from 6 to 12 l s -1 . (b) Normalized accelerations records over the full event duration for the two sequences at station 4550. (c) Zoom on the waveforms at two different time intervals. All records of a sequence are superimposed and are very similar.
depth intervals (Fig. 11). Immediately below this depth interval, in the depth range (3150-3500 m), repeating events are initiated later with a significant increase of the event size as a function of time. In the part of the reservoir just above the injection zone (2700-2850 m deep), a similar history is obtained. Earthquakes have a small radius (smaller than 5 m) at the beginning of the time series but show an abrupt increase of these event radii when the flow rate is increased from 6 to 12 l s -1 as discussed in the previous section. After this change of the flow rate, earthquakes in these two depth intervals, tend to have higher and broadly distributed radii dominated by the second population type of events. Finally, the last depth interval we consider includes the upwards migration of the events towards shallower depth (2000-2700 m). This migration occurs after the increase of the flow rate to 18 l s -1 . These shallower events are the ones with the largest radius illustrating the second population of repeating events. No statistical change of radius is observed for these events from their first occurrences up to the end of the injection (Figs 10 and 11). As a summary, we conclude that as the injection progressed we clearly evidence an increase of the source dimension of the earthquakes that occurred in the reservoir, but we observe that this increase does not occur everywhere in the reservoir: (i) the depth section around the open hole (2850-3150 m) hosts mostly the first population of small and peaky distributed events over the whole course of the injection; (ii) the upper domain (2000-2700 m) corresponds to the second population of large and broadly distributed events and (iii) the two other domains show a mix between these two populations.

Possible role of the geological complexity
Because of the fluid diffusion during the injection and the associated migration of the seismic front, we could wonder if the observed evolution in time of earthquake radius (Fig. 10), is simply explained by a constant triggering effect of the fluid but entering new areas where exist different quenched mechanical properties (e.g. structural and rheological heterogeneities related to different lithologies) leading to different rupture sizes. This could be a simple explanation for the observed change of source radius over the course of the injection and would not require any direct influence of the fluid pressure on the source parameters. This hypothesis could be easily dismiss because if the earthquake source dimensions were controlled uniquely by the size distribution of pre-existing asperities, one should not observe any variation of the earthquake radius as a function of time in a given area. Indeed, in this case, time (or similarly all time-related quantities such as volume of injected fluid or pressure) should not be linked to any variation of the source size. However we observe for several zones, clear changes of the event radius with time (see Fig. 11). For example, the zone just above the injection level, shows a significant increase of the earthquake radius just at the time of an injection step (12 l s -1 ). Furthermore this hypothesis requires a very specific distribution of the earthquake size in order to guarantee the observation of an increase over the course of the injection, which seems quite unlikely. This simple observation ruled-out the possibility that the inferred earthquake dimensions are purely controlled by the size of pre-existing faults within the reservoir, the quenched disorder.

Size-limited events and aseismic slip on pre-existing interfaces
It was previously evidenced that several fault segments intersecting the borehole at the injection depth were sheared during the 1993 episode at Soultz ). The measured displacement Figure 10. Evolution of the earthquake source sizes (radius of the rupture area) as a function of time and depth (colour circles). Only events within identified sequences are reported. The size and colour of the circles refer to the dimension of the event. The blue vertical lines represent the flow-rate steps as introduced in Fig. 8. The horizontal dashed lines represent the various depth intervals chosen for analysing the seismicity (see Fig. 11). on these interfaces is significant and cannot be explained by the seismic slip of a single large earthquake. It has been interpreted that this observed displacement on these faults was caused by aseismic movements. Bourouis & Bernard (2007) consider the particular case of one of this fault, called 'fault F', intersecting the borehole at 2925 m deep. They showed that repeating events occurred along this fault and that the cumulative displacement on the identified seismic asperities matches the offset on this fault recorded at the borehole. It implies that these repeating events represent shear rupture of seismic asperities embedded in an otherwise creeping interface. Similarly to Bourouis & Bernard (2007), we interpret the events that took place close to the injection well at the depth interval 2850-3150 m [corresponding to the section where 80 per cent of the flow rate is lost ], as associated with slow aseismic movements on pre-existing faults within the reservoir. This is in agreement with the mechanical interpretation of Cornet (2016) that proposes that during the first injection steps, the fluid pressure is sufficient for the stress state in the reservoir to reach the slip onset condition on pre-existing structures. The orientation of the seismic cloud as shown in Fig. 11, is then dominated in this regime by the orientation of pre-existing faults (N150 • E, Cornet et al. 2007). From a frictional point of view, because of the reduced effective normal stress (caused by the water injection), and the high temperature at that location, aseismic movements are promoted (Scholz 1998). The recording of events having only small source radius in this time-space window (see Fig. 11, right panel for 3150 > z > 2850) indicates that slow aseismic slip has a distinct signature in terms of source properties: limiting the earthquake size that translates in peaked distribution of small radii. We argue that because this part of the reservoir seems to deform mostly in a ductile regime, most of the stress is relieved constantly by the aseismic movements. This reduces the opportunity for dynamic ruptures to grow to large size because there is not enough stress to drive the rupture over large distances. We note however that this does not totally rule out the possibility of large rupture in this part of the reservoir as pointed out by Wei et al. (2015) but that such occurrences are more unlikely. Cornet (2016) proposed that the shallower earthquakes during the 1993 injection episode, are related to the generation of new fractures. Indeed, after having reached the Byerlee's friction criterion on pre-existing faults around the open hole domain during the first part of the stimulation, the reservoir is entering a second phase with the increase of the fluid pressure corresponding to failure in the bulk of the reservoir following a Mohr-Coulomb or more probably a Hoek and Brown criterion (Jaeger et al. 2009;Villeneuve et al. 2018). This failure process is evolving in two steps as the effective pressure decreases. First, shear hydraulic fractures are created when the Mohr-circle of the stress field intersects the failure criterion and then tensile hydraulic fractures are generated when the minimum principal stress reaches the tensile strength of the reservoir. The later process is unstable as the effective pressure reduces vertically as the crack grows upwards at constant injection pressure (Cornet 2015). The initiation of fresh fractures (first in shear mode and lately in tensile mode) corresponds to rotation of the seismic cloud orientation with shallower depths (Fig. 11, Cornet et al. 2007;Cornet 2016). It then favours the interpretation that the repeating earthquakes observed at these locations are the results of strain localization during the coalescence phase of small subfractures to obtain a large scale fracture in the reservoir (Lockner et al. 1991;Zang et al. 2000). This second stage of the injection is proposed to be associated to a widening of the event radius distribution. Finally, it suggests that the mechanism at the origin of the seismicity (shear on existing fracture or nucleation of a large fresh fracture) has an influence on the dimension of the rupture size of the seismic events. We note that the injection step that produced the most significant variation of the event's radius is when the flow rate was increased from 6 Figure 11. Left-hand panel: map of the events within sequences (colour circles) and the whole earthquake data set (grey dots), for the different depth intervals. The size and colour of the circles refer to the radius of the inferred rupture area. Middle panel: evolution of the earthquake radius as a function of time for repeating events within the considered depth interval (blue) and within the other depth intervals (grey). Right-hand panel: distribution of event size for the data shown in the middle panel.

Largely distributed event sizes and generation of new fractures
to 12 l s -1 and not the increase from 12 to 18 l s -1 as suggested by Cornet et al. (2007).

Observational limitations
Because our frequency range for the analysis of seismic event is limited between 80 Hz (this lower boundary is fixed due to noise contribution) and 500 Hz (attenuation being too severe at higher frequencies), it potentially may affect our rupture sizes estimation. These two limits indeed set a minimum radius estimation limited to 2.1 m while the upper boundary is set to 13.4 m. These two boundaries may control the shapes of the radius distributions. For example at the beginning of the stimulation we indeed observed a narrow distribution peaked in 3.5 m. The lower bound set at 2.1 m, may hide a wider distribution spanned over smaller sized events but that are inaccessible given our limited higher frequency range. We however note that if a majority of the events were actually having radius outside the range imposed by our limit, then we should have observed a concentration of event's radius at the two boundaries, which is not observed. We also note that our size estimate are in similar range as the one reported in Bourouis & Bernard (2007). Furthermore, as our study mostly focus on the variation of the event properties with time, the interpretation that we propose is not altered by such limitation although the amplitude of the variation might actually be larger than what is captured here. We however acknowledge that our interpretation is based on the assumption of a constant rupture velocity for all the different earthquakes. Changing adequately the rupture velocity for each event might balance the observed variation of event's sizes and associated stress drop. However this would requires significant variation of V r in order to maintain a constant stress drop for all events. A simple uniform change of the rupture velocity for all events will indeed cause a shift of all determined radius estimations but would leave our conclusion unchanged because we still will observe 2 populations of asperity with two different size distribution. However a distributed variation of rupture velocity in the different parts of the reservoir could potentially canceled the observed radius increase. There is however no argument for such a reduction of the rupture velocity (a 50 per cent reduction is required in order to balance the increase of radius) as all events are located in the same part of the reservoir (in granitic rocks).
We interpret here the change of event radius with time (and associated moment) as a consequence of a mechanical phenomenon and as a true feature of the seismicity. It is however possible to explain such an evolution simply considering a progressive lack of detection of small events. In this scenario, small events still exist in the same proportion as before but are just not detected because the detection capability of the network degrades over time. We discard this hypothesis because over the course of the injection, the latter seismic events are migrating to shallower depth, closer to the seismic sensors. It results of this shortening of the wave-path distance that the detection capability should actually increase over time because events will face a lower attenuation. Furthermore, explaining the change of the event's radius based on the network detection capability as observed would require that a jump in the detection occurs coincidentally with the wellhead pressure increases (notably when the flow rate reaches 12 l s -1 ), which is quite unlikely. We therefore favour the interpretation of radii variation through time as related to a mechanism acting on the events sources and not an observational bias.

C O N C L U S I O N
Based on accelerometer records of the 1993 injection experiment in the geothermal reservoir of Soultz-sous Forêts, we estimated source parameters of repeating earthquakes. We show that the rupture area of these events is dependent on the mechanism at the origin of the seismicity. We first evidence that at the asperity level, a deviation from the typical scaling law between moment and radius is observed suggesting that the stress drop on each asperity is significantly fluctuating during a repeating sequence. We also show that the source radius of repeating events is a relevant feature for identifying two populations of events: a first population of small radius events is evidenced close to the injection well, in the zone of the reservoir that experienced aseismic slip on pre-existing faults. As the fluid pressure is increased in the reservoir and the seismicity develops farther away from the injection zone, a second population of repeating events emerges with a wider range of radii and a different recurrence behaviour. We propose to relate this second population to the initiation of new large fractures that are fluid induced. Our results suggest that brittle and ductile mode of deformation promotes populations of repeating earthquakes that have distinct properties.

A C K N O W L E D G E M E N T S
We thank F. Cornet, F. Cappa, P. Dublanchet, Z. Duputel, E. Gaucher, A. Zang, A. Jupe, J. Kinscher, G. Dresen and L. Rivera for fruitful comments and suggestions. We thank F. Grigoli and S. Wei for very useful suggestions that helped to improve the manuscript. All date are available at CDGP (CDGP -Data Center for Deep Geothermal Energy) https://doi.org/10.25577/SSFS1993. This work has been done under the framework of the LABEX ANR-11-LABX-0050-G-EAU-THERMIE-PROFONDE and benefits from a state funding managed by the French National Research Agency (ANR) as part of the 'Investments for the Future' program. It has also been funded by the EGS Alsace Grant from ADEME.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online.
Table S1 Station coordinates after Gaucher (1998). Distances are expressed in meters from the injection well GPK1. Figure S1 Example of three Fourier spectra (blue line) fitted with a Brune's model (black dashed line). The 3 events (04/09/1993 11:47:59; 05/09/1993 09:12:07; 05/09/1993 09:39) belong to a same sequence of repeating earthquakes. The corner frequencies (f c ) obtained with this method are 375 ± 140, 420 ± 180 and 265 ± 90 Hz, for the events 1, 5 and 16 of the sequence, respectively. We apply a grid search to retrieve the values of f c and 0 , and estimate the uncertainties from the shape width of the objective function approximated as a normal distribution. Brune's spectra are also shown (light grey dotted lines) for f c ± df. The noise is estimated averaging the signals over 0.25-s-long windows preceding the events (yellow line). Figure S2 The spectral ratio analysis is applied to the events shown in Fig. S1. The spectra on the upper part of the Figs (a,b,c) are the individual amplitude spectra. At the bottom (d,e,f) we present the associated spectral ratio and the linear fitting of the log 10 (G(f). The corner frequencies differences, and in turn, the variations in radius, are derived with eq. (8) of the manuscript from the slope of log 10 (G(f). Figure S3 Comparison between the estimations of source properties for 19 repeating earthquakes of a given multiplet, derived from the Brune's model (in black) and from the inversion of the relative measurements (in red). On top, the results of the logarithm of Moments with the uncertainties are shown. At the bottom, the source dimensions. The uncertainties are notably reduced after the inversion of the relative measurements. In this particular sequence of repeating earthquakes: the properties are re-estimated from 1.05e8 ± 8.2e7 to 1.2e8 ± 9.9e6 Nm for the moments, and from 4.37 ± 1.62 to 3.87 ± 0.11 m for the source dimensions. The grey triangles highlight the results for the events presented in Figs S1 and S2.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.