We present an overview of induced seismicity due to subsurface engineering in the Netherlands. Our overview includes events induced by gas extraction, underground gas storage, geothermal heat extraction, salt solution mining and post-mining water ingress. Compared to natural seismicity, induced events are usually small (magnitudes ≤ 4.0). However, due to the soft topsoils in combination with shallow hypocentres, in the Netherlands events exceeding magnitude 1.5–2.0 may be felt by the public. These events can potentially damage houses and infrastructure, and undermine public acceptance. Felt events were induced by gas production in the north of the Netherlands and by post-mining water ingress in the south-east. Notorious examples are the earthquakes induced by gas production from the large Groningen gas field with magnitudes up to 3.6. Here, extensive non-structural damage incurred and public support was revoked. As a consequence, production will be terminated in 2022 leaving approximately 800 billion cubic metres of gas unexploited. The magnitudes of the events observed at underground gas storage, geothermal heat production and salt solution mining projects have so far been very limited (magnitudes ≤ 1.7). However, in the future larger events cannot be excluded. Project- or industry-specific risk governance protocols, extensive gathering of subsurface data and adequate seismic monitoring are therefore essential to allow sustainable use of the Dutch subsurface now and over the decades to come.

This paper presents a geoarcheological study on potential canal subsections present in the Roman-age Vecht branch of the Rhine-Meuse delta (the Netherlands).The first Roman canals in this delta were dug around 12 BC by Drusus, but their location has been the subject of debate since the 16th century, with various hypotheses proposed. Based on actual palaeogeographical knowledge of the Rhine-Meuse delta, the Utrechtse Vecht hypothesis is considered the most plausible. Within the study area, in the northern part of the Vecht system, natural sections of this river may alternate with possible artificial reaches, created at the time of Drusus.

Such artificial canals, being part of an otherwise natural channel belt system, can widen and deepen overtime, eroding all or most of the recognizable features associated with their original construction. As study area was chosen a relatively straight section of the Vecht between two former lakes. Two approaches were used. The first approach centred upon mapping channel morphology and recording sediment stratigraphy of the river deposits through detailed auger coring. Results corroborated the hypothesis of an originally straight feature (landform), confirming that it might have started life as a dug course, but not providing preserved archaeological remains of this stage. The second approach was chronological, whereby a programme of 14C dating was undertaken to refine the understanding of the origin and development of this reach of the Vecht, allowing earlier chronological investigations to be further contextualised and reassessed. A significant challenge to understand age control and floodplain evolution is the degradation of the top of the clayey peat that was observed below the levee deposits; this degradation is due to the lowering of groundwater levels and causes the end of peat growth to be dated as older than it actually is.

Using new radiocarbon dates we have reconstructed that the Overmeer-Nigtevecht reach of the Vecht between two former lakes started life as a straight channel. We have constrained its age to be closer to the time of Drusus’ activities (early Roman age). Although we have not found in situ remains of Drusus canal(s), these two new insights make the Vecht option, effectuated by a series of short canals, more likely to be the Drusus canal(s).

The reconstruction of submerged palaeolandscapes and detection of settling surfaces along coastal zones became a major research topic within the last two decades. In this context, the WASA project made use of a multidisciplinary approach for defining the extension and describing the characteristics of the Late Pleistocene to Holocene deposits in the central Wadden Sea region. In addition to sub-bottom transects, more than 140 sediment cores were taken in such an area for stratigraphic reconstruction, making use of multi-proxy analysis. To harmonize these data with the existing regional core database (LBEG archive) a new core catalogue was developed, that allows the identification of the local Late Quaternary sedimentary sequences and their characteristic facies. The WASA core catalogue has been successfully applied for reviewing the published data about the stratigraphic sequence of the Wadden Sea, for a better definition in terms of stratigraphic sequences and spatial extent of the Quaternary geological evolution of the region, and for a detailed reconstruction of the coastal palaeoenvironments.

Discussions on the age and the depositional environments of the Veldhoven Formation and its members are persistent in Belgium and the Netherlands. Uncertainties on stratigraphy and the constructive process of sediment accumulation continue today as a result of lack of data on this succession within the Roer Valley Rift System. The present study provides new information on the bio- and lithostratigraphy and facies from two boreholes based on dinoflagellate cyst taxa. The results were correlated by gamma-ray logs towards other key boreholes in the region and show a good consistency for stratigraphy and depositional environments for the members of the Veldhoven Formation.

After marginal to restricted marine conditions in the latest Rupelian (early Oligocene), the start of deposition of the Veldhoven Formation marked the transition towards a higher sea level, expressed by increased glauconite contents and gamma-ray values. The Voort Member in the lower part of the Veldhoven Formation has an early to late Chattian (Late Oligocene) age and comprises predominantly shallow marine (fluctuating restricted to open marine) conditions. The lithology in the lower part of this unit is often very clayey but is coarsening upward into sands. The superjacent Wintelre Member has a latest Chattian to early Aquitanian (early Miocene) age. This member is distinct by its clayey nature which is expressed by relatively high gamma-ray values. Earlier studies suggest a deeper marine facies for the Wintelre Member compared to the Someren and Voort members. However, the dinoflagellate cyst assemblages in this unit are mostly dominated by a single genus indicating a restricted marine setting, including salinities that deviate from normal marine conditions, most probably due to minor ventilation by narrow or lack of connection to the Atlantic Ocean. A glacio-eustatic sea-level fall around the Oligocene/Miocene boundary limited the sea coverage to the strongest subsiding areas, where deposition of the Wintelre Member is recorded, while non-deposition or erosion occurred in the surrounding highs, hence creating an isolated (sub)basin. The superjacent Someren Member was deposited during the late Aquitanian to middle Burdigalian and consists of shallow to open marine clayey fine sands. Increasing clay contents indicate a gradual development towards a higher sea level, which coincide with upward increasing gamma-ray values.

The biostratigraphic results of this study suggest that no major hiatuses are present in the differentially subsiding blocks of the Roer Valley Rift System during the late Oligocene to early Miocene.

Two new records from the Amersfoort glacial basin are investigated by means of pollen analysis. The cores are situated in the deeper part, close to the original Eemian stratotype Amersfoort 1 (Zagwijn, 1961) and at the margin of the basin. The aim is to reconstruct the Eemian and Early Weichselian vegetation development and to explore the impact of accommodation space, influx of allochthonous pollen and geomorphology on the vegetation composition. The results of the Amersfoort Basin are compared to the current Eemian stratotype in the Amsterdam Basin and other Eemian sites in the Netherlands. An almost complete Eemian to Early Weichselian sequence (E2-EWII) was retrieved from the deeper part of the Amersfoort basin. The late Saalian (LS) to early Eemian transition is not recorded in the Amersfoort basin, in contrast to the deeper Amsterdam Basin. The basin marginal core Den Treek reveals a condensed late Eemian (E5-6) and Early Weichselian (EW I-II) succession showing the importance of accommodation space. The first impact of the Eemian transgression is registered at the E3 to E4a boundary in the Amersfoort and Amsterdam basins, and highest sea level is proposed at the end of pollen zone E5. Upstream in the Eemian delta, in the palaeo-Vecht valley and IJssel Basin, the transgression is recorded later.

The influx of reworked (allochthonous) pollen in clastic sediment units hampers vegetation and climatic reconstructions during the LS and Eemian. The early appearance of Picea in zone E4 and Abies in zone E5 in clastic sediment intervals can be related to long-distance transport by the river Rhine and redistribution in the Eemian delta.

Local vegetation development can complicate regional biostratigraphic correlations. Alnus, considered characteristic for the late Eemian (E5-6), shows large differences over short distances in the Amersfoort Basin, related to local alder growth since Eemian E3. Carpinus, diagnostic for pollen zone E5, shows high values in the basins adjacent to higher, well-drained ice-pushed ridges, but low values in low-relief environments. Salt- to brackish-water marshes were present during high sea level in zone E5 in the Amsterdam and Amersfoort basins, while further upstream in the Rhine delta brackish to fresh-water tidal conditions dominated.

In line with Zagwijn (1961), the E6 to EWI boundary is defined at the first opening of the vegetation cover with Calluna, Poaceae and Artemisia increase, often coinciding with a lithological change from organic to clastic deposition, reflecting increased landscape instability.

The cores from the Amersfoort basin reveal a complete Eemian to Early Weichselian record. It is suggested to define the boundary stratotype for the base of the Weichselian Stage in the Amersfoort Basin. The current stratotype Amsterdam-Terminal reveals a fully developed LS to Eemian transition and contains the boundary stratotype for the base of the Eemian Stage.

Herein, Ryazanian (Berriasian) macrofossils from three well cores in the Central Graben (wells B18-02, L06-02, The Netherlands) and on the Jæren High (well 7/7-2, Norway) in the southern North Sea region are described. Macrofossils are mainly represented by buchiid bivalves (Buchia volgensis) and ammonites (Surites, Lynnia and Praetollia?). The genus Lynnia is recorded for the first time outside its topotypical area, and its systematic position and stratigraphic ranges are discussed. Additionally, the studied core sections yielded coleoid remains and a single limid bivalve. Based on the stratigraphic ranges of key ammonite genera (Lynnia, Surites and Bojarkia), the zonation of the Ryazanian stage is reconsidered. Uppermost Volgian to Ryazanian ammonite faunas are quite consistent and diverse but showing a higher degree of similarity throughout the Panboreal Superrealm as compared to those from rest of the Upper Volgian and the Middle Volgian. Buchia volgensis is the only species known from the southern North Sea and East Anglia, which is in strong contrast to the high diversity of Buchia in East Greenland and the remainder of the Boreal Realm. We hypothesise that such differences in the distribution of ammonites and bivalves in general, and the absence of buchiid species other than Buchia volgensis south of East Greenland in particular, are the result of anoxic bottom water conditions in the southern Viking Strait. The unusually wide geographic range of B. volgensis, which is known from such distant areas as Mexico and the Crimea, suggests a potential higher tolerance of this species to adverse conditions.

We reanalysed the induced seismicity data from the Groningen gas reservoir. We used the well-maintained induced event catalogue of the KNMI. The distributions of seismic moments and interevent times show a power law behaviour over several decades, and we find that upon increasing the magnitude threshold, these distributions remained scale-invariant. Because of this scale-invariance, we can put a constraint on the average loading of the elastic energy within the reservoir, which upon reaching a critical value gives rise to the seismic events. We find that the elastic energy roughly increases proportional to time. We also propose a new machine learning approach for declustering the seismic events, separating correlated from independent events. We find that only few events are truly independent, i.e. exponentially distributed over time. There are also few aftershocks following an Omori-type power law. The bulk of events presents a Gamma distribution for interevent times. This gives us an indication that the rigidity in the reservoir is high, but whether this results in overall correlated events should be settled with physics-based arguments rather than statistical ones.

Pleistocene glacial stages were implemented into a 3D basin and petroleum systems model of the northeastern Netherlands to address the influence of low surface temperatures and the mechanical loading of ice sheets on the subsurface. Two ice sheet thickness scenarios were used based on published data. Overall, Quaternary glacial stages have a substantial impact on the temperature and pressure distribution in the subsurface. Subsurface temperatures are significantly reduced during glacial stages, leading to lowered present-day temperatures and a low geothermal gradient in the shallow subsurface. In deeply buried sedimentary formations, pressures build up with every glacial advance resulting in overpressures at the present day. Glacial stages do not directly influence the petroleum generation of petroleum source rocks in the area, but high pressures during loading might have impacted petroleum expulsion of the early mature Coevorden Formation. Hydrocarbon accumulations in the Lower Saxony Basin were simulated to investigate the possible effects of mechanical ice loading and unloading on hydrocarbon migration. A loss of Coevorden Formation-sourced hydrocarbons to the surface was calculated in the Lower Saxony Basin during the glacial stages, indicating an influence of glacial loading on the Mesozoic petroleum system.

The occurrence of felt earthquakes due to gas production in Groningen has initiated numerous studies and model attempts to understand and quantify induced seismicity in this region. The whole bandwidth of available models spans the range from fully deterministic models to purely empirical and stochastic models. In this article, we summarise the most important model approaches, describing their main achievements and limitations. In addition, we discuss remaining open questions and potential future directions of development.

The mid-Maastrichtian carbon isotope event (MME), dated at ∼69 Ma, reflects a perturbation of the global carbon cycle that, in part, correlates with the enigmatic global extinction of ‘true’ (i.e., non-tegulated) inoceramid bivalves. The mechanisms of this extinction event are still debated.

While both the inoceramid extirpation and MME have been recorded in a variety of deep-sea sites, little is known about their expression in epicontinental chalk seas. In order to study the shallow-marine signature of the MME in this epicontinental shelf sea, we have generated quantitative foraminiferal assemblage data for two quarries (Hallembaye, NE Belgium; ENCI, SE Netherlands) in the Maastrichtian type area, complemented by a species-specific benthic δ13C record. In contrast to deep-sea records, no significant changes in benthic foraminiferal assemblages and benthic foraminiferal accumulation rates are observed across the MME in the type-Maastrichtian area. At the Hallembaye quarry, the otherwise rare endobenthic species Cuneus trigona reaches a transient peak abundance of 33.3% at the onset of the MME, likely caused by a local transient change in organic matter flux to the seafloor. Nevertheless, high and near-constant species evenness shows that neither oxygen nor organic matter flux was limited across the extinction level or during the MME. Benthic foraminiferal data from the uppermost part of the studied section, above the MME, indicate a significant increase in food supply to the seafloor. Decreased amounts of terrigenous elements across this interval document a lesser riverine or aeolian influx, which means that the increased benthic productivity is linked to a different origin. Potentially, the continuous precipitation of chalk under nutrient-poor conditions in the Late Cretaceous chalk sea was enabled by efficient nutrient recycling in the water column. In shallower depositional settings, nutrient recycling took place closer to the seafloor, which allowed more organic matter to reach the bottom. These results provide insights in the importance of nutrient cycling for biological productivity in the NW-European chalk sea.

We address aseismic fault slip and the onset of seismicity resulting from depletion-induced or injection-induced stresses in reservoirs with pre-existing vertical or inclined faults. Building on classic results, we discuss semi-analytical modelling techniques for fault slip including dislocation theory, Cauchy-type singular integral equations and the use of Chebyshev polynomials for their solution and an eigenvalue-based stability analysis. We consider slip patch development during depletion for faults with zero, constant static and slip-weakening friction, and our results confirm earlier findings based on numerical simulation, in particular the aseismic growth of two slip patches that may subsequently merge and/or become unstable resulting in nucleation of seismic slip. New findings include improved approximate expressions for the induced seismic moment per unit strike length and a description of the effect of coupling between the slip patches which affects both forward simulation and eigenvalue computation for high values of the ratio of fault throw to reservoir height. Our implementation based on analytical inversion and semi-analytical integration with Chebyshev polynomials is more efficient and more robust than our numerical integration approach. It is not yet well suited for Monte Carlo simulation, which typically requires sub-second simulation times, but with some further development that option seems to be within reach. Moreover, our results offer a possibility for embedded fault modelling in large-scale numerical simulation tools.

From the middle of the 1st century AD, the Lower Rhine was part of the frontier of the Roman Empire (Limes). However, this Limes was not an impermeable line, but rather an open corridor that served as march area and as a trade and supply route for the adjacent military and civilian settlements as well as the hinterland. This required access to the river and permanent harbours. When planning military camps and towns along the Rhine, the existing topography and the challenge of a dynamic riverine landscape had to be taken into account. The prefered location for forts and towns were the raised edges of the Lower Terraces close to the undercut bank of a meander. For many years, research assumed that the river bend adjacent to a Roman site had to be an oxbow lake. The main argument being that only an oxbow situation would have protected the Roman sites against strong currents and the risk of demolition of land and settlements by the main stream (abandoned channel hypothesis). A re-evaluation of the extensive archaeological, sedimentological and palynological archive and published data from the Rhine plain, as well as of supplementary data from new boreholes, lead to a considerable increase in the previously incomplete knowledge of the fluvial history, especially regarding the timing of palaeomeander infilling. It showed that most of the meanders and river bends of the Rhine with Roman settlements on their banks not silted up before Late Antiquity (from the end of the 3rd century AD on). Before that those meanders were part of the main stream. The advantage of steep undercut banks was that a quay could be built on a location that ensured mooring all year round, even at low water level. However, it was necessary that stabilising bank protections were present in order to avoid shifting the course of the river with subsequent destruction of the infrastructure. Such an antique bank protection construction (a so called “Packwerk”) could be recognised in front of the Colonia Ulpia Traiana (CUT, Xanten). With this knowledge in mind, other excavated structures on the banks of palaeomeanders, previously mostly interpreted as Roman harbour remains, could be interpreted as bank protections.

At Wesel near the strategically important mouth of the Lippe, the archaeological, sedimentological and palynological data showed that a meander that had silted up in Prehistoric times (Bronze Age) had been reactivated in Roman times and silted up again in Late Antiquity. A man-made diversion (perhaps by building a groyne) of the main stream could have been responsible for this.

With the beginning of the Late Antiquity crisis of the Roman Empire at the end of the 3rd century, it presumably became increasingly difficult to maintain these water works. This enabled the river to regain its morphodynamics, by cutting of the meanders that were active during Roman times. We hypothesise that this increase in fluvial activities of the Lower Rhine from the end of the 3rd century onwards is due to an anthropogenic trigger: the collapse of the Roman Empire.

Two closely spaced geothermal doublets were operated in the Californië geothermal field near Venlo, the Netherlands. The geothermal wells target the Dinantian Zeeland formation below 2 km depth. For several years, hot fluid was produced from the Tegelen fault, a regional fault in the Roer Valley rift system, until a felt M1.7 earthquake led to the suspension of geothermal activities. The Californië showcase provides a rare opportunity to retrospectively evaluate the assessment and the management of induced seismicity risks for a geothermal project. A seismic hazard assessment was conducted at several stages of the project, and seismicity was continuously monitored with a local station network.

In this paper, we report on the characteristics of the induced seismicity and evaluate the findings of the seismic hazard assessments conducted prior to the earthquakes. Seismic hazard assessments were based on numerical simulations of subsurface stress changes associated with geothermal operations. A geomechanical analysis indicated that the mapped faults in the subsurface are likely to be critically stressed. The largest hazard was inferred to result from thermo-elastic stresses, originating from cold water injection close to the Tegelen fault.

Subsequent earthquakes predominantly occurred near a production well after stopping or reducing production. We attributed this observation to a thermo-elastic stress load caused by cold water injection close to the Tegelen fault, combined with a counter-acting stabilisation of the fault due to pressure depletion during production. This mechanism was consistent with the dominating mechanism considered in the preceeding seismic hazard assessments. Although geothermal operations have not resumed yet, the geomechanical analysis indicates that re-locating one of the injection wells further away from the Tegelen fault could provide an efficient measure for mitigating induced seismicity risks at Californië.

Long-term exploration of the Groningen gas field in the Netherlands led to induced seismicity. Over the past nine years, an increasingly sophisticated Ground Motion Model (GMM) has been developed to assess the site response and the related seismic hazard. The GMM output strongly depends on the shear-wave velocity (VS), among other input parameters. To date, VS model data from soil profiles (Kruiver et al., Bulletin of Earthquake Engineering, 15(9): 3555–3580, 2017; Netherlands Journal of Geosciences, 96(5): s215–s233, 2017) have been used in the GMM. Recently, new VS profiles above the Groningen gas field were constructed using ambient noise surface wave tomography. These so-called field VS data, even though spatially limited, provide an independent source of VS to check whether the level of spatial variability in the GMM is sufficient. Here, we compared amplification factors (AF) for two sites (Borgsweer and Loppersum) calculated with the model VS and the field VS (Chmiel et al., Geophysical Journal International, 218(3), 1781–1795, 2019 and new data). Our AF results over periods relevant for seismic risk (0.01–1.0 s) show that model and field VS profiles agree within the uncertainty range generally accepted in geo-engineering. In addition, we compared modelled spectral accelerations using either field VS or model VS in Loppersum to the recordings of an earthquake that occurred during the monitoring period (ML 3.4 Zeerijp on 8 January 2018). The modelled spectral accelerations at the surface for both field VS and model VS are coherent with the earthquake data for the resonance periods representative of most buildings in Groningen (T = 0.2 and 0.3 s). These results confirm that the currently used VS model in the GMM captures spatial variability in the site response and represents reliable input for the site response calculations.

Palaeoecological reconstructions in the Netherlands are commonly based on pollen and macrofossil analysis, but can be limited if the preservation of organic material is poor. Phytoliths, biogenic silica, do not have this limitation and preserve in settings where other macro- and microfossils do not. Little is known about how phytolith assemblages preserved in soils and sediments reflect the parent vegetation in north-western European systems, so it is currently difficult to contextualise past environments. Here, we characterise phytolith assemblages for soil samples recovered from three major vegetation types in the Netherlands to provide reference data for future reconstructions of past vegetation change. We collected 42 soil surface samples from forests, wetlands and agricultural fields across the Netherlands and characterised the phytolith assemblages they contained. We identified the different phytolith morphotypes and quantified the percentages and concentrations (#phytoliths/cm3 soil) in each sample. We used non-metric multidimensional scaling to assess the variation in phytolith assemblage composition within, and between, the three vegetation types. The phytolith assemblages analysed from the forests, wetlands and agricultural fields were clearly distinguishable from each other. Agricultural fields were dominated by four phytolith morphotypes of grass silica short cells (GSSCs): rondel (tabular), cross type 1 (>15 µm), rondel (elongated) and disturbance or crop phytoliths. Forests settings had significantly higher amounts of different arboreal phytoliths (large and small spheroid rugose) compared with other vegetation types. Wetlands could be identified by significantly higher amounts of Cyperaceae phytoliths (papillate) and other GSSCs (saddle and bilobates with thick castula). Phytolith assemblages could distinguish different subtypes of vegetation within forest and wetland areas, while differences between agricultural systems could not be identified. Our study demonstrates that phytoliths preserved in soils or sediments can be used to separate major vegetation types across the Netherlands. Thus, these results support the hypothesis that phytoliths can be used to infer past environmental conditions in palaeoecological reconstructions. We suggest that future work should: (1) focus on characterising which phytolith types are produced by the commonest tree, wetland, shrub and herb species in the Netherlands and (2) characterise phytolith assemblages across a wider array of vegetation types in north-western European systems to increase the capability for quantitative reconstructions using phytolith assemblages.

We perform statistical analyses on spatiotemporal patterns in the magnitude distribution of induced earthquakes in the Groningen natural gas field. The seismic catalogue contains 336 earthquakes with (local) magnitudes above $1.45$ , observed in the period between 1 January 1995 and 1 January 2022. An exploratory moving-window analysis of maximum-likelihood b-values in both time and space does not reveal any significant variation in time, but does reveal a spatial variation that exceeds the $0.05$ significance level.

In search for improved understanding of the observed spatial variations in physical terms we test five physical reservoir properties as possible b-value predictors. The predictors include two static (spatial, time-independent) properties: the reservoir layer thickness, and the topographic gradient (a measure of the degree of faulting intensity in the reservoir); and three dynamic (spatiotemporal, time-dependent) properties: the pressure drop due to gas extraction, the resulting reservoir compaction, and a measure for the resulting induced stress. The latter property is the one that is currently used in the seismic source models that feed into the state-of-the-art hazard and risk assessment.

We assess the predictive capabilities of the five properties by statistical evaluation of both moving window analysis, and maximum-likelihood parameter estimation for a number of simple functional forms that express the b-value as a function of the predictor. We find significant linear trends of the b-value for both topographic gradient and induced stress, but even more pronouncedly for reservoir thickness. Also for the moving window analysis and the step function fit, the reservoir thickness provides the most significant results.

We conclude that reservoir thickness is a strong predictor for spatial b-value variations in the Groningen field. We propose to develop a forecasting model for Groningen magnitude distributions conditioned on reservoir thickness, to be used alongside, or as a replacement, for the current models conditioned on induced stress.