03_sota.tex

\section{State-of-the-Art} 

\subsection{THM simulations and open source development}

Process-oriented numerical simulation programmes are necessary for predicting possible environmental impacts as well as for the macroeconomic and safety design of geosystems for underground use with, if necessary, different or even multiple management. These programmes must be able to represent the running processes and their interactions. Already in the mid-eighties of the last century, specific models were developed in the USA and partly implemented in scientific simulation platforms in order to describe THM processes taking place in the geological subsoil which are connected with the thermal use of the subsoil as energy source (geothermics) or energy storage. However, these investigations primarily had a basic character. In addition, the numerical calculation tools are often oriented towards the description of special processes and only partially consider couplings of different physical processes. Geotechnical applications can also be simulated with a number of established commercial program systems. For hydraulic processes such as multiphase flow in porous media are simulators from the oil and gas industry available (e.g. ECLIPSE, STARS), for the description of mechanical processes as well (FLAC3D). All mentioned codes can only cover a part of the necessary process spectrum. Therefore, simulation programs are required which can represent thermal, hydraulic, mechanical, and chemical (THMC) processes coupled, such as TOUGH \cite{Pruess2004738}, HYTEC \cite{vanderLee2002599}, DuMuX \cite{Flemisch20111102} or OpenGeoSys \cite{Kolditz2012613}. In order to be able to represent the foreseeable impact area of underground use in realistic simulation areas, efforts have been made to parallelise these codes (e.g. OpenGeoSys \cite{Wang20152269}, TOUGH \cite{Wu2002243}). In particular, the simulation of systems subject to discontinuity requires high performance computing. A major limitation of commercially available numerical simulation programs is that their source codes are not accessible and therefore not transparent and that a further development of such programs is therefore only possible by the commercial developer. In the research project applied for here, the platforms OpenGeoSys (UFZ (coordinating), BGR, CAU, IfG, TUBAF), mD-LEM (CAU) and pythonSPH (Uni Stuttgart) developed by some of the applicants as open source software can be used, so that the limitations mentioned do not exist. The description of discontinuities with different approaches described in the following as well as their processing in the sense of high-performance computing (HPC) requires targeted program extensions.
\index{HPC High-Performance-Computing}

\subsection{Continuum models (XFEM and variational phase field)}

In recent years, extended \cite{Belytschko1999601} or also known as generalized \cite{Strouboulis200043}  finite element methods (XFEM/GFEM) and phase-field methods \cite{Bourdin2000797} for the description of existing and developing discontinuities and singularities within continuum mechanical approaches have established themselves ahead of all others. 
Both methods differ fundamentally and have their own strengths and weaknesses. 
XFEM locally extends the approach and test function space by formulations that can map the discontinuous course of the solution and introduces corresponding additional local degrees of freedom. 
Usually, this approach is combined with so-called level set methods, which help to localize the discontinuity and thus ultimately determine in which elements the solution space has to be extended. 
This approach allows the approximation of discontinuous solutions on comparatively coarse grids, but requires programmatic infra-structures for the treatment of flexible additional degrees of freedom, level sets and other aspects, which require a considerable implementation effort, especially in branched crack systems. 
In contrast, the variational phase field method was originally proposed as a generalized Griffith criterion by~\cite{Francfort1998} and numerically implemented using a phase-field variable by~\cite{Bourdin2000}.
In the variational phase-field model, cracks are represented by a smoothly varying function (phase-field variable) that transitions from intact material (phase-field variable = 1.0) to fully broken state (phase-field variable = 0.0) using a regularization parameter with the dimension of a length and the energy consumed by the cracks is computed from this diffused representation. 
One of the strengths of this approach is to account for arbitrary numbers of pre-existing or propagating cracks in terms of energy minimization, without any a priori assumption on their geometry or restriction on the growth to specific grid directions.

XFEM \cite{Belytschko1999601} was originally developed for crack propagation problems and was also applied in the geotechnical context, e.g. for multiphase flows \cite{Chessa200310},\cite{Mohammadnejad2013327} and heat transport \cite{Khoei2012701},\cite{Shao2014155}. Current developments of generalized and extended finite element methods in the context of hydraulic stimulation are mainly concerned with the efficient coupling of solid-state and flow-mechanical problems \cite{Yazid20094269},\cite{Watanabe20121010},\cite{Meschke2015438}.

The variational phase-field model of fracture has witnessed wide ranging applicability in from dynamic fracture~\cite{Bourdin2011},\cite{Borden2012},\cite{Li2016}, to ductile fracture~\cite{Ambati2015},\cite{Miehe2015486}, \cite{Alessi2017}, to thermal and drying fracture~\cite{Maurini2013},\cite{Bourdin2014},\cite{Miehe2015_thermo}. 
The first application of the variational phase-field model to hydraulically driven crack propagation has been proposed by~\cite{Bourdin2012} and followed by many others~\cite{Wheeler2014},\cite{Wilson2016},\cite{Heider201738},\cite{Santillan2017},
\cite{Chukwudozie2019957},
\cite{Li201942} 
or for land slide modeling \cite{Wei2020} 
with various formulation and numerical implementation. 
While the reported findings are promising thus far, the method still needs more establishments for practical field scale applications. 
The required efforts may include validation against laboratory/field experiments, approaches to recover explicit properties such as fluid leak-off from smeared crack representation, and more complex physics phenomena such as visco-elasticity.
\index{VPF Variational Phase Field method}

\subsection{Discontinuum models}

Discontinuum models directly map forces of interaction between predefined discrete elements. The latter may themselves be discretized and mapped by continuum mechanics. Decisive for the mapping of developing discontinuities, however, are the pre-defined interfaces subject to certain interface formulations. This type of modelling was applied in geotechnics, for example, to geothermal systems \cite{Zeeb2015264} and has also made a decisive contribution to the simulation of the pressure-driven generation of flow paths in polycrystalline salt rocks, which is bound to the discontinuum-mechanical microstructure of the salt rocks. Polycrystalline salt rocks represent a discontinuum of intergrown salt crystals on the micromechanical level. In contrast to porous media, there is no cross-linked pore space in salt rocks. Only by pressure-driven opening and cross-linking of pathways, i.e. generation of connectivity by opening channels along the grain boundaries of the salt crystals, cross-linked flow paths are created in salt rocks. Fluid pressure-driven percolation is direction-dependent and seeks the path of least resistance along the crystal grain boundaries in the polycrystalline salt rock under the effect of the existing stress field. This mechanism of directional percolation can be simulated in coupled HM models on a discontinuity mechanical basis.

The observations on numerical models of pathogenesis by source and shrinkage processes based on a microscale based analysis must be able to map significant structural changes and discontinuity developments in nonisothermal HM coupled processes, which manifest themselves in progressive fracture or self-healing processes under pressure, saturation and temperature influence. First basics of the modelling of fracture processes on the microscale were published at the end of the 1990s with reference to self-organising fracture processes based on Voronoi discretizations \cite{BolanderJr.1998569}. By combining the approaches of HM modelling in saturated media \cite{Asahina201413} and TM modelling \cite{Rizvietal2016}, the connection for the simulation of self-organizing fracture processes in geomaterials shall be established with consideration of complex TH$^2$M processes. Based on the elasticity theory, linear fracture models following Mode I and Mode II were developed for fragile, largely homogeneous material with few inclusions. For materials with high interference, models based on continuum fracture mechanics \cite{Talreja1991165} were developed which require specific information on material microstructure and fracture behaviour. Thermal conductivity in cemented geomaterials is determined by heat transfer between mineral particles, porosity, fluid and contact quality 
\cite{Woodside19611688},\cite{Bahrami20063691},\cite{Widenfeld200315}.
The Thermal Particle Dynamics method can be used to simulate the transient heat propagation in granular media and the associated thermal expansion \cite{Vargas20011052}. This was considered in a thermal DEM \cite{Vargas20023119},\cite{Vargas2007}, but the calculation effort is enormous and the grain or contact shape is greatly simplified \cite{Zhang2011172}. In contrast to the particle methods, the heat propagation in cemented materials can be determined numerically very effectively by classical FEM, but microlevel information disappears due to the underlying homogenization. This poses a problem for the initiation of discontinuities by thermal processes in THM coupling. 

The hybrid lattice models have been developed to tackle the shortages in continuum based models, such as the simplicity to define the heterogeneity/anisotropy as well as the fracture simulation and stress redistribution during the frack propagation (discontinuities) without the need to re-mesh the domain \cite{Bolanderetal1998, vanMieretal2002}. The lattice model is similar to the finite volume (FVM) or finite difference (FDM) methods, with the difference that the FVM or FDM explicitly discretize the continuum \cite{Rizvietal2018a, Rizvietal2018c}. The simplicity and accuracy of lattice models to simulate the fracking in cemented geomaterials, such as rock and concrete \cite{Liuetal2007, Pradoetal2003, Karihalooetal2003}, are well established. The lattice models in comparison to the continuum methods are time consuming and expensive. Therefore, their applicability and development in real engineering applications or commercial softwares are not well developed. However, with the increase of computational power during past years as well as the implementation of parallel computing or GPU computing methods, the application of lattice models in commercial softwares is imminent.
\index{LEM Lattice Element Method}

The mentioned DEM approaches, as classical discontinuum models, have the disadvantage that additional connections between the particles have to be implemented by beam elements, which contain the fracture-mechanical criteria. Lattice based models - LEM \cite{Chessa200310, Chung199615094} have been developed for modeling of fracture mechanical processes considering discontinuity and crack initiation as well as crack propagation. These include a networking of the existing heterogeneity ranges (Voronoi and Delaunay triangulation) and use simple linear fracture criteria on the microscale. The cross-linked two- and three-dimensional continuum regions are microscopically coupled by 1D elements in the center of gravity of the Voronoi cells. In the simplest case these elements are Hookesche springs with a normal stiffness \cite{Curtin1990535}. In three-dimensional space these simple springs already give a good approximation of the Mode I failure model \cite{Wong2015417}. With the use of Born spring models and an additional tangential degree of freedom, shear behaviour can already be modelled \cite{Jagota19933123}. By extending the spring, for example as a beam element \cite{Schlangen1992435}, displacements, rotations and moments can be transferred to the node in addition to the forces, whereby an additional bending contact can be taken into account \cite{Sahimi1993713}. For the spatial lattice network thus generated, the displacements at each point are determined by generating an equilibrium or by minimizing the energy \cite{Meakin1991226} or dynamic relaxation \cite{Cundall197947}. The LEM combines the advantages of simple implementation with the ability to control particle interaction in the model while simultaneously self-organizing initiation and progression organization of a discontinuity, \cite{Sattarietal2019b, Wuttke201785}.
\index{DEM Discrete Element Method}

Additionally, in contrast to discrete models, the lattice models can be implemented to represent a continuum, as the lattice elements do not necessarily define the particle to particle contact mechanics \cite{Rizvietal2019a}. The hybrid lattice model can represent a continuum or particle-to-particle contacts, depending on the objective of the simulation. In both cases, the domain is discretized into series of spring or beam elements, representing the bonds. The regularization of a lattice model grants the independency of the results from the mesh size and meshing technique \cite{Ostojastarzewski2002}. The lattice models were initially emerged in order to simulate the fracture initiation and propagation in cemented geomaterials. With the time, lattice models have been extended to simulate the wide variety of the thermal \cite{Shresthaetal2019, Rizvietal2019d, Rizvietal2016}, thermo-mechanical \cite{Sattarietal2017, Sattarietal2019b} and hydro-mechanical \cite{Grassl2009} problems in the engineering applications. In the recent years, the hybrid lattice models have been extended to determine the granular, cemented or swelling geomaterials response under the coupled thermo-hydro-mechanical (THM) processes.

\subsection{Smoothed Particle Hydrodynamics}
\index{SPH Smoothed Particle Hydrodynamics}

Smoothed Particle Hyrodynamics (SPH) methods are reticule numerical collocation methods for solving partial differential equations. SPH methods were formulated almost 40 years ago to solve astrophysical problems and have been further developed in recent years to solve a variety of problems and models in fluid and solid state mechanics \cite{Monaghan2011323}. The SPH method is particularly suitable for problems with free surfaces or material interfaces such as discontinuities and cracks: SPH methods are Updated or Total Lagrange methods, i.e. boundary conditions at discontinuities can be described numerically well. In recent years, great progress has been made in the efficiency of SPH formulations, especially for questions with internal interfaces, such as for non-Darcy flows in porous media or in the multiphase fluidics of immiscible fluids in porous media \cite{Morris2000333},\cite{Tartakovsky2005610}. The so-called "Whole Domain Formulation", i.e. a numerical procedure in which the surface conservation equations (mass and momentum), such as the Young-Laplace equation in multiphase fluidics, are "smeared" by means of the kernel function and integrated into the bulk conservation equations (Continuum Surface Force - CSF), can be interpreted here as a "phase field method" which "continuously smears" the physical properties of the discontinuities. Besides the consideration of the SPH-inherent kernel function in the CSF methods and the absence of the need to artificially adduce discontinuities, the net-free SPH methods above all show great efficiency advantages when complex and small-scale (pore) geometries are to be precisely mapped \cite{Sivanesapillai2016212}. In addition to small-scale direct numerical simulations on the pore space scale, SPH methods for coupled HM problems in geomechanics have already been developed \cite{Bui2007339},\cite{Bui20141321}. The two HM-coupled biotubes poroelastic equation sets for the porous solid phase and the viscous pore fluid were formulated in these works with two disjunctive particle sets which can lead to difficulties in impulse interaction modelling. A further development of the SPH method for HM processes, also taking into account propagating discontinuities such as cracks and crack networks, is therefore imperative to establish the SPH method as an efficient and reliable tool for geoscientific problems.

All the approaches described above have proved to be suitable in principle for the physical analysis of the growth of discontinuities. However, in connection with the extension of the methods to coupled THM processes, there is still a fundamental need for development in many areas. This is to be supported by an improved process understanding to be worked out, by building on it some of the numerical methods used here are to be further developed purposefully beyond the state of the art. Applications that go beyond the simulation of laboratory experiments and use the methods for solving practically motivated problems of large-scale geosystems have so far hardly been found in the literature or have not even been developed for certain essential process couplings. There is an urgent need for systematic investigations into the questions of how these methods can be translated into practical applications, what computing resources are required, and in which cases certain methods appear more suitable than others. The aim of this project is to develop such an overall view and a systematic comparison of the methods at defined benchmarks as well as their embedding in proven software, partly with the inclusion of methods of high performance computing.