Radial anisotropy and S-wave velocity depict the internal to external zones transition within the Variscan orogen (NW Iberia)

. The cross-correlation of ambient noise records registered by seismic networks has proven to be a valuable tool to obtain new insights into the crustal structure at different scales. Based on 2-to 14-s-period Rayleigh and Love dispersion data extracted from the seismic ambient noise recorded by 20 three-component broadband stations belonging to two different 10 temporary experiments, we present the first i) upper crustal (1-12 km) high-resolution shear wave velocity and ii) radial anisotropy variation models of the continental crust in NW Iberia. The area of study represents one of the best exposed cross-sections along the Variscan orogen of western Europe, showing the transition between the external eastern zones towards the internal areas in the west. Both the 2-D maps and an E-W transect reveal a close correspondence with the main geological domains of the Variscan orogen. The foreland-fold and thrust-belt of the orogen, the Cantabrian Zone, is revealed 15 by a zone of relatively low shear wave velocities (2.3-3.0 km/s), while the internal zones generally display higher homogeneous velocities (> 3.1km/s). The boundary between both zones is clearly delineated in the models, depicting the arcuate shape of the orogenic belt. The velocity patterns also reveal variations of the elastic properties of the upper crust that can be linked to major Variscan structures, such as the basal detachment of the Cantabrian Zone or the stack of nappes involving pre-Variscan basement; or sedimentary features such as the presence of thick syn-orogenic siliciclastic wedges. 20 Overall, the radial anisotropy magnitude varies between -5 and 15 % and increases with depth. The depth pattern suggests that the alignment of cracks is the main source of anisotropy at < 8 km depths, although the intrinsic anisotropy seems to be significant in the West-Asturian Leonese Zone, the low-grade slate belt adjacent to the Cantabrian Zone. At depths > 8 km, widespread high and positive radial anisotropies are observed, which we attribute to the presence of subhorizontal alignments of grains and minerals in relation to the pre-or syn-orogenic deformation associated with the Variscan 25 orogenesis.

refraction/wide angle surveys were carried out in this area. These programs focused in the unveiling of the deep crustal structure of the CM (Pérez-Estaún et al., 1994;Pulgar et al., 1995Pulgar et al., , 1996Álvarez-Marrón et al., 1996Álvarez-Marrón et al., , 1997Gallastegui et al., 1997Gallastegui et al., , 2002Ayarza et al., 1998;Fernández-Viejo et al., 1998, 2000Fernández-Viejo and Gallastegui, 2005;Fernández-Viejo et al., 2011 and broadened the knowledge about the multiorogenic crustal architecture. To improve imaging resolution at crustal depths, several temporary seismic arrays with different specifications and dimensions were deployed in 70 the CM (e.g. López-Fernández et al., 2012;, the most recent of which is featured for the first time in this study.
Similarly to other passive seismic imaging techniques, the ANI reflects the variation of the seismic velocities of the bulk rock. In orogenic belts, where lithologies with different properties are placed in contact by structural features, ANI can also help to infer indirectly the presence of structures at depth and thus it can be valuable for unravelling the upper crustal structure at regional or even local scale in a cost-effective manner (e.g Zigone et al., 2015;Sammarco et al., 2017;Acevedo 75 et al., 2019;Gu et al., 2019;Olivar-Castaño et al., 2020). In any case, the main advantage of ANI is that, unlike earthquakebased tomography, it can be applied in seismically quiet regions due to the ubiquitous presence of noise sources. Up to the date, there is limited knowledge on the seismic anisotropy signature of NW Iberia. Díaz et al. (2002Díaz et al. ( , 2006Díaz et al. ( , 2015 studied the upper mantle and lower crust anisotropy from distant earthquake shear wave splitting measurements and Acevedo et al. (2020) used local earthquake shear wave splitting and ANI analysis to infer the anisotropic properties of the upper crust in 80 the central part of the CM. Both studies consistently reported an average E-W fast direction, but the rotation of the fast orientations in the vicinity of major Alpine structures suggest that the contribution of the Alpine deformation to the anisotropy in the upper crust can be significant (Acevedo et al., 2020).
Here we make use of the ambient noise recordings from a new temporary seismic network in NW Iberia to extend to the North and the West an earlier tomographic S-wave velocity model of the eastern sector of the CM (Acevedo et al., 2019), 85 which was mostly affected by the Alpine deformation, and to further determine the 3-D distribution of the crustal radial anisotropy. From a subsequent work (Acevedo et al., 2020), we learnt that the crustal fabric controlled the orientation of fastest seismic velocities, sometimes in accordance with Variscan structures, others with Alpine structures. However, there was also local coincidence in the orientation of Variscan and Alpine structures, mostly by reworking during the youngest convergence. The aim of the new seismic array, building on our previous experience, was two-fold. Firstly, to apply the ANI 90 method to learn differences in architecture in relation to two orogenic events with contrasting kinematics (therefore orientation of structures). For this reason, we extended the area of observations to include parts of the Iberian crust with stronger Variscan imprint, incorporating domains with strong tectonic fabric susceptible to produce seismic anisotropy (Cárdenes et al., 2021) and with contrasting orientation to Alpine structures (Figs. 1a and b). Secondly, to constrain in the upper crust the contribution to seismic anisotropy by internally deformed rocks versus tectonically fractured rocks. 95 2 Geological setting

Tectonic history
The crustal architecture of the study area is mainly defined by the Variscan orogeny, which affected the continental crust of western Europe and the northwest of Africa between the Late Devonian and the Early Permian (e.g. Matte et al., 1986;Ribeiro et al., 2007;Pastor-Galán et al., 2013). Subsequently to the denudation of the 100 Variscan orogen, several rifting episodes affected the crust in NW Iberia, first in the Permian-Triassic, then in the Mesozoic where the crust was completely dismembered with the final opening of the Atlantic. The most important extensional episode occurred between the Late Jurassic and the early Cretaceous and was related to the opening of the Bay of Biscay (e.g. Cadenas et al., 2018;Tugend et al., 2015). This event configured the crustal structure of the eastern part of the study area, beyond the Ventaniella fault (Fig. 1a). 105 In Cenozoic times, the dominant geodynamic regime of the region shifted again and caused the North-South convergence between Europe and Africa, producing the onset of the Alpine orogeny. The compression started along the northern edge of the Iberian Peninsula in the Pyrenees and migrated to the West to form the CM (e.g. Teixell et al., 2018).
During the convergence, the transitional crust of the Bay of Biscay was underthrusted southwards beneath the Iberian continental crust (e.g. Boillot et al., 1979;Fernández-Viejo et al., 2012;2017). The bulk of 110 the deformation was mainly accommodated by East-West oriented reverse faults and thrusts, but there are numerous examples of the reworking and retightening of favourably oriented earlier Variscan thrusts and folds, with a different degree of development across the study area (e.g. Alonso et al., 1996). Currently, in the North Iberian crust, Alpine structures are dominant in the eastern half of the Variscan foreland-fold and thrust-belt, but they decrease their importance towards the West (Alonso et al., 1996;Martín-González et al., 2012;Llana-Fúnez and López-Fernández, 2015). 115

Crustal structure
Numerous works have contributed to broaden the knowledge about the Variscan imprint in the crust prior to the rise of the CM (e.g. Julivert, 1971;Pérez-Estaún et al., 1988;Alonso et al., 2009) and, more recently, its evolution during the Alpine convergence (Alonso et al., 1996;Pulgar et al., 1996Pulgar et al., , 1999Gallastegui et al., 2002, Martín-González et al., 2012Pedreira et al., 2015). From the early 90's, extensive seismic surveys were carried out in the area. These programs focused in 120 the unveiling of the deep crustal structure of the CM from deep seismic reflection and refraction/wide angle data (Pérez-Estaún et al., 1994;Pulgar et al., 1995Pulgar et al., , 1996Álvarez-Marrón et al., 1996Álvarez-Marrón et al., , 1997Gallastegui et al., 1997;Ayarza et al., 1998;Fernández-Viejo et al., 1998, 2000Gallastegui et al., 2002;Fernández-Viejo and Gallastegui, 2005;Fernández-Viejo et al., 2011 and led to significant findings that improved decisively the geological and geodynamical models of the area. For example, the E-W oriented profiles ESCIN-1 and ESCIN-3 imaged for the first time the Variscan structures of the 125 crust under the CM and their progressive loss of importance compared to the Alpine structures towards the East. Later, the ESCIN-2 profile was performed with the goal of sampling the Alpine features in the transition between the eastern sector of the CM and the Cenozoic Duero basin. Both profiles revealed the presence of an important crustal root beneath the CM. In order to complement the active geophysical surveys, several temporary seismic arrays with different specifications and coverage were deployed in the CM. Besides the seismotectonic characterization of the region (e.g. López-Fernández et al., 130 2012Llana-Fúnez and López-Fernández, 2015) the recorded passive seismic data allowed the study of the crustal configuration of North Iberia using receiver function analysis (Díaz et al., 2003(Díaz et al., , 2009aMancilla and Díaz, 2015) and validated many of the ESCIN observations. The use of ambient noise interferometry is also beginning to provide additional knowledge about the internal structure of the upper crust in the study area (e.g. Acevedo et al., 2019Acevedo et al., , 2020Olivar-Castaño et al., 2020). 135

Upper crustal domains in the study area
In the central sector of the CM the erosion of the Mesozoic sedimentary cover during the mountain building led to the exposure of the underlying Paleozoic basement. The physical properties of large parts of the upper crust targeted in this study with ambient noise tomography can be predefined according to the nature of the rocks and the structures, based on surface geology. As it has already been mentioned, a major part of the crustal architecture was shaped during the 140 development of the Variscan orogeny and for that reason, the crustal domains that will be defined below follow partially the subdivision established for the Variscan orogeny (Julivert et al., 1971). However, subsequent tectonic events, although did not produce a significant change in the nature of rocks, did have an impact on their structure. In most cases, producing a tendency to upright structures by increasing the dip of cartographic units and structures (e.g. Alonso, 1989).
The first crustal domain to the East of the study area is the foreland-fold-and-thrust belt of the Variscan orogen, the 145 external part of the orogenic belt, corresponding to the Cantabrian Zone (CZ) in the literature (e.g. Alonso et al., 2009). It is composed by sedimentary rocks that bear no tectonic or any other type of fabric orienting rock-forming minerals. However, the arrangement of tectonic units during the Variscan orogeny led in the final stages of the orogeny to the upright attitude of thrust sheets (e.g. Merino-Tomé et al., 2009). Folding of syn-orogenic sediments not involved in earlier Variscan thrusting also led to the upright position of lithostratigraphic units (e.g. Alonso y Pulgar, 1989;Aller and Gallastegui, 1995). 150 To the West of the Variscan foreland-fold-and-thrust belt in the basement lies the hinterland of the orogen, showing an approximate N-S arrangement in the study area (Fig. 1a). In terms of the rock types, their fabric, and their attitude we could distinguish four domains: a linear slate belt wrapping around the external part of the Variscan orogen, a subhorizontal slate belt coinciding partly with the Mondoñedo nappe, a subvertical high grade metamorphic domain intruded by Variscan granitoids, and the allochthonous nappes of high-grade rocks sitting on top of the previous domain. 155 The linear slate belt is a domain that can be defined ahead of the Mondoñedo nappe, over the Navia-Alto Sil domain (Marcos, 1973). From the structural point of view, most of the structures are steeply dipping (Marcos, 1973) and they coincide with an area with a N-S lineation fabric (Matte, 1968). The rocks present a tectonic fabric developed at low-grade metamorphic conditions, which imposes a strong mostly planar mechanical anisotropy to the propagation of seismic waves (Cárdenes et al., 2021). To the West, the subhorizontal slate belt coincides with the mid and rear part of the Mondoñedo 6 nappe ( Fig. 1) and its structure is dominated by large recumbent folds (Matte, 1968;Bastida et al., 1986). This isoclinal folding is accompanied by a tectonic fabric equilibrated in medium to high grade conditions that strengthens the subhorizontal anisotropy imposed by the orientation of the layers. Altogether, the linear and the subhorizontal slate belts constitute the West-Asturian Leonese Zone (WALZ).
The high-grade domain West of the Vivero fault ( Fig. 1a) corresponds to the Central-Iberian Zone (CIZ). It is 165 constituted by a rock sequence similar to that of the WALZ, but at higher grade given the pervasive intrusion of Variscan granitoids throughout the history of the orogenic belt. The recognition of large recumbent folds is hindered by the large amount of granitoids and by the development of upright folds and subvertical shear zones in the late stages of the Variscan orogen (e.g. Llana-Fúnez and Marcos, 2007). The folding is accompanied locally by the development of a strong crenulation fabric (e.g. Bastida et al., 2010). Subvertical corridors of foliated rocks in relation to large shear zones are common and often 170 affect several kilometres on either side of the crustal scale structures (e.g. Llana-Fúnez and Marcos, 2001).
The last domain corresponds to the allochthonous nappes, which are thrust sheets of exotic provenance thrusted over the CIZ (e.g. Martínez-Catalán et al., 1997). The presence of these allochthonous complexes defines the Galicia-Trás-os-Montes Zone (GTOMZ). Their essential features are that they are mostly made of high-grade high-pressure rocks, bearing a strong The upper crustal domains presented have been defined by the types of rocks and structures formed during the Variscan orogeny. Subsequent tectonic events in the study area lack the intensity and the pervasive character shown during the Variscan orogeny. No tectonic fabric at regional scale has been reported in relation to either the several extensional episodes leading to the opening of the Atlantic neither during the more recent Alpine convergence. Faults are also more spaced in 180 between.
Earlier studies have reported in the CZ domain a somehow widespread effect of Alpine deformation over previous rocks and structures, dominated by the tightening of previously steep Variscan structures (e.g. Alonso and Pulgar, 1989). Locally in this domain, there are reports of evidence of post-Variscan fracturing. No systematic study of Alpine fracturing has been carried out in the study area, however, separate observations in post-Permian rocks in the confines of the CZ show the 185 development of subvertical N-S set of joints (e.g. Lepvrier and Martínez García, 1990;Pastor-Galán et al., 2011;Uzkeda et al., 2013).

Seismic data
This study used continuous broadband seismic data from two experiments with non-simultaneous recording periods: (1) 190 Between June 2019 and February 2020, a seismic network (Geocantábrica-Costa Seismic Network, GEOCSN) composed of 7 permanent seismic stations belonging to the Spanish Seismological Network (SSN) of the Spanish Geographical Institute (Fig. 2a). The seismic array served the dual purpose of enhancing the location of the intraplate seismicity in the West of the CM and the close continental shelf and acquiring new passive data for ambient noise analysis. The experiment covered an 195 area of ~12,500 km 2 with interstation distances of ~40 km. Seismic stations were equipped with Nanometrics Taurus data loggers and broad-band Trillium 120 s sensors recording 3-component seismic data at a sample rate of 100 Hz. Most of the stations were installed in isolated chapels and sheds with connection to the electricity network, although solar power supply was required in some locations. A few stations were also mounted directly in the field, buried in ~0.5 m pits. The average data recovery reached 96%, with minimum data loss due to occasional power cuts caused by thunderstorms.
(2) To extend 200 our observations to the West, we also processed 12 months (January-December 2012) of continuous seismic data registered by 7 temporary stations of the northern IberArray seismic network (Díaz et al., 2009b) plus the station EPON of the SSN (Fig. 2a). The seismic recording equipment and their technical configuration were the same as in the GEOCSN, except that all the seismic stations were mounted outdoors. The station spacing was also larger, around 60 km, covering an area of ~10,000 km 2 with a 26% overlap with the GEOCSN. 205

Ambient noise interferometry
The treatment of the seismic ambient noise data to obtain Rayleigh-and Love-wave dispersion velocities was based on the widely used processing scheme described by Bensen et al. (2007), with a few modifications. Our processing began with the splitting of the continuous dataset in 24-hour long segments and the removal of all the traces with gaps. The records were downsampled to 25 Hz and the instrument responses were removed. Next, the mean and the trend were also removed from 210 the traces and a band-pass filter with corner frequencies of 0.01 and 2.0 Hz was applied. The described pre-processing benefited from the use of the Phase Cross-Correlation technique (PCC, Schimmel, 1999) for the calculation of the cross correlations. Given that PCC is amplitude unbiased (Schimmel et al., 2011), we did not perform any temporal or spectral normalization, thus avoiding the potential introduction of artifacts and shortening the processing time. The next step consisted in the cross correlation of the 24-hour traces recorded simultaneously by each pair of stations with the 215 aforementioned PCC procedure. In this way, we obtained the Daily Cross-Correlation Functions (DCCF) for lag times of 100 s. The observation of the cross-correlograms revealed a noticeable asymmetry between the causal and the acausal parts of the signal. Since we ignore which of these parts contain more useful information, the causal and the time reversed acausal parts of each station-pair DCCF were stacked using the Time-Frequency domain Phase-Weighted (tf-PWS) procedure of Schimmel and Gallart (2007). This stacking method enhances the coherent signal and cancels incoherent noise, allowing the 220 extraction of the symmetrical empirical Green's functions (EGF) of the medium between each station pair. In this study, we processed the vertical (Z) and the horizontal components (E, N) to obtain the Z-Z (Fig. 3a), E-E and N-N EGFs. Then, the horizontal components were rotated along the azimuth of the interstation path to provide the radial (R-R,  The measurement of the Rayleigh and Love group velocity dispersion curves was performed by using the multiple filter analysis (MFA) method of Dziewonski et al. (1969). In this technique, a narrow-band Gaussian filter is applied over different periods to isolate the wave package around the central period of each filter (Figs. 2d and S1). The group velocities in their fundamental mode were then manually picked to construct dispersion curves. According to Bensen et al. (2007), dispersion curves are limited at the longer periods by the fact that velocity measurements are only reliable for interstation 230 distances greater than three wavelengths. Given the dimensions of our network, the strict compliance of this criterion results in a loss of measurements that affects the resolution of the group velocity maps at periods > 10 s. However, Luo et al. (2015) demonstrated that consistent measurements can also be obtained using a cut-off of one wavelength.
In order to avoid discarding useful data without relaxing the constraints to a point at which the velocity determinations become unreliable, we adopted a two-wavelength criterion in this work. This cut-off has been proven to be a suitable 235 compromise solution in other studies (e.g. Sammarco et al., 2017;Brandmayr et al., 2016, Shapiro et al., 2005. To ensure the quality of the results, our analysis was limited to velocity measurements within two standard deviations from the mean and with more than 60 measurements at each period. The stability of the dispersion curves was also investigated by stacking sets of 60, 90 and the available DCCFs for all the interstation paths (Fig. S2). Overall, Love wave dispersion curves are slightly less stable than Rayleigh wave dispersion measurements. Waveforms and dispersion curves for the GEOCSN paths 240 stabilize at a lower number of days (~60 days) than IberArray paths (~90 days) but, in both cases, the stabilization is reached at a much lower number of days than the available for each interstation path (Fig. S2).
We decided not to extract phase velocities due to the difficulty of resolving the 2-pi phase ambiguity without long period data, due to the small aperture of our seismic network and the complexity of the CM crust. The estimation of the velocity uncertainties is vital to warrant the trustworthiness of the results. For this task, we have checked the group velocity 245 mismatch between the dispersion curves obtained from 10 stacks made with the randomly selected 75% of all the available DCCF (Fig. S2). The results show, in general, uncertainties well below 1% for the Rayleigh and Love wave group velocity determinations from vertical and transversal components respectively. In relative terms, Love waves display higher velocity uncertainties than Rayleigh waves, so as the IberArray dataset when compared with the GEOCSN measurements. Compared to vertical-component observations, Rayleigh wave velocities from radial-component EGFs display higher uncertainties, up 250 to ~4% at the longest periods (Fig. S2), and a lower signal-to-noise ratio. These facts, combined with the loss of some interstation paths due to dubious radial MFA surfaces, drove us to use exclusively Rayleigh group velocities obtained from the high quality vertical EGFs in further calculations.
On the other hand, the stacking of noise-correlations over long time spans ensures that the results are not affected by a seasonal distribution of the seismic sources. Moreover, the directionality analysis of noise sources in N Iberia performed by 255 Olivar-Castaño et al. (2020) showed that, although the most intense noise comes from the N-NW directions (Bay of Biscay and Atlantic Ocean), a sufficient level of noise emerges from all azimuths. Therefore, the effect of the spatial variation of noise sources is small and does not compromise the reliability of the dispersion measurements.

Group velocity tomography and depth inversion
Reliable dispersion measurements of Rayleigh and Love waves were obtained at periods between 2-14 s (Figs. 2e and 260 2f). The transformation of the dispersion velocities between the stations into continuous Rayleigh and Love group velocity surfaces was performed with the Fast-Marching Surface Tomography scheme (Rawlinson, 2005). The FMST is an iterative non-linear method that is implemented in two steps: a) the forward prediction of travel-times with the fast-marching method (Sethian, 1996;Rawlinson and Sambridge, 2004a;Rawlinson and Sambridge, 2004b) and b) the modification of the model parameters to explain the data observations, for which a subspace inversion scheme is used (Kennett et al., 1988). The 265 resulting models are controlled by a grid of nodes with cubic B-spline interpolation that generates a smooth and continuous velocity medium. We selected a grid size of 0.1º x 0.1º for the inversion of a homogeneous starting model that was created from the average velocity measured for a certain period (Nicolson et al., 2014). After extensive tests, the smoothing and damping parameters were set to 0.001 and 0.1, respectively, but their change only had a minor influence in the resulting models. 270 Based on the set of tomographic group velocity maps, we created new pseudo-dispersion curves by extracting the Rayleigh and Love velocity values for each point of a 0.1º x 0.1º grid. The independent inversion of these curves with the CPS surface inversion codes (Herrmann, 2013) allowed us to reconstruct the variation of the Vsv and the Vsh with depth in each grid node. In order to address the non-uniqueness nature of the inversions, we calculated two initial regional S-wave velocity model to be used as input in each of the subsequent computations. These models were obtained in a previous step 275 from the separate inversion of the average Rayleigh and Love dispersion curves in the whole study area (Figs. 2e and 2f).
We started in all cases from a homogeneous 22-layer model, the first four layers with a thickness of 0.5 km and the rest with a thickness of 1 km. Based on published studies in the area (Fernández-Viejo et al., 2000), we set a S-wave velocity of 3.35 km/s for each layer, with a fixed Vp/Vs ratio of 1.74 and a density of 2,600 kg/m 3 . The mismatch between the S-wave estimations from Rayleigh and Love waves allowed us to assess the intensity and distribution of the RA following Eq. (1): 280 (Yanovskaya et al., 1998;Savage, 1998;Ojo et al., 2017). All the inversions were performed with a damping factor of 10 in the first iterations to avoid an excessive shift from the initial model, whereas a value of 0.5 was used for the remaining iterations.

Rayleigh and Love wave group velocities
In total, we extracted Rayleigh and Love dispersion measurements of 78 station pairs from the GEOCSN network and of Rayleigh and Love group velocity variations can be observed in Fig. 4. Considering the Rayleigh waves (Fig. 4, left panels), the absolute group velocities range between 2.3 and 3.1 km/s. On the other hand, Love waves (Fig. 4, right panels) 290 show slightly faster velocities, varying between 2.3 and 3.4 km/s. At all periods, both Rayleigh and Love derived maps image two clear blocks: a high velocity zone occupying the western half of the study area and a relative low velocity zone in the eastern half. Both anomalies are separated by a narrow C-shaped transition sector which follows the boundary between the CZ and the WALZ and the surface trace of the main structures. The western anomaly presents high and moderately stable group velocities in the entire period range. On the contrary, group velocities in the eastern anomaly zone raise 295 constantly as the period increases. The most noticeable difference between Rayleigh and Love maps lies in the shape of the western anomaly at the highest periods. At those periods, Rayleigh group velocities display a light decrease, and the high velocities tend to concentrate in the South. However, Love group velocities delineate a consistent large high velocity surface in the western sector, while the low velocity anomaly migrates towards the South.
The resolution of the final velocity slices was investigated by performing synthetic checkerboard tests (Figs. S3 and S4). 300 We used as reference three models with small (20 x 20 km), medium (30 x 30 km) and large (40 x 40 km) anomaly sizes and velocity perturbations of ±0.4 km/s. Then, we evaluated the performance of our array for recovering this previously known velocity structure at periods of 4, 8 and 12 s. For both Rayleigh and Love waves, at periods up to 8 s, the resolution is sufficient to resolve anomalies of ~20 km in the central and eastern sectors of the study area, and ~30 km in the westernmost zone. At periods longer than 8 s, the resolution decays due to the progressive decrease of the number of interstation paths 305 and signs of smearing arise in the models. The average size of the recovered features increases to ~30 km. At those long periods, the Rayleigh paths recovered more accurate models than the Love paths. The RMS misfit between observed and predicted travel times ranged between 0.3-0.8 s, and 0.7-1.2 s for the Rayleigh and Love waves respectively, which indicates that the inverted models fit the observations and they are not affected by the selection of the regularization parameters (Nicolson et al., 2012). 310 Figure 5 shows examples of the S-wave inversion results in specific grid points located in the CZ (Fig. 5a), the WALZ ( Fig. 5b) and the CIZ+GTOMZ (Fig. 5c), featuring a good fit between the observed and the inverted data. The models and the anisotropy characteristics from the WALZ and the CIZ+GTOMZ show similar patterns, whereas the CZ presents some peculiarities, such as a constantly increasing Vsv and Vsh velocity with depth and a pronounced negative radial anisotropy 315 anomaly near the surface. Figure 6 shows examples of the variation of Vsv (Fig. 6, left panels) and Vsh (Fig. 6, right panels) at depths of 3, 6, 9 and 12 km. Within that depth range, the average Vsv increases from 3.0 to 3.4 km/s, while average Vsh oscillates between 3.1 and 3.5 km/s. Similarly to the group velocity maps, the presence of a low velocity anomaly to the East of the study area stands out both in Vsv and Vsh maps and its limit runs parallel to the curved boundary between the external and the internal zones, between the CZ and the WALZ. However, this anomaly attenuates at depth and finally disappears 320 below ~6 km. At greater depths (>9 km) the correspondence between the features depicted in the Vsv and the Vsh maps begins to fade away. In the Vsv surfaces, two new high velocity anomalies emerge in the South. However, the Vsh maps depict high velocities in the western/central sector and in the northeastern corner.

Radial anisotropy
The spatial distribution of the radial anisotropy (Fig. 7) was calculated from the discrepancies between Vsv and Vsh in 325 each node of the tomographic grid (Figs. 5a, 5b and 5c). In absolute values, the whole study area shows a significant anisotropy that reaches 10-15% in certain areas. Nevertheless, there is a clear differentiation with depth. At depths shallower than ~8 km the anisotropy strength keeps within the ±6% range. Beneath those depths, the anisotropy starts to rise up to ~15%. Negative anomalies (Vsv > Vsh) are dominant at shallow depths in the eastern half and the southwestern limit of the study area, with a weak positive corridor separating them. These negative anomalies become positive (Vsv < Vsh) at ~6 km 330 and ~8 km depth, respectively.
In order to represent the velocity variations within the context of the Variscan structure of the CM, we constructed an E-W cross section that runs perpendicularly to the main Variscan elements (Fig. 8). The most striking feature of the Vsv (

Discussion
Tomographic maps in Figure 6 depict the first high-resolution S-wave models of the transition between the internal and the external parts of the Variscan orogen in NW Iberia. Considering the sensitivity of the method (Fig. S5), the discussion of the results is restricted to the first 12 km of the crust. In addition, at the investigated periods (<14s) Love waves are primarily 350 sensitive to the uppermost crust, while the sensitivity peak of Rayleigh waves deepens as the period increases. The following discussions are only focused on well-resolved areas from the checkerboard tests (Figs. S3 and S4).

Velocity structure in the upper crust
In previous ambient noise regional studies on the Iberian Peninsula, Macquet et al. (2014) and Palomeras et al. (2017) reported S-wave velocity values at 5 km depths ranging between 2.5 km/s and 3.5 km/s in the study area. Likewise, a recent 355 work centered in the BCB also obtained this velocity range for the upper crust in the eastern termination of the CZ (Olivar-Castaño et al., 2020). Silveira (2013) and Palomeras (2017) identified a first order anomaly which corresponds with the Iberian Massif and recognized some of the tectonic units that constitute it. However, due to resolution limitations, the transition between the internal and external zones was poorly constrained. In contrast, the most striking element in both surface-wave group and shear velocity models presented here is the precise delineation of the boundary between these 360 domains by a sharp velocity change (Figs. 4 and 6).
The contour of this limit depicts a clear convex-to-the-west shape with the CZ located at its core. The contact between the high velocity anomaly, to the West, and the low velocity anomaly, to the East, is parallel to the Narcea antiform (NA), which is the structure that marks the boundary between the hinterland and the external zones of the Variscan orogen ( Fig. 1) (Gutiérrez- Alonso, 1992;Pérez-Estaún et al., 1994). The presence of rocks with completely different physical properties at 365 both sides of the eastern boundary of the NA justifies the observed velocity contrast. The CZ represents the foreland fold and thrust belt of the Variscan orogen, constituted by a pre-orogenic sedimentary sequence deformed in a thin-skinned deformation context, and a syn-orogenic clastic sequence that was deposited during the formation of the orogen. On the other hand, the slate belt to the West, is a relatively homogeneous domain mainly composed by siliciclastic rocks affected by regional deformation fabric formed at low grade metamorphic conditions. In summary, the upper crust in the CZ is 370 essentially made of sedimentary rocks without metamorphism while the upper crust in the NA and the WALZ slate belt is composed by igneous and volcanic rocks in the NA and low grade metamorphics towards the West, generating a clear contrast in the intrinsic physical properties of the rocks that are reflected in the crustal velocity patterns observed.
The sedimentary rocks of the CZ display relatively low S-wave velocities, typically under 3.0 km/s in the Vsv models or 3.2 km/s in the Vsh estimations. The profile of Figure 8a displays a good correspondence between the lower velocities and 375 the location of the siliciclastic wedges of syn-orogenic sediments. In fact, the low velocity anomaly reaches its greater thickness, around 5 km, in the CCB, the largest of these wedges. The depth estimation in the CCB is consistent with previous interpretations (Aller and Gallastegui, 1995). Towards the West, both the Vsv and the Vsh models depict high velocities. Vsv remain constant about 3.3 km/s, while Vsh increases steadily with depth (Fig. 6). The rock sequences in the WALZ are quite homogeneous, although in the cross-section of Fig. 8a some high velocity patches can be identified between 5-9 km depth. 380 The anomaly located close to the orogenic boundary is related with the stack and thrusting of rigid lower Paleozoic and pre-Variscan basement rocks from the NA (Pérez-Estaún et al., 1994) and the exhumation of basement rocks, gneiss and pre-Variscan igneous rocks between the NA and the Allande fault (Marcos y Pulgar, 1982;Rubio-Ordóñez et al., 2015). A smaller high velocity anomaly, in the West, is located near the Vivero fault (Fig. 8a) (López-Sánchez et al., 2015). High velocities were expected in this area due to the presence of mafic rocks and Variscan granitoids. However, the Vsh profile 385 (Fig. 8b) registers a velocity reduction in this domain, maybe caused by the fact that the resolution in the edges of the models decays and velocities are less well constrained. The basement rocks at the surface also show pervasive fracturing and weathering, strongly affecting velocities.

Thickness of the sedimentary layer and basement depth
A significant result from S-wave models is the estimation of the thickness of the sedimentary layer. In Macquet et al. 390 (2014), a S-wave velocity value of 2.9 km/s is proposed to indicate the limit between the basement and the sedimentary cover. In our case, this velocity defines quite well the position of the syn-orogenic sequence in the CZ, the presence of small Permo-Triassic basins in the CZ -BCB transition and the presence of Mesozoic sediments in the BCB but fails in recovering the depth of the basement. The main detachment of the CZ is believed to be located near the Precambrian-Cambrian boundary (Pérez-Estaún et al., 1988). Taking into account that the pre-Variscan basement rocks are exposed in the NA and 395 show velocities around 3.1 km/s in Figure 8a, we considered that this can be a more appropriate value to estimate the depth of the basement-cover surface in our study area to infer the position of the detachment. The 3.1 km/s iso-velocity line initiates in the NA area. It stretches towards the East for 20 km at shallow depths and then it sinks almost vertically until reaching 8 km. This deepening coincides with the maximum eastern extension of the pre-Variscan rocks at depth.
Continuing to the East, the surface dips smoothly to the West and alternates ramps and flats until reaching depths of 1 km to 400 the West of the Ventaniella fault. In the extension of Fig. 8a, based on the Vsv model of Acevedo et al. (2019), the isovelocity line deepens again, showing an apparent dip towards the East. This change is related with the different structure of the CZ-BCB transition, which was mainly configured by the Alpine inversion of Mesozoic extensional faults (e.g. Espina, 1997;Tavani et al., 2011). If we consider that the 3.1 km/s surface depicts the location of the base detachment of the CZ, its greatest depth ranges around 8-9 km. This observation agrees with the studies of Pérez-Estaún et al. (1988,1991,1994), 405 although it is slightly shallower than the estimations of Gallastegui et al. (1997Gallastegui et al. ( , 2000.

Radial anisotropy
From the relation between Vsv and Vsh in Eq. 1 we have obtained the first model of RA in the continental upper crust in NW Iberia (Figs. 7, 8c, 9 and 10). Overall, the RA patterns display a fair correspondence with the main Variscan tectonic units, at least in the shallowest crust (<7 km). At those depths, the distribution of the RA negative anomalies bound the 410 GTOMZ-CIZ and the CZ, whereas in the WALZ a constant positive anisotropy is observed (Fig 8c). The sign of the RA is related with the inclination of the anisotropic features. Thereby a negative RA implies the presence of steeply dipping to subvertical structures (60-90º) while a positive RA suggest gently dipping to horizontal features (0-30º, Dreiling et al., 2018;Xie et al., 2013). Thus, in the CZ, the negative anomaly can be explained by the verticalization of Variscan nappes and other structures, and the retightening of folds caused by the Alpine convergence, which was more intense in this domain than in 415 any other. Moreover, the presence of a N-S oriented pervasive Alpine fracture system with subvertical dip has been recognized in the area and may play a significant role in the development of the anisotropy.
As the depth increases, the thrusts reduce their dip as they approach to their imbrication with the CZ basal detachment (Pérez-Estaún et al., 1991), which is reflected in the progressive shift to a positive anomaly (Fig. 8c). There are two factors that suggest that this change can be provoked by the presence of a basal detachment at those depths. First, the positive RA 420 points towards the presence of smoothly dipping anisotropic features. Second, the RA sign inversion occurs at greater depths beneath the western CZ than the eastern CZ, where the positive anomaly gets close to the surface. This point implies the presence of a possible structure with a general dip towards the West. Laterally, RA sign variations are also present, as in the transition between the CZ and the WALZ. From the Narcea antiform to the Mondoñedo thrust, subvertical structures are dominant, but the existence of a Variscan subhorizontal foliation parallel to the fold axis may compensate their effect and 425 explain the slightly positive anomaly instead of the expected negative anisotropy. To the West of the MT, in the CIZ and the GTOMZ, the structures and foliation are mainly gently dipping, contrasting with the presence of a negative RA anomaly that spreads from the surface to 8 km depth. Strong subvertical deformation fabric related to late Variscan shear zones and folds (Llana-Fúnez and Marcos, 2007) and N-S fracture systems associated to the Alpine convergence have been reported in the area, but the resolution of the models at their edges is not optimal and further studies should be carried out to verify the 430 existence of such negative anomaly.
In general, the magnitude of the RA is low, around 4%, at shallow depths (< 8 km). Considering the entire range of investigated depths, except for the WALZ, where the RA magnitude increases almost steadily, a progressive decrease followed by a rapid increment of the RA anisotropy magnitude is observed in the CIZ/GTOMZ and the CZ (Fig. 8c). The latter domain was also targeted by a previous ambient noise-based study that estimated the azimuthal anisotropy in the area 435 (Acevedo et al., 2020). A comparison between the RA measured in a representative CZ node and the azimuthal anisotropy determinations is displayed in Fig. 9a. Despite some discrepancies in their magnitude at shallow depths (< 3 km), both azimuthal and radial anisotropy depict a decline at depths ranging from ~3 to ~7 km. At greater depths, the RA rises gradually until it reaches magnitudes around 9% at 12 km, while azimuthal anisotropy increases abruptly at 7.5 km depth and stabilizes at 9%. The progressive closure of the open cracks due to the increment of the lithostatic pressure can explain 440 the observed radial and azimuthal anisotropy magnitude reduction at shallow depths (Figs 9a and 9b). In crystalline low porosity rocks, the complete closure of microcracks at grain boundaries is observed in laboratory experiments at about ~200 MPa, which is equivalent to depths around 7 km (e.g. Christensen, 1985;Brown et al., 2009;Guo et al., 2014) and coincides with both the radial and the azimuthal anisotropy minimums (Fig. 9b).
At shallow depths (<5 -10 km), the RA is believed to be caused by the alignment of stress- (Crampin, 1978;Crampin 445 and Peacock, 2005) or structure-induced fluid saturated cracks (Zinke and Zoback, 2000). The fast anisotropic directions obtained from shear-wave splitting in Acevedo et al. (2020), which tend to align parallel to major structures, favor a structure-controlled anisotropy model for the CZ. What is more, although the Alpine deformation did not produce pervasive deformation in NW Iberia, it generated shallow fracture systems that can affect decisively the physical properties of the rocks near the surface, especially in the CZ and the CIZ/GTOMZ. On the other hand, in the case of slates from the WALZ, 450 recent laboratory data have shown that the magnitude of the anisotropy increases with the decrease of confining pressure, as the orientation of the microcracks is strongly controlled by the shape fabric in very well foliated slates (Cárdenes et al., 2021). However, this means that another origin must be invoked to produce the RA magnitude increment at higher depths ( Fig. 9a). As it is depicted in Fig. 1, the crystalline basement rock outcrops in the WALZ, showing an anisotropy that increases from 2% to 15% in depth. These Precambrian rocks, as well as the Lower-Paleozoic pre-Variscan siliciclastic 455 sequences that form the WALZ and constitute the basement of the CZ, are affected by low-grade metamorphism and the development of a widespread tectonic foliation. Anisotropic minerals, such as phyllosilicates, are common in these rocks, so the alignment of minerals and grains (Wenk et al., 2020;Cárdenes et al., 2021) seems to be the most likely mechanism controlling the anisotropy in the study area and at depths > 7.5-8 km. In the CZ, this change in the deformation and tectonic history of the rocks controlling the anisotropy is also reflected in the fast azimuthal directions. At depths < 7.5 km, they 460 remain steady at around 84º, whereas at depths at which the pre-Variscan basement gains influence in the results (> 7.5 km) a progressive NE rotation (~70º) can be noticed (Fig. 9a).
Many studies have tried to establish the linkage between tectonic regimes and RA patterns (e.g. Luo et al., 2013;Dreilling et al., 2018;Wang et al., 2020). Generally, compressional and extensional regimes are associated to the vertical and horizontal alignment of cracks and minerals, respectively. However, the specific geodynamic setting of each location 465 must be considered, given that both regimes can produce very different anisotropy signatures. The presence of negative RA at shallow depths in foreland fold-and-thrust belts has been observed in orogens worldwide. For example, Chen et al. (2009) and Guo et al. (2012) reported negative RAs in the Himalayan fold-thrust belt. Recently, Wang et al. (2020) described a widespread negative anomaly of ~5% below the Foreland belt of the Canadian Cordillera that turns into a positive anomaly in the Omineca Belt, the hinterland of the orogen. These authors suggest that the RA in the foreland fold and thrust belts is 470 caused by the horizontal contraction of the upper crust, which creates subvertical faults and/or preferentially alignments of minerals under convergent regimes. The continuity of the negative anomaly at depths up to 20 km allowed them to infer a "thick-skinned" Cordillera-Craton transition. Contrarily, in our model, the negative anomaly in the CZ is distinctly limited at depth (Fig. 10), pointing to a "thin-skinned" style of deformation, in accordance with previous studies (e.g. Pérez-Estaún et al., 1994). On the other hand, positive anisotropy in the hinterland is usually associated to horizontally sheared fabrics, 475 which in the case of the Variscan orogen could either be caused during pervasive deformation associated to plate duplication in convergence (e.g. Martínez-Catalán et al., 1997), or by syn-or post-collisional localized crustal extension, as reported for the Mondoñedo thrust (Marcos, 2013). Ductile crustal flow has also been invoked as the origin of positive RA anomalies, mostly in the middle and lower crust (Lynner et al., 2018), and this has also been suggested to occur at the late stages of the Varican orogeny for a significant part of the hinterland of the Variscan orogeny in NW Iberia (Llana-Fúnez and Marcos, 480 2007).
In this study, we used group velocity measurements from Rayleigh and Love waves to derive the first S-wave velocity model of the continental crust in NW Iberia. An independent inversion of the Rayleigh and the Love group velocity estimations was done to create a vertically polarized S-wave model and a horizontally polarized S-wave model, respectively. 485 The first model displayed velocities between 2.3 and 3.6 km/s, while the velocities of the latter ranged between 2.6 and 4.2 km/s. Both models clearly depict the contrast in physical properties between the rocks located in the external (CZ) and the internal (WALZ, CIZ and GTOMZ) zones of the Iberian Variscan orogen, delineating the core of the large Western European Variscan belt: the core of the Ibero-Armorican arc at the Cantabrian Zone. Moreover, the high resolution of the models and the good existing geological knowledge of the area allowed us to associate the observed velocity changes to the 490 variation of the bulk properties of the rocks caused by other first-order orogenic structures, such as the basal detachment of the CZ or the presence of a pre-Variscan basement in the NA, or sedimentary features such as the presence of syn-orogenic wedges in the CZ.
The crustal radial anisotropy 3-D model has been proved useful to identify and characterize the different sectors and units that form an orogenic belt and show a good correspondence with the dip of the main structures. The anisotropy pattern 495 can also be associated to the deformation history of the orogen. A negative anomaly is found in the foreland fold-and-thrust belt, which can be related to the stack of thrusts and vertical fabrics caused by compression. In the internal areas, widespread positive anisotropy is accounted. The radial anisotropy magnitude oscillates between -4 and 15 % and increases with depth.
Analyzing the depth pattern, we concluded that, at < 8 km depths, the alignment of Alpine fracture systems controls the anisotropy. At depths > 8 km, high and positive radial anisotropies are observed, caused by the presence of subhorizontal 500 alignments of grains and minerals, possibly in relation with shearing in the basal detachment of the CZ and/or other shear zones formed during the development of the Variscan orogeny.

Data and resources
The datasets presented and used in this study were collected using a seismic network funded by project GEOCANTABRICA-COSTA (https://doi.org/10.7914/SN/YR_2019) and it can be released to the public on demand at 505 GEOCANTABRICA@ftp.geol.uniovi.es. The deployment of the IberArray broadband seismic network (https://doi.org/10.7914/SN/IB) was part of the CONSOLIDER CSD2006-00041 (Geosciences in Iberia: Integrated studies on Topography and 4-D Evolution) grant from the Spanish Ministry of Science and Innovation. This dataset is available online at the ORFEUS Data Center (http://www.orfeus-eu.org/). Data processing used the codes: Corr_stack_v04.1 (Schimmel et al., 2011) and Computer Programs in Seismology v3.30 (Herrmann, 2013). Figures were drafted using Generic 510 Mapping Tools (Wessel and Smith, 1998) and QGIS 3.10.   Fig. 2a, except for paths crossing the CIZ, GTOMZ or both, which are indistinctly marked in pink. Grey dispersion curves represent mixed paths that span through more than one region. Red circles illustrate the average dispersion values used in the computation of the initial regional S-wave velocity models. U: group velocity.