Groundwater extraction could be a crucial driver of subsidence; a phenomenon that has had, and still has, implications for the landscape, infrastructure, and environment. Vitens, a Dutch drinking water company, expects an increase of 10% in drinking water demand in 2040 for the Netherlands. This may intensify Subsidence due to Groundwater Extraction (SGE) in the future. The objective of this research is to quantify the historical and future effects of groundwater extraction on subsidence at Vitens extraction site ”Groenekan” in the Province of Utrecht. In this thesis, a methodology is developed to assess whether Vitens can responsibly expand Groenekan groundwater extraction from 5.0 Mm3/yr to 10.0 Mm3/yr. The groundwater extraction is situated in a rural polder landscape where two subsidence-sensitive soil layers are present: a 1 m thick Holocene peat/clay top layer and a 12 m thick aquitard (Waalre clay) at a depth of -48 m NAP. Under this aquitard, Vitens currently extracts ∼5.0 Mm3/yr of drinking water from the second aquifer. Next, a groundwater model is used to study the lowering of the groundwater level and hydraulic head due to the extraction. Subsequently, the effect of both lowerings has been quantified for four subsidence processes, 1) Biochemical Degradation of Organic Material (BDOM), 2) shrinkage, 3) consolidation, and 4) creep. Groundwater lowering affects all four processes, whereas head lowerings only affect consolidation and creep. Occurred subsidence from the start of the groundwater extraction and for eight future scenarios is calculated with the aid of the groundwater model output, analytical approaches and a 2D subsidence model. The cumulative subsidence of these four processes is labelled as SGE. From 1961 to 2023, the groundwater level at Groenekan is lowered with 33 and 60 cm for 5.0 and 10.0 Mm3/yr extraction discharges respectively. In the second aquifer, hydraulic head lowerings of 1.5 to 2.9 m with a symmetrical influence circle are observed at Groenekan. In a radius of 500 m around Groenekan, significant groundwater level lowerings are found, which enforces compaction processes in the form of BDOM and shrinkage. For the period 1961 to 2023, assuming a 1 m thick clean peat Holocene top layer reveals ∼16 cm of potential BDOM. However, potentially ∼11 or ∼18 cm of shrinkage would occur over the same period assuming a clean clay or organic clay Holocene soil type respectively. Compression of the Holocene and Waalre layers due to consolidation and creep has contributed ∼10 times less to SGE around Groenekan compared to compaction due to BDOM and shrinkage. Within a 500 m radius of Groenekan, SGE is noticeable; beyond that, SGE is negligible. The maximum total subsidence through compression is ∼1.0 cm over a period of 62 years (from the start in 1961 to 2023). Generally, consolidation contributes 75% to the total subsidence due to compression, whereas creep is responsible for 25% of the compression. Sensitivity analyses reveal that the Holocene soil type and the thickness of the soft soils are the most sensitive parameters in calculating the final subsidence. Clean peat soil gives a high potential subsidence rate, whereas sand soils hardly subside. Variations in groundwater level and head lowerings are found to be least sensitive. Within the Groenekan system, the 1 metre Holocene thickness is the limiting factor in the groundwater level lowering, since no additional subsidence effects are found if the groundwater level drops below the Holocene layer. The future scenario of 10.0 Mm3/yr + an extreme climate shows ∼1.2 cm of consolidation and creep in the Holocene and Waalre layers over the period 1961 to 2100 (139 years). The base scenario of 5.0 Mm3/yr gives of ∼0.9 cm of subsidence, which results in ∼0.3 cm of additional subsidence. This extra subsidence value is negligible over a period of 139 years compared to other locations in the Netherlands with dozens of centimeters of subsidence. Though, BDOM and shrinkage can potentially cause centimetres of subsidence in the Holocene on the long-term due to groundwater level lowering. In conclusion, expansion of Groenekan from 5.0 to 10.0 Mm3/yr is considered to be responsible from a subsidence perspective as long as BDOM and shrinkage are exercised with extreme caution. Future studies on the effects of extended droughts and the surface water system on the groundwater level are suggested to better distinguish the total subsidence from SGE. Lastly, it is recommended to Vitens to extract groundwater from deep, confined aquifers with a sandy top layer for maximal SGE mitigation.

Hydrological modeling is used to estimate the states and fluxes in a system given inputs. This allows us to gain an understanding of the water system and make predictions such as flooding or drought. Models use parameters to describe a general process such as infiltration. These parameters can be measured, inferred or calibrated. Parameters represent local aspects such as infiltration into the soil and can differ for a catchment in rural Ohio versus mountainous Colorado. Even with the best local parameters, models will always fall short of reality since they contain simplifications. Data assimilation allows the uncertainty in the prediction of stream flow to be reduced or quantified. This thesis aims to answer the research question ’What do the adjustments in parameters and states imposed by data assimilation say about the deficiencies in the hydrological models of observed streamflow?’.

In this thesis, a data assimilation framework was developed for the eWaterCycle platform. This framework takes care of the intricacies of data assimilation for the user. Any model implemented with the Basic Model Interface can make use of this, models within eWaterCycle are focused on FAIR hydrology. The framework was implemented in such a way that any data assimilation scheme can be used, following the FAIR principles. The framework was verified to work with synthetic observations for both the HBV (Hydrologiska Byråns Vattenbalansavdelning) and the Lorenz model. The framework also makes applying data assimilation within eWaterCycle easier and scalable.

To answer the research question, the five-year hydrological response of 671 catchments in the United States of America was modeled. These catchments are from the CAMELS dataset, which is a collection of hydrological data with characteristics for catchments across the United States of America. The data assimilation scheme particle filtering was applied to the conceptual hydrological model HBV using eWaterCycle. The parameters and states of each experiment were stored and analysed. The data assimilation scheme requires hyperparameters. These hyperparameters were optimised using the first 26 catchments in the dataset. This single combination of hyperparameters was used across the 671 catchments.
Taking the best prediction of stream flow at every time step, the data assimilation experiment outperformed the calibrated model in 461 of these catchments. Part of this is due to the bias introduced through the use of one set of hyperparameters. Analysis shows that catchments of the same characteristics perform similarly. Catchments with a low mean stream flow, have a high spread of predicted observations causing very good predictions when selecting the best. Due to this bias, no real correlations can be drawn about the relation of background characteristics and model deficiency. Analysing three catchments on a catchment scale showed that step changes in parameters often occur around flood peaks. The data assimilation scheme adjusts the working of the model to better capture the observed streamflow.

To further incorporate the framework in eWaterCycle more testing should be done on 2D models. Deficiencies in the HBV model can be further analysed by optimising hyperparameters per catchment type or flow regime.

Using a field classification, I have assessed the likelihood that a gas field caused a seismic event. I have put the focus on all of the minor gas fields located in the Rotliegend sandstone formation in the Netherlands. These gas fields were divided into five geologically distinctly different regions: the North Holland Platform (NHP), the Groningen Platform (GRO), the Lauwerszee Trough High Pressures (LTHP), the Frisian Platform (FP) and the Lauwerszee Trough (LT). I have investigated which reservoir factors are most influential in deciding whether or not seismicity occurs. The stress ratio has proven to be the most critical parameter. Then, for the stress ratio, I have computed the stress ratio range for each region at which failure will occur. Above the stress ratio, seismic activity will occur, and below the stress ratio, no seismic activity will occur. To answer the research question "Can the reservoir properties of small gas fields in the Netherlands predict the occurrence of seismicity?", I have taken the following steps.
First, I have examined the probability of all small Dutch gas fields being responsible for a seismic event. Based on the distance between the gas field and the nearest event, as well as the presence of other fields in the surrounding area, I have derived a classification for the likelihood that the field was associated with induced seismicity. Second, I have run a sensitivity analysis to identify which parameter was most significant. I have accomplished this by implementing a semi-analytical model that computed and depicted depletion-induced stresses and fault slip along an inclined fault. The model calculated the depletion pressure at which seismic slip starts to occur, here called the onset pressure, based on reservoir data and fault characteristics. The reservoir data contains the compaction coefficient, critical slip weakening distance, dip angle, dynamic friction coefficient, initial reservoir pressure, Poisson's ratio, porosity, reservoir depth, reservoir thickness, shear modulus, static friction coefficient and stress ratio while the fault characteristics included the absolute fault offset and whether the fault is bounding or not.
Afterwards, I have plotted the onset pressure versus the relative fault offset and assessed the sensitivity of the onset pressure whilst changing input parameter settings, namely the critical slip weakening distance, dip angle, dynamic friction coefficient, initial reservoir pressure, Poisson's ratio, porosity, reservoir thickness, shear modulus, static friction and stress ratio. The most critical parameter has turned out to be the stress ratio. I have examined the stress ratio ranges for each field in order to assess whether this parameter could predict the occurrence of seismicity for an entire region. Some fields that do not suit the optimal regional stress ratio have been considered anomalies and investigated further. The main explanations for the anomalies have been geological complexity, assumed synthetic fault offsets, over- or underestimated offsets, unregistered seismic events and substantial overpressure.
The answer to the research question is, yes, it is possible to predict whether seismicity will occur in a small Dutch gas field located in the Rotliegend formation based on reservoir data and fault characteristics. A regional stress ratio has been determined for each of the five regions. The regional stress ratio is 0.51 for NHP, 0.58 for Groningen, 0.69 for LTHP, 0.46 and 0.47 for FP and 0.50 and 0.56 for LT. Salt layers have most likely contributed to higher regional stress ratios for GRO, LT and LTHP.

Drought and subsidence are two out of several water-related urban climate adaptation challenges many cities in the Netherlands currently face. Drought is expected to increase in frequency and extent due to climate change. Therefore, drought is likely to further pressurize (subsiding) urban areas in the coming decades. Although the impact of drought on soft soils and hence subsidence is described in academic literature, in-practice analyzes are limited. Given the expectation of increased drought impacts, as well as the current urge for new housing developments to combat housing scarcity, this research focuses on in-practice drought impacts on subsidence in two soft-soil urban areas in the Netherlands. This research’s objective is to gain insight into drought-induced subsidence in soft-soil urban areas in the Netherlands, in order to better understand drought impacts on subsidence rates in future housing developments on comparable soils. With a better understanding more appropriate site preparation strategies can be applied. The objective is endeavoured by an in-practice analysis of two study areas: a 90s neighbourhood in Diemen and a recently finished urban area in Kampen. Firstly, variations in relative surface levels, observed by InSAR (2015-2019), are compared for dry and wet periods. Afterwards, this subsidence data is separated based upon local characteristics (soil structure, pavement type, vegetation percentage and vicinity of surface water) to analyze their influence on drought-induced subsidence. Lastly, an expert questionnaire is conducted on suggestions for both mitigation and adaptation strategies. Subsidence is not a linear process in time: surface level movements vary due to the soil’s wetting and drying. The extent of soil compaction in dry periods is found to increase with intenser droughts. The extent of soil swell in wet periods is approximately similar for varying wetness extents. This research found that the severe drought of 2018 thereby caused a soil compaction to such extent it could not be balanced with subsequent winter swell, and hence resulted in approximately 1 to 1.5 millimeter of drought-induced subsidence. Drought-induced subsidence is thus a net result of the (change in) seasonal surface level fluctuation. This net result is assumed to be consequence of extensive groundwater level drops, although the exact share of processes and subsidence mechanisms could not be estimated. Moreover, surface level movements are found to be prone to a lag and difference in duration (hysteresis) in comparison to the start and duration of dry and wet periods. Various local characteristics are found to influence drought-induced subsidence. In general more compressible soils show slightly larger surface level movement between dry and wet periods. The (sand) cover thickness is found to be influential on which subsidence mechanisms are triggered during a drought: a cover thickness resulting in groundwater levels to drop to present clay/ peat layers causes shrinkage and/ or peat oxidation additional to clinch. Unpaved surfaces are found to fluctuate more extensively than paved surfaces, but this does not necessarily result in more irreversible subsidence. Furthermore, abundant vegetation might result in extra irreversible subsidence due to its extensive water usage in dense urban areas. Lastly, in the analyzed soft-soil areas the surface waters seem to influence groundwater levels and hence surface level movements only on short distance. Suggestions on feasible enhancing site preparation strategies are given based upon research results and experts’ opinions. The suggested mitigation strategies consist of two approaches in order to hamper the variations in effective stresses and hence minimize the seasonal surface level fluctuation. The first approach focuses on preventing extensive groundwater level variations: increasing storage and infiltration of water, reversed drainage, building crawl-space free and choosing vegetation types based on their water usage. Additionally, it is suggested to apply a sufficient cover thickness (at raises) if soft soil layers are near the surface. This prevents that future extensive groundwater drops within these layers, and thereby limits (seasonal) shrinkage and/ or peat oxidation. The second approach focuses on reducing the top soil’s weight, via lightweight materials or self-carrying constructions. Adapting to drought-induced subsidence starts with measuring/ monitoring surface level movements in order to analyze spatial and temporal trends. Additionally, drought is to be considered to greater extent in subsidence modelling in order to improve subsidence estimations. This can be done by applying variable or lower groundwater levels in estimations of the (change in) effective stresses, at calculations of consolidation or creep. Lastly, urban utility management should focus on long-term costs via e.g. Life Cycle Analysis, and on overlapping maintenance cycles of surface level raising and e.g. sewer pipe replacements. The most important conclusion derived on urban drought-induced subsidence is that despite individual drought impact on subsidence is limited to 1 to 1.5 millimeters, its seasonal occurrence continuously affects surface levels. Moreover, due to climate change drought is expected to increasingly impact surface levels in Dutch soft-soil urban areas in the coming decades. The suggested strategies mainly hamper variations in soil stresses, via fluctuating groundwater levels, and hence drought impacts on soft soils. These strategies help to limit future soil movements and hence result in more climate adaptable soft-soil urban areas. The focus of this research is on qualitative analyses of in-practice data such that significant processes have been disclosed, rather than statistically verifying the results. Consequently, the research initiates further specific studies on statistical verification of drought-induced subsidence, and moreover topics on its mechanisms; its spatial and temporal variation; its measurement and modelling; the influences of local characteristics hereupon; and the effectiveness of the suggested enhanced site preparation strategies.

Time-varying network data are essential in several real-world applications, such as temperature forecasting and earthquake classification. Spatial and temporal dependencies characterize these data and, therefore, conventional machine learning tools often fail to learn these joint correlations from data. On the one hand, hybrid models to learn from time-varying network data combine several specialized models able to capture these dependencies separately, thus ignoring their joint spatio-temporal interactions. On the other hand, state-of-the-art approaches for jointly learning time-varying network data do not exploit the useful prior provided by a graph-time product graph. This prior structural knowledge can aid learning and help the models improve their performance. For this reason, we propose a novel neural network architecture to learn from time-varying network data using product graphs and graph convolutions, thus exploiting this prior during learning. In particular, our architecture (i) learns the most suitable graph-time product graph to represent the time-varying network data and model the graph-time interactions; (ii) performs graph convolutions over this product graph to learn the graph-time interactions from data; (iii) employs graph-time pooling to reduce the dimensionality over layers. To the best of our knowledge, no research has yet attempted to capture spatial and temporal dependencies using graph convolutions over product graphs. Results on synthetic and real-world data show that the proposed method is useful in learning from time-varying network data for both regression and classification tasks. For classification, we cure a real-world dataset for earthquake classification and compare the proposed approach with state-of-the-art models. Results indicate that graph-aware models outperform graph-unaware models on this task. For regression, we evaluate the proposed method on two real-world temperature datasets and compare our architecture with graph-aware and graph-unaware models. Results on the smaller dataset show that linear models outperform neural-network models. On the larger dataset, we find that prior knowledge about graph-time interactions seems to be less beneficial in case of abundance of data. For the synthetic experiments, we find that (i) learning the structure of the graph-time product graph from data improves the performance compared to adopting a fixed type of product graph; (ii) learning sparser graph-time product graphs further improves the performance; (iii) the proposed graph-time pooling technique contributes to the model's generalization capabilities.

Predictions of the morphology of coastal areas are used to make decisions on coastal defence and nature conservation. To predict this morphology, simulations made by morphological models are used. To base decisions on this morphology predictions, we want the uncertainties in these predictions to be small. The goal of this research is to get a better insight in the uncertainties in a morphological model. By having a better understanding of these uncertainties, the decisions made using modelled predictions are stronger substantiated. The specific area focused on in this research is the Frisian Inlet, which is located between Ameland and Schiermonnikoog in the Wadden Sea. To achieve the goal of this research, data assimilation is used. Data assimilation combines data with prior knowledge of a model to find the distribution of probabilities of estimates of a true state. In this research, the used data are bed level measurements and the estimates are simulations for the bed level height made by Delft3D. Data assimilation methods make use of a distribution of model outcomes which should cover all possible outcomes. So, to set up data assimilation successfully, we want to use a parameter that induces a significant change in the bathymetry outcomes of Delft3D. Therefore, a sensitivity study is performed which considers the following six parameters: current related roughness, wave related roughness, wave-related suspended load sediment transport factor, wave-related bedload sediment transport factor, the transverse bed slope and the tidal amplitude. For each parameter, ten values are chosen to simulate. Which of the parameters induces the most change in bathymetry is assessed, using the difference between simulation result and the observation and the mean squared error skill score. The transverse bed slope shows the most variation and is further used in the data assimilation method, which is a particle filter. Hundred uniformly distributed values of the transverse bed slope, between 0.5-100, are used to create different Delft3D simulations that give a bed level prediction. This initial distribution is used for three different periods: 1970-1975, 1975-1979 and 1979-1982. These epochs are defined by the availability of bed level measurements. After one iteration of the particle filter, a new distribution of the parameter values is found. This is used in the next iteration in the same epoch. In each epoch, three iterations are made. The methodology as described in this thesis leads to a convergence of the initial distribution. In epoch 1 and 3, the distribution focuses on values of the transverse bed slope in the high range of the initial distribution. The resulting distribution of epoch 2 focuses on values around a transverse bed slope value of 70. This study shows that it is possible to get a better understanding of the distribution of \gls{alfa} values that leads to probable bedlevel predictions by applying a particle filter. The application of this method brings the simulations closer to each other, but the simulations did not get closer to the observations. The contribution to our understanding of data assimilation for morphodynamic models is that it can be used as a calibration tool for a specific parameter.

Research on sustainable landfill management has been studied since 30 years ago in the Netherlands, the principle of which is to reduce the emission of harmful substances from the landfill to the surrounding soil and groundwater. As for this purpose, the active treatment is applied on Wieringermeer landfill, meanwhile, the long-term monitoring of substance concentration is of great importance. The measurement frequency of chemical concentration is twice per month, which costs around 48000 / year. To save money by reducing the measurement frequency, SARIMAX model is studied as a tool of data interpolation. For this analysis, we currently focus on the concentration
measurements of chloride and ammonium. By comparing the SARIMAX interpolated data and the data with reduced size, the results indicate that directly dropping half of the measurements can be regarded as an acceptable way to reduce the measurement frequency, as the data properties are well preserved and the errors in estimating the mass of substances leaching out are in the acceptable range. However, interpolating using SARIMAX model doesn’t have significant improvements in preserving the data properties. Further quartering the data can lead to large deviations in data properties.

Injection and production of fluids into/from the subsurface has been known to trigger earthquakes, referred to as induced seismicity. This seismicity may occur when anthropogenically caused changes in the in-situ stress conditions result in reactivation of pre-existing faults in the subsurface causing slip accompanied by sudden release of energy. Several studies have numerically modelled the induced stresses due to production/injection in reservoirs of various geometries. In this report we present a simplified three-dimensional reservoir with a displaced fault and derive analytical expressions for induced stresses in and outside the reservoir due to production and injection of fluids. We use the calculated stresses for the three-dimensional model and analyse onset of slip across the fault. The research builds upon the analytical two-dimensional plane-strain analysis for induced stresses and slip initiation in \citet{Jansen2019} to which our work contributed. We reaffirm the findings from the plane-strain analysis in \citet{Jansen2019} and conclude that the effects of incorporating third dimension on induced stresses and slip behaviour are limited. We find infinite peaks in resultant shear stresses at the reservoir boundaries and observe a distinctly different pattern in induced stresses and slip behaviour between production and injection scenarios. In case of production, the slip patches are predicted to grow inwards into the reservoir initially until they merge, while for injection the slip patches grow separately into the overburden and underburden. The findings in this report are in agreement with the previous analytical and numerical studies on induced seismicity. In this research we also introduce geometrical complexity in the reservoir in the form of laterally varying height of the reservoir and we observe that the effects of variation in reservoir thickness are also minimal, however the induced stress patterns and slip initiation is significantly impacted by fault throw, initial stress conditions and fault frictional characteristics.

Understanding the permeability of the subsurface is a crucial step to simulate fluid flow in the subsurface. A parameter estimation problem for the flow equations can be solved to find the permeability. The robust identification of material parameters remains a significant challenge. In classical approaches, the non-linear least-squares problem is formulated as a non-linear optimization problem in which the partial differential equation governing the permeability field acts as a constraint. These approaches lead to large scale problems and are, therefore, computationally challenging. This thesis proposes a new approach based on mesh decoupling of state and design variables. This approach allows treating the design variables on various scales of resolution without comprising the accuracy of the state and adjoint solver. The method is implemented on a
one-dimensional and two-dimensional Poisson problem taken from the literature. The first numerical results show that the multilevel approach is able to accelerate the initial stages of the search procedure significantly.

Geological storage of COኼ is a crucial and upcoming technology to reduce anthropogenic greenhouse gas emissions. Due to the buoyant characteristic of injected gas, the security of underground storage is a major concern. To asses the security of COኼ storage, an accurate prediction of the gas plume behaviour is essential. In this study, a fully physical 2D model is developed to describe gas behaviour in a saline aquifer. In particular, we investigate the effect of gas impurities on injectivity, macroscopic dissolution and resulting plume migration. The model includes an investigation of 4-component buoyancy driven multiphase convective flow with miscible convective mixing. For an accurate description of the phase behaviour, a recently developed thermodynamic model based on a combination of cubic Equation of state with activity model has been implemented. The implemented thermodynamic model includes a specific description for the behaviour of water and an activity model describing component behaviour of the aqueous phase making this model more accurate than conventional cubic EoS. The phase behaviour based on this thermodynamic model and a consistent set of physical properties have been implemented in Delft Advance Research Terra Simulator, a new simulation framework developed at TU Delft. The results show that the presence of the gas impurities (tested here) have a negative effect on the solubility of CO2 which in turn reduces the security of gas storage in saline aquifers and increase the risk of leaking. This is the first time when this accurate physical model is applied for large-scale simulation of CO2 sequestration. Additionally a framework regarding the energy cost of the total CO2 sequestration process, from separation to injection, is implemented.

Earthquakes in the past few decades has questioned the safety of people and infrastructure in Groningen region and its surroundings areas. The excessive gas extraction from the subsurface has led to human-induced earthquakes in this region. Liquefaction is a phenomenon that occurs as a result of earthquakes, reducing the soil shear strength to zero and the soil in turn behaving like a liquid. It has been observed in other parts of the world in case of natural and human-induced earthquakes alike. Thus, it is necessary to identify the liquefaction prone regions and taking steps towards creating better designs, equipped to handle earthquake loads along with precautionary measures, to save existing building and infrastructure. In order to do so, it is important to examine the existing the methods, various parameters and factors that influence the liquefaction potential analysis. In this project, the different factors that could influence Liquefaction Potential Analysis are studied based on the existing CPT-based correlations between Fines Content (FC) and Soil Behavior Index(IC). The previously proposed correlation by Boulanger and Idriss (2014) is examined for Groningen Soils to evaluate if it needs any modifications. Depth of the sample, Distance between CPT and borehole, Grain Size Distribution and Geology of the sample are main factors considered for the study. Each of these factors are analysed by collecting laboratory samples of Grain size analysis test from Groningen region. The strength of the correlation between FC and IC is evaluated based on the factors. This study would be beneficial for geo-technical software developers, construction and design engineers who highly depend on correlations and it gives insight of how different the onsite scenario can be from the predicted values using correlations. It contributes to the future research on creating a data base of the all samples in the Netherlands, to indicate high risk regions. It also helps in answering if IC is a good parameter to consider while evaluating liquefaction potential for Groningen Soils.

Geothermal energy is a relatively sustainable energy source of which the essence is to extract heat from hot subsurface rocks. Circulating fluids serve as the transport agent of heat. The contact area between the fluids and the rocks is where the relevant heat transfer occurs, i.e., where the water is heated up. In some geothermal reservoirs, this circulation occurs naturally through porous matrix (mainly sediments) or through heavily fractured formations. Enhanced geothermal systems (EGS) are potentially favorable reservoirs where subsurface permeability is increased by means of artificial stimulation techniques. These stimulation techniques often involve hydraulic fracturing where a limited existing fracture network is expanded or “enhanced” by injecting fluids under high pressure conditions. While the geometry of the generated fractures influences permeability, the effect on heat exchange has received less attention. This thesis discusses the effects of fracture geometry on the heat transfer between solid and fluid. Along with laboratory experiments, numerical simulations were conducted. All investigations were performed on igneous granite rocks. Tensile fractures were generated to allow a fluid flow along the otherwise impermeable rock samples. Several parameters were varied throughout the experiments and simulations including volumetric flow rate, fracture aperture, rock temperature and fracture geometry and surface area in order to investigate their impact on heat transfer processes. Flow rate variations in the experiments have shown that higher flow rates cause the fluid to absorb less heat per unit volume and cause the rock to cool down more extensively, therefore thermal depletion of the reservoir is likely to occur within a shorter time frame. The dependency of exchanged heat on fracture aperture variations (in the range of 0.05 to 0.5 mm) did not yield a clear trend within the experiments, but does so in numerical simulations. Aperture variations in the numerical simulations did not cause notable differences in transferred heat as long as the volumetric flow rate is kept constant. However, as the fluid velocity is kept constant the amount of fluid flushed along the fracture per unit time is affected by varying apertures. This causes a difference in heat transfer as well. Increased fracture surface areas alone (more extensive topology/roughness) have shown a minimal impact on the heat production while a more extensive fracture network (additional branches) has shown notable enhancement in the amount of heat produced. Cooling behavior of the rock has shown correlations with Newton’s law of cooling and suggests a limitation of heat production by the heat conduction occurring within the rock.
Experimental findings cannot directly be compared with natural reservoir conditions. The reason for this is a thermal equilibrium that is achieved at each flow experiment, i.e., the heat withdrawn equals the heat resupplied by a heater. In natural reservoirs this is often not the case where a cold front propagates towards the production well and determines the lifetime of how long heat can efficiently be produced from a certain rock mass. This results in an unsteady heat conduction where the heat withdrawn does not equal the heat resupplied.

The subsurface of Groningen gas field composed of several faults. The continuous production has been resulting in several micro-seismicity activities, particularly for the past decades. One of the reasons for the production-induced-seismicity is fault reactivation at depth. The fault that initially in the non-active state becomes active and starts slipping due to pore fluid extraction. The modelling of fault reactivation induced by production can be simulated by coupling of flow and geomechanics.

The project studies the coupling of flow and geomechanics for an idealistic subsurface model with fault. The flow occurs due to pore pressure depletion and is predicted by solving mass conservation. The corresponding deformation is estimated by momentum balance equation.

The observation will be limited to rock deformation that is quantified by effective stress and displacement on the fault due to a pressure change in the subsurface. The simulation is performed using Stanford's Automatic Differentiation General Purpose Reservoir Simulator (ADGPRS).

In the proposed model, the flow equation is discretised using finite volume method, and poromechanics is discretised by finite-element Galerkin's approach. Both problems share the same grid model, and there is no error associated with the information exchange. The coupled problem is being solved fully implicitly where flow and mechanics formulation are solved simultaneously.

The case study for the simulation is a simplified Groningen subsurface model with a fault zone. The impact of variability in production dynamic on fault reactivation is studied. The variations in production rate and production strategy are applied in different formation configurations. It is found that the correlation between the production rates and the stress state in the subsurface depends on the formation offset. However, the stress field does not depend on continuity of the production scheme. Several aspects including the presence of Gas-Water contact in the model was studied in this project as well.

Process-based method forward stratigraphic modelling provides advantages in reservoir modelling by simulating the geological process mathematically, and the genesis of geologic formations over time (Michael et al., 2010). In spite of its advantages, Miller et al. (2008) have recognised significant challenges in process-based simulation models, one of them is the incapability to condition to subsurface data. The problem in conditioning the data can be addressed by using an alternative method named multiple-point geostatistics (MPS) in modelling the subsurface since its introduction in 1993 by Guardiano and Srivastava. MPS considers the relationship between multiple data points that is different from the conventional geostatistical methods that are commonly limited to using a linear relationship between data (Guardiano and Srivastava, 1993). By using the process-based simulation model as training image (TI) for MPS simulation, MPS should be able to address the conditioning issue in process-based simulation model.

Over the past decade, most research in MPS has emphasized on new algorithms for improving efficiency of MPS (Mariethoz and Caers, 2014; Mariethoz and Lefebvre, 2014), but there are still issues remain for the workflow to be widely used in geosciences. Furthermore, using nonstationary TI such as process-based simulation model in MPS are still an issue because the workflow is always different for specific cases. Until recently, little published works are available in applicating nonstationary TI in MPS.

The principal objective of this project was to determine an optimised methodology that allows the use of nonstationary process-based simulation model for TI input with MPS simulation in the fluvial-dominated delta. The process-based simulation model used in the study is a post-processed data from numerical model done in process-based modelling software Delft3D (Lesser et al., 2004) which is the PhD work of van der Vegt in 2018. There were two different cases utilised in this study that represents the whole delta development: Case A with high repetition in the patterns and Case B with low repetition in the patterns.

In order to achieve the desired outcome, this study links process-based simulation model with MPS using unconditional and conditional 2D MPS simulation with two different approaches: zonation approach and control map approach. The realisations from the unconditional simulation have to be validated until successful unconditional MPS simulation. The conditional MPS simulation were carried out when the unconditional realisations have been acknowledged as the approved results in mimicking the patterns of the Delft3D model. Lastly, the results were evaluated with four methods: connectivity function, E-type models, conditional variance models, and analysis of distance (ANODI).

In all of the MPS simulation results, the use of control map approach with unilateral simulation path proved to deliver better realisations for unconditional and conditional MPS simulations. Also, this study has presented an optimised workflow of 2D MPS simulation on using process-based simulation model in fluvial-dominated delta environment as TI based on different conditions of patterns’ repetition and hard data distribution.

Growing interest has been expressed the last decades regarding the Late Cenozoic, gas-bearing sediments of Southern North Sea Basin. Comprising the shallow subsurface of the Dutch offshore sector, these sediments have
been deposited by Eridanos fluvio-deltaic system, draining the Fennoscandian and Baltic shield through the present Baltic Sea. Three successful producing fields - A12-FA (2007), F02a-B-Pliocene (2009) and B13-FA (2011)- and
five under development have triggered the conduction of several studies in the offshore area. However the deltaic environment has been characterized as highly complex owing to its influence by the onset of Northern Hemisphere
Glaciation and thus the processes that governed the system have not become entirely understood. Enhanced cooling followed by the waxing and waning of the glaciers during warmer periods had an immediate impact on sediment
supply, accommodation space and mineralogical input. Although the sediments have been studied in terms of chrono-biostratigraphy, no systematic investigation with respect to the three main aforementioned factors as well
as the regional sequence stratigraphy and its link to the reservoir deposits has been made up to date. Therefore, this study employed sequence stratigraphy as a method to examine Eridanos conditions of deposition, investigate the
interplay between accommodation space and sediment supply and explore the nature of the gas-bearing reservoir sediments. The adopted methodology is comprised of two basic pillars, observation and interpretation. The first
records as clearly as possible the observations arise from the sequence analysis while the second extracts the meaningful information and interprets it in terms of temporal and spatial concepts. Using 2D seismic and well log data the basic observations were the stratal terminations, stacking patterns, seismic facies and the shoreline trajectory.
Seismic and well log interpretation showed that the delta experienced multiple events of sea level fall which forced the shoreline to regress basinward and caused sediment erosion or non-deposition. These events are bounded
by nine time-significant surfaces of subaerial unconformities which constitute the depositional sequences. Normal regression comprised the dominant depositional trend, combined with aggrading-prograding patterns, leading to
characteristic alternations of highstand and lowstand system tracts. Three main depositional environments which correspond to open marine turbidites (submarine lobes), delta front and delta top were identified from well logs,
core descriptions and the seismic facies analysis. The study suggests that the shallow gas is located in the alternations of silty-sandstones and claystones of the delta plain which comprise the vertical stacking of parasequence
topsets within the highstand- lowstand system tracts. Reconstruction of the relative sea level changes and sedimentation rates was made based on a technique introduced by this study. It uses the average thickness of each seismic
unit in order to interpolate time between the three known absolute ages obtained by the literature. The graphs showed that the accommodation space was generated by a low rate of sea level rise while sedimentation rates were
increasing over time. However the scarcity of time constrains in combination with the uncertainty in the estimation of seismic volumes resulted in a low resolution outcome. A comparison between the findings of this study and
those of the existing literature was made. The overall depositional trends and conditions seem to be in accordance with the other surveys. Nonetheless the proposed interplay between accommodation and supply can be assigned as
local since the examined area is limited compare to areas studied in literature. Limitations are traced in the quality of the 2D seismic data hampering the observation regarding the relation between the reflectors and the surfaces.
Consequently, the study offers an insight into the conditions under which the gas-bearing deltaic sediments were deposited and tries to place them in the established sequence stratigraphic framework. Also it provides information
regarding the distribution of the deltaic environments and identification of the reservoir rocks within each setting. The resulting interpretation can be used for prediction of the reservoir formations since the genetically-related
packages where they were identified are distributed in a predictable manner within a sedimentary basin.

There is an increased need for modelling the Dutch overburden to improve the prediction of earthquakes and subsidence due to the extraction of hydrocarbons. Depending on the surface location, the subsidence model of the NAM can under or over predict the amount of subsidence by 4 to 5 cm. This is around 15 % of the current maximum subsidence of 30 cm. The underlying assumption of the NAM subsidence model is that the reservoir, and the overburden above it, behave the same and thus have the same elastic parameters. The purpose of this research is to obtain a better understanding of the role overburden heterogeneities play in man induced subsidence. The main research question is: To what extent do overburden heterogeneities affect subsidence caused by reservoir depletion? Followed by the hypothesis: The heterogeneities of the overburden affect subsidence and should therefore be incorporated.
The results show that modelling the overburden as multi-layered, rather than as one homogeneous layer, leads to a difference of 0.75 cm. This is greater than the measurement error of the InSAR data, which is used to determine the subsidence of the Groningen field, and is significant when compared to the current discrepancy of 4 to 5 cm. Elastic parameters have been calculated from acoustic well data for geological units around the Groningen area. The calculated values differ from the values used by the NAM. When the calculated values are used for the Upper North Sea and Lower North Sea groups, the difference with the reference model is up to 2 cm for a reservoir radius of 3 km and 0.87 cm for a Groningen scale reservoir. The results suggest that the heterogeneities of the overburden affect subsidence significantly enough to warrant further investigation. Due to the fact that the obtained results are for horizontal layers, there is much potential for more complex overburden geology e.g. salt structures and non-horizontal layers.