On the trail of fluids in the northernmost intracontinental earthquake swarm areas of the Leipzig-Regensburg fault zone, Germany

The Leipzig-Regensburg fault zone is documented as a band of seismic activity extending northwards from the earthquake swarm region NW-Bohemia/Vogtland at the Czech-German border area and is intersected by several Hercynian fault zones. Along the fault zone, there are several earthquake swarm areas, the northernmost of which are Schöneck and Werdau. In this study, we investigate the presumably fluid-induced earthquake swarm activity of the Schöneck and Werdau area. For this purpose, we apply two methods: local earthquake tomography and receiver functions to identify the structural composition of the crust, the areas affected by fluids and the origin of the fluids. We detected potential fluid paths characterised by high Vp/Vs ratios and granite intrusions nearby the swarms characterised by low Vp/Vs anomalies. Receiver function analysis yields the Moho at 25 to 33 km depth and two seismic discontinuities at 55 km and 68 km depth.


Introduction
Earthquake swarms are sequences of seismic earthquakes that occur within a local area without a discernible main shock and are therefore distinct from normal earthquake sequences that include a main earthquake and aftershocks. Earthquake swarms are common in volcanic and hydrothermal areas (e.g. Zaliapin and Ben-Zion 2013;Shelly et al. 2016;White et al. 2019;Bachura et al. 2020) or in relation to the injection of industrial fluids (Horton 2012). Magma and fluid movement are suspected of causing seismicity by decreasing the effective normal stress and thus reducing resistance to slip (Hubbert and Rubey 1959). In contrast, intracontinental tectonic earthquake swarms without volcanic or geothermal activity are rare. In this study, we investigate the source of the northernmost earthquake swarm areas of the Leipzig-Regensburg fault zone (LRZ) and the reason for the absence of earthquake swarms further north. For this purpose, we use two different techniques to investigate earthquake swarms. Local earthquake tomography (LET) is used to image the crustal volume; receiver functions (RFs) are used to image the crust and uppermost mantle. Since the Vp/Vs is very sensitive to fluids (O'Connel and Budiansky 1974;Schmeling 1985;Mavko and Mukerji 1995;Takei 2002), the LET can be used to identify areas affected by fluids, such as the earthquake swarm areas, while RFs provide information about the source of the fluids. Previous attenuation studies along the LRZ have shown differences in attenuation behaviour between the central/northern (van Laaten et al. 2022a) and southern parts of the LRZ (Gaebler et al. 2015;Eulenfeld et al. 2021). The fluid-sensitive intrinsic attenuation in the southern part is larger, by a factor of 2 to 8, than in the central/northern part depending on the frequency, indicating a certain level of heterogeneity along the fault zone (van Laaten et al. 2022a).

Earthquake swarms in NW-Bohemia/Vogtland
The earthquake swarms occur in several known clusters along the LRZ (Korn et al. 2008;Fischer et al. 2014). Some of these have been well studied (Horálek and Fischer 2008;Fischer et al. 2014), such as the earthquake swarms in the NW-Bohemia/Vogtland region ( Fig. 1). Various observations indicate that uprising mantle-derived magmatic fluids trigger the earthquake swarms of NW-Bohemia/Vogtland (Špičak and Horálek 2001;Weise et al. 2001;Hofmann et al. 2003;Bräuer et al. 2009;Mousavi et al. 2015;Hrubcová et al. 2017;Vlček et al. 2022) coming from an updomed Moho below NW-Bohemia/Vogtland (Geissler et al. 2005;Hrubcová et al. 2017). A low-velocity zone in the upper mantle around 65 km in depth, interpreted as partial melts or asthenospheric updoming, may be the initial source of these ascending fluids (Heuer et al. 2006). Another connection between fluids and earthquake swarm activity is the spatial distribution of epicentres, which roughly overlaps with gas escapes (Fischer et al. 2014). In contrast to the well-studied NW-Bohemia/Vogtland region, the northern earthquake swarm areas of the LRZ are not yet well-studied. The region Werdau is known as the northernmost earthquake swarm area in the Saxothuringian Seismotectonic Province with the peculiarity of an increased intracontinental seismicity (Neunhöfer 1998). The main differences between the northern and southern earthquake swarm areas are the frequency of occurrence and the focal depths. The earthquakes in NW-Bohemia/Vogtland are much more frequent and shallower, with depths ranging from 6 to 10 km, compared to the northern earthquake swarm areas, which have depths ranging from 10 to 16 km.

Schöneck and Werdau earthquake swarm areas
The earthquake swarm areas Schöneck and Werdau ( Fig. 1), named after nearby cities, are located in Germany. Most earthquakes in the Werdau area occur at the intersection of the fault zones ( Fig. 1). Outside of the intersection, the seismic activity is largely prevented by aseismic creep along the GJZ (Ellenberg 1992), leading to an accumulation of energy at the intersection (Hemmann et al. 2003). As a result, the released energy is higher compared to other earthquakes in NW-Bohemia/Vogtland, leading to an increased seismic hazard. The largest historical earthquake in the region in the past 700 years occurred near the GJZ/LRZ intersection on March 6, 1872, with an estimated moment magnitude of M w = 5.2, and was felt at a distance of several hundred kilometres (von Seebach 1873).
From the period 1903 to 1999, 73 swarms from the LRZ region were analysed by Neunhöfer and Hemmann (2005). The b-values range from 0.5 to 1.5, with the intensive swarms of 1908, 1962, 1985/86, 1997 and 2000 being very close to 1. This value agrees with the observation of Ibs-von Seht et al. (2008), who determined a b-value of 0.8 for intracontinental earthquake swarms in non-volcanic areas. Recent swarm activities were observed in , 2010, 2012/13, and 2017for Schöneck and 1997/98, 2006 for Werdau, with the 2016 swarm taking place northeast of Werdau rather than south of Werdau. The Schöneck earthquake swarm area can be divided into three smaller clusters ( Fig. 1), the western cluster with focal depths between 13 and 17 km, the central cluster with focal depths between 8 and 13 km and the eastern cluster with focal depths between 8 and 13 km and occasionally depths Major fault zones are shown as black lines, such as the Leipzig-Regensburg fault zone (LRZ), the Gera-Jáchymov fault zone (GJZ) and the Bergen-Klingenthal-Chodov fault zone (BKZ). The tectonic earthquakes are shown as red circles, size representing their magnitude, while the earthquake swarms of the Schöneck area are shown in blue and the earthquake swarms of Werdau as green circles. Degassing locations (Mofette/Spring) are marked as yellow and the diatremes of Weißenbrunn and Ebersbrunn as orange and brown diamonds, respectively. In addition, the focal mechanisms of some earthquakes along the LRZ are shown. The areas of Vogtland (blue outline) and NW-Bohemia (pink outline) as well as the locations of the cities of Schöneck and Werdau are indicated. Earthquake locations after Thüringer Seismologisches Netz (2021). Focal mechanisms were obtained with Grond . Faults after LfULG (2018) and degassing locations after Geissler et al. (2005) between 16 and 20 km. The focal depth of the Werdau earthquake swarms is ∼16 km, except for the 2016 swarm northeast of Werdau, which had a focal depth of ∼12 km. The swarm duration of Werdau and Schöneck was about one to two months, with approximately 40 to 70 located events. The maximum magnitudes of the swarms were in the range of 1.1 to 1.3 for Schöneck and 1.7 to 2.8 for Werdau, respectively. In contrast to the NW-Bohemia/Vogtland region, CO 2 degassing near Schöneck and Werdau is not observed (Fischer et al. 2014;Heinicke et al. 2019). Since the Schöneck and Werdau earthquakes show a clear swarm character (Hemmann et al. 2003;Korn et al. 2008;Fischer et al. 2014), fluids might play an important role in the swarm generation. The origin of these fluids is unknown since there is no detected escape at the surface and, therefore, the fluids cannot be analyzed. Hemmann et al. (2003), who investigated the waveform and relocation of the Werdau earthquake swarm, suggested that the fluids rise along the intersection of LRZ and GJZ within the crust. However, no hydro-chemical or younger volcanic activities like those known for the NW-Bohemia/Vogtland region have been observed for the region of Werdau (Hemmann et al. 2003). Only two smaller diatremes (Schmidt et al. 2013) near Werdau indicate past upper mantle fluid activity. Schüller et al. (2012) performed (U-Th)/He age determination test measurements on apatite and determined an age of ∼113-96 Ma for one of the two diatremes, while the exact age of the other diatreme is still unknown. According to Alexowsky et al. (2007) the second diatrem is assumed to be Tertiary or Cretaceous in age.

Geological setting
The Vogtland area is part of the Saxothuringian Zone of the Variscan orogen (Pietzsch 1963). The area formed in several periods (Cambrian to Permian) and consists primarily of metamorphic rocks. It is bordered to the east by intrusions that were emplaced during the Variscan orogeny (Förster et al. 1999). Within the survey area, several Hercynian fault zones (Pietzsch 1963) also strike NW-SE, such as the GJZ located in the northern part and several parallel faults in the south, such as the Bergen-Klingenthal-Chodov fault zone (BKZ). Furthermore, the LRZ (Grünthal et al. 1985) crosses the survey area and is a N-S striking fault zone between the German cities of Leipzig and Regensburg with a seismically active middle part of about 150 km in length and 40 km in width between the Cheb basin in the Czech Republic and Leipzig (Kämpf et al. 1991;Bankwitz et al. 2003). The LRZ cannot be observed with geological methods as no corresponding geological displacements are visible in the investigated area. Older studies based on photolineations of satellite images claimed the observation of the LRZ at the surface (Bankwitz et al. 1979;Krull 1984;Grünthal et al. 1985;Bankwitz et al. 2003;Pohl et al. 2006), but these results were questioned by recent studies using modern satellite-based data (Dahm et al. 2018;Grünthal et al. 2018). Therefore, the existence of the LRZ is mainly proven by the observation of an N-S trending seismicity (Fig. 1). A bundle of short N-S trending geophysical potential field structures provides further evidence for the LRZ (Becker et al. 1989).

Data
We used recordings from 2001 to 2018 at several seismic stations of the BW, CZ, GR, SX, SXNET and TH networks (Fig. 2). For the inversion, we considered earthquakes from the TH catalogue (Thüringer Seismologisches Netz 2021) between latitudes of 50.35 • N and 51.08 • N and longitudes of 12.0 • E and 12.75 • E. We divided the region into several cubes of 2 km in side length to ensure a balanced distribution of identical sources. We calculated the mean magnitude of all events in each cell and considered only events for the inversion whose magnitude was greater than the mean magnitude. In addition, we used several criteria to filter the events further. Each event must fulfil an azimuthal gap criterion of 160 • , the magnitude of the event was required to be greater than −0.7, a minimum of five arrival times each for P-and S-waves were necessary for an event to be used for the inversion, the depth of the earthquake could not be greater than twice the distance to the nearest station and the maximum epicentral distance must be 110 km. After applying the criteria, the number of events decreased from about 2013 to 419 (Fig. 2).
All events were manually re-picked using Seis-micHandler (Stammler 1993). Based on the scheme proposed by Evans et al. (1994), we used arrival time weightings to improve the quality of the data. Regardless of the wave type, each pick is given a quality class. These classes range from 0 (highest quality) to 2 (lowest quality) to account for the picking accuracy of 0.025, 0.05 and 0.1 s. For P-waves, 77, 22 and 1 per cent were in classes 0, 1 and 2, respectively; for S-waves, the proportions were 68, 30 and 2 per cent. The average uncertainty is 27 ms for P-waves and 30 ms for S-waves. This results in 4406 P-phase and 3910 S-phase picks for the

1-D inversion
VELEST (Kissling et al. 1995) was used to invert for a minimum 1-D model. In order to obtain a reliable starting model and relocate the earthquakes for the 3-D inversion, it is necessary to first invert in 1-D to tackle the non-uniqueness of the non-linearized inverse problem's solution. The inversion depends on the starting model and initial hypocenters, so a trial-and-error procedure is necessary (Kissling et al. 1994). We used a wide range of starting models, beginning with a random velocity value for the surface and a random velocity gradient, resulting in several realistic and unrealistic models, as suggested by Boncio et al. (2004). A total of 1000 randomly generated starting models were produced, with a layer thickness of 1 km and a depth extent to the approximate Moho depth of 30 km. To sufficiently relocate the earthquakes and better narrow the solution space of the inversion, the origin time and location of each event varied realistically and unrealistically with randomly Gaussian distributed error in the range of ±10 km for latitude and longitude, ±0.1 s for arrival time and ±5 km for depth for each inversion. We used station SCHF as a reference station for VELEST due to its central location within the survey area and the high number of event recordings. Based on the Wadati diagram, the Vp/Vs ratio was fixed at 1.71, which is a comparable value to the adjacent region of NW-Bohemia (Málek et al. 2005;Mousavi et al. 2015). The final 1-D velocity model with the lowest RMS residual is displayed as a red line in Fig. 3. The inversion reduced the RMS error from 2.6 s in the starting model to 83 ms in the final model, which is still high compared to the average picking uncertainty. Most of the inverted models are similar to the best-fit model, but there are significant differences in the velocity of the near-surface layers. These layers are difficult to resolve with VELEST because of significant deviations from 1D layering near the surface. In addition to the accumulation of models around the best-fit model, there is another area of higher velocities where inverted models accumulate.
These models result from unrealistically high-velocity starting models (Vp at surface ≥8 km/s), which run into a local minimum during the inversion process. However, the RMS of these models with more than 160 ms is significantly higher than the best-fit model's RMS.
In addition, VELEST calculates the station corrections, indicating lateral variations of data fitting within the minimum 1-D velocity model (Kissling et al. 1995). Large-scale velocity variations and nearsurface heterogeneity cause these variations. The station corrections sign represents the velocity variation compared to the 1-D velocity model at the reference station. Positive values indicate that the true velocities are lower compared to the 1-D velocity model and negative values indicate that the true velocity is higher. Figure 4 shows the calculated station corrections for the stations used in this study relative to the reference station SCHF ranging from −0.50 to +0.48 s. The Pwave corrections show similar values at stations with similar geology, with stations close to the reference station showing lower values than stations farther away. Relative positive station corrections of the P-wave can be observed in the north-west of SCHF, which can be explained by the sediments in that area. On the other VELEST does not provide information about the absolute uncertainty of event locations. Therefore, we used HYPOINVERS (Klein 2002) to relocate 183 of the 419 events with at least 15 P-and 10 S-picks, using the new velocity model. The mean location uncertainty was 0.4 km and 1 km for latitude and longitude and 3 km for depth. The uncertainty in longitude was about twice that of latitude and can be explained by the lack of stations in the east and northeast of the survey area. A detailed presentation of the hypocentre change through relocalization can be found in the Appendix (Fig. 14).

3-D inversion
After inverting for the minimum 1-D velocity model, SIMULPS14 (Thurber 1983(Thurber , 1993Eberhart-Phillips 1993) was used to invert simultaneously for hypocentre locations and 3-D velocity structure. The forward modelling uses the pseudo-bending method to calculate the ray paths. SIMULPS14 uses the least-squares algorithm to update the velocity model and hypocentre locations. The inversion of the velocity model can be performed individually or simultaneously for Vp and Vp/Vs.The Vs model generally tends to have a lower resolution than the Vp model, due to the more difficult determination of the S-wave arrival time. Here we directly invert for the ratio of Vp/Vs since this value is sensitive to fluids. According to Thurber (1993) and Thurber and Eberhart-Phillips (1999) the direct inversion for Vp/Vs is more stable than the inversion for Vs and a subsequent computation of Vp/Vs.
The lowest RMS 1-D velocity model was used as the starting model and the inversion was performed first for the Vp model and then jointly with Vp/Vs. As suggested by Evans et al. (1994), the inversion was initially carried out with a large horizontal grid spacing of 30 km and then successively reduced. We used up to 10 iterations before the grid spacing was reduced to a final grid spacing of 5 km in the investigation area's centre. The outer region has a greater spacing to address the heterogeneous ray coverage. Throughout the inversion process, the vertical nodes are fixed at a spacing of 2 km. The damping values of the inversion were determined by a L-curve analysis (Hansen 1992). Damping of 15 for Vp and 25 for Vp/Vs (Fig. 5) was used in the inversion to provide a balance between data misfit reduction and a low model variance. The damping value of Vp/Vs is slightly higher compared to the Lcurve result (around 6) to achieve a smoother result.

Resolution tests
A quantitative analysis of the quality is the resolution matrix. The matrix's rows describe the dependency of the solution for one model parameter on all the other model parameters. The diagonal elements of the resolution matrix (RDE) indicate how the corresponding model parameter's solution depends on any other parameter, especially on its neighbours (Evans et al. 1994). The RDE ranges between 0 and 1, with 0 indicating no resolution and 1 indicating perfect resolution. However, high damping values and larger numbers of model parameters reduce the RDE (Eberhart-Phillips and Reyners 1997). The RDE values are calculated within SIMULPS14 and are shown in Fig. 6. Good resolution, defined as RDE values greater than 0.15 for Vp and 0.3 for Vp/Vs, is achieved in the central and southern areas of the Vp model, where a large number of earthquakes are located. With increasing depth, parts of the northern model are also well resolved due to the deeper earthquakes in this area. In the case of Vp/Vs, the resolved area is smaller and only includes the southern part of the model. Additional figures of the RDE can be found in the Appendix (Fig. 15).
A second evaluation of the resolution can be obtained by a checkerboard test (Hearn and Ni 1994;Zelt 1998). For the checkerboard test, a defined pattern with a specific size and perturbation is applied to a velocity model. Forward modelling is performed with the perturbed velocity model, generating synthetic travel times with the same source-receiver distribution used in the real data set's inversion. According to the pick accuracy, we added Gaussian noise to the data, and the synthetic travel times are inverted using the same parameterization and starting model used for the real data inversion. The checkerboard test is successful if the inversion can reconstruct the checkerboard pattern and the perturbation. During the inversion process for the synthetic test, the hypocenters were fixed to get the same ray distribution. This may enable a slightly easier recovery of the original model (Mousavi et al. 2015), but due to the good relocalization of the earthquakes, this will not be of significant importance. We tested a 5 per cent perturbation with different cell sizes. The smallest pattern that could be reconstructed, and thus the best result of the checkerboard test, is shown in Fig. 7, with a cell size of 15x15x5 km. The slices from 7 to 11 km depth show a good result, and the Vp inversion can reconstruct the pattern except for the corner areas of the model. At 15 km, the sign of the checkerboard pattern is neutral, which explains the diffuse reconstructed pattern. In contrast, the 3 km depth section pattern shows partial smearing effects in the central and southern areas, indicating a lower resolution. A comparison with the RDE confirms the results of the checkerboard test. Areas of reconstructed patterns have a higher RDE value (≥0.15), and vice versa, areas without reconstructed patterns have a lower RDE value (<0.15). The Vp/Vs checkerboard pattern can be reconstructed in the model's southern and some central areas at depths of 3 to 11 km. Comparing the Vp/Vs checkerboard test to the RDE confirms that the resolu- tion is high in the centre and the southern part, whereas the north has no or low resolution. Additional figures of the checkerboard test can be found in the Appendix (Fig. 16). In some areas of the slices there is a high RDE value but no reconstructed checkerboard pattern.
In these areas, interpretation should be done with caution. Figures 8 and 9 show the inversion results for Vp and Vp/Vs. Vp values range from 5.4 to 6.4 km/s, while Vp/Vs varies from 1.63 to 1.75. Horizontal slices of Vp and Vp/Vs (Fig. 8) are displayed at the nodes of the inversion grid from 3 km to 15 km depth in 4 km steps. Black dots display the earthquakes within 1 km of the corresponding depth slices.

Results
In addition, three vertical sections are displayed ( Fig. 9) along the profiles P1, P2 and P3 crossing the earthquake swarm area Schöneck (P1) and Werdau (P2, P3). For orientation, the vertical slices are displayed in the velocity models (Fig. 8). Earthquakes used in the inversion and occurring near the vertical sections are perpendicularly projected and displayed if the projection distance is less than 5 km. The final RMS error of both models was about 54 ms. The 3-D inversion can better explain the residual between observed and calculated travel times (see Fig. 17) over the entire distance range than the 1-D inversion. In particular, the residuals of the P-wave show a clear reduction.

Receiver function
Receiver functions (RFs) allow us to observe the seismic discontinuities of the Earth within the crust and upper mantle below a station (Vinnik 1977). The method of P-wave RFs is based on the effect that a portion of the incident P-wave is converted into a S-wave at a discontinuity. As a result, the propagation velocity of this converted phase is slower and arrives at the station after the primary phase. This time difference can be used to detect discontinuities and determine their depths.

Data and method
All broadband seismic stations that were openly accessible within the latitude range of 49.95 • N to 51.04 • N and longitude range of 12.0 • E to 12.6 • E were used. A total of 18 stations from different networks with network codes BW, CZ, GR, SX and TH are available in the area, distributed along a N-S aligned strip. The Preceiver functions were calculated with the software rf by Eulenfeld (2020). We selected all available earthquakes of the Geofon catalogue (Quinteros et al. 2021) with a moment magnitude M W larger than 5.5 and an epicentral distance between 30 • and 90 • from 2011 to 2020 (Fig. 10). The waveform data was then filtered between 0.05 Hz and 1 Hz. Waveforms with a signal-tonoise ratio smaller than 2 compared to a time window -20 s to -15 s prior to the P-wave were removed. Finally, the data were trimmed to 0 s to 100 s relative to the Ponset, rotated into the LQT system and subsequently deconvolved using iterative time-domain deconvolution (Ligorría and Ammon 1999). We calculated a total of 6844 RFs. Figure 11a) shows the stacked RFs of the Q-component of each station. Earthquakes in the range of 1 km of the corresponding depth are displayed as black dots. The location of profiles P1, P2 and P3 are given for orientation in the 3 km depth slice. The phyllite (S1), the graywacke (S2), presumably fluid pathways (F) and quartz-rich granitic composition intrusions (I1-I3) mentioned in the discussion section are highlighted

Results
At all stations, a strong positive amplitude is visible between 2.5 and 4.5 s, which can be interpreted as the Ps-conversion of the Moho. The weaker Moho multiples are partly visible at about 12.5 to 13.5 s (positive PpSs) and between 14 and 15 s (PpPs+PsPs). Figure 11b) shows a detailed example of individual RFs of the station SCHF. Noticeable differences between the RFs between 60 and 325 degrees back-azimuth and the RFs outside this range are visible at the station SCHF. This effect is also visible at some other stations. Affected stations are slightly south and north of SCHF (MULD to GRZ1) and the effect is always visible at RFs coming from the same back-azimuth range.
According to Cossette et al. (2016), the difference can be explained by two effects. Either there is a largescale structural heterogeneity below the affected stations, which should be visible in the tomography, or elastic anisotropy in the crust causes the effect.

Moho depth and Vp/Vs ratio
Using the method of Zhu and Kanamori (2000), the delay times of the Moho Ps conversion and their multiples can be used to determine the depth of the Moho and the Vp/Vs-ratio (Table 1). For the grid search, an average crustal P-wave velocity of 6.3 km/s was assumed. We estimated the uncertainty by the 95% margin of the In addition above the models, the intersecting areas of the fault zones are marked. The phyllite (S1), the graywacke (S2), presumably fluid pathways (F) and quartz-rich granitic composition intrusions (I1-I4) mentioned in the discussion section are highlighted maximal stacked amplitude (see Appendix Fig. 18). For the station MGBB it was impossible to determine the depth and Vp/Vs, because the Ps-conversion and the multiples did not join together at one point. The grid search yields an average Moho depth of 30.5±1.5 km and an average Vp/Vs value of 1.73±0.07 for the study area. Overall, the results, within the error bars, are consistent with the results of Geissler et al. (2005) and Heuer et al. (2006).
Beneath the NW-Bohemia area, the Moho domes up to 25.0 to 30.0 km. For the Vogtland area, the Moho depth increases to 30.0 to 32.8 km. North of it, the Moho domes up again to 29.1 to 30.9 km. The Vp/Vs ratios follow a slight north-south gradient, with the exception of ROHR. The Vp/Vs values for the southern area are normal to very high (1.69 -1.92) and, with the exception of HKWD where the error interval is significantly larger, the northern area shows normal Vp/Vs values (1.68 -1.74). A systematic error in determining the Moho is caused by using a constant velocity for the crust. Increasing or decreasing the average velocity by 0.2 km/s, the Moho depth increases or decreases by 1 km on average, whereas the Vp/Vs ratio is almost independent of the P-wave velocity. Since this systematic error is smaller than the error estimate at most stations, we can neglect it.

Migration
To create a depth section from the time traces, a migration (Yuan et al. 2000) was carried out. For the migration, we used a modified version of the IASP91 model (Kennet and Engdahl 1991), with the velocities of the first 20 km depth taken from our 1-D velocity model obtained in Section 2.2. Using the velocity model, the RFs are back-projected along the range of Fresnel zones and stacked within a 2-D plane. To estimate the migration uncertainty, we compare the migration once with the modified and unmodified version of the IASP91 model. Using the modified IASP91 model, the anomalies are 2 km deeper than when using the unmodified IASP91 model. Therefore, we propose the 2 km as an uncertainty for the migration. Figure 12 shows the resulting migrated N-S profile calculated from 6844 A distinct positive phase (R1) is visible between 27 and 31 km depth, another positive phase (R2) at 55 km depth which dips to 65 km depth towards the north, and below a negative phase (R3) at about 68 km depth. The positive phase R2 displaces the negative phase R3 at about 50.7 • N and then disappears at 50.8 • N. The reflectors R2 and R3 are not visible north of 50.8 • N. The R2 and R3 phase can also be observed in the single station RFs. The positive R2 phase is visible between GUNZ and PLN at about 6 to 6.5 s and the negative R3 phase is visible between NKC and PLN at about 7.5 to 8 s.

Crustal structure in the area of investigation
A schematic representation of the interpretation is shown in Fig. 13. The depth slices of the velocity models reveal a three-layered crustal structure (Fig. 9), with an upper heterogeneous crust that extends down to 6 -8 km, a middle crust that extends down to a depth of about 17 km and a lower crust underneath. The velocity anomalies within the upper crust can be related to the geology of the surface. Following Christensen (1996), the low P-velocities of 5.4 to 5.6 km/s and Vp/Vs-values of 1.65 -1.68 in the 3 -7 km depth range can be assigned to a graywacke (S2) and the higher P-velocities of 5.9 -6.1 km/s and high Vp/Vs-values of 1.71 -1.74 represent phyllite and/or granite (S1). Since the seismic properties of phyllite and granite are similar and both are present in the area, we cannot further distinguish. The S1 and S2 anomalies are exposed at the surface and follow the NE-SW trend of the Saxothuringian (Pietzsch 1963), which is also visible in the velocity model at 3 km depth (Fig. 8). The thickness of the upper crust of 6 -8 km also fits with the indicated thickness of the Saxothuringian (Pietzsch 1963).
Below the upper crust is the middle crust at depth of 8 -17 km. The study of Enderle et al. (1998) determined a P-wave velocity of 6.3 km/s for the depth range of 7 -15 km using an active seismic refraction study and proposed a felsic granulite composition for the depth range. For the range from 8 to 17 km depth, our velocity model shows Vp values of 5.85 -6.25 km/s, whereas the Vp/Vs ratio shows values between 1.64 and 1.73. Assuming a lithostatic pressure of around 200 -425 MPa for the depths of 8 -17 km, the Vp and Vp/Vs values would, according to Christensen (1996), correspond to a granite gneiss composition rather than the felsic granulite composition proposed by Enderle et al. (1998).
Below 17 km, the velocity increases to more than 6.3 km/s, representing the transition from the middle to the lower crust. Thus, the crustal structure of the study area corresponds to the typical structure of the Central European crust (Tesauro et al. 2008). Furthermore, the brittle-ductile transition (Fig. 13) can be determined with the earthquake distribution. Since very few to no earthquakes occur in the ductile part of the crust (Ellis and Stöckhert 2004), the boundary is 14 km in the NW-Bohemia area and 24 km in the northern part of the LRZ.

Seismic properties of the Schöneck and Werdau earthquake swarm areas
The two earthquake swarm areas Schöneck and Werdau show similar properties in the tomographic models. Both areas (Fig. 9) are inconspicuous in the Vp model and cannot be distinguished from the respective surrounding. In the Vp/Vs ratio, however, there are low Vp/Vs anomalies (I1, I3, I4) of approximately 1.67 near the earthquake swarm areas. Besides these low Vp/Vs anomalies, there is another low Vp/Vs anomaly (I2) in the northwest of Werdau where the 2016 Werdau swarm occurred.
The velocities and ratios (Vp 5.8 km/s and Vp/Vs of 1.67 for all Schöneck clusters and Vp 6.1 km/s and Vp/Vs of 1.68 for Werdau) can be explained by a quartz-rich granitic composition intrusion in the granite gneiss environment of the middle crust. Since we cannot determine the age of the Werdau intrusions, they could be related to the intrusions that emplaced to the east of the study area during the Variscan orogeny (Förster et al. 1999) or, to the diatreme of Ebersbrunn and Weißenbrunn, which lie only a few kilometres   Geissler et al. (2005) south of Werdau (Schmidt et al. 2013). Further investigations are necessary. Geodynamic modelling has shown for the Vogtland area that the deformation and the shear stress caused by the regional stress field are too low to generate earthquake swarms (Kurz et al. 2003). The regional stress field is only responsible for the prestressing of the crust (Kurz et al. 2003) and an additional physical process in the form of fluids triggers the earthquake swarms (Parotidis et al. 2005). The additional periodic pore pressures and temperature changes caused by ascending magmatic fluids ultimately trigger earthquakes in the Vogtland (Kurz et al. 2003). Since the earthquake swarm areas Schöneck and Werdau are directly adjacent to the Vogtland earthquake swarm area, we assume a similar situation where the regional stress field cannot trigger the earthquake swarms on its own. More-over, the earthquakes at Schöneck and Werdau show swarm character, so fluids must also be present. Numerous studies have shown that high Vp/Vs ratios image fluid pathways (Husen et al. 2004;Kato et al. 2010;Muskin et al. 2013). In the Vogtland area, the potential fluid path passing through the earthquake swarm area could already be imaged by Mousavi et al. (2015) by means of a channel-like increased Vp/Vs structure. A similar structure passing through the Schöneck or Werdau earthquake swarm area is not visible in our study. Instead, we see larger areas (F) of presumably high Vp/Vs ratio (1.72 -1.73) around the earthquake swarm areas and the low Vp/Vs anomalies (1.67 -1.68), indicating the presence of fluids (Fig. 9). Since the checkerboard test could not reconstruct the pattern for the area of the fluids, a clear interpretation about the presence of fluids is not possible. Thus we suggest that fluids rise from the lower crust and presumably emerge diffusely at the surface, since in this area, in contrast to NW-Bohemia/Vogtland, there are no known concentrated degassing sites e.g. in the form of mofettes.
Similar to the Vogtland earthquake swarm (Mousavi et al. 2015), the relocated earthquake swarms are distributed along an axis that runs through the low Vp/Vs anomaly in the case of Schöneck and slightly north of the anomaly of Werdau. We assume that smaller fractures, which we cannot resolve with the tomography, run through the intrusive body and fluids rise along these fractures. The fractures, in turn, could be related to the BKZ and GJZ. Therefore we assume the combination of a regional stress field prestressing the crust, fault zones that accumulate stress in the intersection, magmatic fluids ascending along these fault zones and a change in the material property in the form of a quartzrich granitic intrusiva are the cause of the Schöneck and Werdau earthquake swarms.

Crustal mantle transition and fluid migration
Based on the receiver function analysis (Fig. 12), the deeper area of the crust and the transition to the upper mantle can be analysed. The method of Zhu and Kanamori (2000) and the migration, where a strong positive phase (R1) is visible, yield a Moho depth of 25 -33±2.0 km. We can observe an updoming of the Moho below NW-Bohemia, which previous studies had already observed (e.g. Geissler et al. 2005;Heuer et al. 2006;Hrubcová et al. 2013). The attenuation values of the crust in the southern part of the LRZ (Gaebler et al. 2015;Eulenfeld et al. 2021) are higher than in the northern part (van Laaten et al. 2022a), which already indicates an increased magma/fluid/heat transport in the southern part of the LRZ, which in turn can be attributed to the updoming Moho (Heuer et al. 2006). Geissler et al. (2005) and Heuer et al. (2006) had already observed the positive phase (R2) at a depth of 55±2.0 km below NW-Bohemia. According to Heuer et al. (2006), it is a seismic reflector below NW-Bohemia extending further north and south, while Geissler et al. (2005) suggest that the seismic disconti-nuity is related to the base of a metasomatic uppermost mantle containing a few per cent melt. In our study, we can trace the R2 reflector to the GJZ/LRZ intersection, with the reflector dipping to the north. The negative phase (R3) at 68±2.0 km is more prominent than the R2 anomaly, indicating a stronger velocity contrast. Heuer et al. (2006) also observed the R3 anomaly in their study. According to their modelling, it may be a velocity decrease of 3.7 ± 1.0% at about 65 km depth or a 5 km thick gradual velocity decrease between 65 and 70 km. Heuer et al. (2006) interpret the velocity decrease as a local uplift of the lithosphere-asthenosphere boundary below NW-Bohemia. An upwelling of the lithosphereasthenosphere boundary was also interpreted by Plomerová et al. (2007), based on the results of their teleseismic P-velocity tomography study in the Bohemian Massif. The updoming could account for increased partial melt in the lithospheric mantle (Heuer et al. 2006). This in turn may lead to mantle-induced degassing of CO 2 , which rises and causes the earthquake swarms in the NW-Bohemia area. We can also observe the same anomaly and trace it northward up to the GJZ/LRZ intersection. At the intersection, the R3 anomaly is then bounded by the R2 reflector. A schematic interpretation is shown in Fig. 13. The same melting that causes the NW-Bohemia earthquake swarms, thus causing fluids to rise along the region's fault zones and trigger the Schöneck and Werdau earthquake swarms as well.

Conclusion
This study presents the results of LET and RFs to image the crustal and upper mantle structure. The LET shows the transitions from upper to middle and middle to lower crust and that the near-surface geology matches the near-surface velocity models. Magmatic fluids, characterised by increased Vp/Vs ratios, ascend and diffuse towards the surface. The swarms gather in small areas and are always near low Vp/Vs anomalies which we interpret as individual intrusion bodies. This intrusion body is probably a granite of older origin, as recent volcanic activity has not been observed in the area. In sum, the earthquake swarms are triggered by various factors: the regional stress field, which is responsible for the prestressing of the crust, intersecting fault zones causing an accumulation of stress and ascending mantle fluids. Thus, the earthquake swarms of Schöneck and Werdau are subject to a similar trigger mechanism as the earthquake swarms of NW-Bohemia/Vogtland. The RFs results provide the Moho depth, a seismic reflector at 55 km depth dipping to the north, and a zone of partial melting below 65 km depth as the origin of fluids. Werdau is thus the northernmost swarm earthquake area of the LRZ, as the melt zone does not extend further north. and the GFZ (Quinteros et al. 2021) for providing the waveform and earthquake data. We wish to thank Falk Hänel of TU Freiberg for providing the SXNET-offline waveform data and the LfULG (LfULG 2018) for the free access of geological data. Also we would like to thank Xiaohui Yuan and Felix Schneider for providing the RFs migration code and Thomas Burghardt for his support in the usage of SeismicHandler. SeismicHandler was funded by the Deutsche Forschungsgemeinschaft (DFG) within the framework of SPP-1066. We used the rf code from https:// github.com/trichter/rf (Eulenfeld 2020). Figures were generated with Generic Mapping Tools (Wessel et al. 2013) or Python using ObsPy (Krischer et al. 2015), NumPy (Harris et al. 2020) and Matplotlib (Hunter 2007). The focal mechanisms were obtained with Grond .
Author contributions Marcel van Laaten and Ulrich Wegler contributed to the study conception and design. The data processing was performed by Marcel van Laaten. All authors contributed to the data analysis and interpretation. The first draft of the manuscript was written by Marcel van Laaten. Ulrich Wegler and Tom Eulenfeld commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Funding
The authors did not receive support from any organization for the submitted work. Open Access funding enabled and organized by Projekt DEAL.
Availability of data and material The receiver function profile and velocity models are available in the Zenodo repository https://doi.org/10.5281/zenodo.6511543 (van Laaten et al. 2022b).

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.