Conduction-dominated geothermal systems are essential for decarbonizing the built environment, particularly in densely populated areas with high heating demand. Geothermal development in the West Netherlands Basin (WNB) has accelerated but still remains largely uncoordinated, following a "first-come, first-served" model which results in suboptimal subsurface resource utilization.

This study presents a multi-objective optimization approach for geothermal field development that simultaneously considers economic performance and reservoir longevity. The framework applies the Non-dominated Sorting Genetic Algorithm II (NSGA-II) to identify Pareto-optimal configurations for well placement and operational control. Objective functions, Net Present Value (NPV) and system lifetime, are evaluated through fully coupled geothermal reservoir simulations using the Delft Advanced Research Terra Simulator (DARTS). The framework incorporates constraint-aware optimization with regulatory compliance, enabling simultaneous optimization of spatial configuration and flow rates through capacity-dependent rate allocation.

The approach is developed and validated on a small-scale synthetic model before application to a realistic corner-point geometry model of the WNB incorporating heterogeneous fluvial architecture. Multiple well configurations (10, 12, and 20 doublets) are systematically evaluated across different geological realizations.

Results demonstrate that NSGA-II effectively identifies diverse Pareto-optimal solutions spanning NPV ranges of 0.8-1.6 billion euros and system lifetimes of 35-100 years. The analysis reveals that total injection capacity directly correlates with economic performance, with higher well-count configurations achieving superior NPV through increased heat extraction capacity. The optimization consistently reveals a distinctive spatial strategy where injection wells are positioned in the thickest reservoir regions with high-permeability zones, while producers balance maximizing distance from injectors with targeting high-temperature, high-permeability areas.

This framework provides quantitative evidence that coordinated planning strategies yield superior performance compared to the current "first-come, first-served" strategies. By applying multi-objective optimization to geothermal planning, the study advocates for the move towards coordinated, regional-scale planning strategies that enable more sustainable and economically superior use of subsurface resources.

This thesis presents a systematic experimental and numerical investigation into the thermal dynamics of CO2 injection into porous media, focusing on the Joule-Thomson (JT) cooling effect under conditions relevant to Carbon Capture and Storage (CCS) in depleted reservoirs. The purpose of the study was to deconvolve the mechanisms that govern cooling and to validate numerical models against laboratory data.

Experiments were conducted on two contrasting sandstones: a low-permeability (0.37 mD) Kentucky core, analogous to tight reservoirs, and a high-permeability (1-2 D) Bentheimer core. The cores were instrumented with distributed temperature and pressure sensors to capture transient thermal fronts during CO2 injection. In the low-permeability Kentucky core, the JT cooling was driven by the significant pressure gradient across the porous medium itself. In contrast, the high-permeability Bentheimer core required an artificial inlet pressure drop to initiate cooling; subsequent thermal transport was found to be dominated by the high cooling power of phase change (evaporation) rather than continuous isenthalpic expansion.

Across both types of rock, a systematic phase boundary offset was observed, with phase transitions occurring at pressures 2-4 bar lower (or 2-3 K higher) than predicted by bulk CO2 thermodynamics. Rigorous analysis demonstrates that this significant offset cannot be explained by classical confinement theories such as the Gibbs-Thomson or Kelvin effects, pointing to a more complex interplay of non-equilibrium thermodynamics, capillary phenomena, and rock-fluid interactions.

Numerical simulations using the Delft Advanced Research Terra Simulator (DARTS) successfully reproduced general cooling trends but highlighted critical model requirements. High-resolution Equation of State (EOS) tables (>2000 points) were essential for precision, while the absence of a correct CO2 gas-liquid relative permeability model limited the ability to capture the pressure build-up observed experimentally due to accumulation of the liquid phase.

The findings demonstrate that, while JT cooling is a measurable and critical process, its prediction requires models that incorporate pore-scale physics beyond bulk thermodynamics. The results provide a validated data set and a refined understanding of the coupled thermal-hydraulic processes that govern near-wellbore cooling during CCS operations.

The presence of wells in reservoir simulations can be represented by utilizing a well model within the finite volume framework of the simulator. A commonly applied approach is the Peaceman well inflow model, which places wells at the center of a specified grid cell of the reservoir model. Nonetheless, several studies have investigated the ability of placing wells off-center. Capturing the exact well position can significantly impact the reservoir response of the simulation model, particularly for heterogeneous reservoirs, as spatial heterogeneity plays a more pronounced role when wells are no longer located at the center of the grid cell. Furthermore, exact well position modeling would allow for more realistic well presentation within the reservoir simulator. Although the performed studies have successfully modeled
wells off-center, they are course-grid approaches. Therefore, this study investigates the application of an unstructured grid-based well model to accurately account for both the exact well location within a grid block and the surrounding reservoir heterogeneity in the calculation of a steady-state well index.

To achieve this, an improved well inflow model is established using Gmsh (three-dimensional finite element mesh generator) and the open Delft Advanced Research Terra Simulator (open-DARTS). The model implements an unstructured grid on a locally structured reservoir geometry, where typically five structured blocks are integrated. The model is verified against the Peaceman well inflow model for a center position and values match the Peaceman well model with a 5% accuracy. The model is validated for different well placements, permeability profiles and a considered number of structured reservoir blocks. Additionally, the validation study includes an accuracy assessment of the model against a full unstructured approach. The improved well index model is finally applied to the wells of the Geothermie Delft (GTD) project, enabling well index estimations based on true well trajectories.

Validation of the model shows that both well placement and local reservoir permeability impact the computed well index. A comparison of the proposed well index evaluation with the fully unstructured model indicates that incorporating more surrounding reservoir blocks within the well index calculation leads to closer correspondence with the fully unstructured approach. Additionally, the study reveals that the spatial discretization method within the finite volume framework used in open-DARTS influences well index values and can result in errors, depending on the choice of discretization approach. Finally, the results of the application of the proposed approach to the real trajectory of wells in the Geothermie TU Delft (GTD) project demonstrate that the improved well index model generally yields lower well index estimations than the Peaceman well model. The effect of this difference is further tested by simulating the full reservoir with the newly obtained well index values. These results depict that the well pressure profile for both the injector and producer varies depending on the well model, especially for later simulation times. The improved well model captures more variations in the pressure profile, emphasizing its capability for more accurate well modeling.

This thesis explores the occurrence of salt dry-out and hydrate formation when injecting CO2 into porous media. In large-scale CCS projects, injecting CO2 can potentially lead to salt precipitation or hydrate formation. These processes diminish injectivity and negatively alter reservoir rock properties. To gain deeper insight, experiments were conducted utilizing microfluidic setups, which allow for visual observation of salt-crystal or hydrate formation. Using microfluidic chips, ten salt dry-out experiments were conducted with varying pore sizes and six hydrate experiments were conducted with varying pulse trigger times. For the dry-out experiments, it was shown that salt crystals form mostly at the outlet side and that heterogeneity has a large impact on the precipitation process. A heterogeneous pattern results in a shift in salt distribution to the small pores, with results showing salt saturation at 7% in the small pore section of the medium-small pore chip, exceeding the 3% in the larger pore section. This shows the significant role of capillary action on salt precipitation. Results highlight that higher CO2 flow rates accelerate water evaporation and salt formation, yet final salt precipitation levels remain similar across varied flow rates. For instance, in the small pore size chip, final salt saturation was observed at 7% with a flow rate of 4.38 mm/s, decreasing to 4% at 0.78 mm/s. Additionally, the importance of high water saturation for salt dry-out and the impact of water backflow is shown. For hydrate formation, the importance of temperature and pressure was noted in these experiments. Three different pressure pulses were employed: manual control, a 0.5-second electronic pulse, and a 0.2-second electronic pulse, which all showing great effect on hydrate formation, yet no correlation could be determined between pulse length and hydrate saturation. Specifically, manual control yielded a 15% hydrate saturation with a 9.4% conversion factor, while the 0.5-second pulse achieved a 7% saturation and 9.1% conversion factor. The 0.2-second pulse resulted in 8% saturation and a 5.9% conversion factor. For the dissociation, the large effect of temperature was observed. All experiments showed a stable hydrate concentration, and a dissociation temperature between 5 and 9°C where temperature differences as small as 0.1°C were shown to be the difference between no dissociation and complete disappearance of all hydrates. Next to this kinetics, an interesting observation regarding the hydrate morphology was made. In addition to the five hydrate morphologies found in literature, a sixth, ‘sheet’-like type was observed.

Multiple studies have shown the potential for CO2 plume geothermal (CPG) to be a sustainable, reliable energy source that can be utilized in numerous regions worldwide. Compared to conventional
brine-based systems, a significant benefit is that CO2 allows for direct electricity generation at lower temperatures than brine. It could serve as both a continuous source of energy generation and a dispatchable source when energy demand is high. In addition, it could serve as a pre-carbon capture and sequestration (CCS) phase. Where it could verify the integrity of the reservoir and acquire information to characterize the reservoir and understand its behaviour under CO2 injection. Before a proof-ofconcept site can be chosen, candidate fields should be evaluated to find the optimal environment for CPG. In this work, we investigate which systems, aquifer or gasfield, injection-production scheme and what kind of environment provide the best performance for CPG. We use the Open Delft Advanced Research Terra Simulator (Open-DARTS) to simulate this on a reservoir scale. Open-DARTS uses the Operator-Based Linearization (OBL) approach to model all non-linear physics involved. To get an estimate of the electricity and heat generated by the system, we extend the open-DARTs framework to include a simple wellbore model and surface infrastructure. In our results, we look at the performance of two types of reservoirs: aquifers and gas fields. Where we consider the amount of electricity generated energy and other performance metrics. We find major differences between CPG performance in aquifers and gas fields. The results show that maintaining steady electrical energy generation through CO2 production in an aquifer appears to be much easier than it is for gas fields. With Aquifers consistently having a higher water cut. Furthermore, the inclusion of a plume establishment (PE) phase does boost performance once the CPG stage starts for both the aquifer and gas field types.

Induced seismicity refers to seismic events (earthquakes) triggered by human activities. Such events, even when characterized by relatively modest magnitudes, have the potential to jeopardize the safety of individuals, the surrounding environment, and infrastructure.

Production from hydrocarbon reservoirs can alter the in-situ stress state, leading to induced seismicity. This is reported in the Groningen field, where substantial gas production caused fault reactivation and subsequent earthquakes. Understanding events in the deep subsurface is crucial to proactively mitigate future seismic occurrences. To understand the causes of induced seismicity, the underlying physics are examined and defined in terms of relevant governing equations and models. This reveals the interconnected nature of fluid depletion, rock deformation, and fault slip. The goal of this study is to develop simulation techniques to solve these equations.

Towards this end, firstly, a finite volume embedded-numerical simulation method, called the Smoothed Enhanced Finite Volume method (sEFVM), is developed. This method is revealed to be computationally efficient for reservoir-scale modeling of heavily faulted systems and performed well in comparison to known solutions and other simulators.

However, in settings where analytical solutions indicated noncontinuous shear stress profiles, sEFVM accuracy suffers. Recognizing this limitation, a semi-analytical approach is developed, extending analytical expressions to be solved over the sEFVM mesh. This extension allows for more accurate solutions, accommodating complex reservoir and fault geometries. The semi-analytical method is successfully used to estimate the onset of fault nucleation and the magnitude of the seismic moment resulting from depletion.

The semi-analytical approach is limited to simulating fault slip up to the point of nucleation. To overcome this constraint, a hybrid method is developed. With appropriate assumptions regarding the post-nucleation state and the use of sEFVM to numerically calculate post-nucleation stresses, the hybrid method can effectively model multi-fault systems in the seismic stage assuming quasi-static behavior.

In summary, this research contributes by presenting novel computational frameworks for studying fault reactivation in faulted poroelastic media, offering insights into the complex interactions of the physics at play. The proposed embedded-numerical, semi analytical, and hybrid methods pave the way for further advancements in the field.

The increasing global demand for energy has necessitated the exploration of untapped potential energy resources. Unconventional oil reservoirs present a significant opportunity for development through various enhanced oil recovery (EOR) methods. Among these methods, the huff’n’puff scheme using Carbon Dioxide (CO2) EOR has emerged as a potent technique, considering both productivity and its impact on global warming. Unlike water injection, CO2 EOR exhibits high injectivity potential in tight formations, making it a promising approach.
This thesis focuses on investigating the performance of CO2 EOR in unconventional oil reservoirs using the huff’n’puff scheme through a comprehensive reservoir simulation study based on an open-source geological model. The study evaluates the effects of key parameters on oil production, including injection rate, injection time, soaking time, diffusion coefficient, and well completion.
The simulation results demonstrate that CO2 EOR can effectively enhance oil recovery from unconven- tional oil reservoirs. Specifically, the study reveals that implementing CO2 EOR increases cumulative oil production by 27% after 20 years compared to natural depletion cases. Furthermore, when incorporating a two-year primary recovery period, cumulative oil production is further elevated to 45%. By implementing the findings from the sensitivity study, cumulative oil production can be increased up to 67% in the optimum case.
In terms of the influence of the tested parameters on oil production, the study identifies that higher injection rates and longer injection, soaking, and production times for each cycle improve oil recovery until reaching an optimum value. Additionally, extending injection period by 10 times in each cycle with the same CO2 injected volume improved cumulative oil production by 7% compared to fixed injection rate as a result of longer diffusion and dissolution time. Moreover, implementing bottom CO2 injection targeting bottom layers improves vertical sweep efficiency due to buoyancy forces. The study also highlights the importance of the diffusion coefficient, as higher values facilitate faster CO2 transport, resulting in higher oil recovery and reduced CO2 reproduction.
This thesis provides valuable insights into the design and implementation of CO2 EOR projects in tight oil reservoirs. It emphasizes the need for careful selection of optimal operating conditions based on specific reservoir properties and highlights the significance of the diffusion mechanism in CO2 EOR performance. The findings contribute to the development of more effective CO2 EOR strategies for unconventional oil reservoirs, addressing gaps in previous research. Moreover, the study’s outcomes have practical implications for improving the design and implementation of CO2 EOR projects in the field, ultimately leading to increased oil production.

Borehole operations play a crucial role in managing various subsurface activities related to energy, including energy storage, geothermal energy production, CO2 sequestration, oil and gas extraction, wastewater disposal, and thermal recovery processes. In recent times, intelligent well technologies, such as long deviated multi-lateral wells equipped with advanced inflow control valves, have been employed to enhance field profitability and operational reliability. Additionally, different wellbore designs are utilized to efficiently extract heat from the subsurface in advanced geothermal approaches. To ensure accurate modeling, prediction, and optimization of all these energy production processes, it is essential to employ fully-coupled models that accurately capture thermal multi-phase compositional flow and transport in both the reservoir and the boreholes...

Compositional simulation is computationally intensive for high-fidelity models due to thermodynamic equilibrium relations and the coupling of flow, transport and mass transfer. In this report, two methods for accelerated compositional simulation are outlined and demonstrated for a gas vaporization problem. The first method uses a proxy model that reduces the number of components and the second method reduces the number of grid blocks (i.e. upscaling). Both methods are implemented within the operator-based linearization framework of the Delft Advanced Research Terra Simulator.

Lebesgue integration is applied in the loss function of a neural network allowing the neural network to discover the operator space of the reference model in reduced dimensions. Training is carried out for a one-dimensional homogeneous reservoir and minimizes the misfit of the leading and trailing shocks of a compressible pseudo-binary model with respect to observations of the reference model. The operator space of the pseudo-binary model is initially approximated with the method of multiscale reconstruction of physics, a numerical representation of the method of characteristics. Training is carried out in a two-stage transfer learning scheme to increase computational efficiency. In the first stage, neural networks are trained to approximate the analytical reconstruction. In the second stage, a solver is embedded in the loss function of the neural network and the forward solution is used to calculate the Lebesgue integral. The transfer training scheme minimizes the misfit of the leading and trailing shocks for 10 discrete time steps in a one-dimensional homogeneous reservoir. The misfit of the trained model shows a significant improvement in the location of the trailing shock and a modest improvement in the estimation of the leading shock. The trained proxy is applied to the top and bottom 15 layers of the SPE10 model and the estimation of the first and last breakthrough is assessed in conjunction with the error of the phase-state classification. The phase-state classification is significantly improved through time which is also expressed in improvements of the estimation of breakthrough times. The average difference in breakthrough time for the trained and untrained models with respect to the reference model is 293days versus 570days for the trailing shocks and 15days versus 16days for the leading shock. The established training framework enables the development of proxies with increased complexity.

Rigorous upscaling defines the upscaled operator space with dynamic and non-equilibrium thermodynamic upscaling functions. These functions combined, define the upscaled operator space for the three-dimensional compositional space and are inferred from data points gathered from a limited, characteristic portion of the full-size model. Gathered data points are interpreted with an interpolation function or neural networks to construct structured OBL meshes for implementation within DARTS. This upscaled operator space can effectively be used for different boundary conditions without reevaluating the upscaling functions.

In this work, a novel globally implicit framework for reactive multiphase flow and transport has been implemented using the simulator DARTS and the geochemical package Reaktoro. This framework applies the Operator-Based Linearization (OBL) approach together with a Gibbs Energy Minimization (GEM) scheme to model complex nonlinear reactive transport and flow. Using component- to element-based transport allows for numerically stable equations. This formulation is used to model the dry-out effect and subsequent precipitation of minerals in the near-wellbore region of subsurface aquifers when injecting CO$_2$. A kinetic rate in combination with Reaktoro equilibrium calculations is implemented to achieve more accurate precipitation results. The kinetic products are not considered inside of Reaktoro and also not as an element. The consistency at the phase transition boundary has been corrected by using a projection of a plane within the hypercube of zero charge balance. Highly sensitive reactions such as pH, require a high resolution to produce physical results. The parameterization grid can be implemented such that some species have a high resolution within a small variation of concentration. Reaktoro is highly accurate for the complex multiphase thermodynamic and chemical equilibrium calculation, but at phase transition boundaries, the equilibration of the state can be nonphysical. The results have been benchmarked against the existing geochemical simulator from Shell, showing that more physics needs to be implemented, but gives an overall accurate result.

Carbon dioxide (CO2) capture and sequestration (CCS) can play a significant role in reducing anthropogenic CO2 emissions while allowing society to slowly phase out traditional energy sources. One of the main challenge in CCS is the current absence of a clear revenue model. One idea is to use the sequestered supercritical CO2 as a working fluid for geothermal energy extraction from sedimentary reservoirs (Carbon-dioxide Plume Geothermal or CPG). Due to the large variability in density and mobility of CO2 under different temperatures, reservoir heterogeneities can give a rise to a combination of convective and conductive heat transfer. In this work, qualitative and quantitative descriptions
of the effect of reservoir properties on the performance of CPG in depleted
gas fields are provided using an example realistic depleted gas field. The primary focus is on the behaviour of the CO2 plume with regards to different reservoir properties such as porosity, permeability and thermal properties. The effect of large-scale reservoir structure, such as a presence of an aquifer, net-to-gross ratio and layering is also studied. In order to accurately model these effects, a thermal multi-component multiphase model based on a fugacity-activity Equation of State is built and validated for for pressures from 50-400 bar and temperatures from 35°C to 130°C. The developed thermodynamic model
is implemented into the Delft Advanced Research Terra Simulator. Numerous studies of 2D and 3D ensembles and sensitivity studies are carried out to examine the effects of isolated parameters on CPG performance. Results reveal that increased net-to-gross (N/G) ratio is associated with increased recovery factor. In addition, layering architecture becomes an important factor for the importance of conductive flux only at low N/G. Variations in the required pressure to sustain a production rate is associated with fluctuations
in production temperature and density due to expansive cooling. Varying reservoir properties and state also have a significant effect on brine upconing, which is detrimental to CPG performance. It appears that an increase in injection rate have a positive effect on the performance of CPG, but this should be studied in conjunction with a coupled wellbore and power plant model. Heterogeneous porosity-permeability realizations show a strong decrease in reservoir lifetime compared to their upscaled homogeneous counterparts, which is caused by a combination of preferential flow, reduced conductive flux and lower production BHP associated with the upscaled realizations. It was also found that reducing production rate delays the time of thermal breakthrough due to the combined effect of these factors.

The goal of this study is to test out the accuracy of the upscaling approach (Jinyu Tang, 2021) for a geothermal doublet. This approach tries to simulate the physical thermal dispersion. When modeling in a geothermal reservoir it would be easy if a lot of layers can be upscaled. However, when doing upscaling one can not just take an average of the layer properties. This results in a very inaccurate representation of reality. That is where the upscaling comes in. First, some other points need to be considered. The reservoir is simulated in DARTS (Wang, Voskov, Khait, & Bruhn, 2020). This simulator uses numerical methods to model the reservoir. Before starting on modeling geothermal reservoirs first a 1D case is evaluated to check for numerical dispersion that comes into play. With that evaluation done the 2D cases were simulated. The full reservoir has 91 layers, these can be upscaled into 9 layers eventually. To start building up to that scenario first 2 other scenarios are evaluated. First, an upscaled section of the first block of upscaled layers is evaluated. This was originally 10 layers and now upscaled to one layer. This one layer was first simulated in a one grid simulation and then a single grid simulation. These results are compared to a simulation of the original layer properties and a simulation of an arithmetic average of the reservoir properties. This is also done on a second scenario with 2 upscaled groups and eventually the full 9 layer reservoir (originally 91 layers). With these scenarios evaluated the upscaling shows to give a better simulation of the reservoir compared to using an average for the layer properties. This is all compared with a simulation of the original layer properties.

In this study, we introduce the Delft Advanced Research Terra Simulator (DARTS), which employs the Operator-Based Linearization (OBL) approach to provide efficient and precise solutions for non-linear discretized equations, to simulate the polymer flooding process. This compositional simulation framework is feasible for describing complex physical phenomenon with high nonlinearity because detached functional operators have been developed to release the calculation load while maintaining a high degree of stability and accuracy. The main contributions and innovations are shown below: i. A two-phase and three-component chemical flooding mathematical model with some reasonable assumptions was developed to investigate the enhanced-oil-recovery (EOR) mechanism of polymer flooding. A finite volume discretization scheme and a fully implicit method (FIM) are implemented to approximate the conservation equations. The OBL approach is used to reduce the complexity of the simulation framework. The correctness of the model is verified by comparison with analytical model. ii. In addition to shear thinning effect, polymer viscoelasticity has also been taken into account. This mathematical model is able to handle high shear rates during field displacement as well as low to medium shear rates in laboratory experiments. iii. A new model for porosity reduction caused by polymer adsorption is formulated. The effect of polymer concentration on the effective pore volume is incorporated in numerical simulations based on Langmuir adsorption model.

Energy transition extends the range of geological settings and physical processes to be taken into account in subsurface reservoir modelling. Many of these applications consider essentially anisotropic reservoir or require advanced gridding that can not be resolved consistently by conventionally used Two Point Flux Approximation (TPFA). In this project we present a Nonlinear Two Point Flux Approximation (NTPFA) based on gradient reconstruction and homogenization function. The approximation provides consistent solution for full permeability tensor on various grids. The approach combines flux guesses in a nonlinear way such that the obtained approximation is essentially monotone that guarantees the positivity of solution. We demonstrate the consistency of approach on several examples. We also use the multi-physics capabilities to test the simulator on saturation transport of dead oil when displaced with water. The developed approximation was implemented within Delft Advanced Research Terra Simulator (DARTS). Next we propose a new Nonlinear Two Point Stress Approximation technique which follows the collocated finite volume scheme for mechanical problem. In this section we try to discretize the linear elasticity equation by using nonlinear traction flux at interfaces similar to the setup used in fluid flow problem. This is done by balancing each component of traction individually and using the weighting scheme suggested in flux approximations.

To produce accurate 3D fluvial architecture models, it is important to understand the influence of internally generated autogenic controls versus externally generated allogenic controls. When studying a fluvial system, autogenic controls such as river meandering and avulsions generally dominate over 103-104 years. Allogenic controls such as climate change and tectonics dominate more at time scales of 105-106 years. A dominant control is astronomical climate change that may produce cyclical fluvial sedimentary successions. The resulting alluvial architecture of autogenic and of allogenic forcing is highly different and it is expected to largely influence geothermal flow through fluvial reservoirs. The aim of this research is to depict the impact of allogenic (cyclical) versus autogenic (non-cyclical) fluvial sedimentation on geothermal flow. Process-imitating and stochastic based modelling software Flumy is used to generate fluvial facies models where either hypothetical cyclical forcing or hypothetical non-cyclical forcing was the dominant force. These models are subsequently tested for geothermal flow using DARTS (Delft Advanced Terra Simulator). We find that the Flumy numerical model can be used to produce hypothetical cyclical and hypothetical non-cyclical alluvial architecture. The main difference between the architectures is the shape of the sand bodies. The cyclical model has overall thicker, laterally wider sand bodies. The non-cyclical model has thinner, less wide, but more often connecting sand bodies. Geothermal flow modelling shows that matching N/G in the cyclical and non-cyclical model in the 20-40 % N/G range gives similar pressures at a constant, fixed water rate. This points to the hypothesis that, at comparable N/G, the well connectivity must also be similar. The non-cyclical model breakthrough times in the 20-30 % N/G range are generally equal or slower, compared to the cyclical model. This gives credence to the hypothesis that the cyclical model has overall equal or slower flow paths between the well connections in the geothermal doublet. The difference in flow path hypothesis can be properly tested through visualizing 3D streamlines and is a recommendation for the future.

Acidizing is a widely used technique for intensification of geothermal wells. At certain conditions, acid injection results in a formation of highly conductive channels, e.g. wormholes. The associated rock matrix dissolution is highly nonlinear and may influence the fluid flow as well as rock properties. Wormhole patterns and breakthrough time are mainly controlled by acid strength and injection rate. Due to high complexity and, therefore, high computational cost, existing conventional simulators usually do not include complex chemical reactions or use weak (sequential) coupling of flow with reactive transport. Although the sequential approach is good for various applications, it is limited by the CFL condition and might suffer from poor convergence. In this project, we investigate the wormhole phenomenon to provide a stable fully implicit method (FIM) to numerically solve the reactive flow and transport problem coupled with equilibrium and kinetic reactions associated with the acid injection. The developed framework is aimed to predict a pattern of dissolution as well as wormhole breakthrough time. For an accurate description of species interactions, we directly connect PHREEQC equilibrium computation with FIM framework for kinetic dissolution using the Element Balance approach. This coupling was achieved with an application of a recently introduced Operator-Based Linearization scheme utilized within a Delft Advanced Research Terra Simulator (DARTS).

Evaluating an adequate computer model representing the characteristics of a subsurface reservoir is of great importance in the oil and gas industry. With recent data acquisition developments, high-resolution geological models can be generated for the subsurface reservoirs; however, uncertainties of geological data and too many details in the model have brought up remarkable limitations for application of high-fidelity geo-models in reservoir management. To generate efficient low-resolution models, which are relatively representative of the reservoir dynamic behaviour, data-driven models are proposed. The data-driven models are trained principally based on the production history of the reservoir with primary knowledge of structural information. The Data-driven Discrete Well Affinity (DiWA) model is designed to generate efficient proxy models that require basic geological information and full production history of the reservoir. Since it benefits the computational power of the Delft Advanced Terra Simulator (DARTS) and adjoint gradient-based optimization, the DiWA model has become an efficient framework in terms of accuracy and computational costs. This research shows that the DiWA model can relatively reconstruct the reference reservoir parameters out of an ensemble of priors, and the final results can be improved utilizing closer prior knowledge to the true response. Additionally, it is demonstrated that the performance of the DiWA model in a complex, realistic oil field can be enhanced using a comprehensive objective function, applying appropriate well controls, understanding the underlying driving mechanisms, using sufficient degrees of freedom, and applying consistent boundary conditions.

Water injection is widely used in the petroleum industry for the increase of hydrocarbon recovery or disposal of wastewater. Water production and injection are the primary mechanisms in Geothermal Energy. Both include injecting water into a porous formation under matrix injection conditions. While maintaining water injection is vital in these branches in the industry, so is the occurrence of formation damage (FD) due to suspended contaminants or brine incompatibility. Suspended particles in injection water are retained or deposited due to the creation of an External Filter Cake (EFC) or Internal Filter Cake (IFC) impairing the permeability. Consequently, FD results in Water Injectivity Decline (WID). In most cases, the negative impact on injectivity translates into operational and economic targets not being met. WID, as a result of FD, is highly connected to the Water Quality (WQ) of the injection water. A new approach to water quality (WQ) testing is proposed which suffices as a bridging application of membrane filtration and core flooding. To establish a foundation for this method, tests are conducted by performing particle-laden suspension injection experiments with porous outcrop sandstone 8 mm thin discs utilizing the ‘Con-vergence Hydra’. An experimental study is conducted investigating the effect of water quality (WQ) on formation damage, using dilute (20-100 mg/l) Baracarb2 (CaCO3) particle suspension as a model contaminant. Baracarb2 is tested for particle size distribution, mineral content and stability within synthetic brine (resembling Seawater). Subsequently, suspension flow experiments are conducted on porous thin discs (Bentheimer & Berea sandstone) as well as membrane filters (MF0.45μm). As performing suspension injection test with porous thin disc test utilising the Hydra has not been done before, reproducibility of the experiments is tested. The reproducibility of performing experiments utilising the Hydra is high, with very little difference between the experimental outcomes.Varying WQ within each different porous media type illustrates remarkably similar trends. Moreover, it is demonstrated that by performing a similarity curve collapse, a master curve is obtained for each porous medium type which scales with suspension concentration for Baracarb2. The damage mechanisms explaining this master curve all demonstrate a linear permeability impairment in the early part of the experiments. Subsequently, a linear impedance trend is observed where it is assumed that cake filtration is dominant. The latter is demonstrated by calculating the Modified fouling Index (MFI), which shows a linear dependency with suspension concentration. SEM imaging and Micro-CT scan images substantiate the damage mechanisms hypothesized from the pressure and rate data.The Con-vergence Hydra utilising porous thin discs has great potential for on-site testing which allows fast and reliable results on permeability impairment, i.e. formation damage. Subsequently, monitoring of the water quality can be done by performing MFI analysis during cake filtration.