In order to characterise a rock formation prior to subsurface operations, it is required to find a microscale rock volume for which the homogenised property does not fluctuate when the size of the sample is increased; the Representative Elementary Volume (REV). Its determination usually comes at the cost of a large number of simulations, making it overall a computationally expensive process. Therefore, many scientific studies have been dedicated to optimising the process of finding REV. Using statistical numerical methods, it is shown that the fluctuation of the effective property corresponds overall to a cone-like shape convergence. We suggest determining the generic evolution law of the cone of convergence, which can be used to predict the size of the REV and the effective physical property. This study is based on simulations of Stokes flow through idealised microstructures from which the permeability is upscaled. By tracing and plotting the convergence of permeability for multiple samples, the full cone of convergence appears. The cone shows exponential growth and decay, converging towards the effective permeability of the microstructure. By fitting a log-normal distribution on the collected data points, we show that the generic evolution law of the cone of convergence can always be described with two parameters, independently of the porosity. We show that the determined law of the cone also applies to real microstructures, despite the presence of natural heterogeneities. The new method allows us to reduce the computational costs of finding all characteristics related to REV by simulating several subsamples rather than the full-sized sample, unlocking thereby high-resolution samples which are often too computationally expensive. The use of a statistical model provides quantification of the precision level we can obtain on the REV determination.

Faults can fundamentally change a groundwater flow regime and represent a major source of uncertainty in groundwater studies. Much research has been devoted to uncertainty around their location and their barrier-conduit behavior. However, fault timing is one aspect of fault uncertainty that appears to be somewhat overlooked. Many faulted models feature consistent layer offsets, thereby presuming that block faulting has occurred recently and almost instantaneously. Additionally, barrier and/or conduit behavior is often shown to extend vertically through all layers when a fault may in fact terminate well below-ground surface. In this study, we create three plausible geological interpretations for a transect in the Perth Basin. Adjacent boreholes show stratigraphic offsets and thickening which indicate faulting; however, fault timing is unknown. Flow modeling demonstrates that the model with the most recent faulting shows profoundly different flow patterns due to aquifer juxtaposition. Additionally, multiple realizations with stochastically generated parameter sets for layer, fault core, and fault damage zone conductivity show that fault timing influences flow more than layer or fault zone conductivity. Finally, fault conduit behavior that penetrates aquitards has significant implications for transport, while fault barrier behavior has surprisingly little. This research advocates for adequate data collection where faults may cause breaches in aquitards due to layer offsets or conduit behavior in the damage zone. It also promotes the use of multiple geological models to address structural uncertainty, and highlights some of the hurdles in doing so such as computational expense and the availability of seamless geological-flow modeling workflows.

Hydrogen is a promising energy carrier for a low-carbon future energy system, as it can be stored on a megaton scale (equivalent to TWh of energy) in subsurface reservoirs. However, safe and efficient underground hydrogen storage requires a thorough understanding of the geomechanics of the host rock under fluid pressure fluctuations. In this context, we summarize the current state of knowledge regarding geomechanics relevant to carbon dioxide and natural gas storage in salt caverns and depleted reservoirs. We further elaborate on how this knowledge can be applied to underground hydrogen storage. The primary focus lies on the mechanical response of rocks under cyclic hydrogen injection and production, fault reactivation, the impact of hydrogen on rock properties, and other associated risks and challenges. In addition, we discuss wellbore integrity from the perspective of underground hydrogen storage. The paper provides insights into the history of energy storage, laboratory scale experiments, and analytical and simulation studies at the field scale. We also emphasize the current knowledge gaps and the necessity to enhance our understanding of the geomechanical aspects of hydrogen storage. This involves developing predictive models coupled with laboratory scale and field-scale testing, along with benchmarking methodologies.

The seminal work of Gurson (J Eng Mater Technol 99:2–5, 1977) on a simplified pore structure, a single spherical pore, first provided a theoretical relationship between the yield stress and the porosity. This contribution extends the approach to determine the macroscopic yield of a porous material by taking explicitly into account its internal structure. As the yielding of a porous material is controlled by the geometry of its internal structure, we postulate that it is nearly independent of the constitutive plastic behaviour of the material. Here, we show that the influence of that internal structure on the yield could be retrieved from a finite element computation with just an elastoplastic ideal (J2) material equivalent of the skeleton’s. With some basic knowledge about the skeleton’s mechanical properties, this process allows the determination of the yield stress without requiring the experimental compression of the material. We showcase the predictive power of the method against experimental testing, initially for a unit cell following Gurson, i.e., unique cylindrical void in a 3D printed cylinder sample. Eventually, the applicability of the method is demonstrated on a complex 3D printed rock microstructure, reconstructed from a sandpack’s CT-scan.

The foundation of homogenisation methods rests on the postulate of Hill–Mandel, describing energy consistency throughout the transition of scales. The consideration of this principle is therefore crucial in the discipline of Digital Rock Physics which focuses on the upscaling of rock properties. For this reason, numerous studies have developed numerical schemes for porous media to enforce the Hill–Mandel condition to be respected. The most common method is to impose specific boundary conditions, such as periodic ones. However, these boundary conditions influence both the effective property and the size of the REV. The recent study of Thovert and Mourzenko (2020) has shown that most boundary conditions still result in the same intrinsic effective physical property if the averaging is applied outside the range of the boundary layer. From this discovery, it becomes logical to question the status of Hill–Mandel postulate in porous media when homogenising away from the boundary. In this contribution, we simulate Stokes flow through random packings of spheres and a range of rock microstructures. For each, we plot the evolution of the ratio micro- vs macro-scale of the energy of the fluid transport outside the boundary layer, for a growing subsample size of porous media. Here, we prove that we naturally find energy consistency across scales when reaching the size of the Representative Elementary Volume (REV), which is a known condition for rigorous upscaling. Furthermore, we show that this index for the energy consistency is a more accurate indicator of REV convergence since the mean value is already known to be unitary.

Modelling anisotropic flow in rocks requires their full permeability tensor. While theories derived from the upscaling of Stokes flow to Darcy's law may justify the tensor symmetry, homogenisations from micro-scale rock samples often return a non-zero level of asymmetry. Since most studies dismiss these controversial observations as numerical errors, this contribution looks more closely at the physical possibility of such behaviour. Asymmetry of the permeability tensor, which induces a rotatory flow, is manifested at the micro-scale by tortuous streamlines. Conversely, when considering a larger scale – above the Representative Elementary Volume for permeability – these tortuous paths do not statistically affect the flow direction any longer. At this point, the homogenisation of Stokes flow to Darcy's law reaches its domain of validity. We show that the asymmetry in the permeability tensor vanishes for this scale separation, regardless of the choice of boundary conditions unlike previously thought, if the boundary layer effect is disregarded.

Flow simulations on porous media, reconstructed from Micro-Computerised Tomography (μCT) scans, is becoming a common tool to compute the permeability of rocks. Still, some conditions need to be met to obtain accurate results. Only if the sample size is equal or larger than the Representative Elementary Volume will the computed effective permeability be representative of the rock at a continuum scale. Moreover, the numerical discretisation of the digital rock needs to be fine enough to reach numerical convergence. In the particular case of using Finite Elements (FE) and cartesian meshes, studies have shown that the meshes should be at least two times finer than the original image resolution in order to reach the simulation's mesh convergence. These two conditions and the increased resolution of μCT-scans to observe finer details of the microstructure, can lead to extremely computationally expensive numerical simulations. In order to reduce this cost, we couple a FE numerical model for Stokes flow in porous media with an unfitted boundary method for cartesian meshes, which allows to improve results precision for coarse meshes. Indeed, this method enables to obtain a definition of the pore–grain interface as precise as for a conformal mesh, without a computationally expensive and complex mesh generation for μCT-scans of rocks. From the benchmark of three different rock samples, we observe a clear improvement of the mesh convergence for the permeability value using the unfitted boundary method on cartesian meshes. An accurate permeability value is obtained for a mesh coarser than the initial image resolution. The method is then applied to a large sample of a high-resolution μCT-scan to showcase its advantage.

Sedimentary structures have unique geometries and anisotropic hydraulic conductivity, both of which control groundwater flow. Traditional finite-difference simulators (e.g., MODFLOW) have not been able to correctly represent irregular, dipping and anisotropic structures due their use of a simplified conductivity tensor, causing many modelers to turn toward finite-element codes with their sophisticated meshing capabilities. However, the release of MODFLOW 6 with its flexible discretization and multipoint flux approximation scheme prompts us to revisit its capability to compute flow through complex sedimentary structures. Through the use of a novel benchmark and case study, we show that when versions previous to MODFLOW 6 are applied to dipping structures, modeled fluxes and hence flow through the system, can be significantly over or underestimated. For example, effective conductivity for a 30° dipping layer with a 100:1 conductivity ratio is reduced to only 2% of its inputted value. We show that MODFLOW 6, with its XT3D capability and flexible discretization options is far superior to its predecessors, allowing flow through complex sedimentary structures to be simulated more accurately. However, on vertically offset grids, which have been available in all versions of MODFLOW and are often used in practice, loss of accuracy is still a concern when the vertical offset is large, that is, the dip of the sedimentary layer is steep, particularly if the layer is much more conductive than the surrounding material. The hypothesis that vertically offset grids lack sufficient hydraulic connectivity between adjacent model layers to accurately simulate the steeply dipping, highly heterogeneous case is a topic for further investigation.

Digital Rock Physics has reached a level of maturity on the characterisation of primary properties that depend on the microstructure - such as porosity, permeability or elastic moduli - by numerically solving field equations on μCT scan images of rock. After small deformations or at depth though, most rocks eventually reach their limit of elasticity and the complementary plastic properties are needed to describe the full mechanical behaviour. Currently, determination of a rock's yield surface from its microstructure is often restricted to semi-analytical criteria derived by limit analysis or numerical simulations performed on idealised geometries. Such simplification lacks representativeness, particularly for processes that affect directly the pore-grain interface such as the cementation phenomenon, happening during diagenesis. Eventually, only direct numerical simulation of elasto-plasticity performed on digitalised microstructures can be used to assess the strength of different cemented materials and its evolution with the alteration of the microstructure. In this study, we provide a comprehensive parametric study on the impact of cementation on rock strength for real microstructures of cemented granular materials. Compared to most previous studies, the whole yield surface is determined numerically (using Finite Element Method) in order to assess the influence of cementation for different stress-paths. The previously known tendency of rock to strengthen with increasing cementation volume is verified. New results on the influence of cement property namely Young's modulus, friction and cohesion on the rock's yield surface are explored. The envelopes obtained are compared to the ones obtained by experimental data and existing models. The framework presented in this study showcases the wider possibility of determining any rock's or porous material's yield surface from its microstructure.

Faulted and fractured systems form a critical component of fluid flow, especially within low-permeable reservoirs. Therefore, developing suitable methodologies for acquiring structural data and simulating flow through fractured media is vital to improve efficiency and reduce uncertainties in modelling the subsurface. Outcrop analogues provide excellent areas for the analysis and characterization of fractures within the reservoir rocks where subsurface data are limited. Variation in fracture arrangement, distribution and connectivity can be obtained from 2D fractured cliff sections and pavements. These sections can then be used for efficient discretization and homogenization techniques to obtain reliable predictions on permeability distributions in the geothermal reservoirs. Fracture network anisotropy in the Malm reservoir unit is assessed using detailed structural analysis and numerical homogenization of outcrop analogues from an open pit quarry within the Franconian Basin, Germany. Several events are recorded in the fracture networks from the Late Jurassic the Alpine Orogeny and are observed to be influenced by the Kulmbach Fault nearby with a reverse throw of 800Â m. The fractured outcrops are digitized for fluid flow simulations and homogenization to determine the permeability tensors of the networks. The tensors show differences in fluid transport direction where fracture permeability is controlled by orientation compared to a constant value. As a result, it is observed that the orientation of the tensor is influenced by the Kulmbach Fault, and therefore faults within the reservoirs at depth should be considered as important controls on the fracture flow of the geothermal system.

The design of any manufactured material requires the knowledge of its limit of elasticity, called yield strength. Whilst laboratory experiments are currently necessary to do so, this study is part of initiatives which aim at deriving the strength value from simple and fast numerical simulations. The seminal work of Gurson (1977) on a simplified pore structure, a single spherical pore, provided the first theoretical relationship between a material's strength and its porosity, showing that the presence of pore space is responsible for lowering the yield strength. The complexity of new structures requires however to take explicitly into account the internal geometry, usually using direct numerical simulations. This can be particularly challenging since the yield strength of a structure is actually reached after some of its parts have already entered the plastic regime. Therefore, the mere computation of the structure's yield strength currently necessitates the modelling of the exact plastic behaviour of the skeleton's material. In this contribution, we propose to simplify the numerical modelling necessary to predict the yield strength of a porous material, by postulating that yielding is mostly controlled by the geometry of the internal structure. We show that the influence of that internal geometry on the yield can effectively be retrieved from a finite element computation implementing a simple elasto-plastic model to represent the solid phase of the porous material. We showcase the predictive power of this new method against an experimental testing, initially benchmarked for 3D-printed samples with either a unique spherical void or a grid infill, before demonstrating its applicability on a complex 3D-printed real rock microstructure, reconstructed from segmented micro-Computerised Tomography scans.

Permeability is a critical parameter for geological resources characterization. Its evolution with respect to porosity is particularly interesting and many research initiatives focus on deriving such relationships, to understand some hydraulic impacts of microstructure alteration. Permeability evolution from chemical reactions, for instance, can become complex as flow channels may open during rock dissolution. In this contribution, we show that this phenomenon can lead to irregular porosity-permeability curves and permeability hysteresis after reprecipitation. Current approaches describing permeability as a simple function of porosity can therefore not capture this behavior, and we advocate instead the use of dynamic modeling for such scenarios. We demonstrate our approach by modeling a dissolution/precipitation cycle for a unit cell pore channel and quantify the process at larger scale on three different rock samples, whose microstructures are reconstructed from segmented micro-Computerised Tomography scans.

Fluid injection or production in petroleum reservoirs affects the reservoir stresses such that it can even sometime reactivate dormant faults in the vicinity. In the particular case of deep carbonate reservoirs, faults can also be chemically active and chemical dissolution of the fault core can transform an otherwise impermeable barrier to a flow channel. Due to the scale separation between the fault and the reservoir, the implementation of highly non-linear multiphysics processes for the fault, needed for such phenomenon, is not compatible with simpler poromechanics controlling the reservoir behaviour. This contribution presents a three-scale finite element framework using the REDBACK simulator to account for those multiphysics couplings in faults during fluid production. This approach links the reservoir (km) scale – implementing poromechanics both for the fault interface and its surrounding reservoir – with the fault at the meso-scale (m) – implementing a THMC reactivation model – and the micro-scale (μm) – implementing a hydro-chemical model on meshed μCT-scan images. This model can explain the permeability increase during fault reactivation and successfully replicate fault activation, slip evolution and deactivation features, predicted by common fault reactivation models, yet with continuous transitions between the sequences. The multiscale coupling allows to resolve the heterogeneous propagation of the fault slip which can be uncorrelated from the initial highest slip tendency location. We also demonstrate the advantage of dynamically upscaled laws compared to empirical ones as we show that a hydraulically imperceptible alteration of the microstructure's geometry can lead to different durations of the reactivation event at the macro-scale.

Flow simulators have become increasingly popular to compute permeability on digital porous rocks reconstructed from Computerised Tomography (CT) scan data as part of Digital Rock Physics workflows. Various schemes are being used that focus mainly on numerical efficiency but do not necessarily provide the flexibility to account for more complex couplings between other physical models than the flow itself, like irreversible mechanical deformation for instance. To overcome this limitation, we present a finite element implementation of an Eulerian flow solver that is coupled to a Lagrangian solid mechanics solver, both in a sequential or tight manner. The flow simulator is successfully benchmarked for permeability computation on a digital rock that has been previously used for other code validations. The need for a tightly coupled hydro-mechanical simulator is then highlighted by comparing the loose and tight coupling approaches, revealing that the tight approach is more accurate and more efficient computationally. Finally, the use of the simulator in real cases is illustrated by running a hydro-mechanical simulation on a digital rock sample, where plastic deformation and fluid flow in the porous space are resolved simultaneously.