1 Introduction

Gas production in the Groningen field in the northern Netherlands began in 1963 and since 1991 the consequent compaction of the sandstone formation where the reservoir resides has been generating induced seismicity (van Thienen-Visser et al. 2015; Bourne et al. 2018).The largest earthquake to date was the ML 3.6 (M 3.5, Dost et al. 2019) Huizinge earthquake of August 2012, following which the field operator initiated an extensive data acquisition and modelling programme to quantify the resulting seismic hazard and risk (van Elk et al. 2017). The risk estimation model includes four key components: a seismicity model that defines earthquake magnitudes and recurrence rates as a function of gas production scenarios (Bourne et al. 2017a, b); a ground motion model that estimates response spectral accelerations at the ground surface as a result of each induced earthquake scenario (Bommer et al. 2017b); an exposure model specifying the location and construction characteristics of all of the ~ 160,000 buildings within the gas field and a surrounding 5 km buffer (Arup et al. 2020); and fragility and consequence models that transform the predicted ground-motion amplitudes into damage estimates for each building typology and then estimate the potential for injury to building occupants as a function of this damage (Crowley et al. 2017a,b, 2019a). The risk model is calibrated to estimate the impact of both gas production changes and building strengthening interventions on the risk to inhabitants of the region, providing a rational framework for decision-making with regards to mitigation options (van Elk et al. 2019).

The ground motion model, or GMM, has been developed to reflect the specific conditions of the Groningen gas field, including the geological profile above the reservoir, which includes the Zechstein salt formation (Bommer et al. 2016, 2017a). The GMM predicts ground-motion amplitudes at a rock horizon at a depth of about 800 m and these motions are then transformed to estimates of shaking at the ground surface by non-linear frequency-dependent amplification factors (Rodriguez-Marek et al. 2017). The amplification factors are defined for ~ 160 zones within which the dynamic response characteristics are similar. The basis for the site response model is a field-wide model for the shear-wave velocity, \(V_{S}\), profiles from the ground surface down to the reference rock horizon at the base of the North Sea formation (Kruiver et al. 2017a, b). The near-surface \(V_{S}\) profiles are based on the GeoTOP model that defines the lithological profiles across the Netherlands (Stafleu and Dubelaar 2016); the resulting \(V_{S}\) profiles have been validated by site-specific measurements using both borehole and non-invasive techniques (Noorlandt et al. 2018).

The Groningen GMM has evolved over several years through the incorporation of additional data and the refinement of the model. Enhancements of the ground-motion model have included the incorporation of finite fault rupture simulations (Edwards et al. 2019), Groningen-specific models for component-to-component variability and for spatial correlation (Stafford et al. 2019), scenario-dependence of the linear amplification factors for the soft soils encountered in the field (Stafford et al. 2017), and incorporation of uncertainty in the non-linear soil properties (Bahrampouri et al. 2019). In parallel efforts, the Groningen induced seismic risk model has also been extended to consider the potential impact of geotechnical hazards as well as ground shaking. Simplified approaches to the assessment of liquefaction triggering potential were adapted to the specific conditions of the Groningen field (Green et al. 2019). The Groningen-specific models were then embedded in a fully probabilistic assessment of liquefaction triggering hazard (Green et al. 2020), which demonstrated that this is a very minor hazard and effectively a negligible contribution to the induced seismic risk.

This paper addresses the most recent refinement of the ground-motion model for Groningen, which arises as from rather unique geotechnical structures that exist in the region. Throughout the northern coastal region of the Netherlands, a number of buildings are located on dwelling mounds (Fig. 1) that were constructed several centuries ago. The mounds, called terps in Dutch and wierden in the Groningen dialect, were constructed to raise building foundations as a defence against flooding. Although the GeoTOP model does include some anthropogenic soils, the wierden are not specifically represented (as well as natural soils, the mounds often include a large proportion of animal dung) and therefore their presence has not yet been accounted for in the seismic risk modelling until now. While only a very small proportion of the exposed buildings in the Groningen region are founded on these dwelling mounds, it has been a long-term objective to include the wierden into the risk modelling for completeness. Although the total number of buildings located on wierden is small, these include many structures considered to be part of the cultural heritage of the region. In addition, several old and densely populated village centres are located on wierden.

Fig. 1
figure 1

Example of a wierde showing the circular layout and the elevated position in the flat landscape. The wierde in the image is Niehove, the Netherlands. Source https://www.kwaliteitsgidsgroningen.nl/wierdenland-waddenkust/landschap

In Sect. 2, the history, spatial distribution, dimensions and geotechnical characteristics of the dwelling mounds are described. In order to be able to model the dynamic response of the wierden, \(V_{S}\) p measurements were made at a number of selected mounds across the region (Sect. 3). Site response analyses (Sect. 4) were then performed to estimate the dynamic amplification of ground motions by the dwelling mounds. Section 5 discusses the incorporation of the wierden into the risk calculations, which first required identifying the buildings that are founded on dwelling mounds. For the building types located on these geotechnical structures, the risk model needs to consider both the modification of the ground shaking due to the presence of the mound and the possible modification of the fragility functions by the soft anthropogenic soils from which they are constructed. The paper closes with brief discussions and conclusions in Sect. 6.

2 Location and characteristics of the Wierden

Wierden are man-made dwelling mounds and were created from the early pre-Roman Iron Age until the early Middle Ages in a yet undiked coastal salt marsh environment subject to periodic flooding (Bazelmans et al. 2012; Nieuwhof and Schepers, 2016). They are typical landscape elements for the entire Wadden Sea region of The Netherlands, Germany and Denmark. More than 900 wierden are identified in the Province of Groningen in the northern part of the Netherlands (Fig. 2; Miedema, 1983; Kuijer, 1987; Miedema, 1990; Koomen and Maas, 2004; Province Groningen, 2016). They are dome-shaped, relatively small in size, and have an altitude of several metres above their direct surroundings. They were built of almost any material available in the direct surroundings. For example, sods from the surrounding salt marsh deposits, soil from dug out pits and ditches, and animal dung, mainly from cattle (Miedema, 1983; Knol, 1993; Bazelmans et al. 2012; Nieuwhof 2019). The composition of wierden is therefore different from the surrounding land, especially since vertical accretion continued during and after initial construction of the wierden. Given their relatively small size, varying from individual farms to villages, and the regional scale of the geological model used in the site response calculations, the wierden are not well represented in the GMM. Moreover, general geomechanical properties for anthropogenic material were attributed to wierden material, while these properties were primarily based on modern (sandy) foundation layers of modern housing areas.

Fig. 2
figure 2

Wierden (dwelling mounds) in the northern part of the Netherlands and the location of the eight wierden investigated in the study. The background indicates the major physical geographical regions (adapted after Koomen and Maas, 2004)

Their significance in terms of seismic risk lies in the fact that several populated village centres are situated on wierden. In addition, a substantial part of the (built) cultural heritage of the province is located on or associated with wierden, including twelfth century churches, typical regional housing, and village layouts with high cultural heritage values (Bazelmans et al. 2012). Archaeological and lithological variations between wierden can be high, even between locations in the same region and of the same age (Nieuwhof and Schepers, 2016; Meijles et al. 2016; Nieuwhof et al. 2019). This is due to the spatial variation in salt marsh conditions and historical use of the dwelling mounds. Eight wierden were selected as a representative sample for the investigation (Fig. 2, Table 1). The representative wierden were selected based on the major physical geographical regions within the coastal zone, size, age and expected composition. The physical geographical region defines the characteristics of both the natural sequences and the sediments available for building. Other considerations include the number of buildings on the wierde and the social willingness of the inhabitants to cooperate with the research. The selected wierden form a well-balanced mix.

Table 1 Summary of characteristics of selected wierden

A linear transect of hand auger drillings with a spacing of 4 m and a length of ~ 110 m was drilled across each wierde. Cores were described according to NEN5104 (NEN 1989) and relevant archaeological standards. Lithology was converted to lithoclasses as defined in the geological model GeoTOP of TNO-Geological Survey of the Netherlands (Van der Meulen et al. 2013; Stafleu and Dubelaar, 2016) for further processing in the site response analyses. The target depth of the hand auger drillings was 4 to 5 m below the surface or deeper to cover the transition between the anthropogenic material and the natural soil beneath the wierde. The drillings show that there is a large degree of heterogeneity: successions of thin layers of clay, peat, sandy clay and some fine sand. Animal dung or manure is a specific soil type that is only present in wierden. The layers of manure vary in thickness between a few centimetres to more than one metre. Pie charts of typical wierde composition are shown in Fig. 3. Most of the wierden have comparable proportions of soft material, predominantly organic and clayey material. Exceptions are Groot Maarslag, with a relatively large proportion of fine sand, and Grote Houw, which consists of 70% fine sand.

Fig. 3
figure 3

Schematic representation of wierde composition. Lithoclasses conform GeoTOP. Animal dung is included in the “organic” lithoclass

3 Shear-wave velocity measurements on the wierden

The shear-wave velocity (\(V_{S}\) p) structure of the eight wierden was determined using Multi-channel Analysis of Surface Waves (MASW, Park et al. 1999) along the same transect as the hand auger drillings. The aim was to characterize the shallow 2D \(V_{S}\) p structure (upper 10 to 15 m) and develop 1D \(V_{S}\) p models to depths of 150 to 300 m. Survey design was developed by NAM and data acquisition was performed by Rossingh Geophysics. The first survey on the Groot Maarslag wierde was used to test several sources and sampling intervals for sources and receivers. From the results, a generic design was extracted with a strongly reduced field effort and reduced impact on the communities. This optimized design was then applied on all other wierden. Rayleigh and Love wave MASW acquisition and microtremor arrays were performed. One densely sampled line was chosen to characterize lateral heterogeneity. In addition, a circular array was deployed for the derivation of a 1D velocity profile at lower frequencies, complementing the 2D profile at extended depths. The setup of the fieldwork campaign is summarized in Table 2.

Table 2 Summary of survey design

The 2D MASW data were reduced using the common midpoint cross correlation (CMP-cc) approach (Hayashi and Suzuki 2004) to maximize spatial resolution, followed by a wavefield transform to map data from time versus offset to frequency versus phase velocity domain. This approach has been adapted to generate wider bandwidth dispersion curves, especially with multi-source data, by generating CMP-cc gathers over multiple receiver offset ranges and combining the resulting dispersion curves. Near-field effects were reduced by limiting maximum wavelength as a function of the offset range used for analysis based on Yoon and Rix (2009). Shear-wave velocity models were developed at 4 m intervals along each profile using fundamental or effective mode, local search inversion routines and combined to form a 2D image of VS structure along the profile. The 2D VS models developed from the Rayleigh wave data are used in the current analysis. The vertical resolution decreased with increasing depth. The top model layer averages a thickness of 0.7 m and model layer thicknesses increase to ~ 3.5 at ~ 12 m depth.

The 1D VS models were derived from the combined multiple dispersion data sets from microtremor data with different analysis techniques and MASW data with different sources using the extended spatial autocorrelation (ESAC) technique (Okada, 2003), the high-resolution frequency wavenumber transform (Capon, 1969) and recently developed Rayleigh three component beamformer (Wathelet et al. 2018) as implemented in the Geopsy software package (Wathelet et al. 2020). The 1D VS model information was used to inspire a representative half-space VS at the base of the 2D VS models.

In general, the 2D VS information shows an increase in VS values with depth. The transition from anthropogenic material to natural soil roughly corresponds to the transition from lower to higher VS values. This transition is somewhat smeared out because of the vertical resolution of 0.7 m (in the top) to 3.5 m (near the bottom). Within the general increase in VS with depth, there are several local features identified. Two examples of vertical cross-sections of 2D VS data are shown in Figs. 4 and 5, together with the lithoclasses from the drillings. The section of Amsweer (Fig. 4) shows velocities varying between 50 and 230 m/s, with a general increase with depth. At the right side of the figure, the drillings indicate the presence of a zone of peat in the natural soil. This corresponds well with the very low VS zone, which is typical of peat. The profile shows a stepwise increase in VS with the transition at a depth of 2 to 3 m. The section of Groot Maarslag (Fig. 5) shows a sharper transition from anthropogenic material to natural soil compared to Amsweer, especially at the eastern side of the road. There is also relatively more animal dung (grey bars) present in this wierde. The presence of animal dung corresponds to a zone of relatively low VS.

Fig. 4
figure 4

Cross-section showing the VS across Amsweer from 2D MASW data of Rayleigh waves. The grey dotted line denotes the maximum reliable depth of the result. GeoTOP lithoclasses of the drillings are shown in the profile

Fig. 5
figure 5

Cross-section showing the VS across Groot Maarslag from 2D MASW data of Rayleigh waves. The grey dotted line denotes the maximum reliable depth of the result. GeoTOP lithoclasses of the drillings are shown in the profile

The statistical analysis of the 2D VS information of the eight wierden is summarized in Table 3. The full statistical analysis results are included in the Supplementary Material (Table A and B). All investigated wierden show significantly lower VS values in the anthropogenic wierde body than the natural deposits below the wierde, based on analysis of variance (ANOVA) testing. The average VS value of the anthropogenic deposits is 85 m/s, ranging from 68 m/s in Amsweer to 101 m/s in Grote Houw. The average VS of the natural subsurface is 102 m/s, which is significantly higher. The absolute differences in VS between natural subsurface and wierde deposits ranges from 5 m/s in Biessum to 38 m/s in Beswerd.

Table 3 Summary statistics of VS values of anthropogenic and natural lithoclasses of the selected wierden

The transition from anthropogenic material to natural soil roughly corresponds to the transition from lower to higher VS values. In all wierden and their natural subsurface, the Pearson correlation coefficient r between VS and depth is significant, albeit weak in the case of Biessum. In most wierden, VS gradually increases with depth, whereas in three cases (Amsweer, Groot Maarslag and Fransum) the increase is stepwise, with the increment in Vs values usually around the wierde base level.

In the ANOVA analysis (Table 4), lithoclasses show significantly different shear wave velocities (p = 0.000), with peat (69 m/s) being the lowest for all samples within the wierde and below the wierde. Clayey sand, clay and animal dung have significantly higher VS (80–87 m/s). Fine sand shows the highest values (102 m/s). Within the wierde body, however, this seems to be slightly different, with sandy clay providing the lowest values, increasing via animal dung and clay to fine sand being the highest. For natural deposits in the immediate area surrounding the wierde, only measured at Groot Maarslag, there is no significant difference in VS between lithoclasses.

Table 4 ANOVA analysis of difference in VS values between lithoclasses

4 Dynamic response analyses for the Wierden

The first objective of the site response analyses is to quantify the difference in site response in terms of spectral amplification at the wierden with the site response predicted by the GMM. The second objective is to compare the response of the different wierden in order to evaluate if these differences are sufficient to require a wierde-specific model. For the comparison between local wierde information and the GMM, site response analyses were conducted for the three wierden within the GMM area for the full soil column of ~ 800 m, i.e. Amsweer, Biessum and Helwerd (Table 1 and Fig. 2). For the comparison between the eight wierden, site response analyses were conducted on much shorter soil columns of 16 m length for which the local information was available for all wierden.

Site response analyses were conducted using the same approach as during the seismic hazard analysis, which is described in Rodriguez-Marek et al. (2017) and summarized here. Calculations were performed using STRATA (Rathje and Ozbey 2006; Kottke and Rathje, 2008) for one-dimensional vertically propagating shear-waves. The magnitudes covered in the hazard analysis range from 2.0 to 7.25. For earthquakes larger than observed (ML = 3.6), non-linear soil behaviour is to be expected as a result of the soft material and low VS. Therefore, the equivalent linear approach as implemented in STRATA was chosen, combined with Random Vibration Theory (RVT) in frequency domain. The input for the STRATA calculations consists of soil profiles, input motions and soil properties dictating the shear modulus reduction and damping as a function of strain. The three different types of input are described below.

The soil columns were constructed using the local information at the wierden and the GMM. The local information consists of drillings to a depth of 2.3 to 5.8 m and 2D VS profiles to a depth of 16 to 19 m (Fig. 6). The drillings are shorter than the VS profiles. The stratigraphy and lithoclass information for depths larger than the maximum depth of the drillings is derived from the GeoTOP model v1.4 (https://www.dinoloket.nl/en/subsurface-models).Three of the eight wierden fall within the area for which the Groningen GMM was derived, i.e. the outline of the gas field including a 5 km buffer (Fig. 2). The other five wierden fall outside the GMM area, but are also considered representative of wierden in the region. For the wierden within the GMM area, full soil profiles from surface to the GMM reference baserock depth of ~ 800 m are available. These soil profiles consist of stratigraphy, lithoclass descriptions and VS profiles (Kruiver et al. 2017a, b). This enables one-to-one comparison between the local data and the GMM model. In order to derive “wierde response”, the top section of the GMM soil column is replaced by the local VS information from the 2D Rayleigh VS results and the drilling stratigraphy and lithoclasses (Fig. 7). The GMM soil columns are defined on a 100 m × 100 m grid. The local VS data covers a line of 110 m in length, which typically corresponds to two grid cells of the GMM. For the five wierden outside the GMM area, the available information is limited to the depth of the VS data (~ 16–19 m), the drillings (2.3 to 5.8 m) and GeoTOP for stratigraphy and lithoclasses. For the comparison of site response between the wierden, a common length of soil profiles is required. Therefore, the site response analysis calculations for the comparison of all wierden start at the half-space at a depth of 16 m.

Fig. 6
figure 6

Construction of the local soil column from the drilling information, GeoTOP and local VS information

Fig. 7
figure 7

Schematic representation (not to scale) of the soil columns and input motions for site response calculations. The local part is visualized in Fig. 6. Left: used for comparing local data to GMM model for the three GMM wierden. Middle: used for the 16 m columns of all eight wierden. Right: derivation of FAS motions at 16 m depth to be used as input motions for the soil columns of the middle panel. All motions are defined as outcrop motions

The next input for the site response calculations is the base motions. The GMM calculations are based on Frequency Amplitude Spectrum (FAS) motions for the version 6 (V6) version of the GMM (Bommer et al. 2017a, 2019). The FAS motions at the reference baserock horizon (base of the North Sea Supergroup) were obtained from finite rupture EXSIM stochastic simulations (Motazedian & Atkinson, 2005; Boore, 2009) over a range of magnitudes (1.5 to 7.5) and rupture distances (3 to 60 km) and four branches for stress drop. In the GMM, the motions are randomly chosen from the set of 3600 motions. For the wierden site response analysis, a fixed set of 10 motions was selected, enabling direct comparisons of results from different wierden and from the GMM. The selection of motions was such that it covers the full range of peak ground acceleration (PGA) values and representative motions dominating the seismic hazard (epicentral distance of 4–8 km and M 4.2–4.5).

For the wierden within the GMM development area, these motions were applied at the base of the soil columns, which corresponds to the reference baserock horizon (Fig. 7, left). For the wierden outside the GMM area, a representative set of input motions needed to be defined at the alternative reference horizon of 16 m (Fig. 7, middle). This is an uncommonly short soil column in site response analysis, because part of the amplification is expected to have a deeper origin. Indeed, the PGA profiles of the GMM wierden show an increase in PGA generally starting at ~ 50 m depth for Amsweer and Biessum and at ~ 20 m depth for Helwerd. The effect of including only a short column in the site response calculations is that only the shorter periods are included. The effects from periods longer than those corresponding to the quarter wavelength of the input motion will not be captured in the analysis. The purpose of this exercise, however, is not to derive the amplification factors for each wierde, but comparing the response between the wierden. Using the same assumptions for each wierde (uniform reference horizon, same input motions) enables a clean comparison. The approach for deriving the alternative set of input motions is explained in Fig. 7 (right). The selected motions were applied to one full-length standard soil profile from the Helwerd wierde. The FAS motions were then extracted at 16 m and regarded as input motion for the 16 m soil columns. The validity of this approach was checked by using this 16 m FAS motion as input to the top 16 m of the original full column. The FAS motion at the surface obtained in this way was compared to the surface FAS motion resulting from site response calculation of the full column. The FAS motions at the surface from both approaches corresponded well.

The last type of input is the description of the non-linear soil behaviour. This behaviour is defined by the Modulus Reduction and Damping (MRD) curves. In the GMM, the Darendeli (2001) curves are used for the lithoclasses “clay” and “sandy clay and clayey sand”. The Menq (2003) curves are used for sands of various grain sizes. Because of the abundance of peat in the Groningen region, site specific curves were derived for peat (Zwanenburg et al. 2020). For each soil type, described by the combination of stratigraphic unit and lithoclass, appropriate parameters were estimated (Kruiver et al. 2018). For animal dung, there are no parameters available. As a result of loading by the material on top of the dung layers, these layers have compacted. Based on the composition and appearance in the cored material, the properties of Basal Peat are assumed to be representative for dung.

The spectral amplification factors (AF) are defined as the ratio of the spectral acceleration at the surface over the spectral acceleration at the reference baserock horizon (i.e. either 800 m for full columns or 16 m for short columns). The period range which is relevant for risk ranges from 0.01 to 1.0 s. The AF for the three GMM wierden for a representative motion is shown in Fig. 8. The general appearance of the response between the wierden is similar, but the absolute values are different. Amsweer and Biessum have the highest AF, whereas Helwerd shows generally lower AF. The variation in AF between the three wierden is larger than within each wierde.

Fig. 8
figure 8

Amplification Factors (AF) for all GMM wierde profiles for one of the selected motions. Each wierde contains ~ 24 soil profiles

For each of the local data coordinates, the AF from the GMM model VS and soil column (Kruiver et al. 2017a, b) was calculated as well. The local 2D VS data with a horizontal length of 110 m fall within two GeoTOP grid cells. This results in two different GMM soil profiles which were appended below the local data (Fig. 7, left), corresponding with the GMM grid cell of the local VS profile. The AF is calculated for those two soil columns at each coordinate for all 10 motions. For each input motion, spectral AF results from the model profile can be compared to the spectral AF from the soil profile with the top replaced by local data. The relative difference in AF is shown in Fig. 9 for Amsweer as an example. Each dot represents one coordinate on the 2D line and one motion. The average and the standard deviation are represented by the red lines. The average difference over all periods is 8% for Amsweer and 18% for both Biessum and Helwerd. The periods relevant for risk assessment are T = 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.85 and 1.0 s. The average difference in AF for the risk-relevant periods is 7%, 17% and 28% for Amsweer, Biessum and Helwerd respectively. These numbers generally surpass the commonly acceptable uncertainties in geo-engineering. The effect of the wierden on site response therefore is significant. Moreover, the AF from the GMM model is lower on average than the AF from local data. This means that the model underestimates the AF for all three GMM wierden.

Fig. 9
figure 9

Comparison of AF results per soil column and motion for GMM wierde Amsweer (dots). The solid red line represents the average AF for each period, the dashed red lines indicate plus or minus one standard deviation

The average spectra AF behaviour of the three GMM wierden is shown in Fig. 10. The dotted line is based on the GMM soil columns on the location of the wierden and the solid line is based on the soil columns with the top replaced by local data. There is a difference in AF between the local data and the GMM model, not only in maximum value, but also in a shift in the peak period. This means that in the risk analyses, if a building is situated on a wierde, the AF for the specific site needs to be corrected using a penalty function. The difference suggests that a penalty function for buildings situated on wierden should be period- or frequency-dependent. The penalty function is derived in Sect. 5.

Fig. 10
figure 10

Mean Amplification Factors for the local data (solid line) compared to the model data for the grid cells corresponding to the wierde locations (dotted line), calculated for the four motions dominating the seismic hazard

The only wierde where local data are available both on and off the wierde was Groot Maarslag. The 16 m soil columns were used to compare the site response for these two situations. Figure 11 shows that the AF on the wierde (lines 1000 and 3000) are significantly different from the natural ground some distance away from the wierde (line 6000). There is a shift in dominant period and the wierde AF show a significant peak around a period of 0.2 s.

Fig. 11
figure 11

AF results for Groot Maarslag for one representative motion showing the difference in behaviour on wierde (line 1000 and 3000 in grey) and next to wierde (line 6000 in purple). Each line represents the AF curve for one soil column of local VS data

The second objective of the site response analysis is to assess the difference in AF among the eight wierden. This analysis serves to show if an average wierde can be defined and if so, whether the GMM wierden are representative of all investigated wierden. A penalty function can only be derived when considering the full soil column as is done for the GMM. For the comparison between the individual wierden, all columns have been truncated at a common reference depth of 16 m, which is the minimum depth range of the 2D VS profiles (Figs. 6 and 7). The average AF curves for the 16 m columns are shown in Fig. 12. All wierden, except for Grote Houw, show similar AF behaviour among the periods relevant for risk (0.1–1.0 s). Grote Houw shows substantially lower AF compared to the other wierden. This is probably related to the fact that this wierde is much sandier than the other wierden (Fig. 3) and therefore less prone to amplification. The AFs are calculated over relatively short soil columns. The time-averaged VS over the top 16 m on the wierden varies between 114 and 154 m/s. Following the quarter wavelength formula, this would result in a theoretical maximum period of 0.4–0.5 s. Information from deeper layers is included in the input motions at 16 m, but would be identical for all 16 m columns of the different wierden columns by design (Fig. 7). This suggests that the AFs shown for periods longer than T = 0.5 s in Fig. 12 likely do not capture the differences among the wierden. However, it is unlikely that differences in the near-surface deposits affect these low frequencies. Moreover, the overall response of the wierden for shorter periods (T < 0.5 s) do not vary much from wierde to wierde, apart from Grote Houw. Based on the fact that 7 out of 8 wierden show similar AF behaviour we conclude that one average wierde response is sufficient. In addition, the GMM wierden seem to fall on the high side of the average curve. This means that using the average or median AF for these three wierden for the full soil column (Fig. 10) is a conservative choice.

Fig. 12
figure 12

Mean Amplification Factors result for the 16 m soil columns of all wierden

5 Incorporating Wierden into seismic risk calculations

As noted in the Introduction, the influence of the wierden has so far not been accounted for in the estimation of seismic risk due to the induced earthquakes in the Groningen field. The work that has now been carried out, as presented in the preceding sections, makes it possible to explicitly include the dwelling mounds into the risk calculations. To understand the specific steps taken to incorporate the wierden into the risk calculations, it is useful to first summarize how the risk engine works. Figure 13 presents a very simplified schematic for the steps involved in generating the earthquakes and resulting ground motion field in the Monte Carlo simulations. The scheme illustrates the calculation for a single building typology within one grid under one realisation of the ground-motion field for a single earthquake in one realisation of the earthquake catalogue. In the Monte Carlo simulations for a given future gas production scenario, multiple earthquake catalogues are generated and for each earthquake, ground-motion fields are generated by randomly sampling from the between-event and within-event components of variability while imposing an appropriate spatial correlation (Stafford et al. 2019). In order to capture epistemic uncertainty, a logic-tree is defined with alternative models and parameter values for each component of the model, and these are also sampled in proportion to the branch weights in the Monte Carlo simulations (van Elk et al. 2019). Sufficient numbers of simulations are generated to obtain stable estimates of the resulting risk at the annual exceedance frequencies specified in Dutch safety regulations. The red rectangles in Fig. 13 represent the simple scheme that was devised to incorporate the dynamic effects of the wierden into the risk calculations.

Fig. 13
figure 13

Schematic illustration of the steps in the seismic risk calculations; in practice, all buildings across the field are considered and the calculations repeated for every event in every realization of the earthquake catalogue. The red boxes indicate the modified or additional steps for those buildings identified as being located on a dwelling mound

In the following sections, the modifications to this risk calculation scheme made to account for the presence of wierden are described individually, starting with the penalty function for the local site amplification. We then discuss the number and types of buildings encountered on these dwelling mounds in the Groningen field. Finally, the effect of the soft wierden soils on the structural response and the way this could modify the building fragility functions is explored.

5.1 Ground motion model: penalty function for wierden sites

As indicated in Fig. 13, the GMM combines estimates of the response spectral accelerations at the baserock rock horizon with nonlinear, frequency-dependent AFs assigned to each of the 162 site response zones defined for the field to obtain estimates of the ground motions at the ground surface. This is the motion that is then used as input to the risk calculations, which is appropriate since most structures in the region are on rather shallow foundations. For those buildings located on wierden, the input motions need to be defined at the top of dwelling mound, so an additional adjustment—which we choose to call a Penalty Factor—is required to transform the estimated motions at the ground surface into motions at the top of the wierde. To this effect, we use the results of the dynamic response analyses presented in Sect. 4. In order to define a penalty function, AF was calculated for all 75 soil columns of the three GMM wierden for randomly-sampled input motions (from V7 of the GMM) at the reference baserock horizon of ~ 800 m, covering the full range from weak to strong motions (PGAbaserock range 10–8- 1 g). For each column, two calculations were made, each differing on whether the top 18 to 19 m consisted of either the local soil data or the generic soil layers used to develop the GMM model (Figs. 6 and 7). Using the resulting AF values, the residuals, defined as the difference between the AF with local data and the AF for the generic layers were computed in natural log space (i.e., \({\Delta } = \ln \left( {AF_{local} } \right) - \ln \left( {AF_{GMM} } \right)\), where \({\Delta }\) are the residuals). The computed residuals are shown in Fig. 14 (top).

Fig. 14
figure 14

Top: Residuals [ \(\ln \left( {AF_{local} } \right) - \ln \left( {AF_{GMM} } \right)\)] for the three wierden analysed. The solid line is the mean residual for the three wierden and the dotted line is the 75th percentile considering the uncertainty in the mean. Bottom: proposed penalty function for buildings on wierden

To derive the wierden Penalty Factors, several considerations were taken into account. First, the GMM already incorporates a strong degree of epistemic uncertainty. The 95% confidence interval (CI) of the epistemic uncertainty (shown in Fig. 14 for the zones where the three wierden under analysis) generally envelopes the residuals, as shown in Fig. 14 (top). For this reason, it was deemed appropriate to apply a simple modification, consisting of a unique Penalty Factor for all wierden. The Penalty Factor is based on the average of the mean residuals of the three wierden. To account for the fact that only three wierden are used, the standard error of the mean is used to compute the 75th percentile of the mean residual. Finally, the Penalty Factor is defined via a piecewise linear function that envelopes the 75th percentile (Fig. 14, bottom). The corrected amplification factor AFbuilding on wierden is given by:

$$ln\left( {AF_{building \, on \, wierde} } \right) = ln\left( {AF_{GMM zone} } \right) + Penalty Factor$$
(1)

with the Penalty Factor in ln units.

5.2 Exposure model: buildings founded on wierden

As described in van Elk et al. (2019), a probabilistic model of the buildings found within a 5-km buffer around the Groningen gas field outline has been produced, starting from public and private data sets containing information such as e.g. geographical coordinates, year of construction, footprint layout, building height, usage, etc.. Next, using inference rules and engineering-judgement, the structural configuration that is most likely to constitute the load-resisting system for each unique building in the database is defined. Then for each possible classification, a seismic vulnerability class is assigned. Different buildings with different structural systems may be assigned to the same vulnerability class, if their seismic response and performance is similar (further details in Crowley et al. 2017b; Crowley and Pinho, 2020). Finally, it is noted that, out of the more than 260,000 buildings featured in the V7 Groningen exposure model (Arup et al. 2020), approximately 39% refer to structures such as sheds and garages that are not regularly occupied, and which are thus not considered in the risk analyses.

Once the locations of the wierden were identified (Sect. 2), it was possible to spatially link these data with the geographical coordinates of each one of the 157,956 occupied buildings defined in the exposure model, and thus arrive at an estimate of 2862 structures constructed on terrain featuring the presence of wierden soil layers. As summarized in Table 5, nine vulnerability classes characterize the seismic response of the majority (> 88%) of the buildings founded on wierden, which consist of unreinforced masonry (URM) structures founded on shallow foundations, mainly of a detached/terraced houses type (see Fig. 15), with the exception of vulnerability class URM1F_B, which is associated to farmhouse barns (i.e. the barn structure shown together with the vulnerability class building URM1F_HC in Fig. 15). The latter, however, owing to their structural and foundation characteristics (absence of a base slab, light roof loads distributed along lengthy flexible walls), are subjected to no (or negligible) soil-structure-interaction (SSI) effects, for which reason they will not be considered in the subsequent Sect. 5.3, given that their fragility functions are derived without consideration of SSI.

Table 5 Vulnerability classes for which there is a minimum quantity of 1% of the buildings founded on wierden
Fig. 15
figure 15

Screenshots of LS-Dyna numerical structural models (Arup 2017, 2019) used to characterize the seismic response of buildings representative of the vulnerability classes considered in this work

5.3 Fragility model: influence of wierden on building response

The collapse fragility functions (which describe the probability of reaching the collapse state, conditional on the intensity of input ground motion) that have been previously developed for the Groningen region (Crowley et al. 2019b; Crowley and Pinho 2020) have duly considered the effects of SSI through the employment of the methodology proposed and described in Cavalieri et al. (2020a, b). The latter considered geomechanical parameters and soil stratigraphy from Kruiver et al. (2017a, b) and Rodriguez-Marek et al. (2017) that did not take into account the presence of wierden soil layers. In the current endeavour, therefore, collapse fragility functions for the vulnerability classes indicated in Sect. 5.2 above were re-derived considering a representative soil profile (66_101) that features two wierden layers with a total thickness of 1.4 m (Table 6).

Table 6 Geomechanical properties and stratigraphy of the selected representative soil profile (66_101) considered in SSI analyses with wierden. γ is unit weight, PI is plasticity index, OCR is overconsolidation ratio, Cu is coefficient of uniformity for grain sizes, D50 is median grain size, k0 is coefficient of earth pressure and Su is undrained shear strength

The development of the collapse fragility functions followed the procedure described in Crowley et al. (2017a, b) and Cavalieri et al. (2020a, b). Index buildings for each of the vulnerability classes in Table 5, which are typically all constructed with shallow strip foundations, have been used to develop the fragility functions. Fixed-base MDOF models for each index building were produced in LS-DYNA (LSTC, 2013) (Fig. 15) and were subjected to nonlinear dynamic analyses using 11 training records (see Arup, 2017, 2019 for further details). The maximum attic displacement of a given MDOF model under each training record was converted to the equivalent SDOF displacement (see Crowley et al. 2017a, b) and then compared with the displacement obtained under the same records for a fixed-base SDOF model in SeismoStruct. This SDOF system has been modelled in SeismoStruct (Seismosoft, 2020) with the multi_lin model (Sivaselvan and Reinhorn 1999), which is characterised by a polygonal hysteresis loop and can simulate the deteriorating behaviour of strength and stiffness. The soil-structure interactions analyses were carried out using the SeismoStruct software (Seismosoft, 2020), where the non-linear SSI macro-element model by Correia and Paolucci (2021), capable of capturing the nonlinear soil behaviour at the near-field and the ground substratum dynamic characteristics at the far-field, as well as the interaction with the seismic response of the structure, is implemented. The input parameters of such SSI macro-element were defined using the model calibration procedure described in Cavalieri et al. (2020a,b), which involves the analysis of 3D finite element models of the buildings and their predominant foundation system, as well as considering the geomechanical properties indicated in Table 6. The complete set of input parameters defined for each of the structural models considered in this work can be found in the report by Mosayk (2020).

Hazard-compatible records were selected from a large database including European (Akkar et al. 2014) and NGA-West records (Chiou et al. 2008), through disaggregation of Groningen seismic hazard at four return periods (Tr = 500, 2500, 10 k and 100 k years), using the mean magnitude and distance from the disaggregation together with the 2017 ground motion and 5–75% significant duration prediction equations for the Groningen field (Bommer et al. 2018). As illustrated in Fig. 16, the records were selected to match both response spectra and 5–75% significant durations conditioned on four different levels of average spectra acceleration (AvgSa—Baker and Cornell, 2006), corresponding to the four return periods, using the ground motion selection procedure proposed by Baker and Lee (2018), namely the Conditional Spectrum. It is noted that AvgSa was adopted as the intensity measure (IM) in this study not only because it has been shown to be sufficient (Bianchini et al. 2009; Eads et al. 2015; Kohrangi et al. 2017), but also because, unlike spectral acceleration at the period of vibration of the structure (which can also constitute a sufficient intensity measure), it allows a comparison between the fragility functions obtained for the different structural systems considered (each of which has a different period of vibration).

Fig. 16
figure 16

Response spectra of selected records and the conditional spectra (herein represented with the mean and ± 2σ) to which they have been matched (Cavalieri et al. 2020b)

The fragility functions are herein represented by lognormal cumulative distribution functions (CDF) of a scalar IM (i.e. AvgSa). As such, they can be described in terms of two parameters, namely the median of AvgSa (in units of g) corresponding to the collapse limit state and its dispersion (β). Table 7 reports these lognormal distribution parameters for the eight vulnerability classes considered in this study.

Table 7 Lognormal collapse fragility function parameters for the eight vulnerability classes

Figure 17 shows a comparison between the collapse fragility functions obtained considering SSI with and without wierden soil layers. The latter were derived in Mosayk (2019) considering field-wide average soil properties. It can be readily observed that the consideration of wierden leads to little or no impact in the fragility, especially for those vulnerability classes with a higher number/percentage of buildings on wierden soil (i.e. detached/farm houses—URM6L, URM8L, URM7L, URM1F_HA, URM1F_HC). A slightly more noticeable, even if always very modest, change (decrease) in fragility is observed for terraced house vulnerability classes (URM3L, URM4L), as these are weaker structures, which thus benefit more from the energy dissipation that takes place in the soil when it responds in the non-linear range. However, as shown in Table 5, only 0.35% of such buildings are founded on wierden soil and consequently the impact of this reduction in the fragility for this particular building class when located on dwelling mounds on the computed risk estimates over the whole field would be very minor.

Fig. 17
figure 17

Fragility functions developed considering SSI with and without wierden layers. It is noted that whilst these fragility models have been plotted and compared up to a value of AvgSa = 1.5 g, this corresponds to a return period greater than 100,000 years, and the risk is mainly driven by levels of AvgSa less than 0.5 g where the impact of the wierden is even more negligible

Since it was concluded that exactly the same fragility functions can be used for all of these building classes whether founded on wierden or on natural ground, either because the soft soils of the wierden only modify the SSI effects to a very small degree or because the modification is slightly larger for very few cases, the numbers of buildings to which this would apply is too small to warrant the modification. The classification of the building types as being located on wierden then becomes simply a flag in the exposure model to indicate that the penalty functions defined in Sect. 5.1 need to be applied to the surface motions before estimating the building response in terms of damage.

6 Discussion and conclusions

A bespoke seismic risk model has been developed as a tool for the management of induced earthquakes in the Groningen gas field in the Netherlands. The model has been calibrated to local conditions through extensive data gathering efforts, ranging from large networks of accelerographs and shear-wave velocity measurements to full-scale shake table testing of representative buildings constructed from local materials (e.g. Brunesi et al. 2019; Graziotti et al. 2019; Tomassetti et al. 2019). In parallel with the development of this risk model, which can be used to estimate the impact of both changes in gas production and building strengthening interventions on the risk due to ground shaking, studies have also quantified the risk due to liquefaction triggering.

One element that had not been incorporated into the risk model previously was the influence of dwelling mounds, of which some ~ 900 are encountered in the northern Netherlands. These ancient mounds, constructed for flood defence, represent another unique feature of the Groningen field and their explicit inclusion is considered important in order to complete the risk model. Of the 158,000 occupied buildings in the Groningen exposure database, only about 2,900 are situated on wierden (as the dwelling mounds are called locally) but these include some dense village centres and some important cultural heritage as well. From the exposure model, 8 building classes were identified among the houses founded on wierden, which account for almost 90% of the buildings on these mounds. These buildings are now flagged in the exposure database as being located on dwelling mounds and the work described in this paper has been undertaken to enable the risk calculations to explicitly account for the influence of these unusual foundation conditions.

The first step was to characterize the dynamic response of the wierden in order to include in the risk calculations the additional amplification of the shaking due to the presence of the additional layers of soft material above the natural ground surface. For this purpose, 8 wierden were selected as being representative; six of these were predominantly comprised of clayey materials and dung, another was found to be more sandy in composition, although all of them were found to be internally very heterogeneous. Shear-wave velocity measurements were made across these selected wierden using MASW and in general the VS values of the materials comprising the mounds were found to be appreciably lower than those of the natural materials underlying the wierden. Site response analyses were performed to determine the relative increase in ground motion amplification at the top of the mounds relative to the ground surface. The additional amplification of shaking by the sandy dwelling mound was found to be lower than that at the more clayey structures; there were also appreciable differences in the degree of amplification estimated for each of the clay-dung dominated mounds. However, to characterize the composition and dynamic properties of all the wierden in the region covered by the risk model would be a very major undertaking and a disproportionate effort in view of the small proportion of the total building stock founded on these structures. A generic penalty function, quantifying the additional frequency-dependent amplification of response spectral accelerations due to the presence of a wierde, was defined for application at all dwelling mounds. This is a slightly conservative approach but consistent with the geographical spread and purpose of the risk model, which is not to estimate risk to specific individual buildings but rather to quantify the field-wide risk in terms of potential damage and injury. We also acknowledge that the model is a simplification in terms of being based on 1D site response analyses and not taking account of potential topographical amplification effects, especially close to the edges of the mounds. However, such refinements would be challenging to incorporate and probably unwarranted. The risk model already accounts for spatial variations in near-surface geology, over a large region, and non-linear site response conditioned on actual realisations of motion at the reference rock horizon (and in some regards is therefore more sophisticated than many site-specific models); some degree of approximation in the modelling of these geotechnical structures on which less than 2% of the total are founded is entirely appropriate and acceptable.

In addition to modelling the increase of the shaking hazard due to the presence of the wierden, work was also undertaken to explore whether the soft soils of which these structures are composed would require modification of the fragility functions for the building typologies identified as being founded on dwelling mounds. The SSI analyses performed using representative foundation profiles for buildings on wierden found that in most cases the fragility functions were essentially unchanged. For two building classes, the soft wierden soils were found to slightly improve the fragility functions but the changes were modest, and in both cases would apply to less than 0.4% of the total number of buildings in that class. Consequently, it was concluded that there was no justification for defining new building classes with alternative fragility functions. The influence of the wierden is therefore accounted for only by identifying the buildings located on dwelling mounds and applying the penalty function to the calculated ground surface motions. This is an important refinement to this region-specific induced seismic risk model and means that the risk calculations will now take account of the presence of these uniquely Dutch features of the build environment. As with many aspects of the Groningen seismic hazard and risk modelling project, the work of incorporating the wierden has benefited from multi-disciplinary collaborations, including not only seismologists, geotechnical and structural engineers, but in this particular case also geographers and archaeologists.