Crustal S-Wave Velocity Structure Beneath the Northwestern Bohemian Massif, Central Europe, Revealed by the Inversion of Multimodal Ambient Noise Dispersion Curves

The northwestern Bohemian Massif and adjacent areas are a tectonically active region associated with complex geodynamic activities, that manifest as Quaternary volcanism, earthquake swarms in the upper and middle crust, degassing of CO2, and crustal fluid migration. The intricate tectonic evolution and activities of this region reflect the complexity of the crustal structure therein. However, the crustal models derived from previous studies in this area offer different, even contradictory information regarding the existence of a mid-crustal low-velocity zone (LVZ). In this study, we apply the frequency-Bessel transform (F-J) method to extract the fundamental-mode and up to five higher-mode Rayleigh wave dispersion curves from ambient seismic noise data recorded in the study area and perform multimodal ambient noise dispersion curves inversion. The addition of higher-mode dispersion curves enhances the vertical resolution of the velocity structure inversion results. Our models support the view that the general S-wave velocity level of the crust is high within the study area. We detect two S-wave LVZs beneath the study area that are distributed mainly in the middle crust rather than the lower crust, and these LVZs are separated by a high-velocity zone. Considering the results of previous studies in the area, we infer that these S-wave LVZs may be the consequence of crustal fluids, plastic deformation and even partial melting of the felsic middle crust at relatively high crustal temperatures. Furthermore, these S-wave LVZs could be responsible for the origin and foci depth distribution of earthquake swarms. S-wave low-velocity anomalies are also observed in the uppermost mantle beneath the study area. These S-wave models based on the joint inversion of multimodal dispersion curves can provide new references for understanding the tectonic activity and geodynamic evolution of the northwestern Bohemian Massif and adjacent areas.


INTRODUCTION
The Bohemian Massif, which consists of four geological units, namely, the Saxo-Thuringian (ST), Teplá-Barrandian (TB), Moldanubian and Sudetes zones, is one of the largest stable outcrops of basement rocks in Central Europe Karousová et al., 2012;Plomerová et al., 2012). Tectonically, the Bohemian Massif forms the easternmost rim of the Variscan belt that developed between approximately 480 and 290 Ma during the collision between Laurussia and Gondwana (Matte, 2001;Heuer et al., 2006;Hrubcová et al., 2008;Valentová et al., 2017). The area studied herein is the key tectonic area of the Bohemian Massif, namely, the convergence area of the four geological units mentioned above (as shown in Figure 1), including the northwestern part of the Bohemian Massif and adjacent areas. This region, which is presently tectonically active and exhibits complex geological features (Kolínský and Brokešová, 2007;Mousavi et al., 2017), composes a part of the European Cenozoic Rift System and is characterized by a relic Devonian oceanic suture (Karousová et al., 2012;Plomerová et al., 2007). The ongoing geodynamic activities manifest as Quaternary volcanism, earthquake swarms in the upper and middle crust, degassing of CO 2 , and crustal fluid migration (Fischer and Horálek., 2003;Horálek and Fischer., 2008;Babuška et al., 2016). The signatures of the intricate regional tectonic evolution and ongoing tectonic activities are recorded in and reflect the complexity of the crustal structure in this area.
Previous studies on the crust of the northwestern Bohemian Massif and adjacent areas have relied mainly on seismic sounding experiments, that principally provided P-wave velocity models, such as the NW-SE-trending CEL09 profile (Růžek et al., 2007;Hrubcová et al., 2005;Novotný, 2012) and S04 profile (Růžek et al., 2007;Hrubcová et al., 2010) and the NE-SW-trending S01 profile (Růžek et al., 2007;Grad et al., 2008) and 95-B profile (Enderle et al., 1998). These seismic sounding experiments revealed that the middle crust is generally laterally homogeneous, while the lower crust is laterally inhomogeneous. Along profile 95-B, some regions exhibit high P-wave velocities in the upper and middle crust (depths of~5 and 16 km, respectively), suggesting the presence of a low-velocity zone (LVZ) in the depth range of~5-16 km. The crustal structure is more complicated under the Eger Rift along profile S01; P-wave high-velocity bodies (HVBs) exist mainly in the upper crust (~2-10 km) and feature deep roots extending into the middle crust (~18 km).
Moreover, some studies have recently applied passive seismic techniques in this area. By analyzing teleseismic records with the receiver function technique (e.g., the receiver functions at stations MOX, WET, BRG, and NKC shown in Figure 1), Wilde-Piórko et al. (2005) revealed the existence of an S-wave LVZ in the middle crust (depth range of 10-15 km) of the northwestern Bohemian Massif. With the help of Love wave phase velocity dispersion curves, Kolínský et al. (2011) discovered an S-wave LVZ in the middle and lower crust beneath the TB zone and a gradually increasing S-wave velocity structure in the crust of the ST zone. By applying ambient seismic noise interferometric surface wave tomography (ASNT) to the recordings of broadband seismic stations, Růžek et al. (2016) inverted group and phase dispersion curves to obtain crustal velocity models for the Bohemian Massif that monotonically increase with depth; they concluded that the differences among different tectonic units are small and that the most homogeneous part among them generally being the middle crust. Some researchers utilized the ASNT method with data from additional seismic stations in Europe to extract phase velocity dispersion curves (e.g., Kästle et al., 2018) and group velocity dispersion curves (e.g., Lu et al., 2018) for the inversion and derived high-resolution S-wave velocity models that present the almost smoothly increasing velocity structure of the crust in the study area. More recently, Kvapil et al. (2021) reported that the velocity-drop interface (negative velocity gradient) in the lower part of the crust of the Bohemian Massif (depth of 18-30 km); however, the general group velocity level of dispersion curves is lower than the level presented in Lu et al. (2018).
Nevertheless, the crustal models derived from previous studies for the Bohemian Massif present different, even contradictory results regarding the existence of the LVZ in the middle crust (e.g., Enderle et al., 1998;Wilde-Piórko et al., 2005;Grad et al., 2008;Kolínský et al., 2011;Růžek et al., 2016;Kästle et al., 2018;Lu et al., 2018;Kvapil et al., 2021). To address this inconsistency, in this study, we apply our newly developed multimodal ambient noise dispersion curve tomography method, denoted the frequency-Bessel transform (F-J) method , to investigate the crustal structure beneath the study area by using available ambient seismic noise datasets.
In the following, we describe how to process the ambient noise data used in this study. Then, we briefly summarize the principle of the F-J method and apply it to extract the fundamental-mode and higher-mode Rayleigh wave phase velocity dispersion curves from ambient seismic noise data, after which we carry out the joint inversion of these multimodal dispersion curves to obtain the S-wave velocity model for the crust and uppermost mantle with a higher vertical resolution. Finally, we compare our models with the results of previous studies conducted in our region of interest.

DATA AND PREPROCESSING
The continuous ambient seismic noise data analyzed in this study are derived from two independent datasets. The first dataset is from the Bohema digital seismic network (FDSN code: ZV), which is a part of the Bohemian Massif Anisotropy and Heterogeneity (BOHEMA) project (Babuška et al., 2005;Plomerová et al., 2003), from which we select continuous broadband vertical-component seismic noise data (channel code: BHZ/HHZ, sampling rate: 20-100 Hz) from January 2001 to December 2005 recorded by 49 stations. The second dataset is from the Thuringer Seismisches Netz (TSN) network (FDSN code: TH; Jena, 2009), from which we select continuous broadband vertical-component seismic noise data (channel code: HHZ, sampling rate: 100 Hz) from January 2015 to December 2017 recorded by 23 stations. All stations are located mainly in the northwestern Bohemian Massif and adjacent areas and span the area of approximately 6°× 3°, as shown in Figure 1. We take the central longitude and latitude of the area where the seismic array is located as the center of the array (green pentagrams in Figure 1).
The procedures for preprocessing the data are similar to those described in Bensen et al. (2007), including downsampling, tapering, detrending, removing the mean, removing the instrumental response, time-domain normalization and spectral whitening. Finally, we split the whole records from each station into 1-h segments. After processing the data from all stations, we use linear stacking method to compute 1-h stacked noise cross-correlation functions (NCFs) of all station pairs in the two datasets. Supplementary Figure S1 shows the crosscorrelation functions of different arrays (period band 1-50 s).

METHODOLOGY F-J Method
We recently proposed the F-J method , which can effectively extract multimodal dispersion curves from ambient seismic noise data (e.g., Wu et al., 2020;Zhan et al., 2020). Having obtained the stacked NCFs in the preprocessing step, we then extract the multimodal dispersion curves of the areas beneath each dataset by using the F-J method. The main procedures of the F-J method are briefly summarized as follows.
After preprocessing the ambient seismic noise data, we can obtain the frequency spectraC(r, ω) of the NCFs of a series of station pairs, where r is the interstation distance. If the station pair cross-correlation functions are continuously distributed, and the number is infinite, we can calculate the F-J spectrogram of C(r, ω) as follows: Where J 0 is Bessel function of order 0. We have theoretically proven that when c and ω satisfy the dispersion relationship, the value of I(ω, c) tends to infinity. When scanning with limited pixels, the F-J spectrogram I(ω, c) shows significant maximum values in the narrow neighborhood of the dispersion curve of each mode, and thus, we can identify the dispersion curves of the fundamental mode and higher modes from the I(ω, c) diagram. However, in practical applications, the number of stations is finite; hence, we cannot directly apply Eq. 1 to obtain I(ω, c). We can approximate the infinite integral in Eq. 1 by the following truncation: where 0 r 0 < r 1 < r 2 < . . . < r j . . . < r N . Within the interval [r j−1 , r j ], theC(r, ω) can be approximated by a linear function as follows: with a j C r j−1 , ω − ⎡ ⎣C r j , ω −C r j−1 , ω r j − r j−1 ⎤ ⎦ r j−1 and b j C r j , ω −C r j−1 , ω r j − r j−1 .
( 4 ) Substituting Eq. 3 into Eq. 2, and with the following formulas ( J 1 is Bessel function of order 1): We can get the following approximation formula which can be used when scanning to obtain the F-J spectrogram : where c and r j denote the phase velocity and the interstation distance of station pair j, respectively. B 0 (x) x 0 J 0 (η)dη.

Identification of Multimodal Dispersion Curves
To accurately extract higher-mode dispersion curves, we extract the multimodal dispersion curves in two steps. The workflow is illustrated in Figure 2. First, we extract the fundamental-mode dispersion curve from the F-J spectrogram and invert the dispersion curve to obtain a preliminary S-wave velocity model (hereinafter referred to as the fundamental-mode velocity model). Then, we calculate the theoretical higher-mode dispersion curves of the model using the generalized reflection-transmission coefficient method (Chen, 1993) and project them onto the F-J spectrogram. Finally, taking these theoretical higher-mode dispersion curves as a reference, we identify the highlighted areas on the F-J spectrogram and extract the higher-mode dispersion curves.

Inversion Method
Having obtained the dispersion curves, we carry out the inversion by using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (e.g., Byrd et al., 1995). The detailed procedures are described in Pan et al. (2019) and Zhan et al. (2020); a brief summary is given as follows. First, an objective function of multimodal dispersion curves is defined by: where i and k are the indices for the sampled frequency and dispersion curve mode, respectively; c S ik is the phase velocity of the synthetic dispersion curve at the ith frequency and kth mode; c O ik is the observed velocity; m is the number of modes of the dispersion curves used for the inversion; n k is the number of sampled data points for the kth mode; and a k is a weight factor for the kth mode. In this study, we set the weight factor of each higher-mode dispersion curve to 1, and the weight factor of the fundamental mode dispersion curve is equal to the number of all higher-mode dispersion curves. Through this simple strategy, it is possible to ensure that the fundamental mode makes the main contribution to the inversion, as well as the improvement of the inversion by the higher modes. The second term is the smoothing regularization. V S is the shear velocity model;D e where z i and z j are the depths at the top of the ith and jth layers; and d is a smoothing distance (Haney and Tsai, 2017). The smoothing factor γ is near the maximum curvature of the L-curve (Hansen, 2001) and the value is between 3e-3 and 3e-2. In the inversion, we set the layer thickness to 2 km in the depth range of 0-68 km.
In the process of iteratively solving the nonlinear inversion problem, the P-wave velocity v p and density ρ of each iteration are converted by the following empirical formulas: v p 1.67v s and ρ 0.77 + 0.32v p . For the first step in inverting the fundamental-mode dispersion curve, 200 initial models are randomly generated in the range of ±0.4 km/s with the Eurasian 1D average reference model (Marone et al., 2004) as the intermediate value for the inversion, and the best-fitting model is taken as the fundamental-mode velocity model. After obtaining the higher-mode dispersion curves, the fundamentalmode velocity model is used to randomly generate 200 initial models within the range of ±0.4 km/s for the multimodal dispersion curve inversion. The model that minimizes the objective function is taken as the final model for the inversion of multimodal dispersion curves.

Identification of Multimodal Dispersion Curves and Inversion Results
The distribution area of stations in seismic array "TH" is adjacent to the northwestern Bohemian Massif (yellow triangles in Figure 1). We apply the F-J method to this array and the results are shown in Figure 3. Figure 3A clearly shows the fundamental-mode dispersion curve, as well as the possible higher-mode dispersion curves. The fundamental-mode velocity model is obtained by inverting the fundamental-mode dispersion curve. Figure 3B shows the first ten well-fitting inversion results; the red line represents the model with the smallest objective function, namely, the fundamental-mode velocity model, and its fundamental-mode dispersion curve fitting is shown in Figure 3C. According to the fundamentalmode velocity model, we calculate the theoretical higher-mode dispersion curves (yellow solid lines in Figure 3A) and project them onto the F-J spectrogram. With this projection as a reference, we extract the fundamental-mode and five highermode dispersion curves (white dotted lines in Figure 3A).
To invert the multimodal dispersion curves in the area of seismic array "TH", we use the fundamental-mode velocity model (red line in Figure 3B) as the intermediate value within the range of the fundamental-mode velocity model (±0.4 km/s, gray dashed lines in Figure 4A) to randomly generate 200 initial models for the multimodal dispersion curve inversion. Figure 4A shows the first ten well-fitting inversion results (the standard deviations of the ten models are shown in Supplementary Figure S2), while the red line represents the best-fitting model with the smallest objective function value, that is, the final model of the multimodal dispersion curve inversion. The multimodal dispersion curve fitting result of this model is shown in Figure 4B. The model depicts an S-wave LVZ in the middle crust (depths of 12-22 km), and the low-velocity anomalies in the uppermost mantle (depths of 42-62 km).
The stations of seismic array "ZV" are distributed mainly in the northwestern Bohemian Massif (red triangles in Figure 1). Figure 5A shows the picked points of the multimodal dispersion curves and the theoretical multimodal dispersion curves corresponding to the fundamental-mode velocity model (blue line in Figure 5B). The inversion and fitting results of the multimodal dispersion curves are shown in Figures 5B,C. The final model (red line in Figure 5B) of the multimodal dispersion curve inversion similarly reveals an S-wave LVZ in the middle crust (depths of 8-12 km) and low-velocity anomalies in the uppermost mantle (below 50 km).

Sensitivity Kernel Analysis
To illustrate the influence of the higher-mode Rayleigh wave dispersion curves on the velocity structure inversion results, we employ the method of Pan et al. (2019) to calculate both the depth and the frequency distributions of the S-wave sensitivity kernel function of the final inversion model in each of the above two regions (as shown in Figure 6 and Figure 7). The distributions of these sensitivity kernel functions show that the fundamentalmode dispersion curve offers constraint on the entire crust and uppermost mantle in the study area, and the constraint is the strongest at the surface and weakens with depth. Furthermore, the higher-mode dispersion curves strongly constrain both the entire crust and the structure near the crust-mantle boundary  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 838751 6 (depth of approximately 40 km). We also employ the final models of the two arrays as the true models to test the improvement effect of higher modes on the inversion. The Supplementary Figure S3 and Supplementary Figure S4 clearly show that when only the fundamental mode is used for the inversion, the inversion models are considerably different from the true models; with the addition of higher modes to the inversion, the inversion models are much closer to the true models at depths of 0-40 km, and the inversion accuracy at depths of 40-70 km is also improved. Therefore, the addition of higher-mode dispersion curves to the inversion directly and significantly improves the inversion accuracy in the depth range of 0-40 km and thus improves the inversion accuracy in the whole inversion depth range. Introducing highermode dispersion curves on the basis of the fundamental-mode dispersion curve can provide more constraints on inversion, which can help mitigate the non-uniqueness problem.

Subregion Division
Obvious S-wave LVZs are detected in the middle crust from the velocity structure inversion results beneath the two seismic arrays. To further study the distributions of the S-wave LVZs, we divide the stations on the west side of the study area into two subarrays, as shown in Figure 8. Although some of the stations of the two subarrays overlap, the geometric centers of the subarrays are different (the green pentagrams in Figure 8) and are approximately 60 km apart. Figure 9 shows the F-J spectrogram of subarray 1 and the picked points of the dispersion curves, the final inversion model, and the fitting results of the dispersion curves. The inversion result in this area reveals a LVZ at depths of 12-16 km in the middle crust; furthermore, the velocity gradient at depths of 22-26 km is very small, and the velocity in the uppermost mantle (below 42 km) decreases with depth. Likewise, the F-J spectrogram, the picked points of the dispersion curves, and the inversion results for subarray 2 are shown in Figure 10. The final inversion model of this area shows two S-wave LVZs in the crust at depths of 8-12 km and 18-24 km; in addition, low-velocity features appear in the uppermost mantle at depths of 36-44 km and below 50 km.
The current sub-division method considers the quality of the F-J spectrogram and the regional structural characteristics as much as possible. If we try to divide into smaller subregions, the quality of the F-J spectrogram will be lower. For example, based on Figure 8, we remove the stations on the northwest edge of the two subarrays; the   Figure S6. It is clear that the quality of the F-J spectrograms decreases, especially for the higher modes.

The Robustness of the LVZs in the Middle Crust
To verify the robustness of the LVZs in our models, we carry out the random inversions while prohibiting a decrease in velocity with depth. For the joint inversion of multimodal dispersion curves, we obtain 200 inversion models without LVZ for each array. Figure 11 shows the first ten well-fitting models of each array. We calculate the objective function values of the initial models and the inversion models of each array when the velocity decrease is prohibited, and compare them with those when the velocity is allowed to decrease (shown in Figure 12). According to the statistical significance tests and the comparisons shown in Figure 12, the objective function values of the initial models in two cases have little difference (p-value >0.05); furthermore, compared to the initial models, the inversion models in each case are significantly improved (p-value< 0.01), and the inversion models obtained by allowing the velocity to decrease are obviously better than those when the velocity is prohibited from decreasing with depth (p-value <0.01); that is, the models in which the velocity is allowed to decrease are more reasonable. These outcomes confirm that the LVZs in the inversion results are robust. In addition, based on 3D gravity modeling, geological data, seismic refraction (CEL09) and reflection (9HR), Guy et al. (2011) suggested the lower-density middle crust of the ST zone (~10-25 km).

Crustal S-Wave Velocity Models in the Northwestern Bohemian Massif
Recently, some researchers have used the traditional ASNT method to extract phase velocity dispersion curves (e.g., Růžek et al., 2016;Kästle et al., 2018) and group velocity dispersion curves (e.g., Růžek et al., 2016;Lu et al., 2018;Kvapil et al., 2021) to study the velocity structure in the study region and obtained S-wave velocity models. Based on the velocity models of the previous studies, we obtain average 1D S-wave velocity models ( Figures 13A,B; Supplementary Figure S7A velocity models of the grid points nearest to the stations to calculate the averaged 1D models). Furthermore, the average 1D velocity model under the area of arrays "ZV" and subarray 2 (yellow line in Figure 13B and Supplementary Figure S7B) is obtained according to the 1D models of Růžek et al. (2016) in the ST and TB units. In comparison, the S-wave velocities in our models are higher than those in the other models at depths of 30-40 km. In addition, the S-wave velocities in the models from Kvapil et al. (2021) are obviously lower than those in other models (shown in Figure 13, Supplementary Figure S7 and Supplementary Figure S8), even the models of Lu et al. (2018), who also used group velocity dispersion curves for inversion. The general group velocity level of the dispersion curves in Kvapil et al. (2021) is lower than the level presented in Lu et al. (2018), which may be the main reason for the models' discrepancy.
For a comparison with the picked data, we calculate the theoretical dispersion curves of these 1D velocity models  Figure S7D), the results of which clearly demonstrate an unsatisfactory match between the theoretical dispersion curves of these models and the F-J spectrograms (especially at the higher modes). Specifically, the theoretical dispersion curves corresponding to the 1D models under the four arrays extracted from the models of Lu et al. (2018) and Kvapil et al. (2021) are overall lower than picked data; the two 1D models using the phase velocity dispersion curves for the inversion (Kästle et al., 2018;Růžek et al., 2016), especially the models from Kästle et al. (2018) present better fitting results for fundamental mode, but large deviations at the higher modes are still observed for the two models. Our models obtained via the joint inversion of the multimodal dispersion curves match well with the fundamental mode and higher modes for the four arrays. Therefore, introducing higher-mode dispersion curves on the basis of the fundamental-mode dispersion curve is crucial to constrain the structure of the crust in the study area.

S-Wave LVZs in the Middle Crust of Northwestern Bohemian Massif
The models obtained via the joint inversion of multimodal dispersion curves based on the F-J method reveal obvious lowvelocity characteristics in the crust and uppermost mantle beneath the study area. Other studies similarly provided evidence for the existence of an LVZ in this region (e.g., Wilde-Piórko et al., 2005;Kolínský et al., 2011;Kvapil et al., 2021). Based on our inversion results, we construct a simple S-wave velocity profile of the crust and uppermost mantle beneath the western side of the study area ( Figure 14B). The surface line corresponding to this profile in the study area trends NW-SE (the direction of the black arrow in Figure 14A) and mainly covers two geological tectonic units: ST and TB. Along the profile, there are two obvious LVZs (depth range of 8-24 km) in the middle crust separated by the high-velocity zone (HVZ, depth range of 12-18 km), which is different from the crustal S-wave velocity models of previous investigations (e.g., Wilde-Piórko et al., 2005;Růžek et al., 2016;Kästle et al., 2018;Lu et al., 2018;Kvapil et al., 2021). In the middle crust, the uppermost LVZ is thick in the west and thin in the east, and the velocity contrast of the LVZ is strong in the west and weak in the east; the thickness of the lowermost LVZ is relatively uniform, and the velocity contrast of the LVZ is weak in the west and strong in the east. The negligible HVZ between the upper and lower LVZs in the middle crust may be caused by the high-velocity crust beneath the Eger Rift.
In addition to this mid-crustal S-wave LVZs, S-wave lowvelocity anomalies are also discovered in the uppermost mantle beneath the study area. The thickness of uppermost-mantle lowvelocity anomalies is about~8-20 km. These S-wave low-velocity anomalies may be related to the partial melting in the mantle near the Eger Rift or the upwelling of materials from the lithosphereasthenosphere transition zone (e.g., Plomerová et al., 2007;Grad et al., 2008).

Causes of the Mid-crustal S-Wave LVZs in the Study Area
The mid-crustal S-wave LVZs in the study area and the low-velocity anomalies in the uppermost mantle are not directly connected; instead, they are separated by the high-velocity lower crust. Due to the influences of temperature and pressure, the middle crust may be relatively plastic or even partially molten. The surface heat flow is an important parameter for understanding geothermal activity, which is related to regional and global tectonic activities. The average heat flow in the Bohemian Massif is 67.9 mW/m 2 ; low heat flow values are observed in the southern and central parts, while high heat flow values are measured in the northwestern Bohemian Massif (Čermák, 1976). Previous studies have shown that the middle crust below the study area is predominantly felsic. For example, Förster and Förster (2000) calculated the heat budget based on the surface heat flow and radiogenic heat production of the ST unit, and inferred that the middle crust is relatively felsic and the lower crust may be relatively mafic and less felsic. Based on 3D gravity modeling, geological data, and seismic refraction (CEL09) and reflection (9HR), Guy et al. (2011) suggested that the lowerdensity middle crust of the ST zone (~10-25 km) is felsic; similarly, they indicated that the lower-density crust underneath the TB area is also felsic. A higher surface heat flow value usually corresponds to a higher crustal temperature. Hence, the felsic middle crust beneath the northwestern Bohemian Massif could have undergone plastic deformation and even partial melting at relatively higher crustal temperatures, thereby forming the S-wave LVZs in the middle crust. Moreover, the gases and fluids originating from activities involving the upper mantle and the lower crust caused a series of geodynamic activities in the crust of the study area. Microearthquake activities, including single events and earthquake swarms, frequently occur in western Bohemia/Vogtland, with the depths of the microearthquake hypocenters varying between 3 and 23 km . The Novy Kostel focal zone (near seismic station NKC shown in Figures 14A,B) is dominant throughout the whole region of the West Bohemia/Vogtland earthquake swarms.
The foci of the microearthquakes in the Novy Kostel focal zone occur at depths of 6-11 km . These focal depths may be related to the uppermost mid-crustal S-wave LVZ beneath the study area. The crustal fluids in the Western Bohemia/Vogtland play a key role in bringing the faults from the subcritical to the critical state, which triggers the earthquake swarms in this region (Horálek and Fischer, 2008). Fluids in the crust can lower the S-wave velocity, which may be one explanation for these mid-crustal S-wave LVZs. Špičák et al. (1999) suggested that the earthquake swarms in the western Bohemian Massif may be caused by the magma intrusions and related fluid and gas release at depths of~10 km. The plasticity and presence of partial melting in the fluid-rich mid-crust resulting in FIGURE 13 | Velocity models of previous studies and their fitting results with the picked data in the study area. (A) S-wave velocity models for array "TH". Different colors denote the different average 1D S-wave velocity models beneath the stations of array "TH" (Kästle et al., 2018;Lu et al., 2018;Kvapil et al., 2021), and the red solid line denotes our model obtained from the inversion of multimodal dispersion curves. The black solid line denotes the model obtained from the inversion using only fundamental mode. (B) S-wave velocity models for array "ZV". The yellow line is the average of the 1D models below the ST and TB units in the study area from Růžek et al. (2016). (C) The theoretical dispersion curves of the different 1D velocity models for array "TH". The different colors denote the theoretical dispersion curves of the different average 1D S-wave velocity models in (A); for additional explanations, see Figure 4B. (D) The theoretical dispersion curves of different 1D velocity models for array "ZV". The different colors denote the theoretical dispersion curves of the different average 1D S-wave velocity models in (B); for other explanations, see Figure 5C.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 838751 the S-wave LVZs, which could be responsible for the origin and foci depth distribution of earthquake swarms in the study area.

CONCLUSION
In the northwestern Bohemian Massif and adjacent areas, the F-J method was employed to extract up to five higher-mode dispersion curves in addition to the fundamental-mode dispersion curve. The joint inversion of these fundamentalmode and higher-mode dispersion curves improved the vertical resolution of the velocity structure inversion results, allowing us to obtain high-resolution S-wave velocity models of the crust and uppermost mantle beneath the study area. Based on our models, we report the following novel insights: 1. Introducing higher-mode dispersion curves on the basis of the fundamental-mode dispersion curve is crucial to constrain the structure of the entire crust in the study area; 2. The general S-wave velocity level of the crust in the study area is higher than that in the models from Kvapil et al. (2021), which supports the views of previous studies Kästle et al., 2018;Lu et al., 2018); 3. The S-wave velocity of the crust in the study area is relatively high at depths of~30-40 km; 4. The S-wave LVZs are distributed mainly in the middle crust of the study area (~10-20 km) rather than the lower crust (e.g., Kvapil et al., 2021); 5. On the western side of the study area, there are two obvious LVZs in the middle crust which are separated by the HVZ, which is different from the crustal S-wave velocity models of previous studies (e.g., Wilde-Piórko et al., 2005;Růžek et al., 2016;Kästle et al., 2018;Lu et al., 2018;Kvapil et al., 2021).
The mid-crustal S-wave LVZs in the northwestern Bohemian Massif and its adjacent areas may be the consequence of crustal fluids, plastic deformation and even partial melting of the felsic middle crust at relatively high crustal temperatures. Furthermore, these S-wave LVZs could be responsible for the origin and foci depth distribution of earthquake swarms in the study area. In addition, we observed S-wave low-velocity anomalies in the uppermost mantle, especially near the Eger Rift, which may be related to partial melting in the mantle or the upwelling of materials from the lithosphere-asthenosphere transition zone. These S-wave models based on the joint inversion of multimodal dispersion curves can provide new references for understanding the tectonic activity and geodynamic evolution of the northwestern Bohemian Massif and adjacent areas. Moreover, considering the recent discovery of a widespread mid-crustal low-velocity layer beneath Northeast China (Zhan et al., 2020), the existence of a mid-crustal LVZ may be a common feature in tectonically active areas. Considerable in-depth research is needed to confirm this speculation.

AUTHOR CONTRIBUTIONS
The specific contributions of each author are as follows. QM: Conceptualization, Software, Writing-Original Draft, Validation, Data Curation, Investigation, Visualization; LP: Methodology, Software, Validation, Investigation, Visualization; J-nW: Methodology, Software; ZY: Software, Visualization; XC: Conceptualization, Methodology, Writing-Review; Editing, Resources, Supervision, Project administration. All authors read and approved the final manuscript. (B) A simple schematic diagram depicting the distribution of the S-wave LVZs in the middle crust and the S-wave low-velocity anomalies in the uppermost mantle along the profile. The earthquake swarms beneath the station NKC occur at depths of 6-11 km.