In reservoir simulation, it is important to understand the mechanical behaviour of fractured rock and the effect of shear and tensile displacements of fractures on their aperture. Tensile opening directly enhances the fracture aperture, whereas shear of a preexisting rough-walled fracture creates aperture changes dependent on the local stress state. Since fracture dilatation increases reservoir permeability, both processes must be included in a realistic and consistent manner into the mechanical reservoir simulation model.
Here, we use the extended finite volume method (XFVM) to conduct flow and geomechanics simulations. In XFVM, fractures are embedded in a poroelastic matrix and are modelled with discontinuous basis functions. On each fracture segment the tractions and compressive forces are calculated, and one extra degree of freedom is added for the shear and tensile displacement, respectively. To embed multiple intersecting fractures into non-conforming Cartesian grids of a given resolution, a merging technique is utilized which locally adapts fracture geometry and topology. Merging critical segments by retaining only the weakest links allows us to embed any kind of intersecting and branching fractures into the grid.
In this particular XFVM implementation, we assume that linear elasticity and steady state fluid pressure adequately constrain the effective stress. This allows us to integrate the stress equation to obtain force balances and then solve for the grid deformation, and the fracture shear and tensile displacements. In this paper, shear dilatation is not calculated a posteriori, but it enters the equations such that aperture changes directly affect the stress state. This is accomplished by adding shear dilatation to the displacement gradients and, therefore ascertains a consistent representation in the stress-strain relations and force balances. We illustrate and discuss the influence of this newly added term in simple test cases and in a realistic layer-restricted two-dimensional fracture network subjected to plausible in situ stress and pore pressure conditions.

Shallow geothermal reservoirs are excellent candidates for low-enthalpy energy storage, and can serve as heat batteries providing constant discharge of base heat, as well as rapid discharge of heat in periods of high demand. Recharging can be done by pumping down hot water, heated using excess heat from e.g., waste incinerators. In addition to having a very low carbon footprint, such systems also require limited surface infrastructure, and can easily be placed near or under the end-user, such as residential buildings.

The geological setting is typically complex, with horizons, faults, and intertwined patterns of natural fractures, and the nearwell region is often hydraulically fractured to enhance well communication. Therefore, in order to fully utilize the potential of shallow geothermal heat storage, numerical simulations are imperative. In this work, we show how to practically model such systems, including generation of computational grids with a large number of wells and fractures, numerical discretizations with discrete fractures (DFN), and complex storage strategies with multiple wells working together under common group targets. We also discuss how adjoint-based methods can be used to tune model parameters (e.g., well injectivities, rock properties, and hydraulic fracture conductivities) so that the model fits observed data, and to find well controls (e.g., rates and temperatures) that optimize storage operations.

The methodology is demonstrated using a set of real geothermal heat storage projects currently under development, and we highlight important challenges and our suggested solutions related to each of them.

Accurate simulation of hydraulic fracturing is required to understand the response of an enhanced geothermal system or a carbon storage reservoir to fluid injection. Current state-of-the-art approaches have involved the use of analytical solutions, termed asymptotes, to guide the evolution of the hydraulic fracture.
However, such asymptote-based approaches are difficult to derive in complicated geometries, with poro-elastic, thermo-mechanical or inertial forces, and with material inhomogeneities or inelastic deformation. Alternatively, energy-based methods, like the J-integral, have had a large success in simulating dry fractures under such general settings, which are more relevant to geological reservoirs.
However, in hydraulic fracturing, the J-integral has been reported to accumulate errors when simulating viscous-dominated propagation. Consequently, this article assesses a recent modification of the J-integral, the JHFM−integral, to accurately extract energy release rates throughout the toughness-viscous propagation regime. Such a JHFM does not rely on the use of asymptotes to predict the energy-release rate, which renders it extensible to the more general settings mentioned above, in principle. In order to implicitly evolve the fracture geometry, the JHFM-integral is used in combination with a linear extrapolation scheme of the stress intensity factor to predict and iterate on the location of the fracture front. Simulations are performed for the classical Khristianovich-Geertsma-de Klerk problem for comparison against analytical solutions. Furthermore, the JHFM-integral is compared against the use of the multi-scale tip asymptote to extract the propagation velocity and converge on the fracture front. A finite element discretisation using linear elements is shown to excellently predict the
location of the fracture front. Furthermore, less than 3 iterations are required for the extrapolation scheme to converge on the correct fracture length. This is achieved by setting a variable convergence criteria that depends on the local dimensionless toughness, which reduces unnecessary iterations in the viscous-dominated regime. Afterwards, the JHFM-integral is assessed in a mixed propagation regime, and found to preserve its excellent accuracy and convergence speed. Overall, the JHFM-integral combined with the linear extrapolation technique, appears to offer improvements in accuracy and speed compared to the multi-scale tip asymptote. Such results show that the JHFM-integral holds great
promise for simulating hydraulic fracturing in realistic, geological settings.