Articles | Volume 22, issue 1
Nat. Hazards Earth Syst. Sci., 22, 41–63, 2022

Special issue: Earthquake-induced hazards: ground motion amplification and...

Nat. Hazards Earth Syst. Sci., 22, 41–63, 2022

Research article 06 Jan 2022

Research article | 06 Jan 2022

Development of a seismic site-response zonation map for the Netherlands

Development of a seismic site-response zonation map for the Netherlands
Janneke van Ginkel1,2, Elmer Ruigrok2,3, Jan Stafleu4, and Rien Herber1 Janneke van Ginkel et al.
  • 1Energy and Sustainability Research Institute Groningen, University of Groningen, Nijenborgh 6, 9747 AG Groningen, the Netherlands
  • 2R&D Seismology and Acoustics, Royal Netherlands Meteorological Institute, Utrechtseweg 297, 3731 GA De Bilt, the Netherlands
  • 3Department of Earth Sciences, Utrecht University, Princetonlaan 8a, 3584 CB Utrecht, the Netherlands
  • 4TNO – Geological Survey of the Netherlands, Princetonlaan 6, 3584 CB Utrecht, the Netherlands

Correspondence: Janneke van Ginkel (


Earthquake site response is an essential part of seismic hazard assessment, especially in densely populated areas. The shallow geology of the Netherlands consists of a very heterogeneous soft sediment cover, which has a strong effect on the amplitude of ground shaking. Even though the Netherlands is a low- to moderate-seismicity area, the seismic risk cannot be neglected, in particular, because shallow induced earthquakes occur. The aim of this study is to establish a nationwide site-response zonation by combining 3D lithostratigraphic models and earthquake and ambient vibration recordings.

As a first step, we constrain the parameters (velocity contrast and shear-wave velocity) that are indicative of ground motion amplification in the Groningen area. For this, we compare ambient vibration and earthquake recordings using the horizontal-to-vertical spectral ratio (HVSR) method, borehole empirical transfer functions (ETFs), and amplification factors (AFs). This enables us to define an empirical relationship between the amplification measured from earthquakes by using the ETF and AF and the amplification estimated from ambient vibrations by using the HVSR. With this, we show that the HVSR can be used as a first proxy for site response. Subsequently, HVSR curves throughout the Netherlands are estimated. The HVSR amplitude characteristics largely coincide with the in situ lithostratigraphic sequences and the presence of a strong velocity contrast in the near surface. Next, sediment profiles representing the Dutch shallow subsurface are categorised into five classes, where each class represents a level of expected amplification. The mean amplification for each class, and its variability, is quantified using 66 sites with measured earthquake amplification (ETF and AF) and 115 sites with HVSR curves.

The site-response (amplification) zonation map for the Netherlands is designed by transforming geological 3D grid cell models into the five classes, and an AF is assigned to most of the classes. This site-response assessment, presented on a nationwide scale, is important for a first identification of regions with increased seismic hazard potential, for example at locations with mining or geothermal energy activities.

1 Introduction

Site-response estimation is a key parameter for seismic hazard assessment and risk mitigation, since local lithostratigraphic conditions can strongly influence the level of ground motion amplification during an earthquake (e.g. Bard1998; Bonnefoy-Claudet et al.2006b, 2009; Borcherdt1970; Bradley2012). In particular, near-surface low-velocity sediments overlying stiffer bedrock modify earthquake ground motions in terms of amplitude and frequency content, as for instance observed after the Mexico City earthquake in 1985 (Bard et al.1988) as well as more recent ones (e.g. L'Aquila, Italy, 2009; Tokyo, Japan, 2011; Darfield, New Zealand, 2012). Site-response estimations require detailed geological and geotechnical information of the subsurface. This can be retrieved from in situ investigations; however, this is a costly procedure. Because of the time and costs involved, there is a lack of site-response investigations covering large areas, while the availability of detailed and uniform ground motion amplification maps is fundamental for preliminary estimates of damage on buildings (e.g. Falcone et al.2021; Gallipoli et al.2020; Bonnefoy-Claudet et al.2009; Weatherill et al.2020).

Empirical seismic site response is widely investigated by the use of microtremor horizontal-to-vertical spectral ratios (HVSRs; e.g. Fäh et al.2001; Lachetl and Bard1994; Bonnefoy-Claudet et al.2006a; Albarello and Lunedei2013; Molnar et al.2018; Lunedei and Malischewsky2015). The HVSR is obtained by taking the ratio between the Fourier amplitude spectra of the horizontal and the vertical components of a seismic recording. When a shallow velocity contrast is present, the peak in the HVSR curve is closely related to the shear-wave resonance frequency for that site. However, the HVSR peak amplitude cannot be treated as the actual site amplification factor but rather serves as a qualitative estimate (Field and Jacob1995; Lachetl and Bard1994; Lermo and Chavez-Garcia1993).

The Netherlands experiences tectonically related seismic activity in the southern part of the country, with magnitudes up to 5.8 measured so far (Camelbeeck and Van Eck1994; Houtgast and Van Balen2000; Paulssen et al.1992). Additionally, gas extraction in the northern part of the Netherlands regularly causes shallow (3 km), low-magnitude (Mw≤3.6 thus far) induced earthquakes (Dost et al.2017). Over the last decades, an increasing number of these induced seismic events have stimulated the research on earthquake site response in the Netherlands. Various studies (van Ginkel et al.2019; Kruiver et al.2017a, b; Bommer et al.2017; Noorlandt et al.2018) undertaken in the Groningen area (northeastern part of the Netherlands) concluded that the heterogeneous unconsolidated sediments are responsible for significant amplification of seismic waves over a range of frequencies pertinent to engineering interest. Although the local earthquake magnitudes are relatively small, the damage to the houses can be significant. Hence multiple studies (e.g. Rodriguez-Marek et al.2017; Bommer et al.2017; Kruiver et al.2017a; Noorlandt et al.2018) were performed on ground motion modelling including the site amplification factor for the Groningen region.

Groningen forms an excellent study area due to the presence of the permanently operating borehole seismic network (G-network). Local earthquake recordings over the Groningen borehole show that the largest amplification develops in the top 50 m of the sedimentary cover (van Ginkel et al.2019), although the entire sediment layer has a thickness of around 800 m in this region. Furthermore, van Ginkel et al. (2019) showed the existence of a correlation between the spatial distribution of microtremor horizontal-to vertical spectral ratio (HVSR) peak amplitudes and the measured earthquake amplification. This observation is in accordance with those of Pilz et al. (2009), Perron et al. (2018), and Panzera et al. (2021), who show a comparison of site-response techniques using earthquake data and ambient seismic noise analysis.

The aim of this work is to design a site-response zonation map for the Netherlands, which is both detailed and spatially extensive. Rather than using ground motion prediction equations with generic site amplification factors conditioned on Vs30, we propose a novel approach for the development of a nationwide zonation of amplification factors. To this end, we combine multiple seismological records, geophysical data, and detailed 3D lithostratigraphic models in order to estimate and interpret site response. We first select the Groningen borehole network where detailed information on subsurface lithology, numerous earthquake ground motion recordings, and ambient seismic noise recordings is available since their deployment in 2015. From this, we extract empirical relationships between seismic wave amplification and different lithostratigraphic conditions, building upon the proxies defined in van Ginkel et al. (2019).

Next, the ambient vibration measurements of the seismic network across the Netherlands are used, necessary to calibrate the amplification (via HVSR) with the local lithostratigraphic conditions. By combining the detailed 3D geological subsurface models GeoTOP (Stafleu et al.2011, 2021) and NL3D (van der Meulen et al.2013), with a derived classification scheme, a zonation map for the Netherlands is constructed.

The presented site-response zonation map for the Netherlands is especially designed for seismically quiet regions where tectonic seismicity is absent, but with a potential risk of induced seismicity, for example due to mining or geothermal energy activity (Majer et al.2007; Mena et al.2013; Mignan et al.2015). As a result, this map can be implemented in seismic hazard analysis.

2 Geological setting and regional seismicity

The Netherlands is positioned at the southeastern margin of the Cenozoic North Sea basin. The onshore basin infill is characterised by Paleogene, Neogene, and Quaternary sediments reaching a maximum thickness of ∼1800 m. Minimum onshore thicknesses are reached along the basin flanks in the eastern and southern Netherlands and locally at uplifted blocks like the Peel Block. The main tectonic feature of the country is the Roer Valley Graben, bounded by the Peel Boundary Fault in the northeast and the Rijen, Veldhoven, and Feldbiss faults in the south and southwest (Fig. 1).

Figure 1Surface geological map of the Netherlands. Older deposits comprise unconsolidated Neogene and Paleogene deposits as well as Mesozoic and Carboniferous limestones, sandstones, and shales. Modified after Schokker (2010). RF: Rijen Fault; VF: Veldhoven Fault; PBF: Peel Boundary Fault; FBF: Feldbiss Fault; RVG: Roer Valley Graben.

The Paleogene and Neogene sediments are dominated by marine clays and sands that were primarily deposited in shallow marine environments. The Quaternary sediments, reaching a maximum onshore thickness of ∼600 m, reflect a transition from shallow marine to fluvio-deltaic and fluvial depositional environments in the early Quaternary to a complex alternation of shallow marine, estuarine, and fluvial sediments in the younger periods (Zagwijn1989; Rondeel et al.1996; De Gans2007; Busschers et al.2007; Peeters et al.2016). Imprints of glacial conditions are recorded in the upper part of the basin fill by, among others, deep erosive structures (subglacial valleys, tongue basins) and glacial till (Van den Berg and Beets1987).

2.1 Surface geology

The surface geology (Fig. 1) is mainly characterised by a Holocene coastal barrier and coastal plain in the west and north and an interior zone with Pleistocene deposits cut by a Holocene fluvial system (Rondeel et al.1996; Beets and van der Spek2000). The coastal barrier consists of sandy shoreface and dune deposits and is up to 10 km wide. It is intersected in the south by the estuary of the Rhine, Meuse, and Scheldt and in the north by the tidal inlets of the Wadden Sea. The coastal plain is formed by mainly marine clay as well as peat. Although much of the peat has disappeared because of mining and drainage, thick sequences of peat (> 6 m) still occur. The Holocene fluvial deposits of the rivers Rhine and Meuse are characterised by a complex of sandy channel belt systems embedded in flood basin clays (Gouw and Erkens2007). The fluvial channel belts pass downstream into sandy tidal channel systems in the coastal plain (Hijma et al.2009).

The Pleistocene interior of the country mainly consists of glacial, aeolian, and fluvial deposits (Rondeel et al.1996; Van den Berg and Beets1987). Glacial deposits include coarsely grained meltwater sands and tills. Ice-pushed ridges, with heights up to 100 m, occur in the middle and east of the country. Aeolian deposits mainly consist of cover sands and are locally made up by drift sand and inland dunes. In the south and east of the country, sandy channel and clayey flood basin deposits of small rivers occur.

Neogene and older deposits are only exposed in the eastern- and southernmost areas of the country. In the east, these sediments include unconsolidated Paleogene formations as well as Mesozoic limestones, sandstones, and shales. In the south, the older deposits comprise unconsolidated Neogene and Paleogene sands and clays, as well as Cretaceous limestones (chalk), sandstones and shales, and Carboniferous sandstones and shales.

Figure 2Map of the Netherlands depicting epicentres of all induced (Mw 0.5–3.6, orange) and tectonic (Mw 0.5–5.8, yellow) earthquakes from 1910–2020. The diameter of the circles indicates the relative earthquake magnitude. The triangles represent the surface location of the borehole stations (blue), borehole stations with SCPT measurements (purple), and single surface seismometers (green). The triangles with red outlines depict the locations of example HVSR curves presented in Fig. 8. The inset in the northeast depicts the location of the Groningen borehole network (G-network). Here, the red triangle depicts the location of borehole G24. The 19 (Mw≥2) induced earthquakes in this panel are used for the AF and ETF computations. Coordinates are shown within the Dutch national triangulation grid (Rijksdriehoekstelsel or RD), and lat–long coordinates are added in the corners for international referencing.

2.2 Regional seismicity

The Netherlands experiences tectonically related earthquake activity in the southeast, and induced earthquakes occur in the north at shallow depths due to exploitation of gas fields. Tectonic seismicity occurs mainly in the Roer Valley Graben (yellow circles, Fig. 2), which is part of a larger basin and range system in western Europe, the Rhine Graben Rift System. At the beginning of the Quaternary, the subsidence rate in the Roer Valley Graben did increase significantly (Geluk et al.1995; Houtgast and Van Balen2000), and the rift system still shows active extension (Hinzen et al.2020). The largest earthquake recorded (Mw=5.8) in the Netherlands was in Roermond in 1992, due to extensional activity along the Peel Boundary Fault (Paulssen et al.1992). Gariel et al. (1995) quantified the near-surface amplification based on spectral ratios of aftershocks from the 1992 earthquake in Roermond. They observed great variety in ground motion amplitudes over different stations, which is the site effect of shallow sedimentary deposits.

Most induced earthquakes in the Netherlands (orange circles, Fig. 2) have their epicentre in the Groningen region due to production of the gas field. Here, reservoir compaction due to pressure depletion has reactivated the existing normal fault system that traverses the reservoir layer throughout the whole field. (Buijze et al.2017; Bourne et al.2014). Even though the magnitudes are relatively low (van Thienen-Visser and Breunese2015), the damage to buildings in the area is substantial due to shallow hypocentres and amplification on the soft near-surface soils (Bommer et al.2017; Kruiver et al.2017a).

3 Dataset

For this study, we use the seismic network of the Royal Netherlands Meteorological Institute (KNMI1993), consisting of borehole and surface seismometers distributed over the Netherlands (Fig. 2). The network includes 88 locations with vertical borehole arrays where each station is equipped with three-component 4.5 Hz seismometers at 50 m depth intervals (50, 100, 150, 200 m) and an accelerometer at the surface. The southernmost station is at Terziet (TERZ). This station consists of a borehole seismometer at 250 m depth and a surface seismometer. In Groningen, multiple boreholes have seismic cone penetration tests (SCPTs) taken directly adjacent to the borehole. The remaining triangles represent 29 locations of single surface seismometers (accelerometers and broad bands). All seismometers have three components and are continuously recording the ambient seismic field and, when present, local earthquakes. In Sect. 4 we use 19 local earthquakes of M≥2, recorded in the Groningen borehole network between May 2015 and May 2019. All data are available through the KNMI data portal (, last access: 1 June 2021).

In the construction of the site-response map, we have made extensive use of the detailed 3D geological subsurface models GeoTOP and NL3D. Both models are developed and maintained by TNO – Geological Survey of the Netherlands (van der Meulen et al.2013). GeoTOP schematises the shallow subsurface of the Netherlands in a regular grid of rectangular blocks (voxels, tiles, or 3D grid cells), each measuring 100 m by 100 m by 0.5 m (x ,y, z) up to a depth of 50 m below ordnance datum (Stafleu et al.2011, 2021). Each voxel contains multiple properties that describe the geometry of lithostratigraphic units (formations, members, and beds), the spatial variation in lithology, and sand grain size within these units as well as measures of model uncertainty. GeoTOP is publicly available from the TNO's web portal: (last access: 1 August 2021). To date, the GeoTOP model covers about 70 % of the country (including inland waters such as the Wadden Sea). For the missing areas we have used the lower-resolution voxel model NL3D, which is available for the entire country (van der Meulen et al.2013). NL3D models lithology and sand grain-size classes within the geological units of the layer-based subsurface model DGM (Gunnink et al.2013) in voxels measuring 250 m by 250 m by 1 m (x, y, z) up to a depth of 50 m below ordnance datum. To determine the depth of bedrock in the shallow subsurface, we consulted the layer-based subsurface models DGM and DGM-deep. These models are also available from the web portal mentioned above. More details on the models GeoTOP and NL3D are given in Appendix B.

4 Empirical relationships from the Groningen borehole network

The dense Groningen borehole network (G-network) provides the opportunity to derive empirical relationships between measured amplification in the time and frequency domain, estimated amplification from ambient vibrations, and the local lithostratigraphic conditions. First, we present results of three different methods to assess amplification. Next, we compare subsurface parameters with the measured amplification in order to evaluate which parameters are most critical.

4.1 Definition of amplification

The majority of site-response studies define the level of soft sediment amplification with respect to the surface seismic response of a nearby outcropping hard rock. In the Netherlands, no representative seismic response of outcropping bedrock is available. As an alternative, we define reference conditions as found in Groningen at 200 m depth, where we have many seismic recordings. These same reference conditions can be found at this depth over most locations in the Netherlands, though it can sometimes be deeper or shallower than 200 m. The corresponding elastic properties (shear-wave velocity and density) form the basis from which the ground motion amplification effect of any site with respect to the reference can be predicted. This approach is akin to the reference velocity profile that is used in Switzerland (Poggi et al.2011).

For a recording at the Earth's surface, the upward and downward waves are simultaneously recorded, which leads to an amplitude doubling. This is called the free-surface effect. If the amplification is defined with respect to a surface (hard-rock) site, both the site of investigation and the reference site experience the same free-surface effect, and there is no need to remove it in order to isolate the relative amplification. We have a reference site at depth, where no free-surface effect takes place. Hence, we need to remove the free-surface effect first from the surface recording before isolating the relative amplification.

The G-network forms a representative resource for the definition of the reference horizon parameters. Hofman et al. (2017) and Kruiver et al. (2017a) determined shear-wave velocities at all borehole seismometer levels in Groningen. From their velocities found at 200 m depth, the average is taken, resulting in a reference shear-wave velocity of 500 m/s. At this depth, the sediment density is on average 2.0 kg/m3. At 200 m depth, 95 % of the Dutch subsurface is composed of laterally extensive Pleistocene and Pliocene sediments; hence the estimated site-response and corresponding amplification factors can be applied to a large part of the country. The remaining 5 % consists of shallow (<100 m) Triassic, Cretaceous, and locally Carboniferous bedrock, and therefore these locations need to be evaluated separately.

Frequency bandwidth

Data processing is applied to a frequency band-pass filter for 1–10 Hz, covering the periods of interest from an earthquake engineering point of view. Moreover, for these frequencies, the used instrumentation (broadband seismometers, accelerometers, and geophones) has high sensitivity for ground motion.

Figure 3(a) Particle velocity Fourier amplitude spectrum measured on the radial component for the 200 m (grey) and surface seismometer (blue) for borehole G24 for a 20 s time window of the 8 January 2018 Zeerijp M 3.4 earthquake. Borehole G24 has an epicentral distance of 10 km. (b) Shear-wave velocity profile for G24 with the interval velocities from Hofman et al. (2017) in blue and the shear-wave velocity from the adjacent SCPT over the top 30 m (black line). The corresponding lithological profile is derived from GeoTOP (, last access: 15 October 2021).


Since most of the amplification is occurring in the top sedimentary layer (van Ginkel et al.2019), the corresponding resonance frequencies are covered in the used frequency filter as well. Above 10 Hz, the amplitude increase due to soil softening and resonance is counteracted by an elastic attenuation and 3D scattering. Furthermore, what exactly happens above 10 Hz is of little interest since most of the local earthquake energy is contained in the frequency band between 1 and 10 Hz. This is illustrated in Fig. 3, which presents the particle velocity Fourier amplitude spectrum of a local earthquake recorded in borehole G24. The location of this borehole is presented as a red triangle in Fig. 2.

4.2 Amplification factors

We compute an overall amplification factor from the G-network earthquake recordings by taking the ratio of the maximum amplitudes recorded at the surface and the 200 m deep seismometer, following the procedure described in van Ginkel et al. (2019). This ratio is taken for both the radial (R) and transverse (T) components, and the results are averaged. The amplitude at the surface is divided by a factor of 2 in order to remove the effect of the upward and downward waves recorded at the same time. Earthquake records are processed in the time domain for a 20 s time window and frequency band of 1–10 Hz. Next, the AF for each borehole is obtained by repeating the above procedure for 19 available M≥2.0 local induced earthquakes and subsequently averaging the values. It is determined in the time domain and therewith provides an average amplification over the applied frequency band.

Many studies (e.g. Bommer et al.2017; Rodriguez-Marek et al.2017; Borcherdt1994) model site-response amplification factors (AFs) for different periods of spectral accelerations, which tailors to engineering structures with different resonance frequencies. In this study we aim to provide an average amplification level in a broad frequency band. The choice for the specific AF frequency band was supported by evaluating amplification over a range of frequency bands. Figure 4 shows the AFs over the Groningen network for several bands. AF values are highest in the band 1–5 Hz due to the strong resonances of the Holocene infill (van Ginkel et al.2019). In the 1–10 Hz band the AFs are lower but contain considerable earthquake amplitudes (Fig. 3). When the frequency band is extended beyond 10 Hz (Fig. 4c and d), the AF values do not change significantly. We decided to pick the frequency band of 1–10 Hz as the representative AF since it is a better overall representation of amplification by the soft sediments than the higher values of 1–5 Hz, as shown in Wassing et al. (2012).

Figure 4Amplification factors (AFs) for (a) frequency bands 1–5 Hz, (b) 1–10 Hz, (c) 1–15 Hz, and (d) 1–20 Hz and corresponding standard deviations (error bars) associated with the averaged AFs over 19 M≥2.0 local earthquakes.


4.3 Empirical transfer functions

Empirical transfer functions (ETFs) represent shear-wave amplification in the frequency domain. ETFs are defined as a division of the Fourier amplitude spectra at two different depth levels. With a reference horizon at 200 m depth and the level of interest at the Earth's surface, the transfer function has both a causal and acausal part. The causal part maps upward-propagating waves, from the reference level to the surface. The acausal part maps downward-propagating waves back to the free surface (Nakata et al.2013). When describing amplification, we are only interested in the causal part. We select this causal part of the estimated transfer function and compute its Fourier amplitude spectrum to obtain a measure of frequency-dependent amplification. We use 20 s long time windows for particle velocity recordings on the radial component of the G-network seismometers. Subsequently, we average over 19 local earthquakes with magnitudes ≥2.0. This can be seen as an implementation of seismic interferometry by deconvolution (Wapenaar et al.2010). Lastly, the deconvolution results are stacked to enhance stationary contributions.

Figure 5Probability density function of 1 month of daily HVSRs and the mean HVSR (solid black line) and ETF for the borehole seismometer interval of 0–50 m (grey dashed line) and 0–200 m (grey solid line). The selected borehole sites exhibit differences in curve characteristics with (a) G24 illustrating clearly peaked curves (>4), (b) G16 illustrating medium (2–4) peak amplitudes, and (c) G51 with no pronounced peak.


From the ETF between the 200 m depth and surface seismometer (ETF200), peak amplitudes are identified. Some example ETF curves are plotted in Fig. 5. Additionally, ETF curves for the top 50 m are calculated (ETF50). The ETF curves for the different intervals show very similar peak characteristics and peak amplitudes, demonstrating that most amplifications are developed in the top 50 m of the sediment cover. Furthermore, the borehole sites with the highest ETF peak amplitudes can be linked to the local geology, as presented in van Ginkel et al. (2019, Fig. 9). Here, highest peak amplitudes are measured at sites with unconsolidated Holocene alternations of clay and peat, overlying consolidated Pleistocene sediments.

4.4H/V spectral ratios from the ambient seismic field

The Groningen surface seismometers are continuously recording the ambient seismic field (ASF), and these data are used to estimate horizontal-to-vertical spectral ratios (HVSR). In this study we focus on the second peak ( 1 Hz) in the HVSR curve, which represents shear-wave resonances at the shallowest interface of soft sediments on top of more consolidated sediments. Resonances of the complete sediment layer display a peak at lower frequencies (van Ginkel et al.2020). Above 1 Hz the noise field is dominated by anthropogenic sources. The details of the method to obtain stable HVSR curves from the ASF in the Groningen borehole network can be found in van Ginkel et al. (2020). In summary, from power spectral densities for each component, the H/V division is performed for each day. Subsequently, a probability density function is computed over one month of H/V ratios and the mean is extracted. This yields a stable HVSR curve that is not much affected by transients like nearby footsteps or traffic.

Based on the HVSR curve and peak characteristics, different criteria are defined conformable to the SESAME consensus (Bard2002): (1) clear peaked curves (HVSR amplitude >4) related to a sharp velocity contrast in the shallow subsurface. (2) HVSR peak amplitude between 2–4, associated with a weak velocity contrast. (3) No distinguishable peak and a flat curve indicate the absence of a velocity contrast in the shallow subsurface. Example HVSRs for these three criteria are plotted in Fig. 5. Its associated peak amplitudes are derived from the mean HVSR curve and presented in van Ginkel et al. (2019). The correlation between peaks on the HVSR curves and the presence of a velocity contrast at some depth is stressed in studies from Bonnefoy-Claudet et al. (2008), Konno and Ohmachi (1998) and Lermo and Chavez-Garcia (1993) and this contrast is very likely to amplify the ground motion.

Figure 6Relation between the HVSR peak amplitudes (A0) and the ETF peak amplitudes (grey) and the HVSR peak amplitudes and the AF (blue). The solid line represents the fitting function (Eqs. 1 and 2) between the HVSR peak amplitudes and the measured amplification from AF and ETF, respectively. Rsq (R-squared) represents the coefficient of determination of the fitting. Note the log-scale of the y axis.


4.5 Amplification parameters

Across the Netherlands, the ASF is continuously recorded on all seismometers, while many locations lack recordings of local earthquakes. Therefore we further investigate whether the HVSR can be used as a proxy for amplification as measured by the earthquake-derived ETF and AF. Figure 6 displays the correlation between the peak amplitudes of the HVSR and ETF200 as well as HVSR and AF. Secondly, we evaluate the subsurface seismic parameters enhancing amplification (Fig. 7).

Figure 7Each panel depicts the fitted function and the coefficient of determination (Rsq) of the AF (1–10 Hz) per G-network borehole location and the corresponding subsurface parameter. (a) Vs10, (b) Vs20, (c) Vs30, (d) Vs50, (e) depth of the velocity contrast, and (f) size of the velocity contrast.


Based on these data points, relationships are defined in order to be able to estimate amplification factors using HVSR peak amplitudes (A0), amplification factors:

(1) AF = 1.49 + 0.87 log ( 1.12 A 0 ) ,

and maximum amplification as measured by the ETF:

(2) ETF = 1.08 + 6.89 log ( 1.09 A 0 ) .

The relationship between the AF derived from the Groningen borehole locations and local site conditions is investigated in the following.

Many ground motion prediction equations which include site response consider the shear-wave velocity for the top 30 m (Vs30) as the main parameter affecting amplification (Akkar et al.2014; Bindi et al.2014; Kruiver et al.2017b; Wills et al.2000), as well as Eurocode 8 (Comité Européen de Normalisation2004). However, studies (Castellaro et al.2008; Kokusho and Sato2008; Lee and Trifunac2010) have drawn attention to the fact that using only Vs30 as a proxy for site response is inadequate, because it does not uniquely correlate with amplification. They show that amplification is defined by several parameters, including the depth and degree of the seismic impedance contrast. Furthermore Joyner and Boore (1981) introduce the shear-wave velocity ratio between the top and base layers as a proxy for site amplification, and this is further explored by Boore (2003).

In order to assess the impact of different parameters, firstly, the AF at each borehole location is fitted against several subsurface parameters. AF=a+becVs as a functional form, where a, b, and c are unknown coefficients to be fitted. This fitting is applied with averaged shear-wave velocities over various depth intervals. Vs10, Vs20, and Vs30-values are derived from SCPTs, acquired directly adjacent to 53 borehole sites in Groningen (Fig. 2). Hofman et al. (2017) derived interval velocities determined from the G-network, using seismic interferometry applied to local induced events. The velocities from this reference are used to determine Vs50. Secondly, from the SCPT data we derive the depth and size of the velocity contrast. The contrast is computed by the division of the two different velocity values bounding each 1 m interval. This division is done for each 1 m interval over the 30 m of SCPT records. The largest division value is defined as the (main) velocity contrast (VC), and corresponding depth is the depth of the contrast (zVC). Thereafter the VC values and their depth are fitted with the AF using AF=a+blog(cVC) as a functional form. This procedure is also performed for the relation between the subsurface parameters (Vs, VC) and the HVSR. The results are given in Appendix A.

Figure 7 presents the fit between the AF and the six relevant subsurface parameters. Best fit (Rsq=0.47, Fig. 7a and b) is observed between the AF and Vs10 and Vs20. This is in line with findings of Gallipoli and Mucciarelli (2009), who use the Vs10 as the main amplification parameter instead of the more common Vs30. In Groningen, the low-velocity and unconsolidated Holocene sediments have a thickness of 1 to 25 m, and below these depths the velocities increase in the more compacted Pleistocene sediments. The reduced fitting quality of Vs30 and Vs50 arises since the amplification develops mainly in the Holocene sediments (van Ginkel et al.2019). Although the fit is relatively poor, a relationship is observable between a large VC and an elevated AF (Fig. 7e). By contrast, Fig. A1 presents a good correlation between VC and HVSR. A large VC value leads to resonance in the near surface, which is expressed in high-amplitude peaks of the HVSR.

Figure 8Each panel depicts a probability density function from ambient noise HVSR curves and sediment classification profile (Sect. 6) for 16 stations of the NL network. The black line represents the mean HVSR, and the red line in the panels of T06, T010, and TERZ represents the ETF calculated from 10 local earthquakes. The colour bar in the lower right displays the HVSR probability range that is valid for all panels.


5 HVSR estimations throughout the Netherlands

Based on the good relationship between Groningen HVSR peak amplitudes (A0), ETF peak amplitudes, and AF (Fig. 6), we conclude that the HVSR can be further used as proxy for amplification. Therefore, for all surface seismometers in the Netherlands seismic network, the HVSR curves are estimated following the method described in Sect. 4.4. Figure 8 displays a selection of representative examples of HVSR curves. It also includes the sediment classification profile presented in Sect. 6 and Fig. 9. In addition, for boreholes T010, T060, and TERZ, the ETF curve is added, calculated from local earthquakes similar to the approach described in Sect. 4.3. These 16 HVSR curves illustrate the diversity in curve characteristics throughout the Netherlands. Here, we can distinguish the three types of curves as described in Sect. 4.4. The flat curves with no distinguishable peak (FR040, DRA, T060, T010, WTSB, HRKB, ROLD, and BING) are recorded at seismometers on top of outcropping Pleistocene sands (Holocene–Pleistocene aeolian and fluvial deposits in Fig. 1). Also, ALK2 exhibits no peak amplitude since this seismometer is positioned on dune sands (Holocene coastal deposits in Fig. 1). These locations are characterised by absence of a strong velocity contrast in the shallow subsurface. Selected examples of HVSR curves exhibiting clear peak amplitudes (NH01,J01, ZH030, FR010, EETW) are located at sites with a distinct velocity contrast between unconsolidated Holocene marine sediments overlying Pleistocene sediments.

In the southernmost part of the Netherlands, Cretaceous bedrock is either outcropping or situated less than 100 m deep. MAME and TERZ are examples of locations with soft sediments overlying hard rock, and the HVSR curves exhibit a clear peak amplitude.

6 Zonation map for the Netherlands

The effect of local site response on earthquake ground motion is included in the Eurocode 8 seismic design of buildings (Comité Européen de Normalisation2004, EC8). In order to estimate the risk of enhanced site response, EC8 presents five soil types based on shear-wave velocities and stratigraphic profiles. Soil type E is essentially characterised by a sharp contrast of a soft layer overlying a stiffer one. However, in our opinion, this single classification for soft sediments is rather limited, especially concerning the wide variety of lithostratigraphic conditions throughout the Netherlands. Therefore, we present an alternative, or extended, classification for ground characteristics designed to specify the large heterogeneity in site conditions than exists within Eurocode 8 ground type E.

In this section, in a few steps, the site-response zonation map for the Netherlands is derived. For this, the country is subdivided into grid cells. As a result, about 95 % of the grid cells are populated with a site-response class with a corresponding AF.

6.1 Classification scheme

The borehole ETFs confirm that most of the amplification is developed in the top 50 m (Fig. 5) of the sedimentary cover, corresponding to the findings in van Ginkel et al. (2019). Beyond 50 m depth, the Quaternary deposits mainly consist of more compacted marine and fluvial sediments. Therefore the sediment classification presented in this section uses the top 50 m with a special focus on the top 10 m. Also, the presence of a velocity contrast is used in the classification, as it was shown to have a (albeit weaker) link with amplification (Figs. 7 and A1).

Table 1Details of the NL classification. The Vs10 and velocity contrast (VC) values assigned to each class are based on the amplification relationships presented in Sect. 4 and Appendix A. For class V there are no empirical data available relating Vs10 and VC with A0 (HVSR peak amplitude), hence not determined (n.d).

Download Print Version | Download XLSX

Following Convertito et al. (2010) and from the studies by Kruiver et al. (2017a) and van Ginkel et al. (2019), the Eurocode 8 classification requires modification, caused by the heterogeneous shallow sediment composition and bedrock depth of the Dutch subsurface. Table 1 lists the criteria for the classification division defined for the Netherlands (NL classification). The NL classification is divided into five classes based on the top 50 m lithostratigraphic composition, the velocity contrast (VC), and average shear-wave velocity for the top 10 m.

For setting up the classification, we initially use A0, the peak amplitude from HVSR. We only have measured AF in Groningen, whereas we have measured A0 for many sites throughout the Netherlands (all stations depicted in Fig. 2). Moreover, we found a clear relationship between A0 and AF (Eq. 1). The relationships between Vs10, VC, and A0 are estimated from lithological conditions in Groningen, where the sites are assigned to Classes II, III, and IV. For Classes I and V we have insufficient empirical constraint on A0 and AF. Only sites with bedrock at depths shallower than 100 m fall into Class V, for which the resonance over the complete unconsolidated cover can reach frequencies larger than 1 Hz. Therewith, these sites become distinct from Classes II, III, and IV, where such resonances occur at frequencies <1 Hz. At these lower frequencies, there is no match with resonance periods of most building types in the Netherlands.

Figure 9Sediment profiles corresponding to the classification presented in Table 1, where the different columns are typical examples of the top 50 m of the Netherlands. The division in classes is based on the shallow subsurface composition related to the expected level of wave amplification during a seismic event.


Figure 10General outline of the vertical voxel-stack analysis used to assign the appropriate sediment class into each grid cell of the GeoTOP and NL3D geological models in the construction of the site-response zonation map.


The short lithological description in Table 1 is not sufficient to classify each site. To further aid the classification, representative sediment profiles are obtained (Fig. 9) based on the lithologic class sequences of the GeoTOP and NL3D. By grouping the main sediment profiles into the classes, we link the lithostratigraphic conditions to the expected amplification behaviour of the shallow subsurface. The classification is tested and optimised using all the sites with an estimated HVSR curve. The next step is to attribute a class to each lithostratigraphic profile per grid cell in the GeoTOP and NL3D models.

6.2 Lithology-based classification

Based on the site-response amplification estimated with the HVSR peak amplitudes at 115 sites, we have categorised each sedimentary profile (Fig. 9) into a class. The next step is to substitute GeoTOP and NL3D into these five classes. This geologically based method allows the determination of site response on regional and nationwide scales. Figure 10 gives a general outline of the procedure used to assign the appropriate sediment class to each of the voxel stacks in GeoTOP and NL3D. A voxel stack is the vertical sequence of voxels at a particular (x,y) location in GeoTOP or NL3D. Details on each of the processing steps are given in Appendix C. The next step is to attribute an amplification factor to each class.

6.3 Amplification factors for the Netherlands

For shake-map implementations or seismic hazard analysis, amplification factors (AFs) are usually derived from the Vs30 (e.g. Borcherdt1994). In this study, we estimate AFs by substituting the HVSR peak amplitudes (A0) for 115 stations throughout the Netherlands into Eq. (1). This allows the calculation of nationwide applicable AF values (AFNL) assigned to each of the classes presented in Fig. 9.

Figure 11Comparison of calculated AF distribution in terms of manual classification (a) and automatic classification by GeoTOP (b) and by NL3D (c). The locations where the empirical AF relationship is not valid are eliminated (Classes I and V). The red central mark indicates the median; the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points, and the outliers (1.5× away from the interquartile range) are plotted individually as red circles.


In order to obtain an AFNL for each class, the 115 calculated AFs are plotted against their site sediment class in Fig. 11a. For these 115 locations, the sediment classes are manually assigned based on the geological models, SCPT, or other geological data available. From the AF distributions, the mean AF values (AFNL) and corresponding standard deviation (σAF) are calculated for each class (Table 2). In Class II there are a number of sites with exactly the same AF of 1.6. These are sites with no distinguishable peak, where A0 is set to 1, which yields, after filling out in Eq. (1), AF = 1.6.

The AFNL values are valid on a national scale for a frequency range of 1–10 Hz and for reference conditions of Vs=500 m/s (Sect. 4.1). There are no AF values for sites in the farthest south and east of the Netherlands, so these areas fall into Classes I and V. There are too few data to calibrate the corresponding amplification factor.

By applying the workflow that we introduced in Sect. 6.2, automatic classification for the 115 sites is performed based on GeoTOP and NL3D and plotted against the AF (Fig. 11b and c). Due to uncertainties in the models (Appendix C), these distributions deviate from the manual classification (Fig. 11a). Note that for the manual classification SCPT information could be used at 53 sites for which local information is not included in GeoTOP and NL3D. We therefore distinguish two types of uncertainty.

  1. σAF. This is the variability that originates from the classification. Within the classification, a number of different sites are binned into the same class (Fig. 9), although in reality there is still a range of amplification behaviour. This variability is approximated with the outcome of the manual classification (Fig. 11a), which could be done in great detail.

  2. σmod. The geological models are geostatistical models where not all grid cells contain individual lithological data. Hence, there is an uncertainty of the actual lithological succession at each grid cell. The total uncertainty σtot (derived from Fig. 11b and c) can be written as σAF2+σmod2. By additionally averaging over the classes (labelled with subscript i), we find the model uncertainty σmod:

    (3) σ mod = 1 n i = 1 n σ tot , i 2 - σ AF , i 2 .

Table 2 lists the mean AF values, the uncertainty in AF (σAF), and the uncertainty (σ) for the GeoTOP and NL3D models.

6.4 Site-response zonation map

The workflow presented in Fig. 10 results in a class category assigned to each grid cell of the GeoTOP and NL3D models. As a result, we present the national site-response zonation map (Fig. 12), where each class characterises a certain level of expected site-response amplification. Additionally, each class has an AFNL assigned (Table 1). Figure 13 presents four zoom-in panels of the map, each depicting a region of particular interest.

Table 2Amplification factors and standard deviations (σ) for the NL classification. σAF is the uncertainty when a local (HVSR) recording is available. σ GeoTOP and σ NL3D represent the additional uncertainty associated with the GeoTOP and NL3D models.

Download Print Version | Download XLSX

Figure 12Seismic site-response zonation map for the Netherlands designed for low-magnitude induced earthquakes. The GeoTOP model coverage is highlighted in black in the small inset. For the remaining part of the Netherlands, the NL3D model is used as the foundation for the classification. The white spots are water bodies. The amplification factors and related uncertainties are presented in Table 2.

Figure 13Panels highlighting different regions in the site-response zonation map, including the seismometer locations. (a) North Holland: with a heterogeneous pattern between Classes III and IV. (b) The densely urbanised area of South Holland: the red line indicates the S–N cross section through the GeoTOP voxel model (Fig. C1). (c) Limburg is in the north, a quite homogeneous zone of Class II, while the south is dominated by Classes I and V due to shallow and outcropping bedrock. (d) Northeast Groningen is added as a comparison to other studies performed in that region. No seismometer locations are plotted here because of the high density covering the map.

Some areas show a large scatter in classes, which is derived from a large heterogeneity in the near surface as represented in the lithostratigraphic models. Typically, at these places there is large model uncertainty, for example in northeast North Holland (Fig. 13a). Here, the Holocene lithological successions are very heterogeneous in terms of clay, peat, and clayish sand. This region also exhibits discrepancies between the model's lithological successions and HVSR curve characteristics, for instance with seismometers J01 (Fig. 8) and J02. The geological model at these locations presents large portions of clayish sand, resulting in Class III, while the HVSR curves exhibit distinctive, high-amplitude peaks, demonstrating local conditions related to Class IV.

For larger sedimentary bodies, like the dune area, there is less model uncertainty. Dune sand is identified as Class II, and here, the HVSR of the seismometers (e.g. ALK2, Fig. 8) lacks any peak due to the absence of a velocity contrast in the near surface.

Figure 13b covers the “Randstad” region, the most densely urbanised part of the Netherlands, where the class is mainly determined as IV. Figure 13c shows the southeastern part. Most of the northern part of this region is Class II due to Pleistocene sands reaching the surface. Most of the southern part of this region falls into Class V since the bedrock occurs at a depth less than 100 m. A few places with bedrock outcrops fall into Class I.

Since Groningen has been studied in much detail, we also present the site-response zonation for this region (Fig. 13d) and discuss this in Sect. 7.

7 Discussion

The seismic site-response zonation map presented in Fig. 12 distinguishes five classes, each of which defines the potential of occurrence of the related site-response amplification. Here, the lithological conditions are collated into zonations (classes) using the classification as shown in Fig. 9. In the development of the lithostratigraphically based classification, we used (i) HVSR peak amplitudes, (ii) the presence of a velocity contrast at depth, (iii) shear-wave velocities. Amplification factors are assigned to each class. In the following paragraphs we discuss the validity and uncertainties of the classification, the AF distributions, and the usage and limitations of the presented map.

Since the ambient noise sources in the frequency band of interest (1–10 Hz) partly have an anthropogenic origin, one should be careful about contamination by local strong noise because it may seriously affect the amplitude of the HVSR as shown in Guillier et al. (2007) and Molnar et al. (2018). We resolved this problem by using large portions (30 d) of noise data to create stable HVSR curves (van Ginkel et al.2020). It is important to mention the qualitative character of the microtremor HVSR peak amplitudes, which in themselves do not directly relate to the amplification of a signal at the surface during an earthquake. However, the microtremor HVSR curve characteristics show major similarities with the measured amplification from earthquake ETFs (Fig. 5), but not in terms of absolute values. Therefore an additional fitting relationship (Eq. 1) has been defined, suitable for use of the microtremor HVSR peak amplitudes as proxies for amplification. HVSR measurements have proven to be very informative for site-response estimation and remain a valuable input for seismic site-response zonation (Molnar et al.2018; Bonnefoy-Claudet et al.2009).

Considering the difficulty in observing sufficient numbers of earthquake ground motions in areas that are not seismically active, or where no large seismic networks are available, we resorted to deriving and calibrating a lithology-based classification scheme. We took advantage of the detailed models of Cenozoic lithostratigraphy which are available in the Netherlands. As a consequence, the site-response map (Fig. 12) exhibits a regional pattern which is rather similar to the geological map (Fig. 1). We showed that the use of these models yields additional uncertainty in the determination of the AF (Table 2). This uncertainty of the actual lithostratigraphic profile at a site can be circumvented by a local recording. This may be an HVSR to obtain more certainty on the site effect (Table 1), a cone penetration test (CPT) to obtain constraints on the lithology, or, better still, a seismic cone penetration test (SCPT) to get a local shear-wave velocity profile.

Rodriguez-Marek et al. (2017) defined a site-response model including magnitude- and distance-dependent linear amplification factors (AFGr) for several period intervals (0.01–1.0 s) for the Groningen region as input for ground motion prediction equations by Bommer et al. (2017). This site-response model starts from a reference horizon at the interface between the unconsolidated sediments and the stiffer Chalk formation below at around 800–1000 m depth. However, this contrast is variable in both depth and value throughout the Netherlands and therefore not easily applicable as a reference horizon for the purpose of our study. The class-dependent AFNL presented in this paper is defined against a reference rock with a velocity at 500 m/s (which in Groningen is situated at 200 m depth). Therefore the AFGr cannot be directly be quantitatively correlated to the AFNL; this requires a correction which includes the transmission coefficient calculated at the base of the North Sea Group and a damping model. By ignoring the absolute values and comparing both AFs qualitatively, the overall spatial distribution of AFNL in the Groningen region (Fig. 13d, in a frequency band 1–10 Hz) corresponds best with AFGr at a spectral period of 0.01 s (Fig. 10; Rodriguez-Marek et al.2017). This is in line with or findings that AFs do not change much anymore when frequencies above 10 Hz are included (Fig. 4).

Usage of the site-response zonation map

The map presented in Fig. 12 enables a prediction of site response after a local earthquake as recommended in the following. It is very important to note that lithological information from geological voxel models is based on spatial interpolation and aimed at interpretations on regional scale. As a consequence, the presented site-response zonation map is also designed for regional interpretation, and not on an individual grid cell scale. Furthermore, at locations with large subsurface heterogeneity, the interpretation should be handled with care. Additional local investigations like SCPT measurements should be performed at sites of interest in order to assess the site response in detail. For the map presented, the uncertainties to keep in mind are first the AF distribution along the classes (Fig. 11a) and secondly the uncertainty of the geological model used (σ GeoTOP and σ NL3D, Table 2). The AFNL is designed to be added to an input seismic signal at a reference horizon with a shear-wave velocity of 500 m/s. This AFNL is class dependent and covers only frequencies of 1–10 Hz. Furthermore, the AFNL including the σAF does not reflect the maximum amplification that might occur within a smaller frequency band.

The frequency content of large tectonic-related earthquakes differs from induced tremors. The national AF is based on low-magnitude induced earthquakes and incorporates a frequency range of 1–10 Hz. In the case of a strong tectonic earthquake, frequencies below 1 Hz start to play a role, and resonances with deeper velocity contrasts (>100 m) which are not reflected in the current AFNL might become important. Also, for very strong ground motions, which would occur in the epicentral area of large-magnitude tectonic events, non-linearity and distance dependence could become important (Bazzurro and Cornell2004; Kwok et al.2008). Both effects have not been included in the derivation of the AFNL. Moreover, in the country's southern regions, a topographic effect may influence the site response. It is important to mention that for now these areas are aggregated in Class V and require additional detailed site investigations for site-response assessment.

8 Conclusions

In this paper we presented a workflow to create a nationwide site-response zonation, using lithological sequences as proxies for seismic site response. To that end, we first analysed the observed earthquake and ambient seismic field recorded at 69 stations of the Groningen borehole network in order to obtain empirical relationships for amplification. Based on the shallow subsurface resonance frequencies and earthquake amplitude spectra, the earthquake and ambient noise frequency band-pass filtering was applied in the range 1–10 Hz. Derived from the Groningen empirical relationships, we showed that the horizontal-to-vertical spectral ratio (HVSR) approach provides a simple means of determining the amplification potential for most subsurface conditions in the Netherlands. In a second stage, we determined the HVSR curves for additional 46 surface seismometers throughout the Netherlands and calculated the subsequent peak amplitudes. These peak amplitude distributions were related to specific lithological profiles and amplification factors. With the accrued knowledge of amplification potential of different lithological sequences, a classification scheme was designed. This turned out to be a useful tool for translation of the grid cells of the geological models into five classes and therewith establishing a national site-response zonation map. Most classes have an AFNL assigned, to which values can be added to input seismic responses adhering to the reference seismic bedrock conditions.

Class I shows sites with a hard rock setting. These sites can only be found in the very south and east of the Netherlands. An amplification factor (AF) of 1, meaning no amplification, is assigned to these locations. Class II is associated with sites with stiff sands or Pleistocene clays without strong impedance contrasts in the near surface. One may expect only small amplification at these sites. Class III shows sites with relatively soft sediments (clays, sandy clays, loess) overlying stiffer sands, resulting in impedance contrasts in the near surface. Class IV is related mostly to very soft and unconsolidated Holocene clay and peat successions overlying stiffer sands, forming a strong impedance contrast. At these sites, the largest amplification occurs. Class V shows sites at which the bedrock occurs shallower than 100 m, which is not very common in the Netherlands. For these sites there were insufficient data to assign an amplification factor.

Some limitations exist in this study. The method and map proposed are not applicable to regions with strongly deviating lithological sequences, or for earthquakes with very strong low-frequency (f<1 Hz) shaking.

Finally, it is worth noting that the proposed map could be improved by (i) adding new site geotechnical data like SCPTs, (ii) including updates and extensions of GeoTOP, (iii) including amplification factors derived from new KNMI stations, and (iv) adding new records of earthquake motions to constrain amplification factors for Class V.

Appendix A: HVSR amplification parameters

In this Appendix, HVSR peak amplitudes (A0) are fitted with the six parameters that influence ground motion site response (Fig. A1). The best fit (Rsq=0.39) is observed between A0 and Vs10. Hence, the Vs10 is used for further correlation purposes instead of the more common Vs30, supporting the findings of Gallipoli and Mucciarelli (2009) by using the Vs10 as the main amplification parameter. The depth of the first strong velocity contrast (VC, which is defined within the top 50 m) has a poor relation with A0 (Fig. A1e). The size of the velocity contrast, however, does have a strong relation with A0 (Fig. A1f).

Figure A1Each panel depicts the fitted function and the coefficient of determination (Rsq) of the HVSR peak amplitude A0 per G-network borehole location and the corresponding subsurface parameter. The (a) Vs10, (b) Vs20, (c) Vs30, (d) Vs50, (e) depth of the velocity contrast, and (f) size of the velocity contrast.


Figure A2A 3D plot of the two main parameters that define amplification: the velocity contrast and Vs10. The pale blue surface depicts the fitting function between the parameters and divides the data points where the red points are above the surface and black points are below. Blue arrows indicate the difference (error) between the surface and data point.


Compared to individual 1D correlations, a 2D correlation (Fig. A2) using both the VC and the Vs10 results in an improved correlation (Rsq=0.53) and allows us to define an empirical relationship for HVSR peak amplitudes (A0) based on these two parameters:

(A1) A 0 = - 1.29 log ( 0.01 V s 10 ) + 0.99 VC + 1.94 .

Furthermore, this equation supports the hypothesis of Joyner and Boore (1981) and Boore (2003) that A0 also depends on the VC. The motivation for Eq. (A1) is to achieve an amplification equation based on subsurface parameters only. Using Eq. (A1), an estimate of A0. is obtained Subsequently, Eq. (1) can be used to obtain an estimate of the amplification factor.

Appendix B: Geological models


GeoTOP schematises the shallow subsurface of the Netherlands in voxels measuring 100 m by 100 m by 0.5 m (x, y, z) up to a depth of 50 m below ordnance datum (Stafleu et al.2011, 2021). Each voxel contains estimates of the lithostratigraphic unit the voxel belongs to and the lithologic class (including a sand grain-size class) that is representative for the voxel. GeoTOP is publicly available from the web portal of TNO – Geological Survey of the Netherlands (GDN;, last access: 1 August 2021). GeoTOP is constructed using some 275 000 borehole descriptions from DINO, the national Dutch subsurface database operated by GDN, complemented with some 125 000 borehole logs from Utrecht University in the central Rhine–Meuse river area. The modelling procedure involves four steps: first, the borehole descriptions are interpreted into standardised lithostratigraphic units with uniform sediment characteristics. Given the large number of boreholes, automated lithostratigraphic interpretation routines (Python scripts) were developed. These routines combine digital maps, stratigraphic rules (e.g. superposition), and lithologic criteria (e.g. main lithology, admixtures, grain size, and shell content, amongst other criteria) to determine the depth of the top and base of the lithostratigraphic units in each of the borehole descriptions. Next, 2D interpolation techniques are used to construct surfaces bounding the bases of the lithostratigraphic units as observed in the boreholes. The interpolation algorithm allows for the calculation of a mean depth estimate of each surface and its standard deviation. Subsequently, all surfaces are stacked according to their stratigraphical position, resulting in a consistent layer-based model with estimates of top and base of each lithostratigraphic unit. Top surfaces are derived from the bases of the overlying units. The surfaces are then used to place each voxel in the model within the correct lithostratigraphic unit. In the third step, the borehole descriptions are revisited and classified into six different lithologic classes (“peat”, “clay”, “clayey sand & sandy clay”, “fine sand”, “medium sand”, and “coarse sand and gravel”). In the last modelling step, a 3D stochastic simulation is performed for each lithostratigraphic unit separately. The simulation results in 100 equiprobable realisations of lithologic and grain-size class for each voxel. Post-processing of the realisations results in probabilities of occurrence as well as a “most likely” estimate of lithologic and grain-size class. This most likely estimate is used in the construction of the seismic site-response zonation map (Appendix C).


To date, the GeoTOP model covers about 70 % of the country (including inland waters such as the Wadden Sea). For the missing areas we have used the lower-resolution voxel model NL3D, which is available for the entire country (van der Meulen et al.2013). NL3D models lithology and sand grain-size classes within the geological units of the layer-based subsurface model DGM (Gunnink et al.2013) in voxels measuring 250 m by 250 m by 1 m (x, y, z) up to a depth of 50 m below ordnance datum. NL3D uses a much simpler modelling procedure than GeoTOP: first, the borehole descriptions are interpreted by intersecting each borehole with the top and base raster layers from the DGM model. The resulting stratigraphical interpretations are geometrically consistent with the DGM model, but not necessarily consistent with the borehole descriptions (e.g. a borehole interval describing “sand” may erroneously fall within a unit that is characterised by clay deposits). Second, the surfaces of the DGM model are used to place each voxel in the model within the correct lithostratigraphic unit. DGM is a layer-based model using a smaller dataset of some 26 500 manually interpreted borehole descriptions from the DINO database. Consequently, it is less refined than GeoTOP. For instance, DGM combines all Holocene formations into a single unit, whereas GeoTOP features some 25 different Holocene formations, members, and beds. The third and fourth steps are identical to the ones described for GeoTOP. The resulting NL3D model has a similar most likely estimate of lithologic and grain-size class which is used in the construction of the seismic site-response zonation map (Appendix C).

B3 Model uncertainty

The current version of GeoTOP covers about 28 605 km2 using some 400 000 boreholes. This implies that only about 7 % of the voxels at land surface contain a borehole. Moreover, this number rapidly decreases with depth because many boreholes are quite shallow. Therefore, the lithostratigraphic unit and the lithologic class of almost all voxels are estimated on the basis of nearby borehole descriptions. As a rule of thumb, the limited amount of data available deeper than 30 m below land surface strongly reduces the quality of the lithologic class estimates of GeoTOP (Stafleu et al.2021). For NL3D, this number is 15 m.

B4 Applicability

GeoTOP and NL3D model the subsurface at a regional to subregional scale that is suitable for applications at the levels of province, municipality, and district. The models are not suited for applications that require a finer scale at the level of streets or individual buildings.

Appendix C: Workflow site-response map

Figure C1Cross section through the GeoTOP voxel model: (a) lithostratigraphy, (b) lithologic class, and (c) corresponding sediment site-response class. The cross section runs S–N through the city of Alphen aan den Rijn, situated on a sandy Holocene channel belt of the river Oude Rijn (“Old Rhine”). For location see Fig. 13b. Classes III and IV appear where soft, Holocene sediments (clay and peat) are overlying stiff Pleistocene deposits (sand). However, where Holocene sediments are sandy, such as in the channel belt in the centre of the cross section, Class II occurs.


The steps below describe the procedure used to assign the appropriate sediment (site-response) class to each of the voxel stacks in GeoTOP and NL3D, as exemplified in Fig. C1. A voxel stack is the vertical sequence of voxels at a particular (x,y) location in GeoTOP or NL3D. At each voxel there is an estimate of the lithostratigraphic unit and the lithologic class (Appendix B).

  • A.1. Calculate the cumulative thickness for each of the lithologic classes (peat, clay, clayey sand & sandy clay, and sand) in the models for multiple depth intervals (5, 10, 20, and 50 m). The thicknesses of the lithologic classes fine sand, medium sand, and coarse sand & gravel have been added together in the superclass sand.

  • A.2. Calculate the depth of the top of the first consecutive sequence of sand with a minimum thickness of 1.5 m (GeoTOP) or 2 m (NL3D). This depth is further referred to as “top sand”. In general, “thick” sequences of sand represent the stiffer Pleistocene sediments. In other cases, they may represent Holocene sediments of, for example, the fluvial channel belt systems of the Rhine and Meuse, or the coastal dunes. These sands form the contrast with the overlying soft sediments (peat, clay, and clayey sand and sandy clay). Voxel stacks containing a continuous Pleistocene clay sequence (Elsterian tunnel valleys) are included in the depth of the first sand (top sand), since no amplification is estimated here with the HVSR of site N02.

  • A.3. Calculate the percentage of each lithologic class above the top sand. These percentages play an important role in assigning sediment site-response classes as described in steps A.8 and A.9.

  • A.4. If anthropogenic deposits reach up to depths larger than 2 m, no sediment class is assigned. Anthropogenic activities have modified the near-surface composition at many locations in the urbanised areas of the Netherlands. The lithologic class of these sediments is unknown. Therefore, we are not able to assign a sediment site-response class to those locations.

  • A.5. If bedrock outcrops or occurs at a depth smaller than 2 m, the site is assigned to Class I. The depth criterion is set at a maximum of 2 m since a deeper top bedrock would lead to a top layer with a possible resonance in the 1–10 Hz frequency band and hence a different site-response class. The top of the bedrock is determined from the DGM model (Gunnink et al.2013) (top surfaces of the Houthem, Maastricht, and Gulpen formations).

  • A.6. If bedrock in the eastern and southern parts of the country occurs at a depth of more than 100 m, the sediment site-response class is set to V. These are sites where the layer on top of the bedrock could yield a resonance in the 1–10 Hz band, for which resonance has not sufficiently been calibrated to assign an AFNL. Class V thus corresponds to sites with a currently unknown amplification potential. The top of the bedrock is determined from the DGM-deep model (Gunnink et al.2013) (top surfaces of the Rijnland and Chalk groups).

  • A.7. If top sand is less than 2 m, the site-response class is set at II. Examples of HVSR curves with top sand less than 2 m do not exhibit any peak amplitude due to the absence of a resonating soft layer on top of a stiffer one.

  • A.8. If top sand is between 2 and 3 m, the lithologic distribution of the overlying soft sediments determine if the sediment site-response class will be II, II, or IV. Examples of HVSR curves with top sand between 2 and 3 m show peaks for certain lithological successions, forming a resonating layer. Class II is assigned if the overlying sediments contain more than 60 % sand. Class III is assigned if the overlying sediments are mainly composed of clayey sand & sandy clay and Class IV if clay and peat dominate. We do not elaborate on the exact percentages to decide between Classes III and IV. While testing the different criteria, this step appeared to be quite sensitive and needed the implementation of several exceptions to the general rule.

  • A.9. If top sand is larger than 3 m, the approach is basically the same as in A.8. Class II is assigned if the overlying sediments contain more than 60 % sand. However, the exact criteria to decide between Classes III and IV differ from those in A.8.

Code and data availability

The data for the compilation of the site-response zonation map can be downloaded from the KNMI Data Platform: (last access: 1 December 2021; van Ginkel et al.2021).

Author contributions

JvG performed the method development and geophysical and geological data analysis, developed the site-response criteria, and constructed the site-response zonation map. She wrote the manuscript with input from all co-authors. ER acted as a daily advisor and provided input on the data analysis method and results and performed text input. JS provided the GeoTOP and NL3D data, supported and quality-checked the voxel-stack analysis, and advised on the criteria for the construction of the site-response zonation map, as well as performed text input. RH acted as the PhD promotor and advisor and was the initiator of this research project. He contributed to the method development and acted as a final text editor.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Earthquake-induced hazards: ground motion amplification and ground failures”. It is not associated with a conference.


The authors thank the two reviewers for their valuable suggestions. This work is funded by EPI Kenniscentrum. Ambient noise and earthquake recordings were provided by KNMI and are publicly available through the website (last access: 1 June 2021). Figures are produced in MATLAB. We would like to thank Deltares for the use of the SCPT data and lithological interpretations and TNO for the use of the 3D geological models and maps.

Review statement

This paper was edited by Hans-Balder Havenith and reviewed by two anonymous referees.


Akkar, S., Sandıkkaya, M., and Bommer, J. J.: Empirical ground-motion models for point-and extended-source crustal earthquake scenarios in Europe and the Middle East, B. Earthq. Eng., 12, 359–387, 2014. a

Albarello, D. and Lunedei, E.: Combining horizontal ambient vibration components for H/V spectral ratio estimates, Geophys. J. Int., 194, 936–951, 2013. a

Bard, P.-Y.: Microtremor measurements: a tool for site effect estimation, in: Proceeding of the Second International Symposium on the Effects of Surface Geology on Seismic Motion, AA Balkema, Rotterdam, 3, 1251–1279, 1998. a

Bard, P.-Y.: Extracting information from ambient seismic noise: the SESAME project (Site EffectS assessment using AMbient Excitations), European Project EVG1-CT-2000-00026 SESAME, 2002. a

Bard, P.-Y., Campillo, M., Chavez-Garcia, F., and Sanchez-Sesma, F.: The Mexico earthquake of September 19, 1985 – A theoretical investigation of large-and small-scale amplification effects in the Mexico City Valley, Earthq. Spectra, 4, 609–633, 1988. a

Bazzurro, P. and Cornell, C. A.: Nonlinear soil-site effects in probabilistic seismic-hazard analysis, B. Seismol. Soc. Am., 94, 2110–2123, 2004. a

Beets, D. J. and van der Spek, A. J.: The Holocene evolution of the barrier and the back-barrier basins of Belgium and the Netherlands as a function of late Weichselian morphology, relative sea-level rise and sediment supply, Netherlands J. Geosci., 79, 3–16, 2000. a

Bindi, D., Massa, M., Luzi, L., Ameri, G., Pacor, F., Puglia, R., and Augliera, P.: Pan-European ground-motion prediction equations for the average horizontal component of PGA, PGV, and 5 %-damped PSA at spectral periods up to 3.0 s using the RESORCE dataset, B. Earthq. Eng., 12, 391–430, 2014. a

Bommer, J. J., Stafford, P. J., Edwards, B., Dost, B., van Dedem, E., Rodriguez-Marek, A., Kruiver, P., van Elk, J., Doornhof, D., and Ntinalexis, M.: Framework for a ground-motion model for induced seismic hazard and risk analysis in the Groningen gas field, the Netherlands, Earthq. Spectra, 33, 481–498, 2017. a, b, c, d, e

Bonnefoy-Claudet, S., Cornou, C., Bard, P.-Y., Cotton, F., Moczo, P., Kristek, J., and Donat, F.: H/V ratio: a tool for site effects evaluation. Results from 1-D noise simulations, Geophys. J. Int., 167, 827–837, 2006a. a

Bonnefoy-Claudet, S., Cotton, F., and Bard, P.-Y.: The nature of noise wavefield and its applications for site effects studies: A literature review, Earth-Sci. Rev., 79, 205–227, 2006b. a

Bonnefoy-Claudet, S., Köhler, A., Cornou, C., Wathelet, M., and Bard, P.-Y.: Effects of Love waves on microtremor H/V ratio, B. Seismol. Soc. Am., 98, 288–300, 2008. a

Bonnefoy-Claudet, S., Baize, S., Bonilla, L. F., Berge-Thierry, C., Pasten, C., Campos, J., Volant, P., and Verdugo, R.: Site effect evaluation in the basin of Santiago de Chile using ambient noise measurements, Geophys. J. Int., 176, 925–937, 2009. a, b, c

Boore, D. M.: Simulation of ground motion using the stochastic method, Pure Appl. Geophys., 160, 635–676, 2003. a, b

Borcherdt, R. D.: Effects of local geology on ground motion near San Francisco Bay, B. Seismol. Soc. Am., 60, 29–61, 1970. a

Borcherdt, R. D.: Estimates of site-dependent response spectra for design (methodology and justification), Earthq. Spectra, 10, 617–653, 1994. a, b

Bourne, S., Oates, S., van Elk, J., and Doornhof, D.: A seismological model for earthquakes induced by fluid extraction from a subsurface reservoir, J. Geophys. Res.-Sol. Ea., 119, 8991–9015, 2014. a

Bradley, B. A.: Strong ground motion characteristics observed in the 4 September 2010 Darfield, New Zealand earthquake, Soil Dyn. Earthq. Eng., 42, 32–46, 2012. a

Buijze, L., Van Den Bogert, P. A., Wassing, B. B., Orlic, B., and Ten Veen, J.: Fault reactivation mechanisms and dynamic rupture modelling of depletion-induced seismic events in a Rotliegend gas reservoir, Netherlands J. Geosci., 96, s131–s148, 2017. a

Busschers, F., Kasse, C., Van Balen, R., Vandenberghe, J., Cohen, K., Weerts, H., Wallinga, J., Johns, C., Cleveringa, P., and Bunnik, F.: Late Pleistocene evolution of the Rhine-Meuse system in the southern North Sea basin: imprints of climate change, sea-level oscillation and glacio-isostacy, Quaternary Sci. Rev., 26, 3216–3248, 2007. a

Camelbeeck, T. and Van Eck, T.: The 1992 Roermond earthquake, the Netherlands, and its aftershocks, edited by: Camelbeeck, T., van Eck, T., Pelzing, R., Ahorner, L., Loohuis, J., Haak, H. W., Hoang-Trong, P., and Hollnack, D., Geologie en Mijnbouw, 73, 181–197, 1994. a

Castellaro, S., Mulargia, F., and Rossi, P. L.: VS30: Proxy for seismic amplification?, Seismol. Res. Lett., 79, 540–543, 2008. a

Comité Européen de Normalisation (CEN): Eurocode 8: design of structures for earthquake resistance – Part 1: general rules, seismic actions and rules for buildings, Comité Européen de Normalisation, Brussels, 2004. a, b

Convertito, V., De Matteis, R., Cantore, L., Zollo, A., Iannaccone, G., and Caccavale, M.: Rapid estimation of ground-shaking maps for seismic emergency management in the Campania Region of southern Italy, Nat. Hazards, 52, 97–115, 2010. a

De Gans, W.: Quaternary Geology of the Netherlands, in: Geology of the Netherlands, edited by: Wong, T. E., Batjes, D. A. J., and De Jager, J., Royal Netherlands Academy of Arts and Sciences, Amsterdam, 173–195, 2007. a

Dost, B., Ruigrok, E., and Spetzler, J.: Development of seismicity and probabilistic hazard assessment for the Groningen gas field, Netherlands J. Geosci., 96, s235–s245, 2017. a

Fäh, D., Kind, F., and Giardini, D.: A theoretical investigation of average H/V ratios, Geophys. J. Int., 145, 535–549, 2001. a

Falcone, G., Acunzo, G., Mendicelli, A., Mori, F., Naso, G., Peronace, E., Porchia, A., Romagnoli, G., Tarquini, E., and Moscatelli, M.: Seismic amplification maps of Italy based on site-specific microzonation dataset and one-dimensional numerical approach, Eng. Geol., 289, 106170,, 2021. a

Field, E. H. and Jacob, K. H.: A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent, B. Seismol. Soc. Am., 85, 1127–1143, 1995. a

Gallipoli, M. R., Calamita, G., Tragni, N., Pisapia, D., Lupo, M., Mucciarelli, M., Stabile, T. A., Perrone, A., Amato, L., Izzi, F., La Scaleia, G., Maio, D., Salvia, V.: Evaluation of soil-building resonance effect in the urban area of the city of Matera (Italy), Eng. Geol., 272, 105645,, 2020. a

Gallipoli, M. R. and Mucciarelli, M.: Comparison of site classification from VS 30, VS 10, and HVSR in Italy, B. Seismol. Soc. Am., 99, 340–351, 2009. a, b

Gariel, J., Horrent, C., Jongmans, D., and Camelbeeck, T.: Strong ground motion computation of the 1992 Roermond earthquake, the Netherlands, from linear methods using locally recorded aftershocks, Geologie en Mijnbouw, 73, 315–315, 1995. a

Geluk, M., Duin, E. T., Dusar, M., Rijkers, R., Van den Berg, M., and Van Rooijen, P.: Stratigraphy and tectonics of the Roer Valley Graben, Geologie en Mijnbouw, 73, 129–129, 1995. a

Gouw, M. and Erkens, G.: Architecture of the Holocene Rhine-Meuse delta (the Netherlands) – a result of changing external controls, Netherlands J. Geosci., 86, 23–54, 2007. a

Guillier, B., Chatelain, J.-L., Bonnefoy-Claudet, S., and Haghshenas, E.: Use of ambient noise: From spectral amplitude variability to H/V stability, J. Earthq. Eng., 11, 925–942, 2007. a

Gunnink, J., Maljers, D., Van Gessel, S., Menkovic, A., and Hummelman, H.: Digital Geological Model (DGM): a 3D raster model of the subsurface of the Netherlands, Netherlands J. Geosci., 92, 33–46, 2013. a, b, c, d

Hijma, M. P., Cohen, K. M., Hoffmann, G., Van der Spek, A. J., and Stouthamer, E.: From river valley to estuary: the evolution of the Rhine mouth in the early to middle Holocene (western Netherlands, Rhine-Meuse delta), Netherlands J. Geosci., 88, 13–53, 2009. a

Hinzen, K.-G., Reamer, S. K., and Fleischer, C.: Seismicity in the Northern Rhine Area (1995–2018), J. Seismol., 25.2, 1–17, 2020. a

Hofman, L., Ruigrok, E., Dost, B., and Paulssen, H.: A shallow seismic velocity model for the Groningen area in the Netherlands, J. Geophys. Res.-Sol. Ea., 122, 8035–8050, 2017. a, b, c

Houtgast, R. and Van Balen, R.: Neotectonics of the Roer Valley rift system, the Netherlands, Global Planet. Change, 27, 131–146, 2000. a, b

Joyner, W. B. and Boore, D. M.: Peak horizontal acceleration and velocity from strong-motion records including records from the 1979 Imperial Valley, California, earthquake, B. Seismol. Soc. Am., 71, 2011–2038, 1981. a, b

KNMI: Netherlands Seismic and Acoustic Network, Royal Netherlands Meteorological Institute (KNMI), Other/Seismic Network,, 1993. a

Kokusho, T. and Sato, K.: Surface-to-base amplification evaluated from KiK-net vertical array strong motion records, Soil Dyn. Earthq. Eng., 28, 707–716, 2008. a

Konno, K. and Ohmachi, T.: Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor, B. Seismol. Soc. Am., 88, 228–241, 1998. a

Kruiver, P. P., van Dedem, E., Romijn, R., de Lange, G., Korff, M., Stafleu, J., Gunnink, J. L., Rodriguez-Marek, A., Bommer, J. J., van Elk, J., and Doornhof, D.: An integrated shear-wave velocity model for the Groningen gas field, the Netherlands, B. Earthq. Eng., 15, 1–26,, 2017a. a, b, c, d, e

Kruiver, P. P., Wiersma, A., Kloosterman, F. H., de Lange, G., Korff, M., Stafleu, J., Busschers, F. S., Harting, R., Gunnink, J. L., Green, R. A., van Elk, J., and Doornhof, D.: Characterisation of the Groningen subsurface for seismic hazard and risk modelling, Netherlands J. Geosci., 96, s215–s233, 2017b. a, b

Kwok, A. O., Stewart, J. P., and Hashash, Y. M.: Nonlinear ground-response analysis of Turkey Flat shallow stiff-soil site to strong ground motion, B. Seismol. Soc. Am., 98, 331–343, 2008. a

Lachetl, C. and Bard, P.-Y.: Numerical and theoretical investigations on the possibilities and limitations of Nakamura's technique, J. Phys. Earth, 42, 377–397, 1994. a, b

Lee, V. W. and Trifunac, M. D.: Should average shear-wave velocity in the top 30 m of soil be used to describe seismic amplification?, Soil Dyn. Earthq. Eng., 30, 1250–1258, 2010. a

Lermo, J. and Chavez-Garcia, F. J.: Site effect evaluation using spectral ratios with only one station, B. Seismol. Soc. Am., 83, 1574–1594, 1993. a, b

Lunedei, E. and Malischewsky, P.: A review and some new issues on the theory of the H/V technique for ambient vibrations, in: Perspectives on European Earthquake Engineering and Seismology, Springer International Publishing AG, Switzerland, 371–394,, 2015. a

Majer, E. L., Baria, R., Stark, M., Oates, S., Bommer, J., Smith, B., and Asanuma, H.: Induced seismicity associated with enhanced geothermal systems, Geothermics, 36, 185–222, 2007. a

Mena, B., Wiemer, S., and Bachmann, C.: Building robust models to forecast the induced seismicity related to geothermal reservoir enhancement, B. Seismol. Soc. Am., 103, 383–393, 2013. a

Mignan, A., Landtwing, D., Kästli, P., Mena, B., and Wiemer, S.: Induced seismicity risk analysis of the 2006 Basel, Switzerland, Enhanced Geothermal System project: Influence of uncertainties on risk mitigation, Geothermics, 53, 133–146, 2015. a

Molnar, S., Cassidy, J., Castellaro, S., Cornou, C., Crow, H., Hunter, J., Matsushima, S., Sánchez-Sesma, F., and Yong, A.: Application of microtremor horizontal-to-vertical spectral ratio (MHVSR) analysis for site characterization: State of the art, Surv. Geophys., 39, 613–631, 2018. a, b, c

Nakata, N., Snieder, R., Kuroda, S., Ito, S., Aizawa, T., and Kunimi, T.: Monitoring a building using deconvolution interferometry. I: Earthquake-data analysis, B. Seismol. Soc. Am., 103, 1662–1678, 2013. a

Noorlandt, R. P., Kruiver, P. P., de Kleine, M. P. E., Karaoulis, M., de Lange, G., Di Matteo, A., von Ketelhodt, J., Ruigrok, E., Edwards, B., Rodriguez­Marek, A., Bommer, J. J., van Elk, J., and Doornhof, D.: Characterisation of ground motion recording stations in the Groningen gas field, J. Seismol., 22, 605–623, 2018. a, b

Panzera, F., Bergamo, P., and Fäh, D.: Canonical Correlation Analysis Based on Site-Response Proxies to Predict Site-Specific Amplification 695 Functions in Switzerland, B. Seismol. Soc. Am., 111, 1905–1920,, 2021. a

Paulssen, H., Dost, B., and Van Eck, T.: The April 13, 1992 earthquake of Roermond (The Netherlands); first interpretation of the NARS seismograms, Geologie en Mijnbouw, 71, 91–91, 1992. a, b

Peeters, J., Busschers, F., Stouthamer, E., Bosch, J., Van den Berg, M., Wallinga, J., Versendaal, A., Bunnik, F., and Middelkoop, H.: Sedimentary architecture and chronostratigraphy of a late Quaternary incised-valley fill: a case study of the late Middle and Late Pleistocene Rhine system in the Netherlands, Quaternary Sci. Rev., 131, 211–236, 2016. a

Perron, V., Gélis, C., Froment, B., Hollender, F., Bard, P.-Y., Cultrera, G., and Cushing, E. M.: Can broad-band earthquake site responses be predicted by the ambient noise spectral ratio? Insight from observations at two sedimentary basins, Geophys. J. Int., 215, 1442–1454, 2018. a

Pilz, M., Parolai, S., Leyton, F., Campos, J., and Zschau, J.: A comparison of site response techniques using earthquake data and ambient seismic noise analysis in the large urban areas of Santiago de Chile, Geophys. J. Int., 178, 713–728, 2009. a

Poggi, V., Edwards, B., and Fäh, D.: Derivation of a reference shear-wave velocity model from empirical site amplification, B. Seismol. Soc. Am., 101, 258–274, 2011. a

Rodriguez-Marek, A., Kruiver, P. P., Meijers, P., Bommer, J. J., Dost, B., van Elk, J., and Doornhof, D.: A Regional Site-Response Model for the Groningen Gas Field, B. Seismol. Soc. Am., 107, 2067–2077, 2017. a, b, c, d

Rondeel, H., Batjes, D., and Nieuwenhuijs, W.: Synopsis: Petroleum geology of the Netherlands – 1993, Geology of gas and oil under the Netherlands: The Royal Geological and Mining Society of The Netherlands, Utrecht, S3–S20, 1996. a, b, c

Schokker, J.: Geologische overzichtskaart van Nederland, TNO Bouw en Ondergrond, Utrecht, 2010. a

Stafleu, J., Maljers, D., Gunnink, J., Menkovic, A., and Busschers, F.: 3D modelling of the shallow subsurface of Zeeland, the Netherlands, Netherlands J. Geosci., 90, 293–310, 2011. a, b, c

Stafleu, J., Maljers, D., Busschers, F. S., Schokker, J., Gunnink, J. L., and Dambrink, R. M.: Models Created as 3-D Cellular Voxel Arrays, Applied Multidimensional Geological Modeling: Informing sustainable human interactions with the shallow subsurface, John Wiley & Sons, Ltd, 247–271,, 2021. a, b, c, d

Van den Berg, M. and Beets, D.: Saalian glacial deposits and morphology in The Netherlands, Tills and Glaciotectonics, Balkema, Rotterdam, 235–251, 1987. a, b

van der Meulen, M. J., Doornenbal, J. C., Gunnink, J. L., Stafleu, J., Schokker, J., Vernes, R. W., van Geer, F. C., van Gessel, S. F., van Heteren, S., van Leeuwen, R. J. W., Bakker, M. A. J., Bogaard, P. J. F., Busschers, F. S., Griffioen, J., Gruijters, S. H. L. L., Kiden, P., Schroot, B. M., Simmelink, H. J., van Berkel, W. O., van der Krogt, R. A. A., Westerhoff, W. E., and van Daalen, T. M.: 3D geology in a 2D country: perspectives for geological surveying in the Netherlands, Netherlands J. Geosci., 92, 217–241, 2013. a, b, c, d

van Ginkel, J., Ruigrok, E., and Herber, R.: Assessing soil amplifications in Groningen, the Netherlands, First Break, 37, 33–38, 2019. a, b, c, d, e, f, g, h, i, j, k, l

van Ginkel, J., Ruigrok, E., and Herber, R.: Using horizontal-to-vertical spectral ratios to construct shear-wave velocity profiles, Solid Earth, 11, 2015–2030,, 2020. a, b, c

van Ginkel, J., Ruigrok, E., Stafleu, J., and Herber, R.: seismic-site-response-zonation-map, KNMI Data Platform [data set], available at:, last access: 1 December 2021. a

van Thienen-Visser, K. and Breunese, J.: Induced seismicity of the Groningen gas field: history and recent developments, The Leading Edge, 34, 664–671, 2015. a

Wapenaar, K., Slob, E., Snieder, R., and Curtis, A.: Tutorial on seismic interferometry: Part 2 – Underlying theory and new advances, Geophysics, 75, 75A211–75A227, 2010.  a

Wassing, B., Dost, B., and Vermillion, N. P.: Seismisch hazard van geinduceerde aardbevingen: integratie van deelstudies, TNO/KNMI, Utrecht/de Bilt, 2012. a

Weatherill, G., Kotha, S. R., and Cotton, F.: Re-thinking site amplification in regional seismic risk assessment, Earthq. Spectra, 36, 274–297, 2020. a

Wills, C., Petersen, M., Bryant, W., Reichle, M., Saucedo, G., Tan, S., Taylor, G., and Treiman, J.: A site-conditions map for California based on geology and shear-wave velocity, B. Seismol. Soc. Am., 90, S187–S208, 2000. a

Zagwijn, W.: The Netherlands during the Tertiary and the Quaternary: a case history of coastal lowland evolution, in: Coastal Lowlands, Springer, Dordrecht, 107–120, 1989. a

Short summary
A soft, shallow subsurface composition has the tendency to amplify earthquake waves, resulting in increased ground shaking. Therefore, this paper presents a workflow in order to obtain a map classifying the response of the subsurface based on local geology, earthquake signals, and background noise recordings for the Netherlands. The resulting map can be used as a first assessment in regions with earthquake hazard potential by mining or geothermal energy activities, for example.
Final-revised paper