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 stability of rock engineering projects is tightly related to the mechanical behaviour of rock discontinuities. Although the mechanical behaviour of a discontinuity is often associated with shearing (sliding), experimental data have shown that it can be accompanied by more complex phenomena such as dilation and post-peak strength reduction. Numerous researchers have contributed to the understanding and modelling of discontinuities behaviour. However, many of the constitutive models available in the literature are questionable when applied in practice either because they are highly empirical with parameters that are hardly determined or because they oversimplify the examined behaviour. Moreover, the models are often proposed for a certain range of stresses and specific stress paths, which makes numerical implementation difficult for large-scale engineering applications. This dissertation aims to evaluate the capabilities and the limitations of different constitutive models for rock discontinuity in the context of both numerical implementation and simulation of the mechanical behaviour of rock discontinuity. The first part of the research project investigates the main features of the existing models in the literature. From the investigated models, the two models highly adopted in research and engineering practice, namely the Coulomb’s model and Barton-Bandis’s model, are extensively investigated. Some enhancements and modifications are made to these two models to improve their modelling capabilities and ensure the numerical stability of numerical implementation. Regarding the Coulomb model, the adopted modifications include the reformulation of the model within the framework of strain softening providing a rigorous implemented version that describes the post-peak behaviour of a discontinuity adopting a linear reduction of the strength. Additionally, the employed modifications to the Barton-Bandis model provide a robust version of the model applying reformulations to the original yield surface that increase its validity in the whole range of the τ-σn space. Furthermore, a simplified definition of the post-peak behaviour, which aligns with the original formulation of the model but at the same time allows for a straightforward numerical implementation, is proposed. To validate the implementation of these models in PLAXIS (implementation done by the PLAXIS research team), the models are implemented in Python scripts for Constant Normal Load (CNL) shear test configuration. Concretely the implemented models are calibrated with experimental data to simulate CNL tests using a PLAXIS 2D Finite Element (FE) model and the obtained results are compared with both Python theoretical simulation and experimental results to verify the FE implementation. The results of these simulations validate the numerical implementation in PLAXIS and prove the applicability of the enhanced models to reproduce with adequate accuracy the mechanical behaviour of a rock discontinuity on a lab scale. Finally, the implemented constitutive laws are employed to perform a FE analysis of a large-scale application of a deep underground excavation in a discontinuous rock layer using PLAXIS 2D. To facilitate the creation of this complex geometry, an automatic discontinuity network generator is developed and improved using PLAXIS Python scripting API. The implemented discontinuity laws are then applied to the randomly generated discontinuity sets to simulate the behaviour of the rock mass. Stress and failure analyses are performed for the most critical discontinuities and wedges formed around the excavation to validate the numerical implementation and analyze the applicability of the constitutive models. The analysis of the boundary value problem confirms both the reliability of the numerical implementation and the applicability of the enhanced constitutive laws to simulate the analyzed large-scale problem.

Following the 2015 Paris Agreement, many countries are switching from fossil fuel to a more renewable energy supply. However, the CO$_2$ concentration in the atmosphere keeps increasing which makes difficult to accelerate the transition to a net-zero emission. In 2020, a new record of CO$_2$ concentration at 412.5 ppm was achieved, the highest concentration seen in the last 800,000 years. For this reason, several technologies have been in the focus to help reducing the amount of CO$_2$ in the atmosphere as well as keeping it from increasing.

One effective technique in the mitigation of these emissions is called Carbon Capture and Storage (CCS). It is not a new technology, in fact, it has been used in the Oil \& Gas industry since the early 1970s for the purpose of Enhanced Oil Recovery (EOR). In very general terms, when CO$_2$ is mixed with petroleum in the subsurface, its viscosity decreases making it easier to extract more oil. In addition, some of the CO$_2$ gets trapped in the rock which introduced the idea to use the mechanism to reduce CO$_2$ concentrations in the atmosphere.

However, this technique is not absent of risks and it is one of the main drawbacks experienced with activities that involved the subsurface, the surrounding uncertainties are too great sometimes. Some of the risks that could be encountered when producing, injecting or drilling a well in the subsurface are: under or over pressurization, leakage, uplift or subsidence, induced seismicity, and fault reactivation. In order to reduce the mentioned risks, geomechanical studies help with the understanding of the rock behavior and their response to applied stresses.

In this thesis, UCS and Triaxial laboratory experiments were carried out to test the geomechanical behavior of different rocks coming from Norway, Denmark and Germany. For the triaxial tests, five different cycles at different and increasing $P_c$ were performed. Doing so, some main insights were obtained: the modulus of elasticity increases with depth and confining pressure more significantly with $P_c$ lower than 30 MPa, and plateauing at higher $P_c$; the Young's modulus decreases with increasing porosity and increases with increasing cohesion; the P and S waves velocities change when the yield point is surpassed; with low porosity and high cohesion the acoustic velocities are not dependent on the $P_c$, and vice versa; the SYM and DYM follow similar trends; and the elastic modulus is negatively dependent on the Poisson's ratio.

The information and data obtained from this geomechanical study will be of use to model predictions involving CCS projects to avoid reaching the risks factors mentioned. An example of the lack of knowledge and understanding about rock behavior could be the induced seismicity caused by the extraction of gas from the Groningen gas field. These seismic events acted as a turning point, sparking a necessary shift in perspective and paving the way for a more informed and balanced perception of activities involving the use of the subsurface.

In the geothermal exploration well, NLW-GT-01, in the West Netherlands Basin, a very tight reservoir with low porosity (≤5.0%) and matrix permeability (≤ 0.1mD) has been encountered in the Main Buntsandstein Subgroup. The impact of natural fractures, which can significantly impact fluid flow in tight reservoirs, has been studied. However, there is a lack of consensus in different studies on the fracture classification.

In this study, a re-evaluation is performed of the well data of the NLW-GT-01 and VAL-01 wells. An analysis comparing core, geophysical, and image logs is performed to document the drivers and characteristics of natural fracture distribution in the WNB and investigate the geological and geomechanical parameters influencing the development of fractures in the Main Buntsandstein Subgroup.

On the core of NLW-GT-01 and VAL-01, five types of fractures are classified veins, joints, Mode II fractures, stylolites and drilling-induced fractures. The natural fractures have a dominant NW-SE strike orientation. The drilling-induced fractures and the borehole breakouts indicate a NE-SW oriented in-situ minimum stress and an NW-SE oriented maximum stress. The natural fractures are favourably oriented to be open based on their orientation compared to the in-situ stress orientation.

The fractures are associated with large-scale tectonic events. During the burial of the formation related to the extensional phase of the WNB during the Mesozoic, veins, joints, conjugate fractures, and stylolites could be formed. During the inversion of the basin in the Late Mesozoic-Cenozoic, tectonic stylolites could have been formed, and the previously formed conjugate fractures could have been reactivated.

Natural fractures are concentrated in more heterolithic intervals of the VAL-01 and NLW-GT-01 wells in both image logs and cores. These heterolithic sections display alternations of medium sandstones with silt- and claystones. The identification of more fractures in VAL-01 compared to NLW-GT-01 can be explained by the difference in basin location. VAL-01 is located in the centre of the basin where distal playa environments produced more fine-grained material alternating with coarser sands. The more proximal NLW-GT-01 is dominated by fluvial sands. The lithological variability produces Young’s and Shear modulus variability that seems to be driving increased fracture density rather than the
absolute value of the Young’s modulus causes this.

The upper Carboniferous in the southern North Sea, especially the Westphalian B and lower Westphalian C stage is characterized by the deposited cyclothems that are a series of coarsening-upward fluvial sediments cycles. These cyclothems are often coal-bearing and no more than 15 m. Studying the fluvial cyclicity has a good potential for improved stratigraphy because the cyclothems show cyclic signals and provide a possible way of high-resolution correlation. An integrated approach using multiple approaches to come to a solid outcome also enables improving the stratigraphy. The previous studies mainly focused on large-scale stratigraphy with a vertical resolution of a few hundred meters. This study has analyzed the fluvial cyclicity and constructed an improved stratigraphic framework of Westphalian B and lower Westphalian C to better understand the sedimentary system. An integrated approach including lithostratigraphy, biostratigraphy and Integrated Predictive Error Filter Analysis (INPEFA) is applied to the series. Six intermediate-scale biozones with a thickness of tens of meters to around 150 m are applied. Six intermediate-scale units named stratigraphic packages were recognized by INPEFA curves of the gamma-ray log. The correlations of the biozones and the stratigraphic packages are consistent with each other, and the offset between the boundaries of these two zonations is generally between 2-15 m. Individual intermediate-scale units have a highly variable thickness and the standard deviation of one unit can be up to 42 m. Constrained by the biozones and INPEFA stratigraphic packages, four main coal seams are correlated and their lateral extent can reach 15 km. The small-scale cyclothems that range from 7 to 15 m can be extracted by INPEFA of potassium content log, referring to the small-scale INPEFA cyclicity. The cyclothems and the corresponding INPEFA cycles are thicker and better developed in the southern area. The average thickness of small-scale cycles is 12 m and these cycles are related to obliquity with a duration of approximately 35 kyr. The duration of Westphalian B calculated by the number and the time span of obliquity-scale cycles is nearly 1.2 Myr. No long-period, 100 & 400 kyr eccentricity cycles were recognized in this study.