The crustal structure beneath The Netherlands derived from ambient seismic noise

This work presents the ﬁ rst comprehensive 3-D model of the crust beneath The Netherlands. To obtain this model, we designed the NARS-Netherlands project, a dense deployment of broadband stations in the area. Rayleigh and Love wave group velocity dispersion was measured from ambient noise cross-correlations. Azimuthally anisotropic group velocity maps were then constructed and the isotropic part was used to determine a shear wave speed model that includes the e ﬀ ects of radial anisotropy. Employing the Neighbourhood Algorithm for the depth inversion, we obtained probabilistic estimates of the radially anisotropic model parameters. We found that the variations in the thickness of the top layer largely match the transition from sediments of Permian age to those of Carboniferous age. Regions of high faulting density such as the West Netherlands Basin and Roer Valley Graben are recognized in our model by their negative radial anisotropy ( V SH − V SV <0). The model has a mid-crustal discontinuity at a depth of around 13 km and the average Moho depth is 33 km, with most of its depth variations within 2 km. Speci ﬁ cally, a localized Moho uplift to a depth of 29 km is found within Roer Valley Graben, in the Campine region in Belgium. Furthermore, our Rayleigh and Love wave group velocity data at periods of around 20 s show evidence for azimuthal anisotropy with a NW-SE fast direction. This anisotropy is likely related to NW-SE rock fabric in the lower crust thought to originate from the Caledonian


Introduction
The sedimentary structure of The Netherlands, down to about 5 km depth, is well mapped based on ample seismic reflection and borehole data (TNO-NITG, 2004;van Dalfsen et al., 2006;Duin et al., 2006;van Dalfsen et al., 2007). In contrast, the deeper crustal structure is only poorly constrained. Limited seismic evidence is obtained from deep seismic reflection data that show a transparent upper crust and a reflective lower crust (Remmelts and Duin, 1990;Duin et al., 1995). Improved knowledge on the deep crust is, however, essential to better understand the tectonic evolution and the processes of crustal deformation that shaped the Dutch subsurface. For instance, it is unclear whether faults that are recognized in the sedimentary sequence continue down into the basement, and, if they do, to what depth. Such knowledge is essential for seismic hazard assessment. Based on deep seismic reflection data (Duin et al., 1995), it has been suggested that the Moho is shallowing beneath basins and deepening beneath structural highs, but this has not been confirmed by other data.
The low level of seismicity across The Netherlands does not allow for imaging the crustal structure of The Netherlands with local earthquake data. Instead, we used the technique of ambient noise tomography to obtain a 3-D shear wave speed model of the crust. We deployed a temporary array of seismometers, the NARS-Netherlands project (http://www.geo.uu.nl/Research/Seismology/nars.html). The network consists of 19 broadband stations deployed from 2008 to 2012, complemented by other seismic stations in and around The Netherlands (see Fig. 1). Interstation Rayleigh and Love wave group velocity dispersion curves were measured from ambient noise cross-correlations, which in turn were used to construct group velocity maps. A 3-D, radially anisotropic shear wave speed model of the Dutch crust was then obtained by inversion of the group velocity maps using the Neighbourhood Algorithm (Sambridge, 1999a,b), thus providing probabilistic estimates of the final model parameters.
In the following we will first present a brief overview of the geological evolution of The Netherlands. Then we will explain how the 3-D shear wave speed model was obtained from the group velocity measurements using the ambient noise data, the construction of the group velocity maps and the probabilistic depth inversion. We will finally show that our model is generally consistent with previous findings but also provides new ones, especially for the deeper parts of the crust and the anisotropy.

Geological evolution of The Netherlands
The flat Dutch surface topography hides a complex structure of folded and faulted sedimentary layers. These sedimentary layers have been mapped from the surface to depths of 3-6 km based on seismic surveys and borehole data acquired by oil companies (TNO-NITG, 2004;Duin et al., 2006). Here, we present a brief summary of the geology and the main geologic events since the Late Carboniferous. A more extensive overview can be found in Wong et al. (2007).
The Phanerozoic geology of The Netherlands is strongly influenced by the plate tectonic events related to the formation and break-up of the Pangea supercontinent. The area itself, however, was not directly affected by the Caledonian, Variscan and Alpine orogenies.
From the Cambrian to the Silurian, The Netherlands was part of Avalonia, a micro-continent that was rifted off from the supercontinent Gondwana in the southern hemisphere (Ziegler, 1990;Pharaoh, 1999). While moving northwards, Avalonia first collided with Baltica and then with Laurentia, and together they formed the new large Laurussian continent after the Caledonian orogeny. It has been suggested that the present dominant trend of NW-SE faults in The Netherlands dates back to the suturing of Avalonia and Baltica, and that these faults were reactivated during later tectonic events (Ziegler, 1990;de Jager, 2007). The Caledonian collision was followed by a phase of extension that lasted from Middle Devonian to Early Carboniferous. The post-Caledonian sediments are hardly deformed, but the NW-SE faults that were generated during this phase have later been reactivated, producing the main structural pattern present in the Dutch subsurface.
The collision of the Laurussia and Gondwana continents, leading to the formation of the Pangea supercontinent, caused the Variscan orogeny during the Carboniferous. The Variscan Deformation Front moved northward but it did not reach The Netherlands. Instead, a foreland basin developed with Carboniferous sediments showing only minor effects of tectonism (de Jager, 2007;TNO-NITG, 2004).
The break-up of the Pangea supercontinent formed the next major plate tectonic phase, characterized by rifting and tectonic activity from Early Triassic to Early Cretaceous. The opening of the North Atlantic and crustal extension in the North Sea affected the "Dutch" foreland basin north of the Variscan Mountains. While basin sedimentation continued during the Triassic and Jurassic, smaller rifts developed with varying subsidence rates during several episodes of rift tectonics (TNO-NITG, 2004;de Jager, 2007). After the North Atlantic had opened, rift tectonics ceased in the Early Cretaceous and was followed by regional subsidence.
The collision of Africa and Europe at the end of the Cretaceous produced a change to a compressive regime during the Alpine orogeny. In The Netherlands, this resulted in several compressional phases from the Late Cretaceous to the Neogene. Pre-existing faults were reactivated in a reversed sense such that former basins were elevated and former highs started to subside (TNO-NITG, 2004). During the Neogene, a period of relative tectonic quiescence, The Netherlands experienced successive periods of rising and falling sea level. The Roer Valley Graben is the only structure that showed continuous subsidence from Early Miocene to the present. It was, and continues to be, seismically active.

Group velocity measurements
We used ambient noise to measure interstation group velocity dispersion curves within the crust beneath The Netherlands. This technique relies on seismic interferometry which shows that the cross-correlation of ambient seismic noise recorded at two receivers produces an empirical Green's function (EGF), effectively turning one receiver into a virtual source. Seismic interferometry is widely used and reviews may be found in Wapenaar et al. (2010a,b), Galetti and Curtis (2012) and Snieder and Larose (2013). Here we used data from 18 stations of the temporary NARS-Netherlands network which were complemented with data from 9 stations of The Netherlands Seismic Network, 2 stations of the Belgium Seismic Network and 2 stations of the German Regional Seismic Network (see Fig. 1). The data of our study span the period from October 2008 until the middle of 2012. The EGFs were constructed closely following the procedure set out by Bensen et al. (2007). Each continuous station record was cut into daily time series. The instrument response was removed to obtain ground velocity and filtered with a band-pass between 5 and 50 s. We used one-bit temporal normalization to reduce the effect of earthquakes and we normalized the spectra in the frequency domain. For each pair of stations we rotated the time-series of the two horizontal components to the transverse and radial directions of the inter-station azimuth, cross-correlated and stacked the time series. The Rayleigh wave data were obtained from the vertical component cross-correlations, the Love wave data from the transverse components.
From all possible inter-station pairs, 142 pairs were left over after discarding the ones with interstation distances smaller than 80 km (∼3 wavelengths at a period of 10 s) and after discarding the EGFs of poor quality identified visually. The bandpass filtered (between 10 and 35 s) EGFs for Rayleigh and Love waves as a function of interstation distance are shown in Fig. 2 (a) and (b), respectively. For this period range, the dispersive fundamental mode surface wave trains are clearly visible. We found that, in general, the positive and negative parts of the crosscorrelations are not symmetric. This is to be expected due to the localization of the noise sources around The Netherlands (Kimman et al., 2012). For measuring group velocity curves, we selected the part of the EGF with the largest amplitude.
Frequency-time analysis (FTAN) was used to determine inter-station group velocities (Dziewonski et al., 1972;Levshin et al., 1972;Herrin and Goforth, 1977). Here we used the semi-automatic approach developed by Meier et al. (2004) with a time-variable filter (Landisman et al., 1969) to reduce artefacts in the dispersion measurements. First, group velocity measurements were made without a time variable filter. From all these measurements the average group velocity curve for The Netherlands was obtained. This reference group velocity was then used to perform time-variable filtering to the EGFs before measuring the final group velocities. Examples comparing the performance of FTAN without and with time-variable filter are shown in Fig. 3. The figure illustrates that incoherent noise or higher mode signal can distort the selected group velocity curve when time-variable filtering is not used.
To determine the uncertainty on these dispersion measurements, we compared the noise measurements to interstation measurements from teleseismic earthquake data for some suitable path azimuths. The agreement was very good and, using the earthquake data as a reference because they have a stronger signal, we established that the average uncertainty for the measurements on the EGFs is 0.03 km/s for Rayleigh wave group velocities and 0.06 km/s for Love wave group velocities.

Group velocity maps
In order to obtain Rayleigh and Love wave group velocity maps for selected periods, the area of study was parametrized into equally sized rectangular blocks. An important question was whether to include azimuthal anisotropy in the description of the maps. Including azimuthal anisotropy increases the number of unknown parameters and as a consequence tuning regularization becomes more difficult. Based on theoretical considerations, Montagner and Nataf (1986) have shown that 2ψ terms, expressing a 180°azimuthal periodicity in wave speed, are dominant for Rayleigh waves, while 4ψ terms, with a 90°azimuthal periodicity, are dominant for Love waves. Trampert and Woodhouse (2003) proposed a statistical test to analyze the necessity of including azimuthal anisotropy. They concluded that global Love wave data didn't require 2ψ components while global Rayleigh data needed both, although the 2ψ component was stronger. Using the same statistical procedure, we found that our Rayleigh and Love wave dispersion measurements both required modest azimuthal anisotropy, but couldn't distinguish between 2ψ and 4ψ terms. We therefore followed the recommendation by Montagner and Nataf (1986) to parameterize the model using only 2ψ terms for Rayleigh and only 4ψ components for Love. The maps were constructed using a standard linearized leastsquares inversion developed by Barmin et al. (2001). To avoid overregularization, a balance between the number of blocks and the anisotropic components needed to be made. From several tests and considering the difference in the number of data for different periods, we found that an optimal solution was obtained by using 14 × 14 blocks (Fig. 4) and using a 4 times higher damping for the azimuthal terms compared to the isotropic ones.
The maps at all periods were used for the depth inversion. Biases would be introduced in the depth inversion if the individual group velocity maps had different lateral resolutions. We therefore chose to fix the number of independent parameters in each of the Love and Rayleigh wave maps (i.e. aim at resolution matrices with a constant trace). Considering the difference in the number of the data at different periods (Fig. 5), we preferred to be conservative and kept the trace of the resolution matrix at around 12 (roughly half the number of data for the least well sampled period) ensuring a uniform lateral resolution across all periods. This means that the actual regularization parameters changed for each map. For shorter period group velocity data, a stable solution could have been achieved using a higher number of independent parameters. We also note that for the chosen small number of independent parameters, azimuthal anisotropy was restricted to an average direction across the region and the isotropic part remained smooth. The average uncertainty of the individual group velocity measurements determined from teleseismic earthquake data allowed us to define a tolerance on the misfit function used for the inference of the group velocity maps. This resulted in average uncertainties of 0.1 km/s for Rayleigh and 0.2 km/s for Love wave maps. Fig. 6 shows the final group velocity maps for Love and Rayleigh waves at different periods. The group velocities are plotted with respect to the average for each of the maps. The magnitudes and fast directions of the anisotropic components are represented by the bars.
Focusing first on the azimuthal anisotropy, we observe that the fast Rayleigh wave directions are predominantly directed NW-SE, whereas fast Love wave directions are mainly N-S and E-W. This is most clearly seen for periods around 20 s. Several studies have shown that a medium with horizontally aligned olivine (fast) a-axes produces fast Rayleigh wave propagation in that direction, whereas Love wave propagation is fast at angles of 45 degrees with that direction (Crampin and Taylor, 1971;Maupin, 1985;Maupin and Park, 2007). The fast directions of the 20-s Rayleigh and Love wave group velocity maps thus both seem to indicate a fast NW-SE direction in the crust. This orientation corresponds to the dominant direction of faulting present in The Netherlands (e.g., TNO-NITG, 2004). Many of the NW-SE faults originated during the Caledonian and were reactivated during later tectonic events (Ziegler, 1990;de Jager, 2007). The anisotropy can therefore be explained by NW-SE oriented fabric generated during past episodes of deformation. On the other hand, the fast direction also agrees with the NW-SE orientation of maximum horizontal compression as inferred from the World Stress Map database (Heidbach et al., 2016). Preferred alignment of oriented cracks in response to the current stress field could therefore be an alternative explanation. However, crack-induced anisotropy is limited to the upper 3-6 km of the crust (Kern, 1990), which is too shallow for the 20-s group velocity data which have their main sensitivity in the lower part of the crust. Hence we suggest that the lower crustal anisotropy inferred from the 20-s group velocity maps is most likely caused by NW-SE fabric in the lower crust. The lattice preferred orientation of the minerals may have been frozen in after the Caledonian deformation, or may have recurred by fault reactivation during later tectonic events.
The isotropic part of the group velocities of Fig. 6 show coherent patterns between Rayleigh and Love waves as well as from period to period. To further interpret these maps, we inverted them simultaneously to obtain a radially anisotropic shear wave speed model of the crust as described in the next section.

Shear wave speed model
We used the Neighbourhood Algorithm (NA) of Sambridge (1999a,b) to jointly invert the Rayleigh and Love wave group velocity data at every grid point of the map. The method is based on a guided random search through the model space to obtain robust information of the chosen parameters. NA consists of two stages: (1) creating an initial model ensemble used for interpolation and (2) the Bayesian inference. In the first stage of NA, at each iteration, a number of random models is drawn, interpolated using Voronoi cells and ranked according to the misfit. In the next iteration, n s new models are drawn from the n r best fitting Voronoi cells. Then a new Voronoi grid is determined and ranking is applied, after which the procedure is repeated. The second stage of NA consists of the appraisal problem which draws inferences from the initial ensemble obtained in stage 1. Using a Gibbs sampler, this initial ensemble is re-sampled to construct a posterior probability density function for the model. The posterior mean model, the model covariance, and the 1-D and 2-D marginals were determined to evaluate the model characteristics.
We used a chi-squared misfit for the Rayleigh and Love wave group velocity data in the 10 to 30 s period range. (1) Fig. 2. Cross-correlations as a function of distance for Rayleigh (a) and Love (b) wave data. The cross-correlations were bandpass filtered between 10 and 35 s.
T. Yudistira et al. Tectonophysics 721 (2017) 361-371 where σ i is the data uncertainty and N is number of the data. An average uncertainty of 0.1 km/s for Rayleigh wave and 0.2 km/s for Love wave dispersion, as determined during the construction of the maps, was used. The Neighbourhood Algorithm requires the forward calculation of group velocity dispersion curves. We used the codes rayleigh.f and love.f (Nolet, 2008), which provides exact dispersion curves for 1-D models. For each point of the grid, the subsurface was represented by four radially anisotropic layers: three crustal layers (a top sedimentary layer, an upper and a lower crustal layer) overlaying a mantle layer with a constant velocity down to 80 km depth where it links to the reference model PREM (Dziewonski and Anderson, 1981). Radial anisotropy in each of the four layers was parameterized by the average (V s ave ) and the radially anisotropic (V s ani ) shear wave speed with V SV and V SH the speeds of the vertically and horizontally polarized S waves. We assumed that there is no difference in wave speed between the horizontally and vertically propagating P waves and that the anisotropic parameter η equals 1 (Panning and Romanowicz, 2006). P wave speed (V P ) and density (ρ) were scaled to the Voigt average of shear wave speed (V S ) using relations derived by Brocher (2005): (4) where speed is given in km/s and density in g/cm 3 . The crustal layers were allowed to vary in thickness with the bottom depth (Z) as unknown parameter. After some exploratory runs, the search boundaries of the model parameters were selected as given in Table 1. NA was run with n s = n r = 100, which gave a good balance between a broad model search and a guided search. With 5000 iterations in total 500,000 model samples were generated for each of the grid points.
The results of the 1-D models for all individual grid points were combined and are summarized in maps presented in Figs. 7-9. The maps on the left show the posterior mean values of the average shear wave speed (Fig. 7), anisotropic shear wave speed (Fig. 8) and the depth (Fig. 9) for each of the layers. The maps on the right show the standard deviations of the 1-D marginals for each of these parameters. The full 1-D marginals and covariance matrices for all grid points may be found in Yudistira (2015).

Interpretation of the 3-D model
The top sedimentary layer of our 3-D model varies in thickness from 2 km in the east to around 4 km in the west and northwest, with a standard deviation of approximately 0.5 km (Fig. 9). The pattern of Frequency-time representation without (middle) and with (right) time-variable filtering. Bottom row: Similar, but then for the Love wave of station pair NE003-WTSB. The white squares indicate the group velocity measurements and the dashed lines represent the reference group velocity curves. Fig. 4. The area of study was parametrized into 14 × 14 equally sized blocks. Shown are the inter-station paths corresponding to Rayleigh wave data for a period of 15 s.
T. Yudistira et al. Tectonophysics 721 (2017) 361-371 thickness variations of our model is in good agreement with the depth of the base of the Rotliegend, which is the oldest (Permian) and deepest sedimentary layer in The Netherlands that is properly mapped (TNO-NITG, 2004). The strong shear wave speed increase between the Permian and the underlying Carboniferous and Devonian sedimentary layers apparently localizes the first interface of our model. Several stratigraphic units of varying thickness and with different wave speeds are combined in the first layer. It is therefore not straightforward to interpret the variations in shear wave speed (Fig. 7). Nevertheless, the Texel-IJsselmeer High in the northwest stands out as a region of relatively high velocity combined with a large sedimentary thickness. The stratigraphy of this region is anomalous because it was an elevated structural high during much of the geological history and, as a consequence, the entire sedimentary sequence from the Permian to the Jurassic is missing (de Jager, 2007). It is therefore likely that the first layer of our model for the Texel-IJsselmeer High not only represents younger than Permian sediments, as elsewhere in the country, but also includes part of the older and higher velocity pre-Permian units.
Radial anisotropy of the first layer is mostly around 0 (Fig. 8). Positive effective radial anisotropy (V SH − V SV > 0) is produced by horizontal layering (e.g., Babuska and Cara, 1991) and is considered to be the more likely type of radial anisotropy. It is striking that two regions with (resolved) negative radial anisotropy correspond to areas characterized by a high faulting density in the sedimentary sequence, the West Netherlands Basin and the northern part of the Roer Valley Graben (see Duin et al., 2006, their Fig. 3b). We suggest that the alignment of near-vertical faults in these regions produces the effective negative radial anisotropy that is identified in our model.
The upper crustal layer, layer 2 of our model, has a fairly constant posterior mean shear wave speed of around 3.1 km/s. Somewhat lower values are found for the Roer Valley Graben and the West Netherlands Basin and higher values for the Texel-IJsselmeer High (Fig. 7) indicating that these structural elements extend to greater depths. Radial anisotropy is approximately zero and appears to become increasingly more positive towards the west, but this is not fully resolved considering its standard deviation (Fig. 8). The bottom depth of the upper crustal layer is nearly everywhere around 13 km. This is in good agreement with the 12 km deep mid-crustal discontinuity of the 1-D model of Paulssen et al. (1992) for The Netherlands.
The shear wave speed of the third, lower crustal layer is around 3.9 km/s (Fig. 7). Radial anisotropy in this layer is found to be mostly positive, although this is only poorly resolved (Fig. 8). Yet, positive radial anisotropy would be consistent with the lower crustal reflectivity that is observed for the deep seismic reflection lines (Remmelts and Duin, 1990;Rijkers and Duin, 1994).
Our model has Moho depths varying between 29 and 34.5 km with a map-averaged value of 32.8 km and an average standard deviation of 1.3 km (Fig. 9). This 33-km average is somewhat larger than the 30 km that Remmelts and Duin (1990) inferred from seismic reflection data. Their study further suggested a thinner crust, down to 27 km, beneath the Roer Valley Graben. This was, however, not confirmed by Sichien et al. (2012) using Moho-reflected waves from local events. Our results are more in agreement with those of Sichien et al. (2012)  We refrain from interpreting the average and anisotropic components of shear wave speed in the upper mantle because the 30-s Love wave group velocity map is poorly resolved due to the relatively low number of Love wave data for that period. The two parameters are strongly correlated for most grid points and cannot be determined independently.
In summary, our 3-D model identifies known features of the sedimentary sequence and matches some patchy, independent observations of the deeper crustal structure. The model reveals that signatures of structures such as the West Netherlands Basin and the Roer Valley Graben can be recognized in the upper crust, but we do not resolve major lateral variations in the lower crust. Most notably, with 33 km our study finds a slightly deeper Moho than previously suggested, with lateral variations not exceeding 4 km.

Conclusions
The 3-D shear wave speed model of the crust beneath The Netherlands presented in this study is obtained from seismic noise recordings using the technique of ambient noise tomography. For this purpose, we designed the NARS-Netherlands project (2008)(2009)(2010)(2011)(2012) where 19 broadband stations were deployed in and around The Netherlands. Together with the distribution of permanent broadband stations, a dense network with an average station spacing of about 50 km was achieved. Group velocity data measured from inter-station ambient noise cross-correlations provided the basis of our 3-D model. The model includes the effects of radial anisotropy through the simultaneous inversion of Rayleigh and Love wave group velocity data. T. Yudistira et al. Tectonophysics 721 (2017) 361-371 Because the depth inversion was performed with the Neighbourhood Algorithm we also have probabilistic estimates of the model parameters.
The first findings were obtained from intermediate results, the group velocity maps which identified fast directions of azimuthal anisotropy. This effect is the strongest for the 20-s Rayleigh and Love wave group velocity data that both require NW-SE fast directions in the lower part of the crust. We suggest that the anisotropy is related to lower crustal fabric with a dominant NW-SE direction, originally imposed by   the Caledonian deformation, but potentially with additional imprints from later tectonic phases. The pattern of thickness variations of the top layer of our 3-D shear wave velocity model is in good agreement with known depth variations of the base of the Rotliegend, the sedimentary layer of Permian age. This implies a strong shear wave speed increase from the sediments of Permian age to those of the underlying units of Carboniferous and Devonian age. Furthermore, our model identifies regions of high faulting density such as the West Netherlands Basin and Roer Valley Graben by their negative radial anisotropy. Faint signatures of these basins can be recognized in the deeper parts of the model. Similarly, the expression of the Texel-IJsselmeer High seems to extend to greater depths. We find a mid-crustal discontinuity at about 13 km and the Moho depth of around 33 km. It is noteworthy that our average Moho depth is somewhat greater than the 30 km previously suggested from deep seismic reflection data. Most of the variations in Moho depth of our model appear to be limited to 2 km, although an exceptionally shallow Moho of 29 km is localized within the Roer Valley Graben, in the Campine region in Belgium.
the NARS-Netherlands stations. The Royal Netherlands Meteorological Institute, the Royal Observatory of Belgium, and the Bundesanstalt für Geowissenschaften are thanked for allowing the use of their seismic network data. We thank two anonymous reviewer for their constructive comments.
The 3-D model and the group velocity maps are available at http:// www.geo.uu.nl/~jeannot/My_web_pages/Downloads.html. Data of the NARS-Netherlands project are available through the ORFEUS Data Center. Fig. 9. Maps of the posterior mean value of the depth Z (left) and its standard deviation (right) for the first three layers. Z1 represents the base of the (sedimentary) top layer, Z2 the interface between the upper and lower crustal layers, and Z3 the Moho. The depth scales correspond to the NA search limits. RVG = Roer Valley Graben, TIJH = Texel-IJsselmeer High, WNB = West Netherlands Basin.