We recently proposed the F-J method (Wang et al., 2019), 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 spectra C ˜ ( 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: I ( ω , c ) = ∫ 0 + ∞ C ˜ ( r , ω ) J 0 ( ω c r ) r d r , (1) 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: I ( ω , c ) = ∫ 0 + ∞ C ˜ ( r , ω ) J 0 ( ω c r ) r d r ≈ ∑ j = 1 N ∫ r j − 1 r j C ˜ ( r , ω ) J 0 ( ω c r ) r d r , (2) where 0 = r 0 < r 1 < r 2 < … < r j … < r N .

Within the interval [ r j − 1 , r j ] , the C ˜ ( r , ω ) can be approximated by a linear function as follows: C ˜ ( r , ω ) ≈ a j + b j r , (3) 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): ∫ x J 0 ( x ) d x = x J 1 ( x ) , (5) ∫ x 2 J 0 ( x ) d x = x 2 J 1 ( x ) + x J 0 ( x ) − ∫ J 0 ( x ) d x . (6)

We can get the following approximation formula which can be used when scanning to obtain the F-J spectrogram (Wang et al., 2019): I ( ω , c ) ≈ ∑ j = 1 N { ( c ω ) C ˜ ( r , ω ) r J 1 ( w c r ) + b j ( c ω ) 3 [ ω c r J 0 ( ω c r ) − B 0 ( ω c r ) ] } | r j − 1 r j , (7) where c and r _j denote the phase velocity and the interstation distance of station pair j, respectively. B 0 ( x ) = ∫ 0 x J 0 ( η ) d η .

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.

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: f ( V s ) = 1 m ∑ k = 0 m a k n k { ∑ i [ c i k s ( V s ) − c i k o ] 2 } + γ ‖ D ^ V s ‖ 2 2 , (8) where i and k are the indices for the sampled frequency and dispersion curve mode, respectively; c i k S is the phase velocity of the synthetic dispersion curve at the ith frequency and kth mode; c i k O 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 − | z i − z j | d , 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.67 v s and ρ = 0.77 + 0.32 v p (Růžek et al., 2016). 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 fundamental-mode 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.

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 fundamental-mode 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 higher-mode 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).

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 fundamental-mode 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 (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 higher-mode 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.

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 two new subarrays are shown in Supplementary Figure S5 and the F-J spectrograms of the two new subarrays are shown in Supplementary Figure S6. It is clear that the quality of the F-J spectrograms decreases, especially for the higher modes.

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).

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, Supplementary Figure S7B) beneath the stations of the four arrays (we select the 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 (Figures 13C,D; Supplementary Figure S7C, Supplementary 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.

The models obtained via the joint inversion of multimodal dispersion curves based on the F-J method reveal obvious low-velocity 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 low-velocity anomalies are also discovered in the uppermost mantle beneath the study area. The thickness of uppermost-mantle low-velocity 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 lithosphere-asthenosphere transition zone (e.g., Plomerová et al., 2007; Grad et al., 2008).

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²; 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 lower-density 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 (Horálek et al., 2000). 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 (Fischer and Horálek, 2000). 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 the S-wave LVZs, which could be responsible for the origin and foci depth distribution of earthquake swarms in the study area.