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.
@en

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.

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. 

Transition to cleaner energy while sustaining the heat demand necessity has always been an emphasized topic in the recent years. With the slowly declining utilization and production of natural gas in Netherlands, geothermal energy shows a promising future in providing sustainable power for heating and various applications. For majority of economically feasible geothermal project, either natural fractures or hydraulic fractures are present. These fractures served as a preferential fluid pathway between geothermal well couplets, which determines crucial parameters such as breakthrough time and project life. In this study, the main focus is the impact of stochastic distribution of fracture aperture with or without pressure dependence and gravity effect on geothermal energy production. The simulation conditions as well as reservoir physical model are collected from analogous enhanced geothermal systems (EGS)that are classified as deep geothermal energy projects. The test cases include base case, case with no pressure or stress dependence, case with half well doublet distance, and case with gravity that has well positioned in the dip or strike direction in relative to the fault orientation. The ensemble simulation of 100 realization is run for each scenario considered. The relevant parameters analyzed are pressure difference between well couplet, production temperature, NPV, and most importantly, net cumulative energy production. It is found out that the geothermal energy production greatly relies on the fracture aperture distribution with or without pressure dependence. For base case, the range of energy production is from 1.50E13 to 3.75E13 KJ with temperature variation of 35 K between the minimum and maximum value. Without pressure dependence, only pressure difference between injector and producer slightly changes, and production temperature as well net energy pro-duction remains nearly identical. For the case with half of original well-spacing, the spread of energy production distribution reduces from 1.50E13 to 3.25E13 KJ. The diminishing energy production results from less total enthalpy in the flow path with shorter distance, yet the outcome of energy production still demonstrates a considerable range. By taking account of gravity effect and populate the reservoir with non-uniform pressure and temperature distribution with corresponding gradient, the net energy production drops to approximately two times less than the base case. However, the ratio of maximum energy production to the minimum over100 realization is similar to that of the base case. Last but not least, two types of well orientation, along the strike and along the dip of fracture plane, are considered to observe the impact of gravity-assisted flow. It is concluded that with the current model setting, the low temperature gradient of injector in the dip configuration case overweighs the effect of gravity assistance, resulting in a smaller net energy production than the strike case. This project concludes that the fracture aperture heterogeneity is an important factor in predicting the outcome of deep geothermal production. The gravity is another crucial component to portray accurate flow behavior between the wells, but it does not result in additional increase in the range of energy production outcome in comparison to the base come with the simulation conditions applied.

Before performing a field production forecast, an inverse problem has to be solved. Resulting in an ensemble of models that include the integration of real data with a complex physical and geological data describing subsurface processes. For large models, this approach can be very time and computationally expensive, therefore we propose an alternative approach for reservoir forecasting. In this work, we develop a physics-based data-driven model that purely relies on production data of the field and does not require any in-depth knowledge of the reservoir geology and physics. In the proposed approach, we utilize Delft Advanced Reservoir Terra Simulator (DARTS) as a base for our reservoir simulations. DARTS uses an Operator-Based Linearization technique for the approximation of exact physics. It allows us to encounter a more realistic interpretation of physics and is computationally efficient. The physics-based data-driven approach uses sequential regression to the data to increase the fidelity of the model forecast and encounter any significant changes in reservoir dynamics and physics over its history. The model was examined and validated for synthetic and real field production models. We demonstrate that the developed approach is capable of providing accurate and reliable production forecast on a daily basis.

The hot water produced from a geothermal doublet possesses energy, which once utilized, the water cools down and is re-introduced back into the same reservoir at a sufficient distance using an injector well. As cold water flows through the reservoir, it acquires thermal energy from surrounding in-situ rocks. This process recurs until a substantial drop in rock temperature occurs and the water is unable to be recharged adequately; as a result, cold water starts to "break-through" into the production well and the doublet soon needs to be abandoned.
This breakthrough time can be predicted using reservoir simulation. Significant work has been done in the past to determine the effect of different parameters on the breakthrough time and accuracy of models have been improved by incorporating real world physics. Despite being capable to predict breakthrough time, accurate high-fidelity 3D models require significant time in uncertainty quantification and data assimilation analysis due to CPU demanding simulations.
In this project, a physics-based proxy model is developed to predict flow and heat transport in low-enthalpy reservoirs. Streamlines that describe flow in a system and mostly controlled by steady-state pressure distribution are traced using Pollock's method (Pollock, 1988) and the reservoir is divided into streamtubes. Rock-heat depletion is modelled by semi-analytic model along streamlines. The objective is to predict geothermal doublet breakthrough time using only a limited number of streamtubes, thus minimizing simulation time and CPU resources.
Comparison with accurate high-fidelity model reveals that results for proxy model are optimistic; the error for pressure and temperature distributions, as well as the breakthrough curves is within the acceptable tolerance. The time required to simulate the proxy-model is less than the high-fidelity model. And as the number of streamtubes (to simulate the proxy model with) decrease, the time required for simulation further decreases but conversely the error between the breakthrough curves increases.

A coupled description of flow and thermal-reactive transport is spanning a wide range of scales in space and time. Subsurface reservoir heterogeneity with complex multi-scale features increases the complexity further. The spatial resolution required during a simulation is dependent on the type and degree of geological heterogeneity, but also on the flowphysics and thewell locations (Karimi-Fard &Durlofsky, 2014). Traditional upscaling or multi-scale techniques are usually focused on the accuracy of the pressure solution and often ignore the transport. Improving the transport solution can however be quite significant for the performance of the simulation, especially in complex applications related to reactive and compositional flow. The use of a method called Adaptive Mesh Refinement (AMR) enables the grid to adapt dynamically during the course of the simulation, which facilitates the efficient use of computational resources (Karimi-Fard & Durlofsky, 2014; Cusini et al., 2016). In this work, the aim was to develop such an AMR framework for general-purpose reservoir simulation. The approach uses a multi-level grid and is applicable to fully unstructured grids. The adaptivity of the grid in the developed AMR procedure is based on a hierarchical grid. The multi-level grid is constructed starting with the static geological model, a fine-scale model that accurately approximates the prospective reservoir. The control volumes present in the model are described by a list of volumes, depths and reservoir properties, and by a connectivity list for all neighbouring blocks, with corresponding transmissibility. The coarser levels are constructed through cell partitioning of the fine-scale model, called level 0. A global flow-based upscaling is applied to redefine the cell properties at the coarser levels. Each coarser level is also described by a list of volumes, depths, reservoir properties and a connectivity list with corresponding transmissibility. An inter-level connectivity list is constructed which includes the connections between control volumes of consecutive levels. The connectivity lists with corresponding transmissibility, along with the cell properties of each control volume are combined to formglobal arrays describing the full hierarchy of levels. The simulation is then started at the coarsest level, while keeping the wells at the finest level of refinement for more accurate representation of the solution. Dynamic adaptivity is performed, based on adaptivity criteria which were developed specific to the used application, resulting in a new grid configuration. The multi-level connectivity list togetherwith the cell properties are redefined at each time step and used as input for the next simulation. The developed AMR framework was implemented in Delft Advanced Research Terra Simulator (DARTS, 2019) which uses the recently developed linearization technique called Operator-Based Linearization, developed by Voskov (2017). The performance of the proposed approach was illustrated for several applications, including geothermal energy extraction and reactive transport. For geothermal applications, several models were tested: a homogeneous reservoir with unstructured grid, a synthetic fluvial model with shale blocks illustrating the role of shale in thermal recharge, a highly heterogeneous fluvial system with low net-to-gross ratio, and several layers of a highly-heterogeneous SPE10 reservoir. The results of these flow examples demonstrate a high level of solution accuracy relative to the reference fine-scale solution. The solution error is considerably decreased as compared to the coarse-scale model. Important features are captured at a higher resolution while keeping the amount of cells to a reasonable number. For the fluvial model, a different approach was adopted for cell aggregation, consisting of grouping the facies together, where shale-regions are kept coarse, and highpermeability facies are kept at a fine level. This significantly improved the solution as compared to the coarse model, where cold front propagation is overestimated. For the application of reactive flow, the dissolution of calcium carbonate through the formation of wormholes was analyzed. Here as well, the solution is considerably improved when using the AMR model. The wormholes are accurately represented, with improved computational performance.