Analysis of the Spatiotemporal Evolution of the Maurienne Swarm (French Alps) Based on Earthquake Clustering

Between August 2017 and March 2019, an intense seismic swarm was recorded in the Maurienne valley in the north of the Belledonne massif (Western French Alps). In order to study the spatiotemporal evolution of the Maurienne swarm, we created a high‐resolution catalog based on template matching, double‐difference relocation, and moment magnitudes. The catalog includes 71,064 events with a maximum moment magnitude of 3.5 and a magnitude of completeness of 0.7. The seismic activity is interpreted as the reactivation of an N80 strike‐slip fault system called Fond de France. Moreover, earthquake relocation reveals the presence of a shallower fault system with the same strike, but opposite dip direction and smaller size. The presence of two distinct fault systems may explain the observed variation of the b‐value with depth. The seismicity migrated asymmetrically in all directions during the course of about 15 months. Shorter migrations lasting 2–3 days are also observed. The different migration patterns suggest that the swarm is driven by several mechanisms, possibly pore‐pressure diffusion, aseismic slip, and earthquake interactions.

The present study deals with a seismic swarm that occurred between 2017 and 2018 in the Maurienne valley, which is located in the North Western French Alps. Before 2015, seismic activity in the study area was limited (Guéguen et al., 2021). As a result, the Maurienne swarm offered the opportunity to study for the first time the local tectonic setting through passive seismic methods. Guéguen et al. (2021) carried out a preliminary analysis on the Maurienne swarm. Starting from a catalog built using SeisComP3 (Weber et al., 2007), they relocated the earthquakes with a double-difference method, calculated focal mechanisms, and provided a preliminary description of the fault system reactivated during the swarm. In this paper, we expand the SeisComP3 catalog and the catalog of Guéguen et al. (2021) by applying template matching. We also relocate the new detections and estimate their moment magnitude. Finally, we use the new catalog to investigate the driving processes of the swarm. The study revolves around the concept of a cluster, which is here defined as a group of events with highly similar waveforms. We show that by classifying the recorded events into clusters, it is possible to include in the analysis events that would normally be excluded because of their low signal-to-noise ratio (SNR). This method is particularly effective in case studies where good SNR events do not provide enough information to properly interpret the spatiotemporal evolution of a seismic sequence and its triggering processes.

Geological Setting
The study area is located in the northern part of the Maurienne valley, which carves in the north-south direction the Belledonne massif ( Figure 1a). The Belledonne massif corresponds to one of the Paleozoic Crystalline Massifs (ECM) of the French western Alps (Figure 1a), which extends over more than 120 km in the N30 direction. This massif is locally higher than 3,000 m and is bounded to the west by the broad post-glacial Isère valley. The basement is made of Paleozoic rocks deformed and metamorphosed during the Variscan orogeny, unconformably covered by Mesozoic sediments. The Paleozoic lithologies are affected by strong Variscan deformation, such as major dextral strike-slip faults trending from N30 to N60 (Guillot et al., 2009). During the Alpine orogeny, the Belledonne ECM was first underthrusted beneath the internal zones in relation with the activation of the Penninic frontal thrust (Figure 1a) since 33 My (Simon-Labric et al., 2009). Then, the westward propagation of crustal shortening toward the European foreland led to the exhumation of the Belledonne ECM, whose internal deformation involves the reactivation of the main Variscan tectonic structures (Dumont et al., 2012). The Belledonne massif is divided into two major tectonic units, the so-called outer and inner units, from west to east. These two blocks are separated by a major late Variscan, N30-trending subvertical fault (Le Roux et al., 2010) named Belledonne middle fault (BMF) (Figure 1a). On both sides of this fault, the lithologies, ages, and tectono-metamorphic histories are different. The outer domain consists of a meta-sedimentary unit of Upper Ordovician age (Fréville et al., 2018) corresponding to micaschists (Série satinée) affected by low-grade polyphase metamorphism, whereas the inner domain consists of a tectonic stack of amphibolitized rocks of Cambrian to Lower Carboniferous age with Variscan granites (Guillot et al., 2009). These structures can also be crosscut by Variscan faults trending N60. Some N30-trending Variscan faults have been reactivated with dextral strike-slip displacement during Alpine deformation (Figure 1a). The Belledonne ECM is presently affected by recurrent active deformation, distributed along a seismic lineament following at depth the western edge of the massif (Thouvenot et al., 2003). This seismic lineament is characterized by low-magnitude earthquakes (M w < 3.5) located at shallow depths (less than 10 km) and distributed over more than 50 km. The focal mechanisms are compatible with the dextral strike-slip regime (Thouvenot et al., 2003). This activity is interpreted as the tectonic signature of the so-called Belledonne border fault (BBF in Figure 1a). However, this active structure has never been evidenced directly at the surface by geological or geomorphological observations. The seismic swarm of the Maurienne valley is located below the western flank of the Lauzière massif on the right bank of the Arc river ( Figure 1b). Most seismic events belong to an N60-trending cluster steeply dipping northward. The focal mechanisms of the largest events of the swarm show a dominant strike-slip activity associated with a normal component (Guéguen et al., 2021). The 3D geometry of the swarm coincides with the outcropping Fond de France Fault (FFF) system (Figure 1b), which is connected toward the southwest to the BMF.

Data Acquisition and Starting Catalog
We used data recorded by six broadband seismic stations from 1 March 2017 to 31 October 2019 ( Figure 1). Station A181A, part of the AlpArray network (Hetényi et al., 2018), was the only recording station before the beginning of the swarm in July 2017. A181A was composed of a Nanometrics Trillium-120QA sensor (period of 120 s) and of a Nanometrics Centaur digitizer sampling at 100 Hz. The remaining five stations were deployed by the Seismological Alpine service SISmalp from the University Grenoble Alpes, hosted by the Institute of Earth Science (ISTerre) (Guéguen et al., 2017). The five temporary stations were installed between 27 and 29 October 2017, soon after the occurrence of the largest event of the seismic sequence, and were equipped with GURALP CMG40 sensors (period of 30 s) and with Nanometrics Taurus digitizers sampling at 200 Hz.
The initial catalog was composed of 5,367 events extracted from the SISmalp catalog. The events are the result of both automatic algorithms provided by SeisComP3 (Weber et al., 2007) and manual revision. The initial locations are shown in Figures 2a and 2b.

Creation of the New Catalog
We applied a set of techniques to improve the original catalog in terms of number of detections, location accuracy, and magnitude estimation (see Figure 3 for a summary of how the new catalog was created).   Guéguen et al. (2021). The locations of the initial catalog (light gray points in a and b) do not allow to distinguish any clear fault structure. After template matching and double-difference relocation, (c) the epicenters align WWS-ENE and (d) the hypocenters, mainly located between 2 and 6 km depth, define a high-angle fault plane dipping north at about 70°. The events are divided into two groups based on the dip direction of the cluster they belong to. Most clusters (blue) have the same dip direction as one of the main fault planes, but some of them (orange) present the opposite dip direction and tend to be concentrated in the same region. Clusters with no clear planarity are shown in gray. Similar structures can be observed in the catalog of Guéguen et al. (2021) (dark gray points in a and b).

Template Matching
The first step was to increase the number of detected events. For this purpose, we employed template matching (Anstey, 1964;Gibbons & Ringdal, 2006;North, 1963), which commonly allows to increase by a factor of 10 the amount of earthquake detected by methods, such as STA/LTA (e.g., Peng & Zhao, 2009). Its effectiveness has been proved in a variety of seismological settings (e.g., Harris, 1991;Peng & Zhao, 2009;Shelly et al., 2007), including the study of seismic swarms (e.g., Shelly & Hill, 2011). Template matching consists in cross-correlating a set of reference events (called templates) with continuous recordings. When the correlation coefficient is above a predefined threshold, the time producing the largest value is saved as an indicator of a detection. To guarantee its good performance, the set of templates must ideally be composed of high-SNR (signal-to-noise ratio) events that are representative of the variety of waveforms generated in the area under investigation.
We applied template matching on a single station, A181A, the only one covering the entire duration of the swarm. We selected the templates starting from the events already present in the SISmalp catalog. We selected earthquakes with an SNR larger than 8, taking care to exclude overlapping events, which could otherwise lead to potential false detections. The SNR is estimated by calculating the ratio between the root mean square of the signal (containing P and S waves) and the noise. The signal corresponds to the samples inside a window starting 0.1 s before the P-wave arrival and ending 3 s after the P-wave arrival. The noise corresponds to the samples inside a window starting 1.1 s before the P-wave arrival and ending 0.1 s before the P-wave arrival.
To reduce the number of templates, we grouped similar events into clusters based on their similarity by using hierarchical clustering (Sokal, 1958). Finally, for each cluster, we selected the event to be used for template matching, which corresponded to the one producing the largest average correlation coefficient among all the other events part of the same cluster. The input for hierarchical clustering is a dissimilarity matrix, that is (1 -CC) values where CC is the correlation coefficient. To build the dissimilarity matrix, we filtered the selected waveforms between 3 and 40 Hz and we cross-correlated them between each other using all three components. We considered the minimum correlation coefficient between the three components to promote the formation of clusters of similar events.
In case of dominant P and S waves, the correlation coefficient can be large even if the rest of the waveforms can differ significantly. To lessen this effect, which would otherwise lead to false cluster assignments, we used waveform equalization before cross-correlating the events, similar to Herrmann et al. (2019). Waveform equalization is effective in reducing the largest amplitudes, thus lowering their influence, without changing the waveform shape ( Figure S1 in Supporting Information S1).
The dissimilarity matrix was then used as an input to link the events into clusters. We adopted the UPGMA (Unweighted Pair Group Method with Arithmetic mean) algorithm (Sokal, 1958). The algorithm runs iteratively, considering at the first iteration each event as a singleton cluster. At each iteration step, the nearest two clusters Figure 3. Summary of the main steps carried out to create the high-resolution catalog. The steps can be grouped into three categories based on what they helped to achieve: new detections, high-resolution locations, and moment magnitudes. The information contained in the catalog was used to study different aspects of the swarm, which are indicated by the red boxes. are combined into a higher-level cluster to form a binary tree. The distance (expressed as a dissimilarity value) between any two clusters is the mean distance between the elements of the two clusters. The formation of clusters is stopped when a specified maximum distance (dissimilarity threshold) is reached. This threshold controls the final number of clusters and we found that a value of 0.25 allowed to achieve a good balance between the number of clusters and the internal similarity of the clusters ( Figure S2 in Supporting Information S1).
In preparation for template matching, we sliced the templates into windows of 3.1 s (0.1 s preceding the P waves, followed by 3 s containing the actual waveform) filtered between 3 and 40 Hz. We selected this frequency band because the earthquakes recorded in the area have most of their energy inside this interval. The correlation coefficient threshold was set to 0.4. The coefficient threshold was defined after testing different values on short periods of time (days) and visually inspecting the ratio between the number of true and false detections. During template matching, we determined the global correlation coefficient for each sample by choosing the average value among the three components.
Waveform equalization does not offer any advantage if applied to continuous recordings. For this reason, during template matching, it is not possible to prevent two waveforms with strong P or S waves from generating inflated correlation coefficients that would cause them to be wrongly assigned to the same cluster. After template matching, we therefore applied waveform equalization to all waveforms and performed a cluster reassignment by cross-correlating each detection with each template. We took advantage of this step to also remove false detections. To this end, we cross-correlated the detections not only with the templates but also with a set of "anti-templates" corresponding to typical noise signals. Detections producing a higher similarity with the anti-template were considered as false events and removed. The anti-templates were selected among the detections with SNR > 8 and CC < 0.5. Indeed, if a detection has a high SNR but low similarity, we supposed that it was either noise or a real event that did not fit into any cluster. Since the impact of real events classified as anti-templates is negligible, for the sake of simplicity, we simply treated them as false detections instead of integrating them in the catalog. During the selection of the anti-templates, we also applied hierarchical clustering similar to what was done for the templates as a way to reduce their number.

Double-Difference Relocation
We used a double-difference location method (Waldhauser & Ellsworth, 2000) to locate the events detected with template matching. If the hypocentral separation between two earthquakes is small compared to the event-station distance and the scale length of the velocity heterogeneity, then waves traveling from the source region to a common station share almost the same path. In this case, the difference in travel times for two events observed at one station can be attributed to the spatial offset between the events with high accuracy. This method has been successfully applied in several studies of swarms (e.g., De Ruhl et al., 2016;Shelly et al., 2016;Yoshida & Hasegawa, 2018).
To perform a double-difference location, it is necessary to assume an initial location and origin time and to know the wave arrivals. Taking advantage of the fact that the events of a cluster should occur close to each other, we derived the initial location, the origin times, and the arrival times of the detections by assuming that each detection occurred at the same location as the reference event of the same cluster.
The reference events we used for the relocation did not correspond to the reference events (templates) we employed for template matching; only the clusters remained unchanged. The reference events (templates) used for template matching were chosen only considering the SNR at station A181A, and this does not guarantee that they are also suitable for relocation. Indeed, a good reference event for relocation must not only have a good SNR but also a good number of picks and an azimuthal gap < 180°.
The earthquakes part of the starting catalog were located with an absolute location method. In order to provide a more precise starting location for the detections, we first attempted to relocate the detections already part of the original catalog. For this step, we used as input differential times estimated with cross correlation and differential times estimated from catalog arrival times. Both P and S waves were considered, and only pairs with at least six observations were accepted.
The cross correlation of event pairs was carried out using windows of 0.4 s for P waves and of 0.5 s for S waves. The differential times were determined with subsample precision by oversampling the cross-correlation function.
During this operation, we only kept measurements with CC > 0.60. We used a 1D P velocity model taken from Potin (2016) The location errors were estimated with bootstrap (e.g., Efron, 1979;Mesimeri et al., 2018;Waldhauser & Ellsworth, 2000), a statistical resampling approach. The first step was to create a synthetic data set whose events have the same location and origin time as the relocated events. Knowing the location and origin time, we also estimated with a ray tracer the theoretical P-and S-wave arrival times of the same phases used to run the original relocation. The theoretical arrival times were then used to derive theoretical differential times. The residuals calculated during the original relocation process were divided into two sets because P and S residuals show different distributions ( Figure S3 in Supporting Information S1). Next, a residual value (randomly sampled from the set of residuals) is added to each differential time obtained from theoretical arrival times. The perturbed differential times are then used to relocate the events with the starting position corresponding to the location calculated during the original relocation. The operation was repeated 200 times, and the cumulative locations were employed to calculate the error ellipsoid, whose axis is enabled to derive the location errors in the three dimensions.
Among the events that could be relocated, we selected for each cluster the event with the greatest number of picks, which becomes the new reference event of that cluster. The locations of reference events derived from hypoDD were then used to determine the initial location of the detections. In particular, we assigned to each detection the same initial location as the reference event of the cluster it belongs to. If detections and reference events are supposed to have the same location, the travel times are also supposed to be equal. Since the detection time given by template matching is linked to the P-wave arrival time at A181A, it is therefore possible to derive the approximate origin time and subsequently the arrival times at all the stations for both P and S waves. This is true whenever the reference event has available picks. By using cross-correlation differential times, we attempted to relocate all the detections that occurred after the installation of the temporary network. We considered only pairs of earthquakes with initial locations separated by less than 1 km using the same time intervals and correlation threshold employed when relocating the initial catalog. Similar to what was done to relocate the reference events, the correlation threshold was set at 0.60, and only pairs with at least 6 observations and an azimuthal gap < 180° were taken into consideration. We point out that all the events that occurred when only one station was active (see Figure 4) cannot be relocated with hypoDD.

Moment Magnitudes
The last step for building the new catalog was estimating the moment magnitude of the detected events.
Moment magnitudes (M w ) were calculated using the equation (Hanks & Kanamori, 1979): The seismic moment M 0 was derived from the equation proposed by Brune (1970): where ρ is the density (2,650 kg/m 3 ), c is the velocity of P or S waves (5,200 m/s for P and P/1.73 for S waves), Ω 0 is the low-frequency amplitude of the displacement spectrum for each component (Z, N, and E), U Φθ is the average radiation pattern (0.73 for S waves and 0.52 for P waves), and R is the hypocentral distance. In order to include more events than those relocated with hypoDD, we calculated R using hypocenters affected by larger uncertainties. For each event, we selected the hypocenter with the best available resolution, following this order: hypocenter associated with clusters relocated with both cross correlation and catalog differential times (77% of the events), hypocenter associated with clusters relocated with just catalog differential times (19% of the events), and hypocenter associated with clusters reported in the original catalog (4% of the events). We derived Ω 0 by fitting the observed displacement spectra with the theoretical source model of Boatwright (1980): where Q is the frequency-independent quality factor, n is the high-frequency falloff rate, f c is the corner frequency, and γ is a constant with a value of 2. γ = 2 produces a sharper corner than the Brune model (Brune, 1970) with γ = 1, and for this reason, it is considered more suitable for shorter hypocentral distances (Abercrombie, 1995). Because the seismic moment is derived from the amplitude of the spectrum at low frequencies, it is not significantly affected by a potential frequency-dependent Q value (e.g., Aki, 1980). We therefore assume that Q is frequency independent.
Q, f c , and n, along with Ω 0 , are the so-called free parameters of the fitting. Q was allowed to vary between 10 and 1,000, f c between 1 and 40 Hz, and n between 1 and 10. We decided not to use a constant Q because we do not have a reliable value for the area under investigation. To facilitate the fitting, we smoothed the observed spectra with adaptive tapering (Prieto et al., 2009). Moreover, to obtain suitable starting values, we employed the differential evolution (Storn & Price, 1997).
We kept an event only if, for all channels, the minimum value of the signal amplitude below 12 Hz was larger than the maximum noise amplitude below 12 Hz. We chose this interval because it represents the most relevant frequency range for the fitting. Figure S4 in Supporting Information S1 shows some examples of spectra calculated from earthquakes of different sizes.
Spectral fitting using S waves provided more stable values compared to P waves. For this reason, we decided to use moment magnitudes obtained from S waves. We computed moment magnitudes for all stations and computed the average value. Only values computed from at least four stations were kept. In order to also consider the events that occurred when only station A181A was active, we fitted a linear relation between the moment magnitude calculated using only station A181A and the mean moment magnitude. We then applied the resulting correction factor to the magnitudes of the events recorded at station A181A ( Figure S5 in Supporting Information S1). Spectral fitting can be applied only if the events have a high enough SNR. As a result, we could estimate M 0 for only a limited part of the events (4,166 events). To calculate M w also for the smaller events, we again exploited the similarity that characterizes events belonging to the same cluster, similar to Herrmann et al. (2019). Indeed, if two or more events occur at a short distance from one another and their corner frequencies are greater than f max (the highest observable corner frequency due to attenuation), then their waveforms are scaled copies of each other (Deichmann, 2017). If this is the case, M 0 can be computed directly from the maximum displacement amplitude of the waveform (A max ) using the relation: The intercept C varies from cluster to cluster, but it can be directly derived by using the M w values calculated during spectral fitting. For this step, we only used A181A because most small events were hidden by noise at the other stations. Also, even if multiple stations would provide more stable values, the use of a single station is justified by the fact that the maximum amplitude values can be easily calculated and thus are affected by small uncertainties.
f max depends on the Q value of the area. In addition, it is important to notice that f c > f max is valid only for events smaller than a specific magnitude. Since we did not have precise information on the Q of the area, we could not precisely calculate f max , and we therefore had to adopt a conservative value (M w > 2.1) derived from the visual observation of the M w −A max fit or of the normalized spectra. As a result, earthquakes with M w > 2.1 were not employed to determine the intercept C.

Frequency-Magnitude Distribution Analysis
The global magnitude of completeness was estimated with the MBASS (Median-Based Analysis of the Segment Slope) method (Amorese, 2010), while the global b-value was estimated with the method of Tinti and Mulargia (1987). The errors for both magnitude and b-value were calculated with 2000 rounds of bootstrap (Efron, 1982).
We also calculated how the b-value varies with time and depth. For the time analysis, we focused on the period of highest seismic activity (September 2017-May 2018), using all the available magnitude values. For the spatial analysis, we only considered the events part of clusters whose reference events could be relocated with hypoDD.
To estimate the change of b-value with time, we used a rolling window containing 1,200 events, shifting it every 30 events. To avoid the effect of short-term variations of magnitude of completeness (Hainzl, 2016;Kagan & Schoenberg, 2001), which would cause an underestimation of the b-value, we adopted the beta-positive method (van der Elst, 2021), which uses differences of magnitudes. The main advantage of this method is that there is no need to calculate the magnitude of completeness.
To compute the variation of b-value with depth, we used a rolling window of 200 m, shifting it every 20 m. For each shift, we employed the maximum curvature method (Wiemer & Wyss, 2000) to calculate the magnitude of completeness and the method of Tinti and Mulargia (Tinti & Mulargia, 1987) to calculate the b-value. The maximum curvature method provides more stable estimates than the MBASS method, and it is generally well suited for short periods of time when the distributions are less rounded. On the other hand, the maximum curvature method tends to underestimate the magnitude of completeness; for this reason, we applied a correction of 0.2 as suggested by Woessner and Wiemer (2005).

Template Matching
Among the 5,367 events part of the SISmalp catalog, we selected 1,330 templates and used them for template matching. Template matching allowed to detect about 100,120 events, including false detections. By employing 411 anti-templates (noise signals) to remove false detections, the total number of events was reduced to 79,503. Compared to the SISmalp catalog, it is a 14-fold increase in the detection rate. The greatest improvement is observed when only one station was active. Moreover, template matching not only improved the detection of small events but also allowed to detect high-SNR events that were missed by SeisComP3. Most false detections consisted in high-frequency noise, probably of anthropic origin since most of them occurred during the day ( Figure S6 in Supporting Information S1). We also observed that false detections were mainly caused by templates with particularly strong P waves and relatively weak S waves ( Figure S7 in Supporting Information S1). Figure 4 shows the time evolution of the seismic activity linked to the swarm. The seismicity increased abruptly at the end of July 2017, followed by gradually more intense peaks of activity until November 2017, when the maximum rate observed during the swarm was reached (about 1,400 events/day). Next, the seismicity gradually decreased, always with isolated bursts, until reaching a stable rate in January 2019. In total, the period of increased seismicity rate lasted for about 17 months.

Event Relocation
We have relocated 2048 events out of 5,367 events of the starting catalog with the double-difference location method. The location errors, estimated with bootstrap, are 55, 125, and 95 m for latitude, longitude, and depth, respectively. The larger error in easting is mainly due to the network geometry, which is stretched along the north-south direction. These errors are smaller than the errors affecting the events relocated by Guéguen et al. (2021) (175, 138, and 164 m). Among the 2048 relocated events, it was possible to select a reference event for 808 out of the previously defined 1,330 clusters. We then employed these 808 reference events to also relocate the detections of template matching. In total, we relocated 12,424 events out of the 79,503 detected events (Figures 2c and 2d). This allowed to significantly improve spatial resolution compared to the initial catalog (Figures 2a and 2b), revealing the geometry of the reactivated fault systems.
The epicenters show a WE alignment, while most hypocenters define a faulting structure dipping north at about 70° between 2 and 6 km below sea level. The position of events matches with the geometry and position of the nearby Fond de France Fault system (Figure 1). The same structure can also be observed by simply relocating the original catalog as shown by Guéguen et al. (2021) (Figures 2a and 2b).
We applied the principal parameter method (Michelini & Bolt, 1986) to determine the geometry of individual clusters and found that they can be divided into two main groups based on their dip direction (Figures 2c and 2d). The biggest group of clusters has a dip angle and a dip direction similar to the global one (355°/68°), while a smaller group of clusters has an almost opposite dip direction and a lower dip angle (170°/60°). These clusters are concentrated in the southern part of the reactivated region and their events are among the shallowest of the seismic sequence. Only clusters with good planarity (>0.8) were considered. The planarity is a value less than 1, which is derived from the spatial covariance matrix.

Migration of Seismicity
One of the main advantages of template matching is the improvement of the catalog completeness, especially before the installation of the temporary network. However, double-difference relocation prevented us from analyzing the initial phase of the swarm because it can only be applied if multiple stations are available. Moreover, double-difference relocation lowered the spatiotemporal resolution because it requires events with a good SNR. In order to take full advantage of the events detected with template matching, we considered all the events that could not be relocated but that were part of a cluster in which at least the reference event could be relocated. We then assigned to each of these events a random location within a radius of 70 m from the centroid of the corresponding cluster. This radius corresponds to two times the mean standard deviation of the distance between the relocated events (Figures 2c and 2d) and the centroid of their corresponding cluster. By doing this, we assumed that the spread observed after the relocation is representative of the spread within a cluster. This method allowed to add 48,237 more events to the 12,424 events that were directly relocated with the double-difference method.
This means that we could analyze the spatiotemporal variations of seismicity using 60,661 of the 79,503 events detected with template matching. Figure 5 shows that the seismic activity started at about 3.8 km depth and expanded in all directions, but mainly in depth, over the course of several months. We point out that the events exclusively relocated with hypoDD, along with the events of the original catalog (Figures 2c and 2d), do not allow to observe this migration pattern because they lack the spatiotemporal resolution. Figure 6 offers a more explicit representation of the time dependence of the locations. Figure 6a shows the variation of the distance from the center of the swarm, which corresponds to the mean location calculated from the hypocenters that occurred at the start of the main seismic sequence (August 2017). The seismicity migrated for about 15 months within 2 km from the center of the swarm. The expansion was discontinuous in time with bursts of seismicity alternating with less active periods, and it cannot be modeled by a simple diffusion law. However, observing shorter periods of time, it is possible to identify migrations that better agree with pore-pressure diffusion and that are faster than the global migration. In the example shown in Figure 6b  To better understand the process controlling the expansion of the seismicity, we also used the natural time (Rundle et al., 2018). This means that the events are plotted based on their occurrence order (Figure 6c), thus eliminating the influence of time. As suggested by Fischer and Hainzl (2021), a discontinuous or nonlinear spreading in the coordinate-time plot and a continuous spreading in a coordinate-event-index plot are evidence of a self-controlled migration, where a rupture promotes the creation of new ruptures. In contrast, a continuous spreading in both coordinate-time and coordinate-event-index plots, as well as a nonlinear spreading in the coordinate-event-index plot, suggests that the migration is externally controlled by an aseismic process, such as pore-pressure diffusion and aseismic slip. In the case of the Maurienne swarm, the migration has a more linear pattern in the distance-event-index-plot (Figure 6c) than in the distance-time plot shown in Figure 6a.

Frequency-Magnitude Distribution
The spectral fitting method described in Section 4.3 allowed to calculate the moment magnitude of 4,166 events. We derived the magnitude-amplitude relation for 1,225 clusters and 71,064 events. Figure 7a shows the magnitude distribution of the entire recorded period, along with the b-value, magnitude of completeness, and associated uncertainties. The estimated b-value is 1.17 and the magnitude of completeness is 0.7 (Figure 7a). Notably, the Gutenberg-Richter law does not provide a good fit for the whole observed distribution since it deviates from linearity by about M w > 1.2. We then attempted to fit the distribution with two different Gutenberg-Richter laws ( Figure 7b) and with the model of Shapiro and Dinske (2013). The model of Shapiro and Dinske (2013) explains the underrepresentation of large-magnitude events by assuming that rupture surfaces and activated volumes are finite, and that pore-pressure diffusion is the main triggering process. To estimate the parameters of the models, we used maximum likelihood estimation. By comparing AIC (Akaike information criterion) values (Akaike, 1974), we found that the double Gutenberg-Richter law (AIC = 1,678) provides a considerably better fit than the pure Gutenberg-Richter law (AIC = 1,767) and the model of Shapiro and Dinske (2013) (AIC = 1,705). The resulting two lines have a slope (b-value) of 1.08 and 1.39 and meet at a breaking point at M w = 1.2.
The FMD obtained from the original catalog ( Figure 7a) shows a similar b-value and a similar deviation from linearity. It should be noted that in the original catalog, magnitudes are expressed in local magnitude (M L ). By directly comparing M w and M L , we observe that the two parameters are linearly related, but several outliers can be identified ( Figure S8 in Supporting Information S1). We manually recomputed the local magnitude of 75 outliers and we observe that especially when an event is small, the automatic algorithm of SeisComP3 tends to select the wrong time window for the calculation of the peak-to-peak amplitude. If the correct time window and peak-to-peak amplitude are chosen, the local magnitudes generally have a lower value that better matches our value of M w . These observations suggest that the moment magnitudes we estimated are more reliable than the local magnitudes of the original catalog.
Temporal variations (both positive and negative) of the b-value are often observed during seismic swarms (e.g., Jenatton et al., 2007;Passarelli et al., 2015;Shelly et al., 2016). Even though we observe an apparent decrease in b-value during periods of high activity, we think that the variability of the b-value observed during the swarm, together with the uncertainty affecting the estimations, does not indicate any clear correlation between b-value and seismicity rate (see Figure S9 in Supporting Information S1). Significant changes in b-value are instead observed when looking at how the b-value changes with depth ( Figure 8). The b-value is found to decrease from approximately 1.2-1.0 between 3 and 5.2 km depth. The largest values thus appear to correspond to the group of clusters whose geometry differs from the general one. See Figure S10 in Supporting Information S1 for a comparison between the FMD of the events part of the clusters dipping southward (b = 1.29 ± 0.08) and the FMD of the events part of the clusters dipping northward (b = 1.16 ± 0.03). The template-matching catalog allows a much more detailed analysis of the spatiotemporal variations of b-value. When using the original catalog, the uncertainty affecting the b-values is much larger because the number of events is significantly smaller than in our catalog.

Discussion
The relocated earthquakes highlight a high-angle fault system dipping northward that agrees with the position and orientation of the Fond de France Fault strike-slip fault system (Figure 1b). The focal mechanisms of the major earthquakes (Guéguen et al., 2021) also indicate the presence of strike-slip faulting during the swarm. Examining clusters of similar events reveals that there is a group of shallow clusters that dip southward instead of northward. This suggests that during the swarm, fault systems with different geometries were reactivated. Previous studies (Guéguen et al., 2021) do not report any focal mechanism for the shallow events, and for this reason, we cannot derive any additional information on their faulting style. Geological studies show the presence of 120°N-striking normal faults, which do not match with the strike of the clusters, but they indicate that in the region different sets of faults are possible, and the possibility of fault geometries yet to be discovered cannot be excluded. While deeper clusters indicate the presence of a single large fault, shallower clusters appear to be part of several smaller fault systems. The observation of faults with different dimensions may explain the change of b-value with depth. Indeed, large faults have more potential to generate larger magnitude events compared to small faults (e.g., Shelly et al., 2016). As a consequence, primary faults would lead to lower b-values than secondary faults in the same area. As faults with different structures likely require different stress states to be reactivated, our interpretation is in agreement with the theory that states that the b-value is stress dependent (Scholz, 1968(Scholz, , 2015. A stress dependent b-value could also mean that the decrease in b-value with depth may be the result of a change in  Figure 2d are also plotted. The shaded area represents the 2σ uncertainty affecting the b-value. confining pressure due to lithospheric loading (Mori & Abercrombie, 1997;Spada et al., 2013). However, we do not favor this hypothesis since a depth variation of 2 km is generally not enough to produce a change from 1.2 to 1.00 (e.g., Spada et al., 2013).
The identification of short-time migrations following a diffusion law (see Figure 6b) suggests that fluids may play an important role in the reactivation of the fault system as shown for many swarms worldwide (e.g., Hainzl et al., 2012;Ross et al., 2020;Shelly et al., 2016;Yoshida & Hasegawa, 2018). Fluid pressure on preexisting faults would reduce the effective normal stress (Hubbert & Rubey, 1959), triggering earthquakes along those faults already near failure. On the other hand, the brief migration of fluids cannot alone explain the slower, discontinuous, and multidirectional migration pattern observed during most of the recording periods. For this reason, we assume that additional driving processes must be taken into account. The generally linear spread of seismicity observed in the coordinate-event index plot (Figure 6c) indicates that a self-driving expansion may be possible (Fischer & Hainzl, 2021). A self-driving expansion means that ruptures promote fluid migration into the reactivated fractures by dilatancy, which in turn facilitates the creation of adjacent fractures. When looking at shorter periods of time, the linearity is less defined, which may hint at transient changes of the dominant process controlling the expansion or to the presence of multiple and independent rupture fronts.
Fluid-triggered seismic swarms can present an upward migration (Hainzl et al., 2012;Ross et al., 2020), a downward migration (Montgomery-Brown et al., 2019), or a combination of both (Hauksson et al., 2019;Li et al., 2021;Shelly et al., 2013Shelly et al., , 2016Yoshida & Hasegawa, 2018). Upward migrations are generally attributed to rising crustal fluids following the hydraulic gradient (Ross et al., 2020). Dominant downward migrations may be caused by groundwater recharge after snowmelt or rainfalls (Montgomery-Brown et al., 2019). Multidirectional migrations may be explained by assuming the presence of fluids whose pressure is greater than hydrostatic (Shelly et al., 2016). The propagation direction may also be due to local heterogeneities affecting permeability (e.g., Sibson, 1996). Figure 6 indicates that the direction of propagation of the Maurienne swarm belongs to the latter case. We therefore propose that the Maurienne swarm was caused by a fluid source located at about 3.8 km depth injecting fluids at supra-hydrostatic pressure.
Aseismic slip is another aseismic process able to trigger seismicity, but there is no compelling evidence indicating that it is a major driving process of the Maurienne swarm. Seismicity driven by aseismic slip usually occurs within a mostly creeping fault or at the edge of a creeping fault (e.g., Linde et al., 1996;Ozawa et al., 2003), but we think both situations are not consistent with the observed seismicity. Indeed, if the seismicity expanded in multiple directions predominantly along a well-defined fault plane ( Figure 5), it is unlikely that it occurred at the edge of a creeping fault or in a mostly creeping fault. An additional piece of evidence against aseismic slip is the lack of clear repeaters, that is, events with high similarity occurring at regular time intervals. However, aseismic slip can also be triggered by fluid flow  and since fluids appear to be necessary to explain the migration, the presence of aseismic slip cannot be completely excluded. More generally, the difficulty in identifying a clear migration pattern and a dominant driving mechanism may be related to the complex tectonic environment of the area. Indeed, both field studies and our results suggest that fault systems with different geometries coexist in the area. The presence of different fault systems may promote migration in preferential directions, thus leading to an asymmetric spreading and different migration velocities.
The observed FMD is not well-described by the GR law and shows that there is an apparent underrepresentation of large events. The observed FMD is best fitted by a double GR law with a breaking point at M w = 1.2 ( Figure 7b). The worse fit provided by models such as the one proposed by Shapiro and Dinske (2013) may derive from the underlying assumption that earthquakes are triggered by one dominant process. Indeed, the Maurienne swarm is likely the result of the interaction between several triggering processes.
The deviation from linearity cannot be explained by the observed decrease in b-value with depth. Indeed, assuming that the FMD is the sum of two GR laws with different b-values, and fitting this resulting global FMD with a simple GR law, the fit would instead underestimate the number of large events. The deviation from the GR law is a known issue that can affect high-resolution catalogs based on template matching (Herrmann & Marzocchi, 2020). The precise cause that leads to the observed deviation is not well understood, and different explanations that may also apply to our case have been proposed. A first explanation is the failure to detect overlapping events because larger events can mask smaller events in their coda waves. This effect is more easily observable during periods of high seismicity rate. Another cause could be selectiveness, that is, the dependence of the results on events present in the original catalog. The exclusion of a number of clusters during the course of the processing further increases selectiveness. However, this effect should be limited since the overall number of clusters remained large and since we managed to estimate M w of about 90% of the detected events. An additional cause that could contribute to the nonlinearity may be the fact that low-magnitude events are less similar to the templates, resulting in an incorrect magnitude estimation. Finally, the low number of large-magnitude events could be explained by the limited size of the activated volume (Shapiro et al., 2007).
The spatiotemporal and the frequency-magnitude distribution analysis of the Maurienne swarm were only possible by taking advantage of the concept of cluster. By employing clusters, we removed the limitations imposed by the relocation process and by the absence of stations during the first part of the swarm. As a result, we were able to consider more events than normally possible and to include in our study the whole period during which the swarm was active.

Conclusions
By employing clustering-based methods, we obtained a high-resolution catalog of the Maurienne swarm. By exploiting the concept of clusters, we were also able to derive the location and moment magnitudes of events that normally would have been excluded either because their SNR was too low or because they occurred when only one station was active. The set of methods we described could also be applied to other local seismicity studies in which the analysis of spatiotemporal variations of seismicity is hindered by the fact that only a small part of the detected events can be located with high accuracy.
Our analysis shows that the Maurienne swarm occurred along a WE-striking strike-slip fault called Fond de France, which progressively ruptured during the course of about 15 months. The reactivation of the Fond de France fault is also associated with shallow secondary fault systems with the dip direction opposite to the main fault system. The presence of secondary shallow faults could also explain the decrease in b-value with depth if the fault size and geometry control the observed b-values. The observed magnitude distribution shows an apparent underestimation of the number of small events and deviates from a pure Gutenberg-Richter law for magnitudes M > 1.2. This deviation hinders the estimation of a representative b-value with important consequences for earthquake hazard assessment.
The seismicity gradually migrated over the course of about 15 months. Overall, the migration appears to be driven by a self-controlled rupture process based on the interaction of seismic slip and fluid flow. However, we also identified short periods of time lasting 2-3 days during which the migration is better explained by pore-pressure diffusion. The Maurienne swarm thus appears to be the result of the combination of different driving processes.