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

. 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 5 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 estab-lish a nationwide site-response zonation by combining 3D lithostratigraphic models and earthquake and ambient vibra- 10 tion recordings. As a ﬁrst step, we constrain the parameters (velocity contrast and shear-wave velocity) that are indicative of ground motion ampliﬁcation in the Groningen area. For this, we compare ambient vibration and earthquake recordings us- 15 ing the horizontal-to-vertical spectral ratio (HVSR) method, borehole empirical transfer functions (ETFs), and ampliﬁcation factors (AFs). This enables us to deﬁne an empirical relationship between the ampliﬁcation measured from earthquakes by using the ETF and AF and the ampliﬁcation es- timated from ambient vibrations by using the HVSR. With this, we show that the HVSR can be used as a ﬁrst proxy for site response. Subsequently, HVSR curves throughout the Netherlands are estimated. The HVSR amplitude characteristics largely coincide with the in situ lithostratigraphic se- quences and the presence of a strong velocity contrast in the near surface. Next, sediment shallow subsurface classes,

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 vibra-10 tion 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 us-15 ing 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 es-20 timated 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 se- 25 quences 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 quan-30 tified 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 35 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. ence 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 50 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 55 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 60 extensive. Rather than using ground motion prediction equations with generic site amplification factors conditioned on V s30 , 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, 65 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 70 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 75 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(Stafleu et al., , 2021 and NL3D (Van der Meulen et al., 2013), with a derived classi-80 fication 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 85 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 90 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 95 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). 100 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 de-105 positional environments in the early Quaternary to a complex alternation of shallow marine, estuarine, and fluvial sediments in the younger periods (Zagwijn, 1989;Rondeel et al., 1996;De Gans, 2007;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 5 (Van den Berg and Beets, 1987).

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 10 a Holocene fluvial system (Rondeel et al., 1996;Beets and van der Spek, 2000). 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.

15
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 20 systems embedded in flood basin clays (Gouw and Erkens, 2007). 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;25 Van den Berg and Beets, 1987). 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 30 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 forma-35 tions 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 2. Map of the Netherlands depicting epicentres of all induced (M w 0.5-3.6, orange) and tectonic (M w 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 (M w ≥ 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.

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 5 (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 Balen, 2000), and the rift system still shows active extension (Hinzen et al., 2020). The largest earthquake recorded (M w = 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 15 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 cir-20 cles, 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 25 the magnitudes are relatively low (van Thienen-Visser and Breunese, 2015), the damage to buildings in the area is substantial due to shallow hypocentres and amplification on the soft near-surface soils Kruiver et al., 2017a). 30

Dataset
For this study, we use the seismic network of the Royal Netherlands Meteorological Institute (KNMI, 1993), consisting of borehole and surface seismometers distributed over the Netherlands (Fig. 2). The network includes 88 locations with 35 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 sur-face. 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 5 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 Gronin-10 gen borehole network between May 2015 and May 2019. All data are available through the KNMI data portal (http: //rdsa.knmi.nl/network/NL/, last access: TS1 ).
In the construction of the site-response map, we have made extensive use of the detailed 3D geological subsurface mod-15 els 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 mea-20 suring 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(Stafleu et al., , 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 25 within these units as well as measures of model uncertainty. GeoTOP is publicly available from the TNO's web portal: https://www.dinoloket.nl/en/subsurface-models (last access: TS2 ). 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) 35 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 layerbased subsurface models DGM and DGM-deep. These models are also available from the web portal mentioned above. 40 More details on the models GeoTOP and NL3D are given in Appendix B.

Empirical relationships from the Groningen borehole network
The dense Groningen borehole network (G-network) pro-45 vides 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,

50
we compare subsurface parameters with the measured am-plification in order to evaluate which parameters are most critical.

Definition of amplification
The majority of site-response studies define the level of soft 55 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 60 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 65 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 70 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 75 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. 80 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 aver-85 age 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/m 3 . At 200 m depth, 95 % of the Dutch subsurface is composed of laterally extensive Pleistocene and Pliocene sediments; hence the estimated site-response and 90 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 seismome-100 ters, accelerometers, and geophones) has high sensitivity for ground motion.
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 attenu-5 ation 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 earth-10 quake recorded in borehole G24. The location of this borehole is presented as a red triangle in Fig. 2.

Amplification factors
We compute an overall amplification factor from the Gnetwork earthquake recordings by taking the ratio of the 15 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 20 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 In this study we aim to provide an average amplification level in a broad frequency band. The 35 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 40 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 45 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).

Empirical transfer functions
Empirical transfer functions (ETFs) represent shear-wave 50 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 55 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 60 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. 65 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.
From the ETF between the 200 m depth and surface seis-70 mometer (ETF 200 ), peak amplitudes are identified. Some example ETF curves are plotted in Fig. 5. Additionally, ETF curves for the top 50 m are calculated (ETF 50 ). The ETF curves for the different intervals show very similar peak characteristics and peak amplitudes, demonstrating that most am-75 plifications 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 80 Holocene alternations of clay and peat, overlying consolidated Pleistocene sediments.

H /V spectral ratios from the ambient seismic field
The Groningen surface seismometers are continuously recording the ambient seismic field (ASF), and these data 85 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 sedi-90 ment 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 95 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 100 nearby footsteps or traffic.
Based on the HVSR curve and peak characteristics, different criteria are defined conformable to the SESAME consensus (Bard, 2002): (1) clear peaked curves (HVSR ampli-  tude > 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 5 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 10 et al. (2008), Konno and Ohmachi (1998) and Lermo and Chavez-Garcia (1993) and this contrast is very likely to amplify the ground motion.  . Relation between the HVSR peak amplitudes (A 0 ) 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. R sq (R-squared) represents the coefficient of determination of the fitting. Note the log-scale of the y axis.

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 5 by the earthquake-derived ETF and AF. Figure 6 displays the correlation between the peak amplitudes of the HVSR and ETF 200 as well as HVSR and AF. Secondly, we evaluate the subsurface seismic parameters enhancing amplification (Fig. 7).

10
Based on these data points, relationships are defined in order to be able to estimate amplification factors using HVSR peak amplitudes (A 0 ), amplification factors: and maximum amplification as measured by the ETF: 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 20 site response consider the shear-wave velocity for the top 30 m (V s30 ) 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 (CEN et al., 2004). However, studies (Castellaro et al., 2008;Kokusho and Sato, 25 2008; Lee and Trifunac, 2010) have drawn attention to the fact that using only V s30 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 30 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, 35 firstly, the AF at each borehole location is fitted against several subsurface parameters. AF = a + b · e c·V s 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. V s10 , V s20 , and V s30 -values are 40 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 45 V s50 . 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 5 AF = a + b · log(c · VC) as a functional form. This procedure is also performed for the relation between the subsurface parameters (V s , VC) and the HVSR. The results are given in Appendix A. Figure 7 presents the fit between the AF and the six rele-10 vant subsurface parameters. Best fit (R sq = 0.47, Fig. 7a and b) is observed between the AF and V s10 and V s20 . This is in line with findings of Gallipoli and Mucciarelli (2009), who use the V s10 as the main amplification parameter instead of the more common V s30 . In Groningen, the low-velocity and 15 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 V s30 and V s50 arises since the amplification develops mainly in the Holocene sediments (van Ginkel et al.,20 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 25 HVSR.

HVSR estimations throughout the Netherlands
Based on the good relationship between Groningen HVSR peak amplitudes (A 0 ), ETF peak amplitudes, and AF ( Fig. 6), we conclude that the HVSR can be further used as proxy 30 for amplification. Therefore, for all surface seismometers in the Netherlands seismic network, the HVSR curves are estimated following the method described in Sect.  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 ex-50 hibiting clear peak amplitudes (NH01,J01, ZH030, FR010, EETW) are located at sites with a distinct velocity contrast Figure 8. Each 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.
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 5 sediments overlying hard rock, and the HVSR curves exhibit a clear peak amplitude.

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 10 (CEN et al., 2004, 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 15 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 20 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 corre-25 sponding AF.

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). Be-30 yond 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 35 was shown to have a (albeit weaker) link with amplification (Figs. 7 and A1).
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 40 Figure 9. CE1 Sediment 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.
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 ve-5 locity contrast (VC), and average shear-wave velocity for the top 10 m.
For setting up the classification, we initially use A 0 , the peak amplitude from HVSR. We only have measured AF in Groningen, whereas we have measured A 0 for many sites 10 throughout the Netherlands (all stations depicted in Fig. 2). Moreover, we found a clear relationship between A 0 and AF (Eq. 1). The relationships between V s10 , VC, and A 0 are estimated from lithological conditions in Groningen, where the sites are assigned to Classes II, III, and IV. For Classes I and 15 V we have insufficient empirical constraint on A 0 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, 20 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.
The short lithological description in Table 1 is not sufficient to classify each site. To further aid the classification, 25 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 classifica-30 tion 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 35 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 40 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 45 on each of the processing steps are given in Appendix C. The next step is to attribute an amplification factor to each class.

Amplification factors for the Netherlands
For shake-map implementations or seismic hazard analysis, amplification factors (AFs) are usually derived from the V s30 50 (e.g. Borcherdt, 1994). In this study, we estimate AFs by substituting the HVSR peak amplitudes (A 0 ) for 115 stations throughout the Netherlands into Eq. (1). This allows the calculation of nationwide applicable AF values (AF NL ) assigned to each of the classes presented in Fig. 9. 55 In order to obtain an AF NL 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, 60 the mean AF values (AF NL ) 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 A 0 is set to 1, which yields, after filling out in Eq. (1), AF = 1.6. 65 The AF NL values are valid on a national scale for a frequency range of 1-10 Hz and for reference conditions of V s = 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 70 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), 75 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 infor- mation 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 5 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 10 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 σ 2 AF + σ 2 mod . By additionally averaging over the 15 classes (labelled with subscript i ), we find the model uncertainty σ mod : Table 2 lists the mean AF values, the uncertainty in AF (σ AF ), and the uncertainty (σ ) for the GeoTOP and NL3D models. 20

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 25 of expected site-response amplification. Additionally, each class has an AF NL assigned (Table 1). Figure 13 presents four zoom-in panels of the map, each depicting a region of particular interest. Some areas show a large scatter in classes, which is de-30 rived 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, 35 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 40 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) 45 deficit CE2 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. 50 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. 55 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.

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 5 (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. Ampli-fication factors are assigned to each class. In the following 10 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-1 TS4 Hz) partly have an anthropogenic origin, one 15 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   Table 2. 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 5 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 amplifi-10 cation. 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 15 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 lithologybased classification scheme. We took advantage of the detailed models of Cenozoic lithostratigraphy which are avail-20 able 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 lithos-25 tratigraphic 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 veloc-30 ity profile. Rodriguez-Marek et al. (2017) defined a site-response model including magnitude-and distance-dependent linear amplification factors (AF Gr ) for several period intervals (0.01-1.0 s) for the Groningen region as input for ground 35 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 through-40 out the Netherlands and therefore not easily applicable as a reference horizon for the purpose of our study. The classdependent AF NL 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 AF Gr cannot 45 Figure 13. Panels 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. be directly be quantitatively correlated to the AF NL ; 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 5 of AF NL in the Groningen region (Fig. 13d, in a frequency band 1-10 Hz) corresponds best with AF Gr 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 inter-15 polation 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 AF NL is designed to be added to an input seismic signal at a reference horizon with a shear-wave velocity of 500 m/s. This AF NL is class dependent and covers 10 only frequencies of 1-10 Hz. Furthermore, the AF NL 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 15 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 AF NL might 20 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 Cornell, 2004;Kwok et al., 2008). Both effects have not been included in the derivation 25 of the AF NL . 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. 30

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 35 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 40 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 addi-45 tional 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 se-50 quences, 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 AF NL assigned, to which values can be added to input 55 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 am-60 plification, 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) 65 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 70 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 devi-75 ating 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, 80 (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 (A 0 ) are fitted with 85 the six parameters that influence ground motion site response (Fig. A1). The best fit (R sq = 0.39) is observed between A 0 and V s10 . Hence, the V s10 is used for further correlation purposes instead of the more common V s30 , supporting the findings of Gallipoli and Mucciarelli (2009) by using the V s10 90 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 A 0 (Fig. A1e). The size of the velocity contrast, however, does have a strong relation with A 0 (Fig. A1f).

95
Compared to individual 1D correlations, a 2D correlation (Fig. A2) using both the VC and the V s10 results in an improved correlation (R sq = 0.53) and allows us to define an empirical relationship for HVSR peak amplitudes (A 0 ) based on these two parameters: 100 A 0 = −1.29 log(0.01V s10 ) + 0.99VC + 1.94.
Furthermore, this equation supports the hypothesis of Joyner and Boore (1981) and Boore (2003) that A 0 also depends on

B1 GeoTOP
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(Stafleu et al., , 2021. Each voxel contains estimates of the 10 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; https://www.dinoloket.nl/en/subsurface-models, last 15 access: TS6 ). GeoTOP is constructed using some 275 000 borehole descriptions from DINO, the national Dutch subsurface database operated by GDN (https://www.dinoloket. nl/en/subsurface-data, last access: 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 lithos-15 tratigraphic 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 esti-20 mate is used in the construction of the seismic site-response zonation map (Appendix C).

B2 NL3D
To date, the GeoTOP model covers about 70 % of the country (including inland waters such as the Wadden Sea). For the 25 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 grainsize classes within the geological units of the layer-based subsurface model DGM  in voxels 30 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 re- 35 sulting 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 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 km 2 using some 400 000 boreholes. This implies that only about 7 % 55 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, 60 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.

65
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.

70
Appendix C: Workflow site-response map 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 75 (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, 80 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 con-85 secutive 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 90 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 95 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. tant 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 5 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 de-15 termined from the DGM model  (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 sed-20 iment 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 AF NL . Class V thus corresponds to sites with a currently unknown 25 amplification potential. The top of the bedrock is determined from the DGM-deep model  (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 30 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. Exam-35 ples 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 40 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 sev-45 eral 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 50 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: https://dataplatform.knmi.nl/dataset/ seismic-site-response-zonation-map-1-0 (last access: TS8 ; van Ginkel et al., 2021).
Author contributions. JvG: method development, data analysis, construction of the map and underlying criteria, writing of manuscript with input from all co-authors. ER: daily advisor, input on data analysis method and results, performed text input. JS: provided the GeoTOP and NL3D data, supported and quality checked 10 the voxel-stack analysis and advised on the criteria for the construction of the site-response zonation map and performed text input. RH: promotor, initiator of this research project, advisor, method development, final text editor. TS9 Competing interests. The contact author has declared that neither 15 they nor their co-authors have any competing interests.
Disclaimer. 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 20 "Earthquake-induced hazards: ground motion amplification and ground failures". It is not associated with a conference. TS10 Acknowledgements. 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 25 and are publicly available through the website http://rdsa.knmi.nl/ dataportal (last access: TS11 ). 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.

CE2
Do you mean "decrease"? It seems you've used "deficit" as a verb here, which is incorrect. Please check and suggest a suitable alternative.
Remarks from the typesetter TS1 Please provide date of last access.

TS2
Please provide date of last access.

TS3
Please check URL.

TS4
Please check numbers.

TS5
Should this be (f)?

TS6
Please provide date of last access.

TS7
Please provide date of last access.

TS8
Please provide date of last access.

TS9
Please provide this section in complete sentences.

TS11
Please provide date of last access.

TS12
Please ensure that any data sets and software codes used in this work are properly cited in the text and included in this reference list. Thereby, please keep our reference style in mind, including creators, titles, publisher/repository, persistent identifier, and publication year. Regarding the publisher/repository, please add "[data set]" or "[code]" to the entry (e.g. Zenodo [code]).

TS13
Please check name.

TS14
Please provide the complete author list.

TS15
Please provide the page range.

TS16
Please provide publisher and place of publication.

TS17
Please provide the volume and page range/DOI.

TS18
Please provide the complete author list.

TS19
Please provide DOI.

TS20
Please provide the volume.

TS21
Please provide the complete author list.

TS22
Please provide the volume.

TS23
Please provide the complete author list.

TS24
Please provide publisher and place of publication.

TS25
Please provide the complete author list.

TS26
Please provide the volume and page range.

TS27
Please provide place of publication.

TS28
Please provide publisher and place of publication.

TS29
Please provide the complete author list.

TS30
Please provide date of last access.

TS32
Please provide publisher and place of publication.

TS33
Please provide place of publication.