This paper presents a machine learning-based approach for detecting ground deformation patterns from time series indicative of potential sinkholes in the mining district of Lim-burg, the Netherlands. Our research uses large datasets of un-labeled, high-dimensional InSAR time series from Sentinel-1 and RADARSAT-2 satellites to identify anomalies that signal ground subsidence risks. Our methodology involves an anomaly detection process using dimensionality reduction, clustering time series data and cluster labeling. One of the challenges we face in the Limburg mining district is that there are no annotated datasets and no target labels that are required for a supervised learning approach, therefore we label the clusters by expert judgement. By adopting unsupervised learning, we aim to uncover unknown patterns in the time series, grouping similar temporal features into clusters for further analysis.@en

In this study, we develop a methodology to evaluate and benchmark the European Ground Motion Service (EGMS) products. Our approach integrates Global Navigation Satellite Systems (GNSS) measurements and measurements at and around corner reflectors (CRs), targeting the three product levels of the EGMS. We employ statistical tests to identify significant differences in the comparisons, focusing on time series, velocity estimation, and geopositioning. The study applies this practical methodology to evaluate the first EGMS product from 2015 and 2021 over various European test sites with significant surface deformation. The results consistently demonstrate a solid time series and velocity correlation between EGMS and GNSS at the tested sites, with, e.g., most of the differences between velocities being below 2 mm/year. This methodology, implemented in Python notebooks, ensures reproducible results and allows for validating the EGMS over other locations.

@en

The European Ground Motion Service (EGMS) is part of the Copernicus Land Monitoring Service (CLMS) managed by the EEA (European Environment Agency) [1]. EGMS is based on the full resolution InSAR processing of ESA Sentinel-1 (S1) acquisitions over Europe (Copernicus Participating states) [2]. The first release or baseline includes ground motion time-series between 2015 and 2020. This open data service has yearly updates.The EGMS employs persistent scatterer (PS) and distributed scatterer (DS) in combination with a Global Navigation Satellite System (GNSS) model to calibrate the ground motion products. This public dataset consists of three products levels (L2a/Basic, L2b/Calibrated and L3/Ortho). The Basic and Calibrated product levels are full resolution (20x5m) Line of sight (Los) velocity maps from ascending/descending orbits. The Ortho product offers horizontal (East-West) and vertical (Up-Down) anchored to the reference geodetic model resampled at 100x100m.The objective of this paper is to describe the independent validation of this continental scale ground motion time-series dataset. The goal is to ensure that the EGMS products are consistent with user requirements and product specifications and they cover the expected range of applications. To evaluate the fitness of the EGMS ground motion data service 7 reproducible validation activities (VA) have been developed gathering validation data from 50 sites across 12 European countries. Last but not least, a software environment has been developed to carry out all the validation activities and describe/store the data coming from the validation sites.

@en

Reclaimed coastal plains often experience significant subsidence as a result of phreatic water level lowering, which induces oxidation of organic material, shrinkage of clay, and sediment compaction. A primary example in the center of the Netherlands, is the 'South Flevopolder' which was reclaimed in 1968 and transformed into an area for residential, industrial and agricultural use. The area subsided over 1.5 m since its reclamation and its surface is still lowering.The city of Almere, with roughly 200.000 inhabitants and a surface area of c. 250 km2, is situated in the South Flevopolder. Most buildings in the city are founded in deeper Pleistocene sand, whilst objects such as parking lots, sport fields and playgrounds are often unfounded and are directly situated on the younger Holocene coastal deposits. Currently, the unfounded objects show subsidence rates as high as 5 mm per year for which the different subsidence rates may be related to subsurface heterogeneities. The upper layers in the area are dominated by clay and sand, up to a few meters in thickness, which overly peat and highly organic layers. The lowering of the phreatic surface results in an erratic pattern of subsidence over the area.We present a workflow to disentangle and parameterize the different contributions of shallow subsidence from Interferometric synthetic-aperture radar (InSAR) measurements. InSAR measurements from founded and unfounded scatterers are separated with a dimensionality reduction technique, t-Distributed Stochastic Neighbor Embedding (t-SNE), followed by an automatic detection of clusters with Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN). We have limited ourselves to structures with a construction date of >10 years with respect to the first InSAR measurement date to reduce the effect of construction-related consolidation and isolate the effect of shallow subsidence related to reclamation and phreatic surface level changes. The filtered dataset represents the surface response of unfounded structures in the urbanized reclaimed coastal area.The subsidence processes are disentangled and parameterized with an Ensemble Smoother with Multiple Data Assimilation (ES-MDA). Additional input data for this method is provided by a phreatic groundwater level model and a voxelized lithological model from the surface towards the top of the Pleistocene sand layers.We show that the automated data selection method prevents bias by selecting unfounded objects and the proposed workflow can be of aid when studying shallow subsidence in urbanized areas, where most objects are founded below the level at which shallow subsidence takes place. The results of this study quantify the rate of the different subsidence processes on a spatiotemporal scale and thus provide insights for tailored decision making to mitigate subsidence.@en

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

@en

In autumn 2017 a network of 14 broadband seismic stations was deployed at the Theistareykir high temperature geothermal field (NE Iceland). This experiment was conducted as part of the current efforts to characterize the field's main structures, and possible short and long term stress variations due to the ongoing fluid injection and extraction operations which started in autumn 2017. In this work, we use two years of continuous seismic records (October 2017–October 2019) to compute a 3D shear wave velocity model of the geothermal field and to detect possible crustal stress changes related to the injection and production activities. From phase cross-correlations of the vertical component recordings, we measure the Rayleigh wave group velocity dispersion curves to obtain 2D group velocity maps between 1 and 5 s. Subsequently, we use a neighborhood algorithm to retrieve the 3D shear wave velocity model of Theistareykir. Mainly, two sets of elongated high and low velocity anomalies can be observed oriented in a NW/WNW direction, parallel to the lineaments of the active Tjörnes fracture zone. Velocity reductions west of Ketilfjall and at Baerjafjall could indicate the location of upflow zones of the magmatic reservoir or hydrothermal system. We analyzed the temporal evolution of phase and amplitude of phase auto-correlations using the stretching technique and discuss their behavior in relation to the geothermal field operations. We notice a slightly stronger long-term velocity decrease in the reservoir region compared to outer regions. This could be related to the mass depletion in that area (higher fluid extraction compared to the water reinjection). In summary, our findings show how a monitoring network can be set up to enable a detailed imaging and monitoring of reservoir behavior in general.

@en

Within the Netherlands, 24 geothermal systems have been realized of which 18 added to production of 3.7 PJ of heat in 2018. The main application domain is horticulture. The ambition of the Dutch government is to increase the production of geothermal energy in the Netherlands to 15 PJ in 2030, subsequently 110 PJ in 2050 and to expand the application domain to district heating and industry. Kramers et al. (2012) published an article on the direct heat resource assessment stating that the recoverable heat amounts to just over 63000 PJ, assuming a reinjection temperature of 35 °C. This resource estimate combines all direct heat projects regardless of their maturity. A state-of-the-art resource assessment can support the Dutch government policy to realize its ambition by classifying the resources based on project maturity, including reference to the obstructing issues in the maturation of projects such as technical, financing, licensing and/or social and environmental complexities. In January 2019, 78 exploration licenses are active or have been applied for. At least one geothermal project is defined in each of these licenses, at different levels of advancement. Outside the existing licenses and licenses applied for, The Netherlands provides ample space for geothermal development by geothermal projects, yet in a conceptual stage. At the moment, there are no tools or practical methodology to integrate all the available information in a coherent and automated manner. Furthermore, a method to quantify the geothermal energy resources of projects ‘not yet realized’ using the United Nations Framework Classification (UNFC) resource classification system, would help policymakers to decide upon the most appropriate measures to relieve the obstruction issues to reach the 2030-2050 geothermal ambitions. In this study, we show the results of a resource assessment and classification system for geothermal energy while incorporating environmental and social issues in the classification, as required by the UNFC-2009. As basis for this resource assessment, we use the webtool ThermoGIS. ThermoGIS which uses information on the Dutch subsurface derived from the regional mapping by the Geological Survey and which is updated by incorporating newly gained insights and data, provides a comprehensive overview of the geothermal potential for a selected set of aquifers. The resource estimate, including the uncertainty range related to these figures, at the status date 1-1-2019 of the Dutch geothermal “project portfolio” is classified using UNFC, resulting in a set of resource figures per resource class according to project maturity: commercial, potentially commercial, non-commercial and exploration projects. Furthermore, some additional development options are introduced like installing heat pumps on relevant projects and stimulation on projects with poor productivity to increase the resource potential. Next level analysis envisages that future projects are given an attractivity attribute reflecting, amongst others, proximity to areas of favourable heat demand, proximity to heat network, which will guide the order of the exploration drilling sequence. An estimate on the geothermal portfolio maturation speed will then reflect on the minimum number of wells to be drilled to reach the geothermal ambition. @en

After decades of oil, gas, and coal exploitation, we have learned about some of the unpleasant aftereffects of subsurface resource exploration. Adverse long-term impacts, some known during exploration periods, others only afterwards, may include induced seismicity, land subsidence, or even sinkholes. While geothermal is currently seen as a sustainable source of energy, seismicity induced by inappropriate operational procedures or lack of knowledge of the subsurface may incite doubt and public sensitivity about its future use. A problem frequently posed before and during geothermal exploration is the cost of geophysical measurements for resource assessment, subsurface characterization during the prospection phase, and monitoring accompanying production. In this study, we investigate and discuss the potential of two economic geophysical measurement techniques for geothermal reservoir characterization and monitoring: passive seismic interferometry for better subsurface characterization through seismic imaging (static model), and satellite-based radar interferometry for geodetic imaging (dynamic model). Seismic imaging using passive seismic techniques allows for subsurface characterization via Ambient Noise Tomography, and supports the assessment of geothermal resources without requiring the use of shooting, which reduces the cost compared to active seismics. Geodetic imaging, by measuring the surface displacements during and after production, allows for the monitoring of the effects of production and constrains reservoir modelling, and can be achieved through the use of (freely available) satellite imagery. We discuss the results of both techniques over two high enthalpy geothermal sites in Iceland: Reykjanes Peninsula and Torfajökull volcano. While the Reykjanes Peninsula has geothermal fields that have been producing for decades, Torfajökull’s geothermal field, despite being the largest in Iceland, is not producing. For the subsurface characterization, we use S-wave velocity tomographic images derived from ambient noise seismic interferometry over the two geothermal sites. Within the tomographic images, low- and high-velocity anomalies are used to characterize subsurface structures, which complement current geological models with information at greater depths. From the monitoring point of view, radar satellite deformation measurements over both areas show displacements (subsidence) due to production (Reykjanes) and due to natural phenomena (Torfajökull). Finally, we summarize the lessons learnt and discuss outcomes on each technique.@en

The coastal plains of the Netherlands are subject to anthropogenic and natural subsidence with rates which are an order of magnitude higher than sea-level rise. Because one-third of the Netherlands lies below mean sea level, subsidence may threaten the country’s subsistence with major socio-economic consequences. Subsidence is a normal natural process but is overprinted and accelerated by anthropogenic activities which drives deep or shallow processes. Deep processes are caused by the extraction of hydrocarbons and salt, whereas shallow processes are primarily caused by lowering of phreatic groundwater levels. At present, the relative contribution of each process to total subsidence (i.e. natural plus anthropogenic) is unclear. Such information is important for stakeholders to support decision making on subsidence mitigation and it should be substantiated with independent scientific studies.

We present the outline of a hybrid Artificial Intelligence (AI) big data and model workflow to disentangle different subsidence forcing and we report on preliminary results for an area covering a gas field in the peat-rich Friesland coastal plain. The proposed workflow is a hybrid approach between a knowledge-based physical model and machine-learning techniques. The big input data comprises a suite of static (structural model) and time-dependent subsurface data (phreatic groundwater level, reservoir pressure), and geodetic measurements. Geomechanical models provide the connection between the drivers (groundwater levels and reservoir pressures) and the surface movement.@en

We investigate ground movements induced by the 8 February 2016, Mw=4.2 earthquake at the Los Humeros Geothermal Field (Mexico) using Sentinel-1 radar interferometry. Previous estimated focal mechanism solution based on seismic data with a hypocentral depth of 1900 m could not resolve the measured coseismic surface deformation pattern. In this study, we applied inverse elastic dislocation models to estimate the source parameters of the seismic event. Our models suggest the reverse reactivation of the Los Humeros normal fault at a shallower depth (<1000 m), with a more significant left lateral component below ~400 m depth. The occurrence of such shallow events at Los Humeros pose increased risks for the neighboring communities and infrastructure. Therefore, continuous monitoring of seismicity and cautious planning of field operations are crucial.

A NNW-SSE striking fault swarm, including the Los Humeros fault, acts as a major boundary of the subsiding area observed by InSAR time-series between February 2016 and May 2019. A potential explanation of the reverse reactivation of the Los Humeros fault and following downward movement of the eastern fault block is the depressurization of the whole hydrothermal system. Such depressurization can occur due to the exploitation of the geothermal field and/or due to natural pressure/temperature changes related to magmatic activity.@en

Surface deformation can be detected with satellite-based geodetic sensors, providing important insights on subsurface geomechanical properties. In this study, we use Interferometric Synthetic Aperture Radar (InSAR) observations to measure ground deformation due to fluid extraction and due to an induced seismic event at the Los Humeros Geothermal Field (Puebla, Mexico). The result of the PS-InSAR (Persistent Scatterer by Synthetic Aperture Radar Interferometry) analysis shows that the subsidence at the LHGF was up to 8 mm/year between April 2003 - March 2007, which is small relative to the produced volume of 5x106 m3/year. The subsidence pattern indicates that the geothermal field is controlled by sealing faults separating the reservoir into several blocks. Our models imply small volume changes in the reservoir and the different nuclei of strain solutions differ only slightly. These findings suggest that the pressure within the reservoir is well supported and that reservoir recharge is taking place. The occurrence of induced seismicity is independent evidence for the development of effective stresses. The inversion of the coseismic deformation due to the 8 February 2016, Mw 4.2 earthquake showed activation of a shallow normal fault in a reverse manner. Even with the limited compaction observed from the subsidence, the stress development was enough to induce seismicity. Pore pressure contrasts over a fault will increase the stress development. However, the fault structure will be more complex than assumed in the present study, because an acceptable joint match of ascending and descending surface movement data was not possible. An integrated study into the complex reservoir behaviour and the associated stress development over complex faults is therefore warranted.@en

Earthquakes carry important information about the current state of stress in the subsurface as well as information on the location of weaknesses. Energy exploitation activities and energy storage are inherently connected to changes in pressure in the subsurface and varying pressure rates are applied depending on the level of activity. Especially rapid changes in pressure are known to lead to induced and triggered earthquakes, and in some cases lead to reactivation of otherwise stable and unknown faults. In some cases, increased small magnitude induced seismicity is an indication of possible larger events to follow.Small earthquakes can be elusive and hard to locate precisely due to low signal-to-noise levels, an insufficient number of seismograph stations as well as over simplified methods and subsurface models. These challenges need to be overcome to be able to more accurately relate microseismicity to anthropogenic activities, and to be able to relate earthquakes to individual faults - both known faults and faults previously unknown. We explore ways to improve the hypocenter locations of earthquakes in exploited areas, principally by studying the effects of improving velocity models from standard 1D models to 3D models. When the geology is relatively uniform over the study area and the number of seismographs is abundant, using a 1D velocity model and a relatively simple method allows for fast and efficient processing yielding useful results. However, a 1D layered velocity model is not the best approximation where anisotropy is high, in areas with localized velocity anomalies, and where substantial velocity jumps in a layered medium control the differential arrival between P- and S-waves resulting in significant depth estimation uncertainties.We present results using the NonLinLoc software for locating small earthquakes in Denmark, Netherlands and Iceland using 3D velocity models. The location method provides good uncertainty estimates on the hypocenters. To test the quality of the hypocenters calculated with 3D velocity models, results from NonLinLoc is compared to existing hypocenter solutions for the same earthquakes. The results are evaluated with respect to deviations in hypocenters, uncertainty estimates provided by the different methods, efficiency and applicability for different geological settings. The focus of the Danish case study is the oil and gas fields in the North Sea, however due to sparse data and low seismicity the entire Danish region is included in the study. For The Netherlands the focus of the case study is the two decommissioned gas fields, Roswinkel and Castricum, where seismicity occurred after the end of production. The Icelandic case study focusses on the Reykjanes geothermal field located at the southwest point of the Reykjanes Peninsula. Most activity is natural but induced activity has been recorded near the production area of the Reykjanes power plant.The GeoERA-HIKE project has received funding from the European Union's Horizon 2020 research and innovation programme under agreement No. 731166 @en

After decades of coal exploitation during the early 20th century in Limburg, lagging aftereffects at the surface, such as sinkholes, are still experienced today. InSAR time-series have shown to contain information on accelerating surface deformation prior to sinkhole formation. In principle this information may be used as an early warning system. However, the signature for sinkhole deformation are often single point-scatterers which are difficult to identify due to the large surface area of the former mining concessions. Here we exploit georeferenced historical mining maps, an upward drilling data set and building locations to significantly reduce the search space for identifiers of future sinkholes. @en

The Netherlands is subject to anthropogenic and natural subsidence with rates for which the relative contribution of each process to total subsidence (i.e. natural plus anthropogenic) is still unclear. Such information is important for stakeholders to support decision making on subsidence mitigation and it should be substantiated with independent scientific studies. We present the outline of a hybrid Artificial Intelligence (AI) and knowledge-based physical model (geomechanical models) workflow to disentangle different subsidence forcing. We use static (structural model) and time-dependent data for both surface (geodetic measurements: GPS, levelling and InSAR) and subsurface (phreatic groundwater level, reservoir pressure). We show preliminary results of the approach for an area covering a gas field in the Friesland coastal plain and an area in the peat-rich central Rhine-Meuse delta plain. @en

Within the Netherlands, 24 geothermal systems have been realized of which 16 are operational, together producing 3.4 PJ of heat in 2018 (status date 1-1-2019). The main application domain is horticulture. The ambition of the Dutch government is to see production of geothermal energy increase in the Netherlands to 50 PJ in 2030, subsequently 200 PJ in 2050 and also to expand the application domain to district heating and industry. A state of the art resource assessment can support the Dutch government policy to realize its ambition. By classifying the resources based on project maturity including reference to the obstructing issues in the maturation of projects such as technical, financing, licensing and or social and environmental complexities targeted policy can be drafted. As of January 2019, 51 exploration licenses have been granted, and 31 more have been applied for (MEA2019). At least one geothermal project is defined in each of these exploration licenses, at different levels of advancement. Outside the boundaries of existing licenses and licenses under application, the Netherlands provides ample space for geothermal development. At the moment, there are no tools or practical methodologies to integrate all available information in a coherent and automated manner. As such, a method to quantify the geothermal energy resources of projects ‘not yet realized’ using the United Nations Framework Classification (UNFC) resource classification system would help policymakers to decide upon the most appropriate measures to remove obstructions for reaching the 2030-2050 geothermal ambitions. In this study, we show the initial results of a resource assessment and classification system for geothermal energy while incorporating environmental and social issues in the classification, as required by the UNFC2009. As basis for this resource assessment, we use the web tool ThermoGIS-v2.1 (Vrijlandt et al. 2019). ThermoGIS-v2.1 using information on the Dutch subsurface derived from the regional mapping by the Geological Survey. It is updated through incorporation of newly gained insights and data and provides a comprehensive overview of the geothermal potential for a selected set of aquifers. The resource estimates, including the uncertainty range related to these figures, at the status date 1-1-2019 of the Dutch geothermal “project portfolio” are classified using UNFC, resulting in a set of resource figures per resource class according to project maturity: commercial, potentially commercial, non-commercial and exploration projects. @en

Torfajökull volcano, Iceland, has not erupted since 1477. However, intense geothermal activity, deformation, and seismicity suggest a long‐lasting magmatic system. In this paper, we use ambient noise tomography to image the magmatic system beneath Torfajökull volcano. One hundred days of ambient noise data from 23 broadband seismometers show the consistent presence of double‐frequency microseism noise with significant power between ∼0.1 and 0.5 Hz. Beamforming results indicate microseism noise with persistent higher energy propagating from west and SE directions and apparent velocities below 3 km/s. We use ambient noise seismic interferometry to retrieve Rayleigh waves, and we introduce a method to estimate the reliability of the retrieved surface waves. We find stable estimation of surface wave phase velocities between 0.16 and 0.38 Hz. Azimuthal velocity variations show a trend of higher velocities in the NE/SW direction, the strike of the rift zone intersecting Torfajökull, and orientation of erupted lavas on a NE‐SW fissure swarm. Tomographic results indicate low‐velocity anomalies beneath the volcano caldera (between −5% and −10%) and even lower velocity variations in the southeast and southwest study area (below −10%), outside the volcano caldera. Low anomalies may indicate the existence of hot material, more prominent outside the caldera outskirts. High‐velocity variations (between 5% and 10%) outline the volcano caldera between 4‐ and 5‐km depth and more pronounced velocities (between 10% and 15%) up to 5‐km depth in the north of the volcano caldera. We interpret the former as possible caldera collapse structure and the latest as solidified intrusive magma from the old preferred magma paths.@en

Tomographic imaging based on ambient seismic noise measurements has shown to be a powerful tool, especially in areas like Iceland, where the microseism illumination is excellent. In this paper, we produce a 3D S-wave tomographic image over the western Reykjanes Peninsula high-enthalpy geothermal fields and evaluate the reliability of the tomographic results for different resolutions through simulated and real data. We use 30 broadband stations operating for approximately one-and-a-half year and apply ambient noise seismic interferometry for each station-pair. This results in empirical Green's functions in which especially the ballistic surface waves (BSW) are well resolved. The retrieved BSW exhibit a high signal-to-noise ratio between 0.1 and 0.5 Hz, and the beamforming analysis indicates an apparent surface-wave velocity of 3 km/s over a broad azimuthal range. For the tomographic inversion, we invert the estimated phase velocities between all station pairs to frequency-dependent phase velocity maps in four different resolutions (1, 2, 3, and 4 km) using a Tikhonov regularisation. With the estimated regularisation parameter per frequency per resolution, we invert simulated data for checkerboard sensitivity tests per frequency for different combinations of velocity anomaly sizes and resolutions.

Finally, after the inversion to depth, we detect S-wave velocity anomalies with variations between −15% and 15% with reference to an estimated average velocity using 1 km and 3 km of lateral resolutions and 1 km of vertical resolution. This study shows the potential of ambient noise tomography as complementary seismological tool for reservoir characterization.@en