Transversely Isotropic Lower Crust of Variscan Central Europe imaged by Ambient Noise Tomography of the Bohemian Massif

. Recent development of ambient noise tomography, in combination with increasing number of permanent seismic stations and dense networks of temporary stations operated during passive seismic experiments, provides a unique 10 opportunity to build the first high-resolution 3-D shear wave velocity (v S ) model of the crust of the Bohemian Massif (BM). The velocity model with a cell size of 22 km is built by conventional two-step inversion approach from Rayleigh wave group velocity dispersion curves measured at more than 400 stations. The shear velocities within the upper crust of the BM are ~0.2 km s -1 higher than those in its surroundings. The highest crustal velocities appear in its southern part, the Moldanubian unit. The model provides compelling evidence for a regional-scale of velocity distribution. The Cadomian part of the region 15 has a thinner crust, while the crust assembled, or tectonically transformed in the Variscan period, is thicker. The sharp Moho discontinuity preserves traces of its dynamic development expressed in remnants of Variscan subductions imprinted in bands of crustal thickenings. A significant feature of the presented model is the velocity-drop interface (VDI) modelled in the lower part of the crust. We explain this feature by anisotropic fabric of the lower crust, which is characterized as vertical transverse isotropy with the low-velocity being the symmetry axis. The VDI is often interrupted around the boundaries of the 20 crustal units, usually above locally increased velocities in the lowermost crust. Due to the NW-SE shortening of the crust and the late-Variscan strike-slip movements along the NE-SW oriented sutures preserved in the BM lithosphere, the anisotropic fabric of the lower crust was partly or fully erased along the boundaries of original microplates. These weakened zones accompanied by a velocity increase above the Moho, which indicate an extrusion of mantle rocks into the lower crust, can represent channels through which portions of subducted and later molten rocks have percolated upwards providing magma to 25 subsequently form granitoid plutons. homogeneous high-resolution 3-D shear velocity model of the whole BM crust, infer Moho depth and map its lateral variations. high resolution of the velocity model data from temporary stations of the large-scale AlpArray passive experiment and several regional passive 55 seismic experiments organized the BM and We introduce a semi-automated procedure for ambient noise tomography, which consists of data gathering, waveform pre-processing, calculation of station-pair cross correlations, the Rayleigh wave extraction, frequency-time analysis (FTAN), discretization of the velocity dispersion curves of individual station pairs and cross-dimensional dispersion curve 60 tomography. The procedure leads to the homogeneous 3-D velocity model constructed from 1-D velocity-depth functions at a dense grid. The stochastic inversion with the use of the Improved Neighbourhood Algorithm (Wathelet, 2008) provides several additional attributes, e.g., misfit, skewness, standard deviation, as well as estimates of uncertainty, allowing us to P velocities of 5.6–6.4 km s -1 and the middle crust with intermediate composition and velocities 6.4–6.8 km s -1 . The lower part of the crust (LPC) denotes the mafic lower crust with 315 velocities 6.8–7.2 km s -1 and the high-velocity lowermost layer with v P higher than 7.2 km s -1 . observed P-velocity lower shear velocities in the LPC in the lower part of the UPC. (VDI) the UPC 320 and LPC rheological the modern view of the Conrad discontinuity

Early Controlled Source Seismic (CSS) research of the BM crust in the 1960s found the thickest crust of about 40 km in the 35 MD unit and the thinnest crust (~27 km) in the western part of the BM, in the ST unit (e.g., Beránek et al., 1975). The initial research was followed by many refraction and reflection seismic experiments which provided additional data for modelling the BM crust. Karousová et al. (2012 and references therein) interpolated the compiled findings of wide-angle refraction and near-vertical reflection seismic profiles into the first three-dimensional P-wave velocity model of the BM. The model maintains the gross features of the individual earlier models of the BM crust. 40 For a long time, the active seismic experiments have represented the only source of data for studies of the BM crust due to the lack of natural earthquakes needed for local earthquake tomography of the massif. Complementary approaches, particularly the receiver function (RF) method (e.g., Hetényi et al., 2018b) and the more recently introduced ambient noise tomography (ANT) have been proved to be useful tools for imaging the crust. The first one analyses delay times of 45 converted phases from teleseismic earthquakes and the second is based on cross-correlation technique (e.g., Bensen et al., 2007). The ANT has the ability to model the Moho and to retrieve shear wave velocities across the Earth crust and of the uppermost mantle when the network aperture is sufficiently large. One of the main advantages of the ambient noise tomography is in its ability to exploit long-term series of seismic ambient noise recorded at seismic stations regardless of earthquakes occurrence and their spatial distribution. Resolution of velocity models depends prevailingly on station density, 50 inter-station distance and azimuth variability and the overall path coverage of the studied region.
In this paper, we present the first homogeneous high-resolution 3-D shear velocity model of the whole BM crust, infer Moho depth and map its lateral variations. The unprecedented high resolution of the velocity model is achieved thanks to data from temporary stations of the large-scale AlpArray passive experiment (Hetényi et al., 2018a) and several regional passive 55 seismic experiments organized in the BM and its surroundings during the last two decades (see below).
We introduce a semi-automated procedure for ambient noise tomography, which consists of data gathering, waveform preprocessing, calculation of station-pair cross correlations, the Rayleigh wave extraction, frequency-time analysis (FTAN), discretization of the velocity dispersion curves of individual station pairs and cross-dimensional dispersion curve 60 tomography. The procedure leads to the homogeneous 3-D velocity model constructed from 1-D velocity-depth functions at a dense grid. The stochastic inversion with the use of the Improved Neighbourhood Algorithm (Wathelet, 2008) provides several additional attributes, e.g., misfit, skewness, standard deviation, as well as estimates of uncertainty, allowing us to https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License. assess reliability of the resulting velocity model. We compare the new 3-D model with several 2-D CSS profiles published by other authors and with estimates of the Moho depths from RF, resolved within this study. Velocity decrease modelled in 65 the lower crust is explained by the vertical transverse isotropic (VTI) fabric with the ‗slow' symmetry axis.

Data
To evaluate the shear velocities within the BM crust, we processed vertical components of daily recordings with sampling rate from 20 to 200 Hz from 410 stations located in the massif and its broader surroundings framed at 7. 1°-21.4° E and 46.2°-56.3° N and operated between 2002 and 2016. We used recordings from permanent observatories, temporary stations 70 of the AlpArray passive seismic experiment (AlpArray Seismic Network, 2015), its target-oriented complementary field experiment EASI (Hetényi et al., 2018b;network doi: 10.12686/alparray/xt_2014), as well as recordings from previous temporary measurements in the BM, e.g., BOHEMA I-IV (Plomerová et al., 2007;; Karousová et al., 2013), PASSEQ (Wilde-Piorko et al., 2008), and EgerRift (Babuška and Plomerová, 2013) (Fig. 2). Thorough analyses of noise sources showed that 90-day-long interval in the summer season (from June to August) is the optimal time period for stacking cross-75 correlation functions in the BM (Fig. 3). Strong microseismic energy arriving from the North Atlantic Ocean, particularly in the winter, strongly disturb the ambient noise wave field in central Europe and thus corrupt the fundamental assumption of the method to have the isotropic distribution of ambient noise sources. The summer season selection of data ensures similarity of the causal and acausal parts of the CCF (Fig. 3c) and the isotropy of the signal-to-noise ratio (Fig. 3d).

80
For our study, we created 83 845 station pairs from all of the available data and kept only those which satisfy conditions as follows: inter-station distance is between 50 km and 600 km, station-pair midpoint is inside the area framed at 11.3°-19.3° E and 47.5°-52.4° N, and station pairs contain at least 61 days of simultaneous recordings. With these selection criteria, we keep 24 258 good quality station-pair data suitable for ambient noise processing.

Methods 85
Theoretical framework for calculation of surface wave velocity dispersion curves from ambient noise was established in many experimental (e.g., Ren et al, 2013;Kästle et al., 2018;Lu et al., 2018;Schippkus et al., 2018;Soergel et al., 2020;and Guerin et al., 2020) and theoretical studies (e.g., Dunkin, 1965;Lobkis and Weaver, 2001;Levshin and Ritzwoller, 2001;Sambridge, 1999;and Yoshizawa and Kennet, 2005). The studies are based on cross correlations of long sequences of broadband seismic recordings and subsequent inversions of dispersion-curves for shear wave velocity model. Among others, 90 Lobkis and Weaver (2001) and Shapiro and Campillo (2004) showed that the cross correlations represent estimates of Green's functions for the crust and uppermost mantle along paths connecting two arbitrary sites. The reliable Green's functions can be retrieved, if the medium is homogeneous and isotropic and if the ambient noise wave field is diffuse or it https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License. arrives from uniformly distributed noise sources. However in seismology, none of these assumptions is fully satisfied.
Therefore an additional processing is needed for modelling structure of the crust. Bensen et al. (2007) describe data 95 processing that delivers reliable surface-wave dispersion measurements, widely used in many recent studies mentioned above, even if the above conditions are not fully satisfied. We have implemented the approach in building the first detailed 3-D velocity model of the BM. Thanks to significant dispersion of surface waves and to their frequency-dependent depth penetration, the surface wave tomographic inversion can be solved in computationally efficient two-step approach. The twostep method integrates an iterative fast marching method resulting in shear velocity maps per frequency bands and a 100 stochastic inversion of dispersion curves for a collection of multi-layered shear velocity models. The stochastic inversion addresses a non-uniqueness of Dunkin's (1965) formulae by a search for the best fit between the modelled and measured dispersion curves.
We build the 3-D shear wave velocity model in three semi-automated phases (Fig. 4) by joining three existing software 105 packages. The MSNoise package of python codes (Lecocq et al., 2014) is applied for data pre-processing, calculation of cross correlations and their stacking in Phase 1. The package is complemented by a new set of Python codes for frequency-time analysis (FTAN) of Levshin and Ritzwoller (2001). This step includes careful velocity-frequency picking of the group velocity dispersion curves. The FMST package of the Fortran codes (Rawlinson, 2005) is used to solve the 2-D surface wave travel time tomography in Phase 2. It employs eikonal equation forward modelling and the iterative subspace 110 method to minimise RMS error between measured and computed travel times (Rawlinson and Sambridge, 2005). The nonlinear stochastic inversion of dispersion curves (Wathelet, 2008) is solved in Phase 3 with the use of the Dinver software package (Wathelet et al., 2020).It employs the Neighbourhood Algorithm (NA) as an improved inversion method to search for the best fitting models (Sambridge, 1999;Wathelet et al., 2004;Wathelet, 2008).

115
In Phase 1, we demean daily seismic records, resample them to 0.1 s interval with anti-aliasing filter, remove instrument responses, normalise amplitudes and apply spectral whitening. Then we compute daily cross-correlation functions (CCFs) for all station pairs by averaging cross-correlations computed for one-hour time windows with 30 minutes overlaps. We extract fundamental mode of the Rayleigh waveforms and apply commonly used frequency-time analysis (FTAN). The FTAN diagrams (Fig 5a) display amplitudes of waveform envelopes as a function of the group velocity and the period of 120 narrow band-pass filters. On the FTAN diagrams we measure group velocities of the Rayleigh waves. The group velocity dispersion curves are discretized efficiently by an automated picking with non-linear one-third-octave band segmentation of the frequency axis in the FTAN diagrams.
The automated picker starts from the third sample before the long-period cut-off point and continues in direction towards 125 both the short and long periods. In the short-period band, the picker prioritizes the low velocity local maxima. This prevents picking the higher modes, as well as skipping and mixing different phases at the shot-period range. The cut-off period is a https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License. function of interstation distance and average group velocity, and relates to the separation of the causal and acausal parts of the Rayleigh wave on the CCF. By visual inspection of the FTAN diagrams, we reject 13 % of station pairs and end up with 21 066 dispersion curve measurements providing velocity-period pairs in a period range from 2 to 74 s (Fig 5b and  In Phase 2, we split whole set of station-pair dispersion curves and cluster the velocities according to the common periods. We convert measured group velocities into travel times and solve the travel time tomography by the Fast Marching Method (Rawlinson, 2005). Residual travel time RMS errors converge to optimal velocity model after 6 iterations. The dense inter-135 station spacing allowed us to construct velocity maps at a grid of 0.1° E-W by 0.1° N-S. In the initial model, we use constant velocity of 3.5 km s -1 and constant radius smoothing in order to avoid introduction of artificial anomalies. Due to the wellknown lower resolution at longer periods we increase the size of the smoothing factor proportionally to each period. We construct Rayleigh wave group velocity maps for period range from 2 to 74 s in an area framed at 7°-15° E-W and 47°-54° N-S. 140 In Phase 3, we assemble the velocity slices computed in Phase 2 and construct dispersion curves at a regular grid of 0.3° E-W by 0.2° N-S. Then we invert the regularised dispersion curves for local 1-D layered shear velocity models with the use of NA in the stochastic inversion program. The reliability of the NA solution depends mainly on the number of randomly generated starting models, number of iterations and the model space. These parameters need to be carefully tested to 145 guarantee the complete sampling of the solution space. The overall quality of the inversion results is assessed by misfit values, which indicate how each model fits the input dispersion curve by the normalised RMS error. We tune the NA parameters by evaluation of the misfit values, the skewness and standard deviations of the group of best-fitting models with the aim of achieving the Gaussian distribution of group of best-fitting solutions without a bias of the NA search constraints.

150
The NA models consist of six-layers over a half-space (Table 1). The top three layers encompass heterogeneous shallowest structures and the three deeper layers represent the upper part and lower part of the crust and the uppermost mantle. The NA searches for the best fitting model in 360 iterations starting from 280 randomly generated initial models. To get a representative model at each grid node, we average 10 % of the best-fitting models from the resulting set of 400 000 models.
We apply two approaches of averaging. First, the conventional one, which averages the v S over the fixed depths, results in a 155 gradient-velocity model (Fig. 6, left, blue curve). The second approach averages over layer indexes (v i, z i ) and results in a layered-velocity model (Fig 6a, left, red curves). This approach allows us to retrieve the velocity boundaries to which the surface waves are insensitive in principle. The dispersion curve, calculated from the layered-velocity model, matches better the input dispersion curve (Fig. 6a, right). The averaging over layer indexes leads to the straightforward implementation of a layer-stripping technique into the stochastic inversion. The layer stripping appears to be efficient tool to enhance the NA search in deeper layers of the model, where the solution suffers from the significant decrease of sensitivity kernels with depth (Yoshizawa and Kennet, 2005). The recursive stripping of the upper layers increases depth sensitivity and leads to better fit at long period range of the dispersion curves (Fig 6b). Figure 6e reveals the linear improvement of average misfits for two selected period intervals, after each of 165 the three steps of the subsequent layer stripping. The misfit reduction is feasible and consistent for all nodes. Similar improvement concerns setting the inversion parameter that allows or does not allow the velocity decrease in the lower part of the crust. Figure 6b,c shows that the models, in which velocity decrease with depth is allowed, fit better the dispersion curve with a local minimum around the 20 s periods. In our data set, a large majority of the resulting models consistently prefer the velocity decrease in the lower crust (see also Fig. 10). 170 We assess the sensitivity of our inversion approach to the sharpness of the Moho as an interface by adding an additional gradient layer into the final 1-D model (Fig.6d). In order to mimic the gradational Moho we insert the layer on top of the uppermost mantle layer (Table 1) and varied its thickness from 2 to 10 km. The test shows that the larger is the thickness of the transitional layer; the worse is the misfit in the dispersion curves, which is detectable from layer thickness larger than 2 175 km. This means that in case of the gradational Moho instead of a sharp discontinuity, the thickness of the transitional layer does not exceed 2 km.

Results
The presented high-resolution model shows shear wave velocities ~0.2 km s -1 higher in almost the entire BM upper crust, on average, than those in its surroundings, with the highest velocities in the southern part of the Moldanubian (MD) unit (Fig 7). 180 The higher velocities clearly contour the BM shape. Very low velocities down to 10 km depth mark the SE rim of the BM beneath the Carpathian Foredeep. At greater depths (15 km), the relatively higher velocities are found in the MD unit.
Locally increased velocities occur in the lower crust at several places, e.g. beneath the Eger Rift at 25 km depth. The upper mantle velocities arise first at 28 km depth slice beneath the Eger Rift (ER), marking thus the massif thinning along the rift. Regardless of the shallowest, not well resolved upper 3 km, the model comprises two well-resolved layers, which we relate to the upper and lower parts of the crust (UPC and LPC in Fig. 8b,c,). In the UPC the average velocities vary from ~3.2 km 195 s -1 in the north and north-west of the massif up to more than 3.5 km s -1 in the southernmost part of the BM. With two small-size exceptions (along BOG and Lugian/Moravo-Silesian contact), the LPC of the BM can be divided into the south-western and north-eastern half, exhibiting relatively high and low average velocities, respectively.
Individual 1-D velocity-depth models resulting from the stochastic inversions and their corresponding dispersion curves 200 exhibit distinct regional variations of the velocity interface between the LPC and UPC, which are characteristic for tectonic units of the BM (Fig. 9a,b). Comparison of 1-D models with the radial Earth model IASP'91 (Kennet and Engdahl, 1991) shows that the shear velocities retrieved in the upper part of the crust are mostly similar, while in the lower part of the crust the ANT model velocities are significantly lower.

205
Most of the massif is characterized by a velocity decrease in the lower crust, which can be quantified as dv S = v S (LPC 1 )v S (UPC 7 ), where indexes denote respective sub-layers (Table 1). The velocity-drop interface (VDI) with negative dv S occurs at variable depths and with different size of the velocity drop (Fig. 10). Within the MD, the interface with velocity drop exceeding -0.1 km s -1 lies at the 26-30 km depth, being the deepest in the south, though the average velocity decrease is not the largest one (14.7° E 49.5° N, depth=28 km, dv S =-0.18 km s -1 ). The MD is surrounded by a band, in which the VDI uplifts 210 into shallower depths. In the southwest continuation of the MD unit out of the BM(see Fig1), the size of the velocity decrease become even stronger locally, but the depth of the interface shallows by ~10km (11.7° E 48.7° N, depth=19 km, dv S =-0.37 km s -1 ). The amplitude of the negative dv S slightly increases also towards the TB unit, but its depth remains deeper than that in the westward continuation of the MD unit out of the BM (13.5° E 50.1° N, depth=23 km, dv S =-0.27 km s -1 ). The SW-NE trend of the Variscan structures, including the ER, is reflected in a separation of the ST, characterized by the 215 weaker velocity decrease across the VDI at shallower depths (12.0° E 50.7° N, depth ~18 km, dv S =-0.21 km s -1 ) in comparison with the rest of the BM. In the south-western part of the region (11.4° E 49.3° N), where the ST unit neighbours the crust of the MD zone (out of the BM), the VDI is absent. The narrow zone without the velocity drop continues north-eastward along the ER (ST/TB boundary). Further to NE, the ER continuation approximately divides the regions with the shallow and deep VDI in the West Sudetes and Lugian Domain, respectively. The VDI is a little bit more expressed 220 around the Elbe Fault Zone (EFZ) and lies deeper to the north-east of the EFZ in comparison with the ST south-west of the EFZ, but it remains shallower in comparison with the MD. Further to the east in the northern BM, the VDI deepens to ~25 km, in the northern Lugian Domain and shallows again to ~20 km in the Silesian (Fig. 10).
Besides the tiny, but systematic differences between average velocities in the upper and lower crust (see Fig. 7b,c,), there are 225 also differences in the velocity gradients in the individual regions of the BM (Fig. 10). The EFZ divides the BM https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.
approximately into two zones. The southern zone of the BM, which includes the MD and TB units, the eastern ST unit and southern part of the Lugian Domain, is characterized by higher velocity gradients below the VDI than that in the layer above the VDI. In the north-eastern rim of the BM, particularly north of the ISF, in the East Sudetes and Silesian Domain, the gradient difference at VDI is negative, i.e., the gradient in the layer above the VDI is higher than that below the interface. 230 The differences between velocity gradients in the UPC and LPC clearly delineate zones of the BM that match with the BM zonation according to the size and depth of the velocity drop at the VDI (Fig. 10).
Distinct lateral variations of the v S velocities, the VDI and the Moho depths, as well as thickness of the LPC are shown along several cross sections in the 3-D view (Fig.11, for more views see S1-S4). The sharp dark blue-light green velocity contrast 235 marks the Moho derived from the ambient noise tomography. The figures show crustal thinning, e.g., beneath the ER (see also Fig.15) and Moho undulations, e.g., around crossing point of Profile 7 and Profile 12 beneath the southern SI (eastern continuation of the EFZ). The shallow Moho beneath the ST unit is also visible along the Profile 14 (south-western part in view from SW).

240
The Ps receiver functions represent an independent method to estimate the Moho depth from teleseismic P-to-S converted phases at the bottom of the crust. We picked the Ps-P delay times for all the stations and calculated the Moho depth with the use of the IASP91 velocity model (for more details see Hetényi et al., 2018b). We interpolate the Moho depths estimates

255
Mostly in the northern and eastern parts of the BM the Moho depths differ more than would correspond to the estimated accuracy of +/-2 km (e.g., Fig. 11, Profiles10, 14 and 15, and S1-4). Exceptional large discrepancies occur in regions southeast of the BM with uncompensated sediments in Moho RF calculations. Lower average velocities in the ANT model shallow the Moho depth estimates from the RF by 1 km, relative to that if the IASP'91 velocities are used, on average.
The most distinct phenomenon of the presented model is the detection of the velocity-drop interface (VDI) in the lower crust with the velocity decrease exceeding -0.1 km s -1 (marked by yellow lines in Fig. 11). Variations of both the interface sharpness (size of the velocity drop) and its depth change laterally. These changes can be related to individual tectonic units of the BM crust (see Figs. 9 and 10). On the other hand, there is no direct correlation between the Moho, or the VDI depths, and the lower crust thickness. Though the thin and thick lower parts of the crust beneath the VDI are visible above the deep 265 and shallow Moho, respectively (e.g. P16, view from the W), the reversed relation, i.e., the thin and thick LPC above the shallower and deeper Moho, respectively, can be seen as well (e.g., P10, view from the W, or P14 from the E). This may reflect the fact that the accuracy of determining both the Moho and VDI depths is not sufficient to quantify the subtle relationship between the thickness of LPC and Moho depth.

Resolution and Sensitivity Tests 270
We use standard synthetic checkerboard tests to assess the resolution of the fast marching tomographic inversions. We search for spatial resolution on checkerboard patterns with +/-3 % velocity perturbation relative to the initial velocity model (3.5 km s -1 ) and with three cell sizes (Fig. 12a). As expected, the spatial resolution relates to the ray-path coverage for each period. In the central and southern parts of the model, we retrieve velocity heterogeneities at size of 44 km by 44 km using a ray-path coverage at period 6 s ( Fig. 12a,b left), 66 km by 66 km at period 19 s (Fig. 12b central) and 88 km by 88 km at 275 period 48 s (Fig. 12b right). The resolution decreases towards the northern part of the model and towards the model margins with sparse ray-path coverage. The ray-path coverage of the region is shown in Fig. 12c by the ray-path hitcount per 0.3° E-W by 0.2° N-S cell size.
Further, we verify reliability of the fast marching inversion at different periods based on travel time residuals computed for 280 all input travel paths through the updated 2-D velocity slices. The travel time residuals document that the FMST inversion returns reliable solution in a period range from 6 to 48 s (Fig. 13b). The computed travel time residuals at short periods (<5 s) exceed significantly the short wavelengths due to increased uncertainty of dispersion curve picking. Therefore, the fast marching inversion does not provide reliable results at short periods. The lower reliability of the fast marching inversion at long periods (60 and 74 s) is affected by reduced number of the travel time inputs into the inversion (Fig. 13a). 285 Finally, we also calculate the depth sensitivity kernels at each cell of the final velocity model and demonstrate their well-known dependence on period (Fig. 13b). The average sensitivity kernels, calculated across cells in the BM, show the highest sensitivity to structures in the upper part of the crust (~3-20 km) for periods between 6-19 s. The lower part of the crust down to Moho is sampled well with periods >19 s. Based on the resolution and sensitivity tests, we infer that the VDI -290 the new phenomenon retrieved in the presented model at depth ~15-30 km -lies in a high-fidelity area for both the sensitivity kernels and minimum travel time residuals. Results of the spatial resolution tests confirm that we can resolve https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.
structures at least as small as ~40 km in the upper part of the crust and ~80 km in the lower crust in the velocity model from the ANT.

Discussion 295
Traditionally, the upper and lower crust are distinguished in velocity models of the Earth crust, by higher velocities in the lower layer and the strong positive velocity contrast at the bottom of the crust (Moho). The mid-crust boundary, called Conrad discontinuity, was originally supposed to separate the granitic and mafic crustal layers. However, suggested Conrad discontinuity is missing in many seismological models with gradient velocities. In this paper, we show that an intra-crust boundary can appear with a negative velocity contrast in the lower crust. 300 The dispersion curves of group velocity computed by FMST in the BM, are characterized by steep slope of velocity increase in short period range <8 s. The curves have flat local maxima between 10 and 15 s periods, followed by a subtle velocity decrease with local minima at ~20 s and moderate velocity increase for periods larger than 20 s. General level of the local maxima of the dispersion curves attains group velocities of ~3.1 km s -1 . This value is lower than velocities from steeply-305 increasing dispersion curves of Lu et al. (2018) in the BM. But, the general shape and velocity level of our dispersion curves are comparable with those of Ren et al. (2013), Shippkus et al. (2020), or Qorbani et al. (2020). Also the ANT models of all these authors show higher v S in the BM upper crust than in its surroundings.
Wide ranges of 1-D layer thicknesses and velocities in our ANT input parameters (see Table 1) results in the v S model of the 310 BM crust with three characteristic depth zones. The uppermost part, thinner than 3 km, is not well resolved due to locally heterogeneous structures including sediments in some places. For the two main parts of the crust we use the same notation as introduced by Artemieva and Thybo (2013) in the European EUNA model derived from a compilation of the CSS models.
The upper part of the crust (UPC) includes felsic material with P velocities of 5.6-6.4 km s -1 and the middle crust with intermediate composition and velocities 6.4-6.8 km s -1 . The lower part of the crust (LPC) denotes the mafic lower crust with 315 velocities 6.8-7.2 km s -1 and the high-velocity lowermost layer with v P higher than 7.2 km s -1 .
However, contrary to the commonly observed P-velocity increase with depth, we retrieved lower shear velocities in the LPC than in the lower part of the UPC. No mid-crust layer or an interface with positive impedance (traditional Conrad discontinuity) is required in our gradient velocity model. The characteristic velocity-drop interface (VDI) between the UPC 320 and LPC can be associated with a thermodynamically controlled interface or rheological boundary representing the modern view of the Conrad discontinuity (Lay and Wallace, 1995).

BM crust images in ANT and CSS
Several CSS profiles criss-cross the BM (Fig. 14) and allow us to compare characteristics of the velocity cross-sections through the 3-D ANT v S model with the 2-D v P models from series of the CSS measurements. For that we superimpose the 325 v S velocity model with main features of the CSS profiles, especially with the Moho discontinuity and inter-crustal boundaries.
Two parallel CSS profiles Sudetes S04 (Hrubcová et al., 2010) and Celebration CEL09   The Moho depths from different data sets and the three independent techniques fit best along the N-S profile ALP01  and the AlpArray-EASI (Hetényi et al., 2018b) running through the BM (Fig. 14e). In general, discrepancies in the Moho depth determination of less than 5 km, which is the most frequent case for the BM, are considered as a good match (Kästle et al., 2018). We found a good fit also with the two-layer model of the crust beneath the KHC station (13.58N 360 49.13N), located near the AlpArray-EASI, ALP01 (~150 km) and 9HR CSS profiles (Ĉervený et al., 1977). The authors applied the spectral ratio method by Phinney (1964) on data of the first European stations providing the continuous broadband digital recordings since 1972 (Plešinger and Horálek, 1976). The preferred two-layer model, with constant velocities, sets the Moho at 35 km, which is in agreement both with CSS results and the ANT model presented here (Table 2). However, the upper/lower crust boundary is 5 km deeper in our model. While the velocity gradient in the upper 365 crust is constant in our model too, we found the 0.2 km s -1 v S decrease at the VDI and relatively steep positive velocity gradient in the lower crust. The velocity remains lower even at its bottom (v S =3.6 km s -1 ). Also the uppermost mantle velocities are lower in the ANT model than those derived from the dispersion curves of teleseismic Rayleigh waves. Waves modelled in Ĉervený et al. (1977) arrived from back-azimuths of ~30°, which is close to the high-velocity directions derived from body-wave propagation of regional earthquakes (Plomerová et al., 1984). On the other hand, the ANT accounts for 370 dispersion curves from all directions and thus results in lower velocities. An intra-crustal boundary in the Alp01 (Fig. 14d) follows the VDI interface similarly to that in the CEL10, where it coincides with the top of the positive high-v P -gradient transitional layer. In this case, the results from the trial and error technique of the 2-D modelling of the wide-angle CSS data differ from results obtained from the 3-D ANT, which model relatively low v S velocities in the LPC and the sharp Moho.

375
When comparing results from the CSS and 3-D ANT, one has to be bear in mind, that the 2-D CSS models are significantly affected by 2-D processing assumptions and narrow azimuth data acquisition, whilst the 3-D ANT model is based on 3-D homogeneous sampling and imaging the crust, but with relatively lower resolution in the lower crust. Therefore, the ANT maps 3-D structures of the crust in a more realistic way. Thanks to the layer-stripping method the ANT reveals the Moho as a sharp discontinuity beneath the whole BM. However, the mantle velocities beneath the Moho remain lower than expected. 380 Nevertheless, the lower velocities are in agreement with ANT of Ren et al.(2013), or results of teleseismic body-wave tomography (e.g., Amaru, 2007;Plomerová et al., 2016 and reference there in), as well as with full wave-form inversion (Fichtner and Villasenor, 2015). The shear velocities in the upper crust are in agreement with the v P CSS models considering vertical directions prevail. Therefore, we can explain the discrepancy between the CSS and our ANT model by anisotropic 390 structure (transverse isotropy with slow vertical symmetry axis) in the LPC (see the text below). Weak VTI anisotropy explaining the velocity drop in the lower crust is in agreement with results of other authors elsewhere (e.g., Luo et al., 2013;Shapiro et al., 2004;and Qiao et al., 2018). The subtle v S decrease in the top of the LPC of the ANT model (see Figs. 8 and 10) can be masked by uncertainties of the trial and error technique of the CSS models (Majdanski and Polkowski, 2014).

395
Some recent ANT models of the European crust (Lu et al., 2020), or the Alpine region (Kästle et al., 2018), represent themselves as candidates to a new reference model of the crust and uppermost mantle for the RF migration purposes. We show (see Section4 and S8) that in regions, where the v S ANT model detects low velocities in the LPC due to anisotropy, the average v P /v S ratio is not valid for the simple v S conversion to v P in the crust.

Dynamic development of the BM crust 400
We aim at improving the current knowledge of tectonic processes that formed the BM by linking the most significant

Traces of Variscan subductions
Crustal thickness varies between ~ 28 and 41 km in the ANT model of the BM. The relatively thin crust (<33 km) forms the north-and south-western part of the massif (see Fig. 8a). The region with the thin crust is bounded by the EFZ in the north 410 and by a 33 km depth band paralleling the ER in the east. The thinnest crust, where the Moho shallows up to 28 km, occurs near the intersection of the ER and its hidden north-eastern continuation east of the EFZ. The strike-slip EFZ is the boundary between the two unrelated parts of the ST unit (see Fig. 1), which obviously played an important role in the BM tectonics.
Transitions between the relatively thinner and thicker crust in the central and south-eastern BM spatially do not coincide 415 with the surface tectonic boundaries of the TB unit. The thinner crust in the north-western part of the TB unit is characteristic for the older (Cadomian) part of the BM. To explain the crust thickening beneath the south-eastern part of the TB unit, we suggest the thrusting of the MD crust and the whole MD lithosphere under the TB unit (Babuška and Plomerová, 2013). The age of the ST/TB collision is estimated at ~380 Ma (Franke et al., 2017). On the other hand, the thicker MD and BV units in the south-eastern part of the BM, were accreted to the TB block during the Variscan orogeny (e.g., Finger et al., 2007). Pitra 420 et al. (1999) has already pointed out that there is no conformity of the TB ductile structures with those of the MD crust. The https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.
Cambrian and Carboniferous tectono-metamorphic histories of the TB and MD, respectively, are separated by nearly 200 Ma. This age difference is reflected also in the difference between thicknesses of the crust.
Three discontinuous belts with thickness exceeding 38 km are mapped in the MD part of the BM (Fig. 8a). We explain the 425 crust thickening most probably by remnants of three Variscan paleo-subductions. Two of them are in the MD adjacent to the two plutons the CBPC and the MDB. The third one relates to the subduction of the Brunovistulian lithosphere in the south-eastern rim of the MD. These processes that enlarged the pre-Variscan assembly of the ST and TB units, started by westward subductions of smaller oceanic plates beneath the TB block (Konopásek et al., 2002;Medaris et al., 2006;Finger et al., 2007). Each of the oceanic subductions was followed by the continent-continent collisions. Thrusting of the two MD 430 domains beneath the TB block was followed by under-thrusting of the BV domain beneath the MD continental lithosphere (see also Fig. 6 in Babuška and Plomerová, 2013).
Similarly, the thick crust in the Lugian Domain (Fig. 8) can represent a remnant of a small microplate that Konopásek et al.  (Burov and Watts, 2006;Heron et al., 2016), the major stress guide in plate tectonic processes. In geologic time scale, the Moho is a dynamic, continuously modified feature. The images that we obtain are snapshots of the Moho depth at the present time, which preserve actions of tectonic processes involved in development of the crust over geological timescales (Carbonell et al., 2013).

Anisotropy of lower crust 460
The dominant feature of the presented 3-D v S model is the prevailingly low velocity in the lower crust (see Figs. 7,9,10 or 5bc), relative to the upper crust velocity. This is in contradiction with the high v P in the lower BM crust modelled by the CSS profiles. Velocities v P as high as 6.90-7.50 km s -1 were modelled at depths 25-35 km in the ST crust and 6.60-6.95 km s -1 in the TB and MD units (Profile CEL 09, Hrubcová et al, 2005). Such high v P velocities are close or even higher than those interpreted in the lower crust of cratons (Mooney et al., 2002). Thybo et al. (2013) note that the top of the high P velocity 465 lower crust might be mistakenly interpreted as the Moho, because (1) the high velocity lower crust may be a -hidden layer‖ for refraction seismic interpretations, (2) it may cause the strongest observed wide-angle reflection in seismic sections, or (3) it may be the strongest convertor in receiver function images. We explain the difference between the high v P in the BM lower crust derived from the CSS and the velocity drop mapped by the v S velocities derived from the Rayleigh waves from ambient noise by anisotropic structure of the lower crust. Refracted and/or wide-angle reflected P waves in the CSS sample 470 prevailingly high velocities, while the Rayleigh waves sample the low velocities in the anisotropic medium-transversely isotropic with vertical ‗low-velocity' symmetry axis.
The anisotropy of the lower crust plays an important role in velocity modelling and may explain some inconsistencies in isotropic models of the crust and namely why the v P derived from the CSS profiles is high and on the other hand, the v S from 475 the Rayleigh waves is low in the same crustal layer. Several v S models incorporating the AlpArray data (Lu et al., 2018(Lu et al., , 2020Kästle et al., 2018;Schippkus et al., 2018;Guerin et al., 2020) have been derived from the ambient noise recorded at European stations. However, their resolution in the lower crust is often low outside of the Alpine region due to sparse station coverage in their ANT studies. Nevertheless, no velocity decrease in the lower crust that we observe in the BM, is mapped either in their well resolved regions. On the other hand, a low-velocity zone in the mid-lower crust, representing a 480 mechanically weak zone channelling material transport, was observed in the SE Tibet (Qiao et al., 2018). Also Luo et al. (2013)  considered as average velocity of isotropic rocks or non-oriented material at the corresponding depth. Then the velocity decrease at the top of the lower crust can be estimated as the velocity difference in the vertical and horizontal directions of 490 the anisotropic lower crust, approximated by the sub-horizontally oriented high velocities of the VTI. Then the minimum strength of the lower crust anisotropy ranges from 3 to 9 % depending on the average v S .
We have looked for rocks that could be the most suitable candidates for composition of the lower crust, namely as to velocities and fabrics expressed in this type of seismic anisotropy. The lower crust represented by the Neukirchen-Kdyně 495 Massif (NK in Fig. 1 ) in the south-western tip of the TB unit, was thrust over the MD unit during their Upper Devonian collision (Buess et al., 2002). The massif, where amphibolite-facies rocks are exhumed (Cháb and Ţáĉek, 1994), can serve as a model for a composition of the lower crust of the BM. Amphibolite, consisting mainly of hornblende and plagioclase, is a most frequent candidate that corresponds best to interpretations in papers focused on seismic and gravity modelling of the continental lower crust. Mafic granulite and gabbro are also considered in literature as potential components. However, both 500 rocks have typically very low v P and v S anisotropy measured in laboratory at high hydrostatic pressure (Almqvist and Mainprice, 2017). Ji et al. (2013) measured v P and v S velocities and anisotropy of 17 amphibole-rich rock samples at hydrostatic pressures up to 650 MPa corresponding to depths of about 20 km. Preferred orientation of hornblende plays decisive role in the formation of 505 seismic anisotropy of amphibolite, whereas other minerals as plagioclase, pyroxene or garnet diminish the anisotropy. The authors found large variability of v P velocities (6.4-7.5 km s -1 ) and concentration of very high velocities along foliation planes (7.3-7.5 km s -1 ). Similar high velocities were locally interpreted on the CSS profiles, e.g., within the lowermost crust of the ST unit (Profile CEL 09, Hrubcová et al., 2005). Shear wave velocities v S measured by Ji et al. (2013) are also very high (3.7-4.25 km s -1 ), the lowest velocities (3.7-3.8 km) being for waves propagating perpendicularly to the foliation of 510 amphibolite samples. Such velocities correspond to the low-velocity propagations of Rayleigh waves in the vertically transverse isotropic BM lower crust, considering Poisson solid and relating phase v S of body and surface wave velocities.
It is generally known that velocities measured in laboratory on fresh rock samples are usually higher than velocities measured in large rock massifs that are heterogeneous and dissected by faults. The lower crust of the NK, taken as a model, 515 consists also of other rocks (meta-sediments, gabbro, see, e.g., Cháb and Ţáĉek, 1994) that can decrease both the v S and anisotropy. On the other hand, an effective anisotropy that is reflected in velocities of Rayleigh waves can be enhanced by metamorphic and deformational processes that produce fabrics such as schistosity and can form extensive regions of foliated and laminated rocks (Okaya et al., 2018). Well documented example of a laminated lowermost crust provides the refraction and wide-angle reflection experiment along profile ALP 01 , parallel to the EASI profile (see Fig. 14c, 520 80-280 km). Several kilometre-thick laminated transition zone above the Moho is interpreted within the BM northward of the contact with the Alps. https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.
Assuming a mafic composition (amphibolite-facies metamorphic rocks), a ductile flow in the lower crust is probably controlled by the deformation of plagioclase feldspar. Tendency of the lower crust to flow (e.g., Hacker et al., 2015) is 525 controlled by its rheology. The onset of plagioclase viscous behaviour corresponds to temperature of ~450°C (Cole et al., 2007). Based on heat-flow measurements in western Bohemia (Šafanda et al., 2007), such temperature can be reached at depths around 20 km. Of course, the present-day depth of the top of the ductile lower crust, which shallows up to ~15 km at some places, was modified by denudation that is estimated at about 7 km in the western BM (de Wall et al., 2019). On the other hand, the top of the lower crust deepens to ~28 km (Fig. 10) in the south-eastern part of the BM core, which probably 530 reflects thickening of the crust by Variscan (first oceanic, later continental) subductions (Medaris et al., 2006;Finger et al., 2007). Consequently, a thermal gradient was lowered due to subductions of cold materials and by the time elapsed since these events took place, and thus the onset of ductile flow might have been deeper in the geological past Decoupling between the upper and lower crust due to rheological contrast between the two thermally different levels was 535 suggested in the French Pyrenees (Denele et al., 2009). Similar model of mechanical decoupling between the brittle upper crust and a -hot deeper crustal level‖ is used to explain the crustal shortening during the India-Asia collision (Van Buer et al., 2015). Many seismic and magneto-telluric experiments within Tibet provide data that require models characterised by a less viscous mid-lower crust in-between more viscous upper crust and mantle lithosphere (e.g., Klemperer, 2006). Such interpretations agree with the widely accepted -jelly sandwich‖ model of Burov and Watts (2006) that can explain best the 540 gross structural styles of plate tectonic collisional systems and also a large-scale flexure of the India Plate during the collision (e.g., Hetényi et al., 2006).

Regional distribution of VDI and late Variscan strike-slip tectonics
Most of the BM crust and its southern and western surroundings are characterized by the abrupt shear velocity decrease in the lower crust exceeding -0.4 km s -1 at some places (see Fig. 10). Though the amplitude of the velocity drop and depth of 545 the drop can be clustered according to individual tectonic units, there is only a weak relation between the depth of the VDI or the thickness of the LPC layer and the crustal thickness. The shallowest VDI (green in Fig. 10) is mostly mapped in regions with thinner crust (Fig. 8). On the other hand, the most significant velocity drop appears both within regions with the thin and thick crust.

550
The VDI horizon, marked at sites where the velocity difference attains at least -0.1 km s -1 , is interrupted at several zones, where the negative v S drop does not reach such value or is absent. Evidently, such zones appear also as regions with no difference between the upper and lower crust velocities (Fig. 10) and often with a local velocity increase in the lowermost crust. Figure 15 depicts a NW-SE cross section through the core of the BM, showing the v S in the crust of the four major tectonic units, former microplates assembled during the Variscan orogeny. The horizon of the VDI is interrupted above mantle sutures and near the boundaries of crustal units, often laterally shifted from their mantle counterparts. The sutures probably served as subduction channels that later became paths, along which subducted rocks were brought to the upper crust and provided magma for large granitoid massifs like the Central Bohemian Plutonic Complex or the Moldanubian Batholith 560 (Babuška and Plomerová, 2013). The interruption of the VDI and relatively increased velocities in the lower crust correlate with the deep sutures also on other cross sections (Fig. 14b-c). No VDI interruption within the BM occurs along the S04 (Fig.14a) which runs parallel to the EFZ and does not cross any distinct suture. This feature supports the importance of the late-Variscan strike-slip tectonics, in which the Elbe shear system played important role (Elter et al., 2020).

565
Contraction of the eastern Variscides caused by the north-westward advance of Gondwana resulted in a collisional shortening of the BM (Franke et al., 2017). The contraction was compensated by Late-Variscan strike-slip zones that considerably modified an original assembly of the BM (Pitra et al., 1999). Field measurements by these authors document dextral strike-slip movements along the CBSZ. This is in agreement with the model of the westward movement of a part of the TB lithosphere along the WSW oriented sutures on both sides of the TB unit (Babuška and Plomerová, 2017). 570 The model we present in Fig. 16 can explain the spatial distribution of VDI, why it is interrupted at places where the velocity decrease is weaker than -0.1 km s -1 or is absent. Most of the interruptions correlate with the deep boundaries of the tectonic units, e.g., ST/TB, TB/MD, and MSD/BV (Babuška and Plomerová, 2013). Late Variscan strike-slip movements preferably took place along such boundaries and most likely modified the sub-horizontal fabric of the lower crust by a system of 575 parallel, deep-reaching faults. Strike-slip zones, like the CBSZ or the suture between the ST and TB units (see Fig. 15), undoubtedly cut the whole crust by sub-vertical faults that continue towards the mantle lithosphere. Movements along these zones disrupted the horizontal fabric and the VTI in the lower crust. Fabrics around the block boundaries might be modified locally either into a more complex symmetry, e.g., the orthorhombic symmetry with ‖fast‖ horizontal axis oriented in direction of movement, or might be entirely overwritten. 580

Conclusions
We present the first high-resolution 3-D shear velocity model of the BM crust and a map of Moho depths retrieved from ambient noise tomography. The model is based on recordings of more than 400 stations operated in several passive seismic experiments organised during the last two decades, including the large-scale AlpArray network (Hetényi et al., 2018a). The shear velocities modelled from dispersion curves of Rayleigh waves are ~0.2 km s -1 higher in the BM upper crust than in its 585 surroundings with the highest velocities in the Moldanubian part. The region north of the Elbe Fault Zone exhibits relatively lower velocities. The presented 3-D model provides compelling evidence of a regional-scale velocity distribution. In general, the Cadomian part of the region has a thinner crust, while the crust assembled, or tectonically transformed in the Variscan period, is thicker. The Moho discontinuity, as the dynamic phenomenon, retains remnants of Variscan subductions imprinted in bands of local crustal thickenings.

595
The dominant feature of the presented 3-D velocity model is relatively low v S in the lower crust that is in contradiction with the high v P modelled in the BM by the Control Source Seismic profiles. We explain this seeming disagreement by anisotropic structure of the lower crust characterized by a transversely isotropic medium with vertical ‗low-velocity' symmetry axis. We estimate strength of anisotropy at 3-9 % in dependence on the average v S . Considering the observed anisotropy and both v P and v S interpreted from different field seismic experiments, amphibolite is the best candidate for the 600 composition of the BM lower crust.
Thanks to the high density of data, thorough processing, enhanced depth sensitivity of trans-dimensional tomographic inversion by layer stripping, layer averaging and allowing for a velocity decrease with depth during the stochastic inversion, we were able to retrieve, for the first time, a distinct intra-crustal interface (VDI) with a velocity decrease as strong as to -0.4 605 km s -1 in the lower part of the crust at depths 18-30 km. Particularly the averaging of the best-fitting models over layers instead of the commonly used averaging over fixed depth, proved to be crucial for the detection of velocity interfaces.
The VDI is interrupted around boundaries of the crustal units, usually above locally increased velocities in the lowermost crust. Boundaries of the tectonic units likely served as subduction channels and partly as strike-slip zones in late-Variscan 610 tectonic processes. Such zones can represent channels through which the high-velocity mantle material intruded into the                https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.       https://doi.org/10.5194/se-2020-176 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.