Hydro-mechanical fault reactivation modeling based on elasto-plasticity with embedded weakness planes

Lu Urpi, Bstin Grupner, Wenqing Wng, Thoms Ngel,d, Antonio P. Rinldi

a Swiss Seismological Service, Swiss Federal Institute of Technology, ETH Zürich, Zürich, Switzerland

b Swiss Federal Nuclear Safety Inspectorate ENSI, Brugg, Switzerland

c Department of Environmental Informatics, Helmholtz Centre for Environmental Research—UFZ, Leipzig, Germany

d Technische Universit?t Bergakademie Freiberg, Freiberg, Germany

ABSTRACT In this paper,an elasto-plastic constitutive model is employed to capture the shear failure that may occur in a rock mass presenting mechanical discontinuities, such as faults, fractures, bedding planes or other planar weak structures. The failure may occur in two modes: a sliding failure on the weak plane or an intrinsic failure of the rock mass. The rock matrix is expected to behave elastically or fail in a brittle manner, being represented by a non-associated Mohr-Coulomb behavior, while the sliding failure is represented by the evaluation of the Coulomb criterion on an explicitly defined plane. Failure may furthermore affect the hydraulic properties of the rock mass: the shearing of the weakness plane may create a transmissive fluid pathway.Verification of the mechanical submodel is conducted by comparison with an analytical solution,while the coupled hydro-mechanical behavior is validated with field data and will be applied within a model and code validation initiative. The work presented here aims at documenting the progress in code development, while accurate match of the field data with the numerical results is current work in progress.

Keywords:Fault reactivation Plane of weakness Finite element Argillaceous material Clay Permeability

1. Introduction

Rock masses may exhibit anisotropy in their shear strength due to the presence of internal discontinuities. External or internal loads(due to pressure,temperature changes,or mass removal)may alter the effective stress state acting on these internal discontinuities. Consequently, plastic failure may occur even under compressive load. In this work, the plane of weakness approach formalized by Jaeger (1960) (derived from theK0model for static earth-pressure) is implemented and validated in the open-source finite element (FE) simulator OpenGeoSys (Kolditz et al., 2012;Bilke et al.,2019).The plane of weakness approach postulates that the rock mass is characterized by an elastic behavior bounded by two possible failure criteria:one related to the rock matrix,and the other to a pre-determined oriented internal structure. While the elastic behavior is completely characterized by the rock matrix characteristics,failure may occur either in the rock matrix or along the internal structure. The method will then be applied to investigating the shearing reactivation of a fracture due to fluid injection,and the permeability change associated with plastic rupture,and to evaluating possible dilatancy effects.In general,the method allows the investigation of the mutual coupling between mechanical and hydraulic processes in the presence of discontinuities such as fractures, faults and bedding planes.

The pore pressure and the stress field are both strongly affected by these internal oriented structures and their interaction with the domain in which they are embedded (Ingebritsen et al., 2006;Hardebol et al.,2015);therefore,it is important to explicitly include their distinguishing features (i.e. flow properties, shear rupture direction, and strength anisotropy) in the solution of poromechanical problem.

The plane of weakness model is in fact suitable when one or several discontinuities with well-defined orientations are present,as noted by Brady and Brown (2006). Discontinuities may be present at different scales in complex materials, being the expression of small structures such as bedding planes or fractures,or they may cut across different formations,manifesting themselves as faults or interfaces between different formations:it is then advantageous to investigate the behavior of such discontinuities by allowing a single mesh element to be subject to two distinct failure criteria, where one criterion is due to the rock matrix(possibly with consideration of the intrinsic anisotropy or other characteristics of the material)and the other is due to the structure or geometry and defined independently from the rock matrix properties. Clearly, the distinction between large- and small-scale is somehow arbitrary,depending on the size of the computational elements themselves,as it can be the case for the bedding planes, which can be individually defined if the element size is small enough(10-3—10-2m),but in that case the computational requirements for a large-scale model (with a size on the order of 103or 102m) are extremely difficult to be satisfied. This approach is therefore useful to investigate idealized situations and speculative scenarios, as well as to improve our understanding of complex experiments, such as the Fault Slip (FS) experiment in the Mont Terri underground rock laboratory (Switzerland) (Guglielmi et al., 2020a). The FS experiment has been conducted to investigate the conditions and processes associated with the reactivation of faults in argillaceous rocks.In particular,it aims at investigating(i)the conditions for slip activation and stability of faults in clay rock,and(ii)the evolution of the coupling between fault slip,pore pressure,and fluid migration with the help of a newly developed probe capable of measuring pressure and displacement at high sampling rates(Guglielmi et al.,2014).

The FS experiment is one of the experiments performed or soon to be performed at an intermediate scale between the laboratory scale and the in situ scale of subsurface technological activities,such as the EGS-COLLAB project(Kneafsey et al.,2018),the fatigue fracture experiment performed at the ?sp? Hard Rock Laboratory in Sweden(Zimmermann et al.,2019),the series of Swiss experiments performed in Grimsel, Switzerland (Amann et al., 2018), and the STIMTEC experiment in Freiberg,Germany(Dresen et al.,2019).All of these experiments aim at understanding the plastic failure along discontinuities (faults or fractures) under in situ conditions,providing viable strategies to upscale results from laboratory-scale experiments to the field scale. These experiments shall provide relevant information on the physical processes that need to be included in forward numerical modeling tools, and to correctly understand and predict the rock mass behavior.

The DECOVALEX initiative (Birkholzer et al., 2019) for the validation of coupled models and the benchmark of numerical codes based on experimental data has been instrumental in the development of advanced thermo-hydro-mechano-chemical (THMC)models and their implementation into various codes,among them OpenGeoSys. Two stages of the FS experiment were subject of model development and validation in the context of“Task B”in the Decovalex-2019 phase (Rutqvist et al., 2020). In this validation attempt,two step rate tests from the FS experiment were modeled.Different codes were employed to provide understanding on the conceptual uncertainty associated with different modeling and code choices. The codes include finite difference methods, among them the TOUGH-FLAC combination with sequential coupling using interface elements (Park et al., 2020), as well as FE techniques(Nguyen et al., 2019). An overview of the different Task B approaches and their results is provided by Graupner et al. (2020).

In this paper,we present the validation of the plane of weakness model developed for the OGS version 5 open-source code and its application to the second stage of the FS experiment. The tools developed here are fully accessible and new rheological models can be implemented in an open and transparent way, allowing to address the fundamental processes involved in the activation of faults and fractures.This is of paramount importance in the context of various geo-energy activities.An approach based on the plane of weakness model has been employed in a range of scenarios: to assess the caprock integrity at CO2geologic sequestration sites(Rutqvist et al.,2016),to characterize the hydraulic stimulation of a geothermal well (Rinaldi and Rutqvist, 2019), to verify the improved stimulation and production from shale gas reservoirs(Lei et al.,2017),to evaluate the occurrence of induced seismicity during gas production(Zbinden et al.,2017),as well as for the assessment of a safe permanent geologic disposal of nuclear waste in argillaceous clay formations(Urpi et al.,2019).

The main concern associated with nuclear waste disposal is the potential creation of permeable flow paths through the initially impermeable or low-permeability host rock barrier. The tool presented here can address whether such a flow path can be created and under what hydro-mechanical conditions it may remain hydraulically conductive, possibly providing a transport path for radionuclides.

Our approach allows for the inclusion of large-scale discontinuities in a representative way, solving the fully coupled poromechanical problem within a single code.For mathematically wellposed problems, the fully coupled method provides unconditional and convergent numerical solutions (Noorishad et al., 1982;Rutqvist and Stephansson,2003).It has been shown that a subset of the iterative, sequentially coupled schemes can be as robust and efficient (in the numerical sense) as the fully coupled approach(Mikeli′c et al., 2014). However, a fully coupled approach avoids possible practical issues arising from the iteration between different codes(such as inconsistencies between input parameters,double mesh/grid realizations, transfer of partial results from different codes), while at the same time simplifying the access to large computational clusters, to successfully solve numerically expensive three-dimensional (3D) geomechanics problems, a welcome feature when the problem belongs to the rock-mechanics“data-limited model” (Starfield and Cundall,1988) class.

2. Numerical method

The FE method-based simulator OpenGeoSys (Kolditz et al.,2012; Bilke et al., 2019) is used to solve the flow of a fluid in a deforming porous medium. The system of partial differential equations with defined initial and boundary values is solved to determine the unknown variables: fluid pressurepand the solid displacement vector(so-calledformulation).Details for the Galerkin formulation applied here can be found in Korsawe et al.(2006).

Given the initial and boundary conditions (that may vary in time), an approximate solution of the problem is obtained by applying the Galerkin FE method for spatial discretization with linear interpolation functions for pressure and quadratic ones for displacement, as well as the backward Euler scheme for temporal discretization (Wang and Kolditz, 2007): the solution to the equation system is computed with the direct solver PARDISO (Schenk and G?rtner,2004)based on the OpenMP parallel implementation.

The deformation process is modeled by the following system of equations:

whereIis the identity tensor,represents the body forces.Eq.(1a)represents the momentum balance equation, taking into account the effective stresses σ′,given by the total stresses σ minus the fluid pressurep(scaled by the Biot’s coefficient) and the effect of body forces(if gravity is acting, this term equates toEq. (1b)defines the constitutive behavior in an incremental form,showing the elastic relation between the effective stress and strain tensor dεe, defined by the fourth-order stiffness tensorCe: in the following, this tensor is assumed to represent isotropic elasticity,therefore only two parameters will be needed to completely define its elements(in 3D).If the elastic guess of the strain tensor leads to a stress state violating the failure criterion, a plastic correction is applied to the strain tensor(see Appendix for its calculation in case of failure along the plane of weakness).

The mass and momentum balance for the fluid phase in conjunction with Darcy’s law are represented by the following equations:

While the equation for the fluid mass balance of the porous media reads

The numerical code OGS solves the discretized equation via the Galerkin approach(Wang and Kolditz,2007).The fluid flux and the effective stresses are derived respectively from the solution fieldspandvia the two constitutive equations,i.e.Hooke’s law(Eq.(1b))and Darcy’s law (Eq. (2b)).

2.1. Plastic behavior

The plastic behavior of the media is evaluated according to the Mohr-Coulomb criterion for the rock matrix and Coulomb’s law for the oriented plane of weakness (Jaeger et al., 2007). Internal cohesionc0and friction angle θ0characterize the maximum allowed shear stress τ for the rock matrix failure at a given normal effective stress:

While failure along the weakness plane is characterized by internal cohesioncpwand friction angle θpw, as well as by stresses calculated on the plane of weakness (shear τpwand normal):

Assuming a two-dimensional (2D) stress state determined by the principal stresses σ1and σ2, where the plane of weakness is at an angle β with the major principal stress,rupture takes place in the matrix if σ1reaches a certain critical value σc:

For a given stress state and a given plane of weakness orientation,rupture would take place in the rock matrix or along the plane of weakness depending on which of the σccalculated in Eqs.(5)and(6) is smaller, i.e. depending on whether the available strength is exceeded first in the matrix or along the weak plane. To test the validity of the implementation,we numerically reproduced a set of uniaxial (i.e. σ2= 0) drained compression tests performed on a cylindrical sample,comparing the critical stress value σc(which can be identified then as the uniaxial compressive strength(UCS)of the sample) as computed by Eqs. (5) and (6) and the UCS values obtained numerically. In Table 1, the parameters characterizing the synthetic rock sample are noted, while in Fig.1, the analytical and numerical UCS values are compared.

The tests are assumed to be conducted in fully drained conditions, and during compression of the sample, no pore pressure is built up. Therefore, stress conditions are homogeneous across the sample and total stresses equal effective stresses(initial pressure is set to be 0).

The numerical results show good agreement with respect to the theory, both in predicting the variation in UCS value and on the failure characteristics, discriminating as expected between plastic behavior along the plane of weakness or due to the rock matrix.

2.2. Representation of cubic-law flow in finite thickness elements

In the following, it is assumed that the flow takes place in a given fracture with hydraulic apertureb, with flow rate directly proportional to the fluid pressure gradient and gravitational force,inversely proportional to the dynamic viscosity μ and with low Reynolds number. Therefore, the fluid flow in the fracture is governed by Darcy’s law, with a cubic dependency between flow rate and hydraulic conducting aperture (Witherspoon et al.,1980). The fluid flow per fracture cross-sectionqcan be expressed as

Table 1Parameter of the synthetic rock sample to test UCS.

Fig.1. Resulting UCS from numerical simulations (dots) and comparison with analytical solution: the failure may occur along the weakness plane only for certain dip angles.

The hydraulic openingbis the results of the sum of three different terms:

Each term represents a different physical process and they are defined as follows:b0is the initial aperture;bpis the aperture induced by the fracture plastic shear dilation.It may be composed of dilation or sudden increase in permeability due to topological changes associated with shear in the intrinsic fault structure.beis the aperture induced by the fracture elastic opening. It will be considered to be affected by change in normal effective stresswith the normal stiffness of the fracturekn:

In the code formulation,the thickness of the element(indicated astel)must be given,to allow for the scaling of the permeability to correctly reproduce the cubic flow rule,therefore the permeability is internally computed as

This value is fed into Eq. (4).

3. Model set-up

The model mimics the presence of a borehole injecting directly into a low-permeable fault structure,embedded in an impermeable and intact rock mass. The borehole is assumed to be connected to an injection chamber having a radius of 0.5 m, where the fluid is accommodated before eventually entering the fault structure. The fault structure is composed of two structures,subsequently defined as minor fault and major fault, having the same strike orientation and dipping at 30°and 60°,respectively.The injection and the fault structure represent the injection experiment performed at Mont Terri (Guglielmi et al., 2020a) in an idealized manner, where the injection and the monitoring of the displacements via a threecomponent borehole sensor allowed to image clearly different processes related to elastic and plastic deformation of the fault material. The scenario described here can be seen in Fig. 2.

The mesh presented here is generated with the meshing tool Gmsh (Geuzaine and Remacle, 2009) and is composed of three material groups,two defining the two faults and one for the intact rock mass.Constitutive parameters are denoted in Table 2.Normal(shear) stiffness is evaluated considering the thickness of the fault and its Young’s (shear) modulus according to simplified relations derived by Singh(1973)and summarized in its practical implication in ITASCA (2017):

Fig. 2. Vertical section showing the set-up of the Step rate test 3. On the left, the monitoring probe partially inserted into the fault zone with dip 60°, and on the right,the injection probe hydraulically connected to a minor structure with dip 30°, expected to encounter the fault at some distance (as extrapolated from structures observed in borehole logs). Figure originally published by Guglielmi et al. (2020a).

whereErepresents the Young’s modulus(an equivalent expression for the shear stiffnesskswould have here the shear modulusGinplace ofE)and ν is the Poisson’s ratio.From the assumption that the shear and normal stiffnesses of the rock mass under analysis are equal,it can be derived that Young’s and shear moduli have the same value.Therefore,Poisson’s ratio must be-0.5 sinceE= 2G(1 + ν).

Table 2Constitutive parameters.

The fault is represented by hexahedra in a structured fashion,since automated computation of the thickness necessary to derive the scaled the correct permeability (see Eq. (10)) of the element is not yet included in the code presented here, and thus elements representing the fault must have the same thickness(with respect to the imposed flow direction along fault).Injection into the fault is defined as pressure imposed at its elements’ nodal points. An overview of the mesh is visible in Fig. 3.

In the field experiment, injected flow rate is controlled, to achieve a pre-determined step-increase in fluid pressure in the injection chamber.In the numerical model,pressure is imposed at grid points defining elements at the injection point, flow rate and pressure are monitored at a monitoring point in the major fault(located at 4.24 m away from the injection point).Displacement of the two discontinuities (red and orange planes in Fig. 3) represented with planes of weakness is monitored both at the injection point and monitoring location in the major fault.The displacements should indicate whether the fault is opening elastically due to the pressure increase, or if sliding takes place due to the shear stress overcoming the shear strength.

Fault opening and/or shearing is driven by its properties and by the in situ stress state.The imposed stress state is derived from the stress state obtained by Martin and Lanyon (2003) as well as the stress determined by analysis of the experimental data (Kakurina et al., 2019; Guglielmi et al., 2020b), and preliminary results from Task B of DECOVALEX2019 (Park et al., 2019).

Fig.3. Overview of the model,with marks at the injection point and monitoring point.Data points can be compared with field data.The x-direction is aligned with strike,the y-direction is aligned with dip (horizontal projection), and the z-direction is aligned with the vertical (positive up).

The assumed stress state is presented in Table 3.

4. Results

The step injection and the comparison between monitored and computed pressure values at the monitoring point, located in the major structure, as well as the measured and numerically computed injected flow rates are presented in Fig. 4. There is no pressure response in the monitoring point until the rupture occurs.Similarly, before the rupture flow rate is limited to the amount needed to pressurize the hydraulically connected section of the minor fault. After the rupture, a sudden pressure increase at the observation point is observed, resulting from the increased permeability due to shearing. The elastic opening of the fault provides additional fluid transmissivity.

Given the assumed normal stiffness of the fault (last row of Table 2),an increase of pressure of some MPa will produce an elastic opening of the fracture/fault on the order of tens of μm. The deformations(in terms of normal and shear displacements)computed at the injection point and monitoring point are presented in Fig.5.

Two types of behavior can be observed in the displacement.First, before 770 s, there is mostly normal opening in the minor fault,due to the increased pressure not propagating into the major fault(although after 600 s the minor fault starts to fail in shearing mode,as denoted by the rapid shear displacement increase).Due to stress and fault orientation, the minor structure behaves like a thrust fault, moving upwards and north/eastward. Second, after 770 s, the major fault is reactivated and the shear of this major structure affects the displacement measured at the injection point as well. Due to the major fault orientation and stress state, the resulting shear is downward and south/westward, therefore the computed shear displacement at the injection point is reduced.

Associated with the flow injection, the Step 3 experiment provides the displacements of the injection point and monitoring point,via the two SIMFIP probes.The comparison of displacements for the injection point is presented in Fig.6(displacements are the relative displacements of the top anchor with respect to the bottom anchor).

Regarding the injection point displacements, the observed behavior is quite different, however, field data are affected by the elastic response of the borehole that is not considered nor correctedhere. Nonetheless, the initial response up to 500 s shows a difference in response: orientation and magnitude of vertical and northern displacements are initially correct, but after the first injection step,the behavior is barely reproduced by the model. Final displacements are nonetheless in good agreement: this may indicate that the borehole elastic effects are dominant in the field data,while the permanent deformation induced by plastic shear deformation is captured by the numerical model.

Table 3Initial stress state. SV, SH, and Sh are respectively the vertical, maximum horizontal and minimum horizontal stress components.

Fig. 4. Flow rate and pressure evolution at the monitoring point. Sudden increase is related to slip in the major fault and subsequent inflow. Dashed line indicates field data, and solid dotted line shows results from numerical simulation.

Fig.5. Computed displacement at the injection point(solid line) and monitoring point (dotted line). Shear and normal displacements calculated considering the dip of the minor fault and the major fault (30° and 60°, respectively). Positive shear means upward displacement for the minor fault, and downward for the major fault.

Fig. 6. Comparison between computed (solid line) and measured (dashed line) displacements for the injection point (a) and for the monitoring point (b).

Regarding the deformations at the monitoring point, although the model shows earlier rupture and the magnitudes are smaller,the general orientation of the deformations is captured, including the recovery of the elastic deformation due to the pressure decrease after 1050 s.The final displacements show discrepancies which are due to both the shear displacement and the different pressure conditions: observed pressure in the field data shows a pressure of approximate 0.01 MPa at the monitoring point,while in the model,there is an established hydraulic connection between the monitoring point and the injection point(where pressure is kept at 1 MPa).

5. Discussion

The step injection and the FS experiment in the Mont Terri underground rock laboratory (Switzerland) were conducted to investigate the conditions and processes associated with the reactivation of faults in argillaceous rocks. Especially, they aim at investigating (i) the conditions for slip activation and stability of clay faults,and(ii)the evolution of the coupling between fault slip,pore pressure and fluid migration.

The numerical simulation presented here is part of the international validation project DECOVALEX and aims at understanding and reproducing the field data measured during part of the FS experiment.The results show that assuming a fault behavior for the two structures represents correctly some general features of the in situ experimental conditions. With respect to pressure evolution(see Fig.4), the model reproduces the sudden increase in pressure after injection pressure reached its peak at 5.3 MPa.A change in the injected flow rate is associated with the rupture: before rupture,flow rate shows only some spikes, needed to raise the pressure to the desired value, while after the rupture, an almost constant injection rate is achieved/needed to maintain pressure.

The numerical results reproduce this behavior, although the spikes have a smaller magnitude and the post-rupture flow rate is smaller by a factor 4: nonetheless, the qualitative agreement with the field data supports a basic validity of the underlying assumptions of the numerical model in capturing the most important physical processes occurring in the field, while the discrepancies can be attributed to the parameters values presented in Table 2.A fit could improve the results.Unfortunately,this was not possible due to time constraints, while achieving a meaningful inclusion of physical processes was deemed satisfactory enough for the validation project.

By comparing the pressure changes and the displacements at the monitoring locations in the major fault,the field data present an increase in pressure before the shear slip takes place, while the numerical results show contemporary increase of pressure with the occurrence of shear slip.The field data have been discussed already by Guglielmi et al.(2020a).They concluded that it is likely that the fault leakage occurs before reactivation and the numerical results presented here tend to favor this interpretation:rupture may occur only after fault leakage. In this study, the fault leakage causes are not yet completely reproduced, since of the two most likely mechanisms, off-fault plasticity and pre-rupture tensile opening/dilatant behavior, only the second is partially implemented.

Our simulations estimate the shear displacement(Fig.5)on the major fault at the monitoring point at 80 μm, overestimating the 10—50 μm shear displacement range computed by Guglielmi et al.(2020a). Final displacements at the monitoring point show differences in the order of tens of μm which arises from the different pressures in the numerical model (determined from the pressure 1 MPa imposed at the injection point)and the pressure measured in the field (around 0.01 MPa).

The complexity of the material response is visible in the comparison between the displacements at the injection point(Fig.6a).The numerical results reproduce only qualitatively the experimental data.Here a large influence is played by the assumption that a fracture/fault-like structure is connecting the injection well with the major fault. With respect to this, it is worthwhile to mention that triaxial experiments performed on Opalinus Clay samples(Gr?sle and Plischke, 2011) and successive independent interpretations from Parisio et al.(2015)and Xu et al.(2013)showed that, in case of plastic failure, the rupture can occur along a plane that is co-located with neither the bedding planes nor the expected direction from a Mohr-Coulomb criterion in the principal stress space. This behavior has been observed also for Tourmaline Clay samples and these experimental results have been investigated analytically by Nguyen and Le (2014). In these studies, the anisotropic elastic stiffness is the main drive for the failure along certain planes and can be correctly captured by the approach proposed by Pietruszczak and Mroz(2001).Regarding our results,they seems to indicate that the assumption of our model on the minor fault as a plane of weakness is not correct and that the minor fault is possibly the expression of failure along the rupture plane arising from anisotropic stiffness and different strength properties.If this is the case, the discrepancies between the numerical displacements at the injection point and the field data may be explained by material having lower cohesion and strain hardening behavior, inducing earlier,smaller rupture.The agreement of final displacement at the injection point between simulation and field may hint at the need of including elastic anisotropy to characterize the elastic displacement,while the residual plastic deformation is adequately captured by the plane of weakness orientation defined in this work. Future work is planned to discern if the minor fault structure effectively exists before the injection experiment (as it is currently assumed)or if it develops due to the fluid injection from successive plastic failures into the anisotropic rock mass in a direction determined by the combination of in situ stress state and material anisotropy(mostly determined by the bedding planes).

6. Conclusions

An extension to the numerical simulator OpenGeoSys was presented in this paper.With this model,we aimed to characterize the pressure evolution associated with the elastoplastic deformation of the fault material.The application of the approach is not limited to a fault in a clay formation, but has a more generic applicability. In addition to the described approach, other methods to represent evolving discontinuities are presented in OpenGeoSys,such as co—dimensional interface elements,phase-field methods,or non-local damage approaches (Yoshioka et al., 2019).

The fault is represented with elastoplastic elements, where the possible failure is considered with an oriented frictional failure criterion.Once the shear stress overcomes the shear strength of the element, plastic deformation takes place, with the strains prescribed to develop according to the flow rule in the direction determined by fault orientation and shear stress, as well as the dilation or compaction of sheared material. Theoretically, the number of faults that can be accounted for is limited only by the hardware constraints on the size of the mesh, given that with the direct solver PARDISO the entire discrete system must reside in the memory (i.e. it is not decomposed in different independent domains as in other solution strategies based iterative solvers).Such alternative solvers can be selected in OpenGeoSys, but their robustness has not been studied here.

The plane of weakness feature has been validated analytically and the code is now participating in the international validation project DECOVALEX,with the aim of correctly representing the field data briefly presented here. The field data show some peculiar characteristics, due to the unique challenges of the in situ conditions. Nonetheless, the code reproduces correctly the shear orientation with respect to the fault orientation and in situ stresses.The discrepancies between the numerical results arise from uncertainties in fault properties (normal and shear stiffnesses, cohesion,strain softening/hardening behavior)and on characterization of the minor fault structure, as well as from neglecting anisotropy in strength and stiffness of the clay matrix.

The approach presented here allows for the monolithic solution of 3D poromechanical problems in presence of large-scale faults and fractures. The use of a single open-source code provides flexibility in problem definition. It can be applied in many geo-energy and geo-engineering fields and it allows the implementation of different rheological models while including the hydraulic and plastic approaches presented here.

Declaration of Competing Interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

DECOVALEX is an international research project comprising participants from industry,government and academia,focusing on development of understanding, models and codes in complex coupled problems in sub-surface geological and engineering applications; DECOVALEX-2019 is the current phase of the project.The authors appreciate and thank the DECOVALEX-2019 funding organisations of Andra, BGR/UFZ, CNSC, US DOE, ENSI, JAEA, IRSN,KAERI,NWMO,RWM,SúRAO,SSM and Taipower for their financial and technical support of the work described in this paper. The statements made are, however, solely those of the authors and do not necessarily reflect those of the funding organisations.

Source code of the OpenGeoSys5 version used here can be obtained from https://github.com/lurpi/ogs5/tree/master.

Input files can be obtained from https://github.com/lurpi/TaskB_input_DECO19/tree/step3.

Figures have been plotted with the aid of Matplotlib (Hunter,2007).

The authors would like to thank Hua Shao,Jobst Maβmann and Gesa Ziefle from the Federal Institute for Geosciences and Natural Resources (BGR) and Wenjie Xu from the MOE Key Laboratory of Soft Soils and Geoenvironmental Engineering for fruitful discussions and troubleshooting with OpenGeoSys.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2020.06.001.