Crustal structure beneath Portugal from teleseismic Rayleigh Wave Ellipticity

Up until now, Portugal lacked a countrywide shear velocity model sampling short length-scale crustal structure, which limits interpretations of seismicity and tectonics, and predictions of strong ground motion. In turn, such interpretations and predictions are important to help mitigate risk of destruction from future large on- and off- shore earthquakes similar to those that Portugal has experienced in the past (e.g. the M w 8.5 – 8.7 tsunamigenic event in 1755). In this study, we measured teleseismic Rayleigh Wave Ellipticity (RWE) from 33 permanent and temporary seismic stations in Portugal with wave periods between 15 s and 60 s, and inverted it for 1-D models of shear wave velocity (V s ) structure beneath each station using a fully non-linear Monte Carlo method. BecauseRWEisstronglysensitiveto theuppermostfewkilometresofthecrust,bothRWEmeasurements andV s models are spatially correlated with surface geology in Portugal. For instance, we ﬁ nd that sedimentary basins produced by rifting that had begun in the Mesozoic such as the Lusitanian Basin (LB) and the Lower Tagus- Sado Basin (LTSB) are characterised by higher RWE (lower V s ). Interestingly, we observe similar RWE (and V s ) values in the interior of the Central Iberian Zone (CIZ), which is a metamorphic belt of Paleozoic age. Together with reduced crustal thickness previously estimated for the same parts of the CIZ, this suggests that the CIZ might have experienced an episode of extension possibly simultaneous to Mesozoic rifting. The Galicia-Tras- os-Montes-Zone (GTMZ) that has undergone polyphased deformation since the Paleozoic is characterised by the lowest RWE (highestV s ) inPortugal.OssaMorenaZoneandtheSouth Portuguese Zoneexhibit intermediate V s values when compared to that of basins and the GTMZ. Our crustal V s model can be used to provide new in- sights into the tectonics, seismicity and strong ground motion in Portugal. the (http:// creativecommons.org/licenses/by/4.0/).


Introduction
Strict constraints on the seismic structure of the crust are important for interpreting tectonics and seismicity accurately. Also, complementary information that can be gathered from accurate 3-D models of crustal structure (Denolle et al., 2014a,b;Bowden and Tsai, 2017) could improve ongoing seismic hazard assessment efforts (e.g. Vilanova and Fonseca, 2007). Although ambient noise and teleseismic surface wave dispersion data can provide accurate images of the subsurface (e.g. Yang et al., 2007;Morelli, 2009, 2011;Zhu et al., 2015), their sensitivity to the uppermost few kilometres of the crust is limited unless periods less than about 5 s are considered (Lebedev et al., 2013).
Mapping crustal structure with such short period dispersion data is challenging due to attenuation effects from intrinsic mechanisms and scattering. Further, because these surface wave dispersion methods provide path-averaged images, the lateral resolution of structures mapped can be in the order of 100 s of kilometres (Bensen et al., 2008), effectively smearing short-scale features. Both these limitations of routine surface wave dispersion methods contribute to model non-uniqueness. These issues can be addressed if complementary measurements that have stronger sensitivity to the uppermost crust are used. Teleseismic Rayleigh Wave Ellipticity (RWE), defined as the horizontal-to-vertical ratio of fundamental mode Rayleigh wave particle motion, depends mainly on the crustal structure immediately beneath a given station (e.g. Malischewsky and Scherbaum, 2004;Tanimoto and Alvizuri, 2006;Ferreira and Woodhouse, 2007a,b;Tanimoto and Rivera, 2008), with strong sensitivity to the uppermost crust (e.g., Tanimoto and Tsuboi, 2009;Berbellini et al., 2016). Since it is an observable based on amplitude data, which have stronger sensitivity to small-scale structure than travel-times, it allows illuminating shorter-scale heterogeneity in the uppermost crust than dispersion data. Moreover, being a single station observable, in principle, RWE is little affected by path-averaged effects. With the expansion of high-quality seismic networks in the past decade, this method has produced robust results (e.g. Yano et al., 2009;Berbellini et al., 2017) as opposed to limitations faced in early studies (Boore and Toksöz, 1969). In particular, RWE has been combined with other data types to reduce model non-uniqueness (e.g. Lin et al., 2012;Li et al., 2016).
Despite being located in a region of slow lithospheric deformation in an intraplate setting, Portugal is exposed to significant seismic risk. Its history is fraught with destruction from large on-and offshore earthquakes, notably the 1755 great Lisbon earthquake, which is the largest magnitude (M w 8.5-8.7) historical earthquake recorded in Europe (Vilanova et al., 2003;Ferrão et al., 2016). The lack of high quality broadband (BB) seismic networks in Portugal has hindered the construction of detailed seismic images of Earth structure in the region, impairing the accuracy of models of seismotectonics and precision of strong ground motion predictions. For instance, the Moho depth model compiled for the Iberian Peninsula proposed by Díaz and Gallart (2009) contains active seismic experiment results largely restricted to the Northernmost and Southernmost regions in Portugal, due to which effects of spatial aliasing of structure are present. Since 2006, however, the permanent Portuguese seismic network has greatly expanded along with dense temporary deployments (Custódio et al., 2014), which resulted in a continuously accumulating large volume of high quality BB data. As a result, new constraints on the crustal structure beneath Portugal have begun to emerge. For instance, recently, Dündar et al. (2016) have used this dataset to map depth to Moho beneath Portugal using P-wave receiver functions, and found that the average depth to Moho is 30 km, with modest lateral variations. Previously, Salah et al. (2011) mapped the crustal structure beneath the lower Tagus valley combining receiver functions and surface wave dispersion data, while Morais et al. (2015) inverted receiver functions for the crustal and mantle structure beneath the permanent station MTE. In addition, Silveira et al. (2013) measured Rayleigh wave dispersion in western Iberia to produce maps of group velocity perturbations. These new measurements and models of Earth structure are bound to improve interpretations of tectonics, characterisations of seismicity, and strong ground motion predictions in Portugal, which is important given its exposure to significant seismic risk.
In this study, we invert teleseismic RWE measurements obtained from 33 permanent and temporary onshore seismic stations for the crustal V s structure beneath Portugal. In the next sections we review the geology and seismicity of Portugal briefly, followed by a description of our data and methods. In Section 5 we present our results, which is followed by a discussion in Section 6 and by the conclusions in Section 7.

Geologic units in Portugal
The present-day geologic makeup of the Western and Central Europe is primarily associated with the Variscan Orogeny of late Paleozoic age (480-290 Ma), which resulted from the collision of the Gondwana supercontinent with Laurussia, and the closure of the Rheic Ocean (Matte, 1986;Matte, 2001;Pous et al., 2011). There is, however, no consensus on the exact model of collision, which varies from a simple twoplate model to a multiple plate model (Scotese and McKerrow, 1990;Bois et al., 1994;Tait et al., 2000;Franke, 2000;Kroner and Romer, 2013). The juxtaposition of high-and low-grade terrains along with heterogeneous distribution of metamorphic belts and sedimentary basins possibly reflects the complexity of the Variscan Orogeny (e.g. Cymerman et al., 1997;Kroner et al., 2007;Hahn et al., 2010).
The Iberian Massif, the westernmost outcrop of the Variscan Orogeny, covers much of Western Iberia (Kroner and Romer, 2013). In Portugal, four out of six distinct tectonostratographic units of the Iberian Massif are exposed (Fig. 1a). In the North, the Galicia-Trás-os-Montes Zone (GTMZ) is an elongated tectonic unit characterised by Biotite and Andalusite bearing meta-igneous and metasedimentary rocks. The presence of Andalusite in the GTMZ indicates deformation due to low pressure-mid temperature metamorphism. The Central Iberian Zone (CIZ) covers much of central Portugal and is a Schist-Greywacke lithotectonic unit. Both GTMZ and CIZ bear evidence of melting of orthogneisses and other metamorphic rocks at crustal level (Castro et al., 1999;Fernández-Suárez et al., 2000;Neiva, 2002;Vegas et al., 2001). The Tomar-Badajoz-Córdoba Shear Zone (TBCSZ) has a sinistral sense of motion and separates the CIZ in the North from the Ossa Morena Zone (OMZ) in the South (Burg et al., 1981;Ribeiro et al., 1990;Azor et al., 1994;Simancas et al., 2001;Gomes et al., 2007). The OMZ is divided into two sub-zones: Estremoz-Brancos to the North and Beja to the South (Fig. 1b). These subzones are separated by the Terena syncline from each other (Ribeiro et al., 2010). Thrusting of allocthonous complexes onto autochthon ones has generated complex structures in the Beja sub-zone with rocks exhibiting variable metamorphic facies (e.g. Pin et al., 2008), whereas the Estremoz-Brancos subzone appears simpler in structure and contains rocks of greenschist facies (Araújo et al., 2006). There are also two Ophiolite series (Internal Ossa-Morena Zone Ophiolite Sequences and Beja-Acebuches Ophiolite Complex) that differ in age, evolution, and emplacement mechanisms, which form the boundary between OMZ and the Southern Portuguese Zone (SPZ [Ribeiro et al., 2010]). This boundary continues westward and links up with the Porto-Tomar-Ferreira-Alentejo Shear Zone (PTFASZ) that exhibits a dextral sense of motion. The southernmost extent of the Iberian Massif is marked by SPZ whose polyphased tectonic evolution is constrained from detailed analysis of structure and petrology (e.g. Onézime et al., 2003).
Geologically more recent large-scale structures are also present in mainland Portugal. The opening of the North Atlantic Ocean associated with more recent extensional tectonics of Mesozoic age (~180 Ma) has produced rift basins to the West (Lusitanian Basin [LB]) and South (Algarve Basin [AB]) of Portugal (Dewey et al., 1989;Wilson et al., 1989;Roest and Srivastava, 1991;Pinheiro et al., 1996). Locked between Mesozoic LB to the West and pre-Mesozoic tectonostratigraphic units of the Iberian Massif to the East, is the Lower Tagus-Sado Basin (LTSB). It is interpreted to be a geologic structure originating from tectonic inversion of LB taking place since the Eocene. Alternatively, Carvalho et al. (2016) hypothesized that the formation of a pull-apart basin in the Palaeogene followed by polyphased-deformation might be responsible for the formation of the LTSB. The origin and evolution of these basins appear to be controlled by pre-existing tectonic structures including systems of faults particularly related to the Variscan Orogeny (Carvalho et al., 2005;Lopes et al., 2006). The relationship between the surface geological features in Portugal and structure at depth, however, is still poorly understood. One of the aims of this study is to fill this knowledge gap.

Seismicity in Portugal
Present-day seismicity in mainland Portugal appears to be related to onshore slow-deformation (~1 mm/yr or less) driven by offshore oblique convergence of the Eurasian and Nubian plates occurring at a rate of about 5 mm/yr (Argus et al., 1989;DeMets et al., 1994;Borges et al., 2001;Kreemer and Holt, 2001;DeMets et al., 2010;Nocquet, 2012). The geometry of the plate boundary between the two plates is more complex than a classical narrow plate boundary (Fig. 2a, Buforn et al., 1988;Jimenez-Munt et al., 2001;Bezzeghoud et al., 2014). From the Azores to the Madeira Tore Rise, a fairly well defined narrow plate boundary characterised by right-lateral strike-slip faulting exists. Movement along the Gloria fault that exists between 24°W and 19°W along this plate boundary segment is suggested to have generated the M8.3 tsunamigenic event in 1941 (Baptista et al., 2016). The plate boundary between the Madeira Tore Rise and the easternmost region of Gulf of Cadiz is poorly defined and appears to be characterised by NW-SE to NNW-ESE oblique convergence. The observed diffuse distribution of seismicity and active structures along this segment suggest that deformation is accommodated across a broad system of fault networks under spatially varying tectonic forces rather than along a single segment of plate boundary sustained by regionally uniform deformation (Sartori et al., 1994;Grácia et al., 2003;Gutscher et al., 2009;Zitellini et al., 2009;Palano et al., 2015). It is hypothesized that deformation on this segment of the plate boundary is responsible for large events including the 1755 Lisbon earthquake (Buforn et al., 2004). Still further east, plate motion is distributed across the Alboran Sea Basin and Betics and Rif Cordilleras.
The rate of deformation varies from about 1 mm/yr in the southernmost regions to b1 mm/yr in Northern Portugal (Negredo et al., 2002;Nocquet, 2012). Despite this low deformation rate, historical records point to numerous large onshore earthquakes that have occurred since the 14th century (Fig. 2a, Vilanova et al., 2014;Stucchi et al., 2013;Ferrão et al., 2016). The 1909 April 23rd M6.3 Benavente earthquake is the most recent significant onshore event. The reverse faulting mechanism estimated for this particular event (Stich et al., 2005) is consistent with the stress regime predicted for the region (Stich et al., 2006;Heidbach et al., 2009). Also, a recent morphometric analysis suggests that the Lower Tagus Valley Fault Zone has the potential of producing an M7.3 earthquake if the entire length of the fault were to fail in a single rupture (Besana-Ostman et al., 2012). In general, however, correlating large historical on-and offshore events with geologic structures has been challenging (Cabral, 2012) due to limited availability of information about sources and Earth structure (e.g. Cabral et al., 2013;Vilanova et al., 2014).
Since the 1960s, the Portuguese seismic network has undergone vast improvements. The volume and the quality of data have improved manifold particularly since 2006. Currently there exist over 30 permanent broadband stations, of which stations used in this study are shown in Fig. 2b along with stations from the temporary WILAS array (Custódio et al., 2014). The expansion of the Portuguese seismic network significantly increased the number of seismic events detected per year. Moreover, the detection threshold has dropped from about NM3 in the early days to b M1 at present (Sousa et al., 1992;Carrilho et al., 2004;Pena et al., 2014). Taking advantage of these developments, Custódio et al. (2015) re-analysed the Portuguese seismic dataset to produce a dense map of seismicity containing nearly 15,000 events mostly in the M0.5-M4 range (Fig. 2a). They found distinct clustering of on-and offshore earthquakes, many of which are possibly associated with unmapped networks of faults. One possibility for this observation can be errors in the 1-D crustal models that were used for relocation. Their map also alludes to the fact that there is heterogeneous spatial distribution of seismicity with the cumulative seismic moment and the number of events decreasing from South to North. This observation correlates well with the spatially varying average displacement across mainland Portugal predicted in plate kinematic models (Negredo et al., 2002;Nocquet, 2012).
Further improvements to this emerging new picture of seismicity and tectonics can be made if more accurate models of Earth structure become available. Because the accuracy of earthquake locations and focal mechanisms is strictly controlled by the quality of constraints on Earth structure (e.g., Ferreira and Woodhouse, 2006;Ferreira et al., 2011;Weston et al., 2014), it is crucial to build precise models of 3-D Earth structure. In addition, more advanced studies of waveform modelling to constrain source physics and better understand the relation between seismicity and tectonics can only be undertaken with the availability of such Earth models. Thus, our countrywide V s crustal structure model has potential wide-ranging applications.

Data and method
We collected seismic waveforms from all available permanent and temporary seismic stations in Portugal for the period 01/Jan/2009/-28/Feb/2015. Our earthquake dataset is limited to the magnitude (M w ) range 6.0-7.8 to optimize the Signal-to-Noise Ratio (SNR). We considered all possible source-receiver azimuths. Seismograms were selected between 50°and 120°epicentral distance range with the intention of minimizing near-source and wave interference effects (Tanimoto and Rivera, 2008). To avoid potential undesirable effects of interference from multiple events, we only considered events at least 90 min apart. We found a total of 454 events that fall within these criteria (Fig. 3). A list of these events is provided in the Supplementary information. These seismograms were windowed between the origin time in the GCMT catalogue ; Ekström et al., 2012; http://www.globalcmt.org/CMTsearch.html) and 5400 s thereon to include the complete surface wave train, after which standard processing steps were applied, i.e., instrument response deconvolution, decimation to 1 s and rotating horizontal components to radial and transverse directions. Finally, using a 4th order bandpass Butterworth filter (also see Section 4.3), the data were filtered between central periods of 15 s and 60 s at an interval of 5 s. The range of wave periods used in this study is determined by the response curves of the seismometers. The response spectra for stations from the permanent networks are generally flat up to T~120 s, whereas they are only flat up to T 60 s for the temporary network 8A (Fig. 2). Because RWE is less sensitive to deep crustal and mantle structure (corresponding to long wave periods) than to shallow crustal structure, and to maintain a uniform approach across all stations, we decided to only use measurements up to T~60 s. Custódio et al. (2014) reported that stations PCVE and PVAQ suffer from significant misalignment of sensor orientation. Upon conducting our own analysis of sensor orientation following the same approach as in Berbellini et al. (2016), we found that these two stations are indeed not favourable to measure RWE and thus we discarded them from the analysis. We also discarded station PBAR because our analysis identified calibration issues with its horizontal and vertical component data.

Measurement of RWE
We used a fully automated time-domain method to measure RWE. Time-domain measurements are shown to be better than frequencydomain ones because effects of wave interference are reduced in time domain (Sexton et al., 1977), and there is no need to deal with spectral holes. In the first step, we windowed the surface wave train using fundamental mode phase velocities predicted for PREM (Dziewonski and Anderson, 1981) applicable to the range of central periods of interest. We then picked a fixed secondary time window around the absolute peak (or trough) of the vertical component found within the time window picked in the first step. The length of this secondary time window was based on synthetic tests and usually corresponds to 2-3 wavelengths at a given central period. We then followed the same strategy as Ferreira and Woodhouse, 2007a (method I), whereby, after Hilberttransforming the radial component and picking the same secondary time window as for the vertical component, we obtained RWE by dividing the sum of the square of each time point of the radial component by that of the vertical component, and taking the square root of the resulting ratio (Eq. (1)). In Eq. (1), RWE is Rayleigh Wave Ellipticity, i = 1,n are the time points in the secondary time window, R is the Hilbert-transformed radial component ground velocity, and V is the vertical component velocity. This window-picking procedure is motivated by the fact that we are specifically searching for records dominated by fundamental mode Rayleigh waves to maximize SNR and reduce contamination by higher mode energy.

Quality control of data
We applied three main quality control steps to ensure that RWE is measured within a time window dominated by fundamental mode Rayleigh waves. Firstly, the phase shift between the vertical component and the Hilbert-transformed radial component must be in the 70°-90°range. Ideally, for fundamental mode Rayleigh waves, this phase shift should be 90° (Tanimoto and Rivera, 2005). As Ferreira and Woodhouse (2007a) have pointed out, however, deviations of phase shift from 90°should be expected due to interaction of waves with small-scale heterogeneities having length-scales shorter than the dominant wavelength of Rayleigh waves and also due to contamination from higher mode energy. Considering the short length-scale nature of geology manifested in the lithologic map of Portugal, and the observation of complex wave propagation effects of short period fundamental mode Rayleigh waves in general (Pedersen et al., 2015), we deemed that a 20°deviation is reasonable for our case.
Secondly, we calculated the cross correlation coefficient between vertical component waveforms and corresponding Hilbert-transformed radial component waveforms within the picked measurement window. The effect of the Hilbert transform is to advance the phase of waveforms by 90°. A characteristic feature of Rayleigh waves is that the cross correlation coefficient between the vertical and Hilbert-transformed radial waveforms approaches 1.0 under ideal conditions (Tanimoto and Rivera, 2008). Thus, we measure RWE only when our data have a cross correlation coefficient equal or N0.8. Fig. 4 shows an illustrative example of measurements for station GGNV both on observed and synthetic waveforms at 25 s period along with the corresponding phase differences, correlation coefficients and RWE values.
Thirdly, we applied lower and upper cut-offs to RWE measurements before finalizing the dataset. We used a lower bound of 0.2 and an upper cut off of 2.0 to our dataset. Both these bounds were picked based on synthetic tests. RWE b 0.2 indicate unrealistic high seismic velocities immediately beneath a given station, whereas RWE N 2.0 indicate anomalously low velocities. Recent maps of short period Rayleigh wave dispersion in western Iberia (Silveira et al., 2013) do not indicate the presence of such unrealistic low or high velocities. The RWE bounds we used are also consistent with those employed in previous studies (e.g. Lin et al., 2012).
Finally, we obtained RWE for a given wave period and station by computing the median of the measurements that passed quality control, which reduces the effect of outliers. By averaging the first and the third percentiles, we obtained our bounds of uncertainty. Fig. 5 shows the distribution of the measurements obtained for three illustrative stations, along with the retrieved RWE values and uncertainties (also see Fig.  S2 in Supplementary information). The measurements show a smooth and stable behaviour as a function of wave period, exhibiting near-normal distributions.

Handling higher mode energy
A potential source of error affecting the RWE measurements is the contamination of the fundamental mode Rayleigh wave by higher mode energy. To test if our measurement scheme was susceptible to such a contamination, we carried out synthetic tests. Since higher modes are more efficiently excited by deep earthquakes rather than by shallow events, we investigated synthetic seismograms computed by normal mode summation (Gilbert, 1970) for the PREM model (Dziewonski and Anderson, 1981) by successively changing the source depth from 10 km to 50 km with an increment of 10 km (see examples of synthetic waveforms for a source's depth of 12 km in Fig. 4). We then used these synthetic seismograms to calculate RWE using our automated measurement scheme. Fig. 6 shows an example of the measurements obtained for the same event-station pair used in Fig. 4, which shows an overall good retrieval of the theoretical ellipticity values predicted for PREM (black curve). In our synthetic tests we clearly observed that, at certain wave periods, higher mode energy dominated the surface wave train. We found, however, that there was a clear separation between fundamental mode and higher mode energy in time. Upon closer examination, we observed that, as expected, there is a tendency for higher mode energy to grow with earthquake depth, where the radial component is more affected by it than the vertical component. Our tests suggested several steps that can be taken to reduce contamination of RWE measurements by higher mode energy: 1. Time windowing based on fundamental mode phase velocities: We recall that our measurement scheme uses a primary window defined by fundamental mode phase velocities, within which a secondary window is picked on the vertical component (which is less affected by higher modes) to measure RWE. Because higher mode energy travels at a different velocity than fundamental mode energy, this is the simplest first precaution that can be taken to avoid measurement bias. 2. Adjusting the bandwidth of filters: Our observations suggest that increasing the bandwidth of filters as the central period increases could minimize the influence of higher mode energy. From our synthetic tests for the source-receiver paths considered in this study, we found that a good strategy was to use a 10% bandwidth about the central period for all measurements below 40 s and a 25% bandwidth for measurements between 40 s and 60 s. 3. Using phase shift as a quality control measure: While the above steps help minimize measurement bias, we found that applying strict controls over phase shift between vertical and radial components is the most effective way of excluding erroneous RWE measurements.
With the margin of phase shift (70°-90°) we applied, inaccurate RWE measurements were minimized because phase shift in time windows dominated by higher mode energy always fell outside the 70°-90°range.
Our synthetic tests suggest that higher modes lead to scatter in the RWE measurements in the order of about ±0.05, which is much smaller than the measured uncertainties (Figs. 5 and 6). Thus, at least in our case, unmodelled wave propagation effects (Ferreira and Woodhouse, 2007a;Pedersen et al., 2015) might be more important than effects of contamination from higher mode energy.  Berbellini et al. (2016Berbellini et al. ( , 2017 show that RWE has good sensitivity to V s in the uppermost 30 km or so, enabling the inversion of RWE for crustal shear wave velocity structure. Since RWE is non-linearly related to Earth structure, we inverted our RWE measurements using the Neighbourhood Algorithm (NA) of Sambridge (1999) -a fully non-linear Monte Carlo technique -for the shear wave velocity structure beneath each station, following a similar approach to Berbellini et al. (2017). This is a direct search technique that avoids calculating partial derivatives as in classical linearized inversion schemes. The NA samples the model space searching for the best-fitting models in a self-adaptive efficient way. It produces an ensemble of models that minimize a cost function, whose variability give empirical information about model uncertainty. We solved the forward problem (theoretical calculation of RWE) using a normal mode formalism (Herrmann, 2013) -a fast and efficient method for calculating normal mode eigenfunctions and eigenvalues numerically in 1-D layered media. We then minimize the following cost function:

Inversion Scheme
where d i obs are the observed data, g i (m) are theoretical RWE values calculated from the normal mode's eigenfunctions for the sampled model m, σ i D is the variance of the measurement, N m is the number of measurements, N i is the number of layers in the parameterization, A is a scaling factor, and v s is the shear-wave velocity.
The first term of the cost function is the misfit between observed data and theoretical values computed using the normal mode formalism for the candidate model. We assume here that the measurements have uncorrelated Gaussian variance σ i D . The second term is a regularization constraint to avoid unrealistic models that better fit the data. We assume here that realistic models have smoother profiles than unrealistic ones. We define a "roughness function" using a second-order derivative operator following the approach of Constable et al. (1987). We multiply the roughness function by a scaling factor A -defined empirically -and by the number of measurements performed. A is obtained by trial and error from synthetic tests. The larger the value of A is, the greater the estimated misfit becomes. Conversely, the smaller the value of A is, the more unrealistic the models are. We used A = 0.0001 because we found a good balance between data-fit and the soundness of models when using this value. By minimizing the cost function, we obtained the best-fitting model, and an ensemble of models is obtained by selecting all models within a 20% threshold from the minimum of the cost function.
The main advantage of using the NA is that it does not need an a priori model; the model space is sampled inside the parameter range chosen and it does not depend strongly on subjective choices. This method provides reliable results especially when there is no detailed a priori knowledge of the structure beneath a given station, which is the case for a majority of stations in Portugal.
The parameterization of the crustal model is one of the most complex aspects because the inversion is strongly non-linear. In the absence of detailed complementary information (e.g. depth of discontinuities, average velocity) beneath a majority of stations in Portugal, we tested different parameterizations (see also Berbellini et al., 2017 for more details) and chose to use the simple 4-layer parameterization proposed by Lin et al. (2012) because it produced consistent and stable results for all stations considered. In this parameterization, we fixed the Moho depth beneath a given station based on the receiver function analysis of Dündar et al. (2016) because the latter study has the best receiver function station coverage for Portugal to-date. Also, as mentioned earlier, synthetic tests carried out by Berbellini et al. (2017) showed that, for the wave periods considered in this study, RWE data have good resolution at crustal depths. Based on this, we parameterized the crust as follows: The thickness of the uppermost layer of the crust was fixed to 3 km, and we accounted for topography by adding station elevation above the Mean Sea Level (MSL), whereby we increased the total thickness of this layer at most to~4 km. The thickness of the second-fromthe-top layer was fixed to 8 km. Finally, the thickness of each of the two bottom layers was determined by: We only inverted RWE for V s and obtained the corresponding Pwave velocity (V p ) and density (ρ) based on the empirical scaling relations of Brocher (2005).

RWE measurements
RWE measurements of stations from permanent networks (PM, LX, GE, and WM; see Fig. 2) are more numerous than those from the temporary WILAS network (8A) because the 8A network operated for b 2 years and was deployed in stages (Dias et al., 2011). Nonetheless, we were able to obtain stable RWE estimates, with a majority of 8A stations recording well over 20 measurements per period, which is statistically robust for RWE measurements (Li et al., 2016). The number of measurements by station retrieved following all our quality control tests is given in Table S1 in the Supplementary information. Although stations PW14 and PW22 consistently recorded b20 measurements across a majority of periods, we did not find anomalous RWE results or uncertainties for these stations compared to all other stations, and thus we preferred to use them for completeness. We consider this a result of the robustness of our measurement procedure. Table S2 given in Supplementary information summarizes our RWE measurements as a function of wave period and geologic unit. The uncertainties we Fig. 6. Synthetic RWE measurements for the same source parameters used in Fig. 4 except for hypocentral depth. Synthetic seismograms were computed using PREM and by incorporating the first 20 higher modes and the fundamental mode. Measurements contaminated by higher mode energy were automatically discarded by our measurement scheme. For example, T = 30 s and 35 s RWEs corresponding to an event depth of 20 km is contaminated by higher mode energy and were automatically discarded. The black solid curve indicates the fundamental mode RWE prediction for PREM. RWE scatter is about ±0.05 about this predicted curve. estimate are consistent with recent studies (e.g. Lin et al., 2012;Li et al., 2016;Chong et al., 2016). As given in this table, the average deviation of phase shift remains b10°from 90°except at 20 s, where the deviation appears to be slightly over 10°. The observed deviation probably results from sensitivity of RWE to short-scale lateral heterogeneity in the vicinity of stations and/or effects of wave interference, which are likely larger for T~20 s than for the other wave periods. Indeed, assuming V s = 3.5 km/s at a depth of 3 km, this would indicate a first Fresnel zone radius of about 10 km for T~20 s Rayleigh waves, which is shorter than the large-scale geological features shown in Fig. 1. We show examples of RWE measurements as a function of wave period for selected stations in Fig. 7 (see also Fig. S3 in Supplementary information). These indicate that the measured RWE curves deviate significantly from those predicted for reference Earth models PREM and AK135 (Kennett et al., 1995), which provides motivation for inverting RWE for more accurate models of Earth structure.
We investigated whether our RWE measurements depend on back azimuth, epicentral distance, source depth, magnitude and focal mechanism (see Figs. S4-S8 in Supplementary information). None of the stations exhibited discernible patterns with respect to these parameters at any period and the observed scatter in RWE is comparable with previous studies (Yano et al., 2009;Lin et al., 2012;Li et al., 2016), being likely related to complex Rayleigh wave propagation effects (Sexton et al., 1977;Pedersen et al., 2015). These observations are compatible with the theoretical prediction that RWE depends only on the structure beneath any given station, which gives us additional confidence that our dataset of RWE is suitable for inversions for V s structure.
In Fig. 8 we show the geographical distribution of RWE measurements at three illustrative wave periods, which exhibit systematic variations. For instance, we observe relatively higher RWE indicating lower seismic velocities at shallow depth in coastal basins, whereas further inland is characterised by lower RWE (higher V s at shallow depths). In addition, RWE tends to increase with period of measurement. A clear correlation of RWE with geological units emerges when the RWE point measurements obtained beneath each station are interpolated on to a 0.1°×0.1°uniform grid with Kriging (Fig. 8b). With this interpolation, we observe consistently higher RWE associated with the Lusitanian Basin (LB) and the Lower Tagus-Sado Basin (LTSB) and lower RWE further inland, particularly beneath the Galicia-Tras-os-Montes Zone (GTMZ). The central Iberian Zone (CIZ) is characterised by relatively higher RWE. Note that the two basins (LB and LTSB) and CIZ have better station coverage relative to other geologic units. Also, one must consider that our station distribution is not uniform, and coverage in the Ossa Morena Zone (OMZ) is poor. Thus, we do not observe sharp RWE variations across tectonic sutures. The segment of the Tagus-Badajoz-Cordoba Shear Zone (TBCSZ) that separates CIZ from OMZ and the ophiolite series separating OMZ from the South Portuguese Zone (SPZ) are considered two prominent sutures in Southwest Iberia. Although smeared, the separation between CIZ and OMZ appears to be somewhat discernible, whereas that between OMZ and SPZ is not clear.

Velocity model
Random sampling of the model space with the neighbourhood algorithm produced a large ensemble of 1-D V s models for each station. We found that solutions converge within about 1000 iterations at most, and in many cases this number is b800. Taking this fact into account and to minimize computational costs, we randomly sampled 4437 models per station. In Fig. 9, we show the results from inversions for six illustrative stations, where all models within a 20% threshold from the minimum of the cost function are shown in green (see also Fig. S9 in Supplementary information). On average, N3700 models occupying a narrow region in the model space fall within the 20% threshold for any given station. This large ensemble of models results from the fast convergence we achieved in our inversion. In Fig. 10, we show examples of RWE curves predicted for our minimum misfit models, which show excellent fit to the observations. We interpolated 1-D shear-wave velocity profiles on to a 0.1°×0.1°u niform grid by applying Kriging. This result is presented in Fig. 11. As expected, the velocity maps follow the spatial characteristics presented in our interpolated RWE map (Fig. 8). The uppermost~3 km of the LB and LTSB is characterised by lower V s of~3 km/s (V p~5 km/s) and it Fig. 9. V s profiles inverted from RWE curves for example stations. An ensemble of models is selected based on a 20% misfit threshold above the minimum misfit value. The minimum misfitmodel is represented by the black line, whereas the ensemble of models is represented by green lines. The station name, minimum misfit value, and the number of models in the ensemble are also shown. See also Supplementary information. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) gradually increases to about~3.75 km/s (V p~6 .25 km/s) in the bottom-most~7-10 km of the crust. On the other hand, the GTMZ is characterised by higher velocities, where V s in the uppermost layer is on the order of 3.75 km/s (V p~6 .5 km/s), and there appears to be a slight drop in velocity to about~3.5 km/s (V p~5 .5 km/s) in the underlying layer. This feature appears to be present in some regions (basin and metamorphic). The V s picks up from the third layer to the bottommost layer, where it is in the order of~4.0 km/s (V p~7 km). Other tectonic units exhibit velocities intermediate to that of the basins and the GTMZ. Note that we have no station coverage in the Algarve Basin (AB) in the southernmost part of our study region, and thus velocity in AB is not constrained. Also, the interior of OMZ and SPZ units is not well sampled. This limited sampling is due to eliminating stations PBAR, PVAQ, and PCVE. These stations had quality issues as we mentioned earlier. With our station distribution, we are only sampling OMZ and SPZ at their margins. Thus, V s remains largely unconstrained in the interior of these units, and corresponding geological interpretation must be made with care.

Discussion
Our interpolation of RWEs and V s profiles correlate well with the geology in Portugal. LB and LTSB Basins (LB, LTSB) whose formation is associated with rifting in the Mesozoic show relatively low seismic velocities particularly west of PTFASZ, where the sediment cover is expected to be about 2 km (Carvalho et al., 2016). We also found that the Central region of CIZ exhibits low seismic velocities at times comparable to those in the basins. This is somewhat surprising given that the CIZ is associated with older tectonostratigraphic units of the Variscan Orogeny and is expected to exhibit higher velocities. While there is a systematic reduction in the number of measurements that were made with temporary stations, there is no systematic variation of RWE or of the uncertainty between permanent and temporary stations. In fact, the temporary stations that we used from this region have the highest number of measurements from similar stations in Portugal. Our quality control steps are consistent and uniform and there is no reason to believe that our method affected measurements in this region in a different way than those in any other region. Thus, we confidently exclude the possibility that the observed signal is due to biased RWEs. Alternatively, Dündar et al. (2016) presented a revealing detail that might explain this anomaly. They found that the crustal thickness reduces by about 2 km in the same region that lower V s is observed. This might imply an episode of extension that this region has undergone in the post-Paleozoic period possibly simultaneous to the opening of the Atlantic in the Mesozoic. Such a tectonic event could in turn modify velocities as observed if the mechanical strength of the crust is changed by geologic processes such as intrusion, faulting, and the introduction of water (e.g. Chen et al., 2016). This is a hypothesis that needs to be tested by other methods including new analysis of surface geologic observations. Indeed, the central CIZ is amongst the least studied regions in Fig. 10. Fundamental mode RWE predictions (black curves) for the best-fitting V s models of stations given in Fig. 9. The median and error bars of observed RWE are also shown. For stations MTE, PESTR, and PMAFR, RWE predictions for V s models inverted from receiver functions by Morais et al. (2015) and Salah et al. (2011), respectively, are also shown (dashed blue curves). The misfit between these curves and ours at shorter periods is clearly observed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Portugal, and thus further geophysical and geological studies of this region are needed to better understand its tectonic evolution. Moho depths obtained by Mancilla and Díaz (2015) for the same region do not show a systematic reduction in crustal thickness, which is perhaps due to the fact that they use a smaller number of stations than Dündar et al. (2016) and average structure over entire western Iberia.
No published countrywide crustal V s model exists for Portugal based on passive seismic data to-date. The EPcrust model (Molinari and Fig. 11. V s at different depths for each station (triangles) interpolated on to a 0.1°×0.1°uniform grid. See main text for specifics of parameterization and interpolation. The colour within a triangle indicates corresponding V s of that particular station. Basins exhibit lower velocities than metamorphic belts. Note that station coverage in OMZ and SPZ is poor and, therefore, the interpolation has significant uncertainties in that region. , which is being used on a regular basis for regional scale analysis, is compiled from interpolating a collection of local models predicted for the European Plate. The only information included in EPcrust describing the structure beneath Portugal is Moho depth based on the work of Díaz and Gallart (2009) that used results from active seismic experiments. This particular Moho depth model itself is an interpolated one, whose source information is restricted to Southernmost and Northernmost Portugal. EPcrust, thus, might not be a representative model of crustal velocities beneath Portugal. This is evident in Fig. 12, where we show V s as given in the EPcrust model. Not only are the seismic velocities in the EPcrust model poorly correlated with the geology, but also the magnitude of velocities in the basins and the metamorphic belts are inverted at lower crustal depths, which seems unrealistic. Similar limitations are expected in other European-scale models (e.g. Ziegler and Dezes, 2006;Tesauro et al., 2008;Grad et al., 2009;Artemieva and Thybo, 2013). More recently, Silveira et al. (2013) produced group velocity perturbation maps for the entire western Iberia at wave periods 8 s, 10 s, 20 s, and 30 s. They seem to resolve metamorphic belts and sedimentary basins in Portugal at T~8 s and 10 s, although resolution at 20 s and 30 s does not seem to be appropriate to make a meaningful comparison with our RWE maps. They resolve a SSE-NNW oriented high velocity elongated body in Northwestern Portugal at 20 s. Because this feature is not observed at other neighbouring periods and it is clear from their ray path density maps (see Fig. 4 in their study) that this region has the poorest ray coverage, interpreting this particular body as a robust feature is questionable. Neither our RWE maps nor our V s models are consistent with such a geologic feature.
The average Moho depth is about 30 km with little lateral variation beneath Portugal. The most significant deviation of Moho depth from this average is observed in LB and LTSB, which is about −4 km (Salah et al., 2011;Díaz et al., 2016;Dündar et al., 2016). For comparison, we also inverted RWE for V s using the Moho depths given in the CRUST1.0 global crustal model (Laske et al., 2013), and the resulting best-fitting models are very similar to those presented in Section 5. This is not surprising given that the difference between Moho depth estimated by Dündar et al. (2016) and CRUST1.0 is less than ±2 km for any given station on average. RWE is not very sensitive to these small variations of crustal thickness.
Amongst the seismic stations used in this study, there are only three stations (MTE, PESTER, PMAFR) for which 1-D seismic wave speed models have been previously published, which can be compared with our results (Fig. 13). Morais et al. (2015) inverted P and S wave receiver functions for structure beneath MTE to a depth of about 300 km. While receiver functions have stronger sensitivity to seismic discontinuities than RWE, RWE is more sensitive to shear wave speed, particularly in the uppermost few kilometres of the subsurface. Moreover, Morais et al.'s work was carried out at the higher end of the frequency range used in our study. Thus, the two types of analysis are largely independent yet complementary. Fig. 13 shows that there is a good agreement between the two types of models when differences in parameterization are considered.
On the other hand, Salah et al. (2011) simultaneously inverted P wave receiver functions and Rayleigh wave phase velocity data for 6 stations in Portugal, of which 5 (PACT, PAZA, PMST, PSSC, PMAFR) were located in the LB and LTSB, with the other (PESTR) located in the OMZ. Unfortunately, we could only use data from 2 of these stations (PMAFR, PESTR; network PM) due to issues with data quality and availability. Below~20 km depth, our V s models for stations PMAFR and PESTR agree well with the results from Salah et al. (2011). For shallower depths, however, there are substantial differences. In particular, their models show sharp jumps in velocity in the uppermost 2 km, which may be due to the fine parameterization used in the inversion (e.g., the uppermost 6 km were represented by 1 km-thick layers). It is unlikely that receiver functions and/or surface wave phase velocities have strong sensitivity to such fine-scale structure (e.g., Lebedev et al., 2013;Chong et al., 2016). Fig. 10 shows clearly that, while the model obtained by Morais et al. (2015) leads to a reasonable fit of the RWE observations at station MTE, the models obtained by Salah et al. (2011) for stations PESTR and PMAFR have a poor fit to short period RWE data.
To further analyse our models beneath stations PMAFR in the LB and PESTR in the OMZ, we compared models based on reflection seismology (Matias, 1996) with ours (Fig. 13). Beneath PMAFR, prediction of reflection seismology is very close to that of our study in the sedimentary column, whereas the reflection model is more consistent with that of Salah et al.'s (2011) at mid-crustal depths. On the other hand, the lower crustal structure in Salah et al. (2011) agrees well with our study, and there are discrepancies with the reflection model. A similar pattern is observed for structure beneath PESTR. These differences could largely be accounted for by considering different sensitivities of the data used to invert for structure and variations in parameterizations. Jointly analysing RWE and other data types in the future could possibly reduce the non-uniqueness in the inversions. In general, we found that the LTSB and LB are characterised by lower V s at all depths than that in regional models such as EPcrust by about 0.25-0.5 km/s. Processes similar to those that might have modified crustal strength in the CIZ such as some combination of temperature and compositional variations, introduction of water, and faulting (e.g. Chen et al., 2016) associated with polyphase rifting (Rasmussen et al., 1998;Carvalho et al., 2016) could explain our observations. Future work including the joint analysis of RWE and other data types will help elucidate these ideas. Díaz and Gallart (2009) provided regionally averaged seismic velocity models along three transects, of which one (see transect C in their Fig. 2) cuts across central Portugal from West to East. Although we do not have stations overlapping this transect representative of laterally averaged structure, there is a good agreement between the seismic structure predicted for it and that obtained for stations in the neighbourhood (e.g. PCAS and MTE). For instance, they predict V p between 5 and 6 km/s in the uppermost 5 km in LB, whereas prediction from station PCAS in our study for the same region and depth range is 5.5 km/s. Also, their prediction of V p = 6.0-6.1 km/s between 5 km and 10 km depth is consistent with that we obtained for station PCAS (V p =~6.1 km/s). Note, however, that Díaz and Gallert's model is based on dense sampling in the Northernmost and Southernmost regions of Portugal with only one transect sampling parts of central Portugal, where our station coverage is greatest. Thus, our data illuminate crustal structure in regions largely ignored by previous active seismic studies.
Our results highlight short length-scale features, which are a direct result of polyphase deformation the crust has undergone. Sexton et al. (1977) hypothesized that RWE is particularly suitable for constraining seismic velocities of thick sedimentary layers, for which excellent support is found in recent studies (Lin et al., 2012;Berbellini et al., 2016Berbellini et al., , 2017Li et al., 2016). The onshore basinal structures in Portugal are not as coherent or extensive as in other regions because they have undergone intense deformation since the Paleozoic. Where the basinal structures are coherent and somewhat simple, RWE measurements exhibit a clear signature (e.g. Southern LB and LTSB).
On the other hand, we do not observe a clear signature of the Porto-Tomar-Ferreira-Alentejo Shear Zone or the Tomar-Badajoz-Córdoba Shear Zone (Fig. 1). While sharp discontinuities might have been smeared due to interpolation of V s , it is also possible that such sharp transitions are not present at depth beneath Portugal. Indeed, intense Paleozoic and Mesozoic deformation may have equilibrated the mechanical strength of metamorphic units. This is consistent with recent estimates of V p /V s ratios of metamorphic terrains that show little-tono variation across GTMZ, CIZ, OMZ, and SPZ (Dündar et al., 2016).

Conclusions
We have built a first countrywide V s model of crustal structure beneath mainland Portugal based on Rayleigh Wave Ellipticity (RWE) data. Our RWE measurements show a good correlation with the geology in Portugal, where recent extensional structures exhibit low RWE as opposed to higher RWE associated with older metamorphic belts. Our high quality RWE measurements enabled us to invert them for crustal shear wave speed using a fully non-linear Monte Carlo inversion scheme. Being a single station technique allowing the characterisation of crustal seismic structure beneath a station, RWE is particularly well suited for regions with low or uneven seismic data coverage, such as Portugal. Our V s model shows low shear wave speeds associated with sedimentary basins, notably with the Lusitanian Basin (LB) and the Lower Tagus-Sado Basin (LTSB). Parts of the Central Iberian Zone also show low V s , which we hypothesize to be related to extensional deformation in the post-Paleozoic period possibly simultaneous to Mesozoic rifting. On the other hand, the Galicia-Tras-os-Montes-Zone (GTMZ), a metamorphic belt, exhibits higher seismic velocities. Despite having sampled Fig. 13. Comparison of V s structure obtained from receiver function studies with our best-fitting models. The model for MTE is from Morais et al. (2015). Models for PESTR and PMAFR are from Salah et al. (2011). Also shown are models based on reflection seismology for PESTR and PMAFR from Matias (1996). The misfit of our best model is shown below the station name.
predominantly at their margins, the Ossa Morena Zone (OMZ) and the South Portuguese Zone (SPZ), on the other hand, appear to have seismic velocities intermediate to those of the basins and the GTMZ. Our new crustal model provides useful information to better understand the complex tectonic evolution of Portugal. Moreover, the new constrains on the shallow crustal structure of the region are valuable for future seismic hazard and ground-shaking simulations. Finally, future deployments of regularly-spaced permanent and temporary seismic stations along with other geophysical and geological analyses (e.g. Gravity, Magnetotelluric, Isotopic analyses) will help obtain more detailed images of the subsurface structure, and test hypotheses presented in our study.

Acknowledgement
This work is a part of the AQUAREL project (PTDC/CTE-GIC/116819/ 2010) funded by the Fundação para a Ciência e a Tecnologia (FCT), Portugal. AMGF also thanks support from NERC grant NE/K005669/1. We thank Instituto Português do Mar e da Atmosfera (IPMA), notably Fernando Carrilho, and Dina Vales for sharing data from the PM network and for fruitful discussions. Thanks is also due to Graça Silveira and Nuno Dias for sharing data from the 8A and WM networks and for many discussions that improved this manuscript. We also thank Luís Pinto and Luis Matias for their assistance with data at various stages of the analysis and for discussions on the geology and tectonics of the region. Susana Custódio and Susana Vilanova shared seismicity data for Portugal, and Iolanda Morais shared her receiver function results, which is very much appreciated. The data for LX stations were downloaded from IRIS-DMC, while the data for GE network was downloaded from GEOFON. Standard data processing was handled with the Seismic Analysis Code (Goldstein and Snoke, 2005). All figures were made using Generic Mapping Tools (Smith and Wessel, 1990;Wessel et al., 2013). RWE predictions for best-fitting models in figures were calculated using the MINEOS package (Masters, Woodhouse, Gilbert, and Backus) obtained from Computational Infrastructure for Geodynamics (CIG) at www.geodynamics.org. We acknowledge discussions and mobility supported by COST Action ES1401-TIDES. This manuscript benefitted from comments made by two anonymous reviewers and JFBD Fonseca. Thanks are also due to them.