Numerical models of pressure-driven fluid percolation in rock salt: nucleation and propagation of flow pathways under variable stress conditions

Success of our ongoing energy transition largely depends on subsurface exploitation. The subsurface can act as a “battery” to store energy dense fluids such as hydrogen, or a “host” to sequester unwanted substances such as carbon dioxide or radioactive waste. On the other hand, these operations cause the subsurface pressure and/or temperature to change and induce various (or cyclical) loadings to the surrounding formations. Their operational safety crucially hinges upon the subsurface integrity. The most imminent risk is nucleation of cracks that can lead to loss of mechanical integrity. Unlike hydraulic fracturing in geoenergy applications where one deliberately initiates cracks at certain targets, we normally design a system to avoid fracturing. At the designing stage, we thus have no prior knowledge of crack initiation locations or propagation paths. And, the computational designing tools should be able to assess the fracturing risk without such prior knowledge. In this study, we compared three computational approaches that do not require prescribed crack geometries—the discrete element method, the lattice element method, and the variational phase-field approach—against percolation experiments on rock salt. The experimental results show different fracture propagation paths depending on the boundary loads. The fracture geometries were reasonably matched by all approaches despite some differences in path irregularities. While the variational phase-field approach predicts relatively regular fracture paths, the paths predicted by the discrete and the lattice element methods are more irregular. These irregularities may seem more comparable to intergrain failure in real rocks, but they are also necessary triggers for fracture initiation in the discrete and the lattice element methods. In contrast, the fracture initiation in the variational phase-field approach is a realization of the energy minimization in the system, and the grain level descriptions are absent in the current formulation. These findings highlight their predictive capabilities and gaps to be bridged between the grain and continuum scales for field-scale applications.


Introduction
It is crucial to assess interacting processes in subsurface applications (Bauer et al. 2017;Martens et al. 2019;Volchko et al. 2020). One needs to select a site based on a long-term subsurface integrity in geological storage of resources such as hydrogen, carbohydrates, or water. The time scale in focus spans beyond hundreds of thousands of years for the storage or sequestration of harmful substances such as chemotoxic or nuclear waste, carbon dioxide, etc.
We rely on geological barriers to isolate the stored or disposed material from the larger hydrogeological cycle during and after the anthropogenic intervention through a low permeability, retention, and retardation mechanisms (McCartney et al. 2016;Martens et al. 2019). The barrier rocks also need to have sufficient strength and ductility to accommodate loads induced by the construction and operation of the storage space. Ideally, the rocks can reverse damage through mechanisms such as sealing or healing.
Several phenomena may challenge the integrity of barrier rocks . The excavation of tunnels, shafts, or caverns may alter the in-situ stress state and exceed the elastic deformation limits, which leads to enhanced permeability and/ or porosity. Temperature changes can also cause inelastic deformation or mineralogical changes that impair the barrier rocks.
This study deals with the pressure-driven percolation of fluids through barrier rocks. Fluids stored under pressure can act as mechanical loads on the barrier rocks, open secondary pathways, and thus compromise the barrier integrity. We particularly look into rock salt, considered an ideal barrier material, because of its ductility, healing potential, high thermal conductivity, and tightness. We use rock salt for storing pressurized fluids in solution-mined caverns, chemotoxic and nuclear waste, archival material, etc.
Undamaged rock salt is almost impermeable as long as the fluid pressure remains below the so-called percolation threshold. This threshold is given by the sum of the minimum principal stress and the tensile strength of the grain boundaries of the polycrystalline material. Above this threshold, the fluid pressure can open up grain boundaries and generate secondary permeability through a connected grain boundary network. It is via this network that stored fluids can escape.
Here, we compare different numerical methods in their ability to reproduce experimental observations of pressuredriven fluid percolation in rock salt. The methods considered are the lattice element method, the discrete element method, and the variational phase-field method.
The paper is structured as follows. Section "Fluid percolation (hydraulic fracturing) experiments on rock salt" describes a classical set of fluid percolation experiments in rock salt. In Sect. "Modeling approaches", we introduce the three different modeling approaches. We apply each approach to the simulation of the experiments described in Sect. "Fluid percolation (hydraulic fracturing) experiments on rock salt". Section "The discrete element method" presents the results of these simulations which are discussed in Sect. "Discussion". The paper closes with concluding remarks in Sect. "Conclusions".

Fluid percolation (hydraulic fracturing) experiments on rock salt
To compare with numerical models, we used published experiments on pressure-driven fluid percolation (hydraulic fracturing) in rock salt (Kamlot 2009). The experiments were conducted in rock salt to study stress-dependent hydraulic fracture propagation. The rock salt studied stems from a mine near Bernburg, Germany, in the Leine formation. Its density is in the order of 2.15-2.2 g cm −3 , and its bulk and shear moduli are 16.7 GPa and 10 GPa, respectively. The shortterm uniaxial compression strength is on the order of 25 to 30 MPa while the long-term uniaxial compression strength is at about half of the short-term values. Table 1 lists the mechanical properties of rock salt used for our simulations.
Cubic samples with 100 mm edge length were prepared with a small borehole in the upper middle boundary (Fig. 1a). The samples were loaded with a true-triaxial apparatus and pressurized fluid was injected through a borehole drilled in the middle of the sample to induce hydraulic fracture (Fig. 1b). The depth and the diameter of the hole are 40 mm and 16 mm, respectively. The drilled hole was cased off to a depth of 30 mm leaving a 10 mm open section in the bottom for fluid entry to the sample.
Two different stress states were applied to the samples. The first case is subjected to the minimum loading from the top surface, mimicking a reverse faulting stress regime 1 (Fig. 2a), while the second case is in a normal faulting regime (Fig. 2b).
The flow rate was increased gradually for stability. The pressure peaked at 10.1 MPa in the first case and 5.8 MPa in the second case. They are both higher than the minimum stress (8 MPa and 4 MPa in the two cases, respectively) but lower than the intermediate stress.
The resulting crack patterns are visible on the sample surfaces (Fig. 2). The first case shows a horizontal crack and the second shows a vertical crack both of which developed on the plane orthogonal to the respective minimum loading direction.

Discrete element method
The discrete or distinct element method (DEM) extends the capabilities of continuum-mechanical approaches by introducing a new level of discretization, which allows it to describe independent deformable bodies that can interact via their contact points and surfaces. The behaviour of these contacts can be modeled using joint constitutive models. This approach is especially suitable for materials with a pronounced grain structure, such as granular materials like loose rocks or polycrystalline materials like salt rock. The constitutive laws employed in this study for both crystal (grain) and contact behaviour, developed by Minkley et al. (2001), Minkley (2004), Minkley and Muhlbauer (2007), and Minkley et al. (2012), are implemented as DLLs (Dynamic Link Libraries) for the program 3DEC of Itasca CG, Inc [12]. The salt grains are described using the elasto-visco-plastic Minkley model. Its rheological description contains four parts: the elastic response, primary and secondary creep, and plastic deformation. The primary creep is implemented using a Kelvin element, and the secondary creep is a modified Maxwell element with a generalized, nonlinear stress-strain rate relation. These two processes do not play a role during the short timescale considered here. The plastic response follows a modified Mohr-Coulomb model with a nonlinear yield function: where D is the uniaxial compression strength, max the maximum effective strength, the curvature parameter, and 3 the minimum principal stress. This yield function was specifically developed to model the plastic deformation of salt. The parameters in Table 2 were determined experimentally in this study to adequately reflect the dependence of the maximum stress a rock can carry on the confining stress.
The interfaces between the distinct elements serve two purposes. On one hand, they allow the explicit description of crack formation, i.e., the mechanical opening of joints between salt crystals. The shear strength of the grain contacts depends on a nonlinear way on the normal stress The friction coefficient is given by Minkley and Muhlbauer (2007) where i 0 is an upslide angle which takes surface roughness into account and is modified by an exponential function to model abrasive effects, Φ R is the residual friction angle, max describes the adhesive friction, and k 1 and k 2 are curvature parameters. As an additional property, a tensile strength can be assigned.
We used the elastic properties in Table 1 for grains, and the plastic and fracture properties are given in Table 2. For the liquid, we used a bulk modulus of 2000 MPa and a viscosity of 1 mPa s.
The second purpose is the description of the hydro-mechanical coupling. The interfaces carry a network of fluid knots (see Fig. 3), each of which contains a set of fluid-related data, such as a knot volume, aperture, and fluid pressure. The fluid flow is then described by Darcy's Law with permeabilities derived from the hydraulic aperture. This is a feature that is directly provided by 3DEC, and it enables a nonlinear pressure dependence of the aperture. If the fluid pressure surpasses the normal stress plus the tensile strength of the interface, the knot is marked as mechanically open, and a crack is formed or propagated. For further details, we refer to references Minkley and Muhlbauer (2007) Minkley et al. (2012), and the 3DEC documentation (Itasca Consulting Group Inc. 2016].

Lattice element method
In the hydro-mechanical lattice model, the domain is discretized into a series of Voronoi cells connected through mechanical lattice elements (mechanical structure) and conduit lattice elements (flow channels). The vectorizable random lattice (VRL) method is applied here to position the nuclei and generate the mesh. The irregularity factor known as the randomness factor ( R ), which varies between 0 and 1, is used in this method to control the regularity of the mesh (Moukarzel and Herrmann 1992). When the randomness factor is 0, the generated mesh is regular, and when it is equal to 1, it reaches the maximum irregularity. After generating a mesh, Voronoi tessellation is deployed and polygonal cells are generated. Then, the mechanical lattice elements are generated based on the Delaunay triangulation connecting the neighbouring Voronoi cells.
The developed hydro-mechanical lattice model is based on the assumption of the dual lattice network (Grassl 2009;Grassl et al. 2013). The mechanical lattice elements transfer the mechanical loads between two adjacent Voronoi nodes (nuclei). In the same manner, the fluid flow is transferred through conduit elements, where the fluid mass is stored within defined physical or artificial cavities (Lisjak et al. 2017). Each conduit node represents an artificial cavity connected through conduit elements, through which the fluid mass is transferred (Fig. 4).

Implementation of a mechanical lattice model
In this study, the mechanical lattice elements are represented by a series of Euler-Bernoulli beam (3DOF) (Fig. 5) elements. Following (Ostoja-Starzewski 2002;Karihaloo et al. 2003), the regularization of the lattice model is performed and a relationship between the continuum and element properties is derived. In this approach, the stored strain energy of the continuum, U ℝ , is considered equal to the stored strain energy in all Voronoi cells, U cell where ℝ are the applied stresses, ℝ are the resultant strains, and V ℝ is the volume of the domain. For a 2D Euler-Bernoulli beam element, the stored strain energy is given by where V is the volume of the Voronoi cell, i,j,k,m is the strain, C ijkm is the stiffness tensor, i,j the curvature strain, and D ij is the curvature stiffness. The stiffness tensor is defined as where N b is the number of lattice elements representing the each Voronoi Cell, n i,j,k,m represent unit vectors of a lattice element, R ′ is the first stiffness coefficient, and R ′′ is the second stiffness coefficient. In each loading step, the minimization of the potential energy of the domain is carried out and the displacements of the cell nodes are determined where, for a single element b, U b t is the stored total strain energy, U b a is the axial strain energy, U b s is the shear strain energy, U b m is the moment strain energy, L b is the length of the lattice element, f x is the axial force, f y is the shear force, The Euler-Bernoulli beam element representing the bond between two cells  is the shear modulus, E b is the Young's modulus, and I b is the moment of inertia. The fracture initiation and propagation in the domain is based on linear elastic fracture mechanics (Rizvi et al. 2019). In this study, only Mode I (tensile) and II (shear) failures of the elements are considered. The failure envelop is defined based on the Mohr-Coulomb tension cut-off model (Bolander and Saito 1998). To simulate the material response of quasi-brittle geomaterials, a bi-linear softening scheme is implemented (Ince et al. 2003). The adjusted E b values after reaching the failure state are calculated based on where p is the peak strain, f is the failure strain, b is current element strain, and f p is the peak load before the stiffness degradation.

Implementation of the coupled hydro-mechanical lattice model
The dual lattice model is adopted here for the simulation of the coupled hydro-mechanical processes in rock material. The conduit nodes are defined as artificial cavities (Fig. 4), which are interconnected through conduit elements. The flow between the defined artificial cavities follows the mass conservation in a discrete time step as follows (Lisjak et al. 2017): where m t+1 f is the fluid mass in the next time step, m t f is the fluid mass in the current time step, Δm f is the total change of fluid mass in an artificial cavity, Sr is the saturation degree of a cavity, V cav is the volume of the cavity, f is the fluid density, P f is the hydraulic pressure, K f is the fluid bulk modulus, f(Sr) is the saturation function which is equal to 0 and 1 in a dry and saturated conditions, respectively, Z i,j is the relative coordinate of the i, j conduit nodes, g is the gravitational acceleration, R h is the hydraulic resistance, Δt is the time step, and Δm f ,ij is the mass of the fluid transported between cavity i and j. R h is calculated from the cubic law where a h is the hydraulic aperture, f is the fluid viscosity, and L ′ b is the length of the conduit element (flow length). The pore pressure inside a cavity develops when the saturation degree is equal to 1. Otherwise, the pore pressure is equal to zero. The pore pressure in each time step is calculated based on the amount of fluid mass flowing into the cavity. In this study, the capillary flow and the capillary rise are neglected The developed pore pressures are then transferred into the mechanical lattice nodes. Through a weak coupling scheme, the pore pressure diffusion and the change of the hydraulic conductivity with the crack opening and closure are simulated. Both pressure-and flow rate-controlled boundary conditions can be considered in the coupled scheme.

Variational phase-field model
A variational approach to fracture has been proposed by Francfort and Marigo (1998) as a recast of Griffith's criterion for fracture propagation and its numerical approximation through a phase-field method by Bourdin et al. (2000). Since its inception, the approach, now often called "phasefield model", has revolutionized the field of computational fracture mechanics. We refer to Bourdin and Francfort (2019); Francfort (2021) for recent reviews. In the following, we briefly describe the "phase-field model" applied in this study.

Mathematical model
We follow the regularization of Francfort and Marigo's energy functional (Francfort and Marigo 1998) by Bourdin et al. (2000), where the total energy is defined as where d is the phase-field order parameter that represents a fully damaged state of the material ( d = 1 ) and the intact state ( d = 0 ) with a continuous function, and is the regularization parameter with a unit of length. The model with n = 1 is typically referred to as AT 1 (Pham et al. 2011) and n = 2 as AT 2 (Bourdin et al. 2014). AT 1 retains an elastic phase before failure, while the damage evolves immediately with the AT 2 model. We used AT 1 in this study. c n is a normalization parameter given by (13) where ⟨⋅⟩ denotes the Macaulay brackets defined as ⟨a⟩ ± = (�a� ± a)∕2.
To account for the hydraulic force, we need to add the work done by the fluid pressure within the crack to the total energy (Bourdin et al. 2012;Chukwudozie et al. 2019) In arriving at Eq. (19), we used the following approximation:

Numerical implementation
The mass conservation is simply given by Bourdin et al. (1 − s) n ds.
Given the low permeability of intact rock salt, we neglect fluid leak-off to the rock formation and assume that the pressure within the fracture is spatially constant (ı.e., no pressure loss) because of the small sample size. Furthermore, we consider no material property dependency of pressure, and thus, the energy functional Eq. (19) (19). We descretize the directional derivative with a Galerkin finite-element method with a first-order shape function. Algorithm 1 shows the computational procedure implemented in OpenGeoSys (Bilke et al. 2019). Since we have the irreversibility constraint in solving for d, we employed a variational inequality solver from PETSc library . Table 3 compares characteristics of the employed approaches. Because the variational phase-field model (VPF) assumes no leak-off with and incompressible and inviscid fluid, no hydraulic parameters are considered. As described above, the discrete element method (DEM) employs Minkley model. Since the number of parameters is big, we refer to Minkley et al. (2012); Minkley (2004) for details.

Results
In this section, we discuss the simulation results from the three models presented.

The discrete element method
To simulate the granular structure, a randomised assembly of polyhedral blocks based on the Voronoi-discretization was built (Fig. 6a). The random distribution and grain sizes of the order of 1 cm were used, which are typical values for natural rock salt. Gases or liquids can move on the interfaces between the blocks (Fig. 6b), if permitted by the local stress and pressure conditions. At the center of the sample, a fluid pressure was applied in terms of a boundary condition. Figure 6c shows a vertical cut through the model, with the initial fluid distribution shown in blue. The central borehole is also visible. For the stress configuration, where the minimal principal stress is oriented in the vertical direction, a horizontal crack is expected. The fluid pressure at the center was increased until the first intergranular surfaces began to open, and then held constant. Figure 6d shows the final distribution of the  (c).] The same picture holds for the outer surface of the cube (Fig. 6e), where the fluids leak out forming a horizontal band matching the experimental observations. For the second stress configuration where the minimum principle stress is horizontal, the fluid is expected to expand in a vertical plane. Using the same model, this is indeed the case (Fig. 7).
On a Desktop-PC, with an Intel i7 CPU and four cores, a single calculation took about 10 h.

The coupled hydro-mechanical lattice model
The dual lattice model was implemented to simulate the fluid-driven percolation in the rock salt samples. The applied hydraulic pressures are transformed into the mechanical model using the weak coupling scheme, and subsequently, the elements' failure and change of hydraulic aperture are determined and transformed back to the hydraulic model. The considered mass conservation law results in the prediction of flow rate and change of reservoirs pressure as well as the flow and fracturing paths, which are then compared to the experimental data.
The total number of conduit and mechanical lattice elements are approximately 853,800 and 225,720, respectively (Fig. 8). The simulations were performed on a Desktop-PC with a Xeon processor (2.10 GHz) with a total number of 16 cores and the computational time for a single simulation is around 12 h. For the lattice element simulations, the required mechanical parameters were taken from Table 1. For the hydraulic parameters, the fluid properties of oil were assigned, where the fluid density is 870 kg/m 3 , the fluid viscosity is 6.5 × 10 −6 m/s 2 , the fluid bulk modulus is 1,4 GPa, and the initial hydraulic aperture is 1 × 10 −9 m. The two experiments were simulated using the dual LEM (Fig. 9). In these simulations, the Young's modulus was assumed to be 25 GPa. The crack propagates on the horizontal plane (visible on the surface fracture path) for the first stress configuration (Fig. 9a) and on the vertical plane for the second configuration (Fig. 9b) similarly to the experimental result (Fig. 2).

The variational phase-field model
Relying on the symmetry of the domain, we simulated a 1/4 of the domain (Fig. 10). We discretized the computation domain with first-order tetrahedral elements (27,917,126 elements with 5,432,325 nodes 3 ) and applied the two different boundary loadings as specified in Fig. 2a and b. Table 1 lists Young's modulus E, Poisson's ratio , and the fracture toughness G c used in the simulations. We used 0.75 mm for the regularization length , which is around 1.5 times of the tetrahedral mesh size used in the computations.
Without prescription of fracture nucleation or propagation, the crack set and the displacement fields are obtained through minimization of the total energy at each quasi-static step. Because of the compression-tension split implemented in this model, the phase field (damage) tends to evolve where For the first stress configuration, a crack propagates mostly horizontally towards the boundary of the sample as the vertical stress is the minimum principal stress (Fig. 11a).
The crack initiates with some angle at the bottom of the injection borehole and gradually turns to align with the horizontal plane which is orthogonal to the maximum loading direction. As expected, for the second configuration, the crack propagates on the vertical plane, but its propagation is hinged by the presence of the casing tube in the borehole (Fig. 11b). Table 4 lists simulated breakdown pressures from each approach.

Discussion
The variational phase model is a continuum-based approach. The governing equations are discretized with a continuous Galerkin finite-element method with linear interpolation in this study. The model represents a crack, which is a discontinuous sharp interface in reality, with a mathematically diffused phase-field variable. Because the mesh does not need to conform to the discontinuities (ı.e., cracks), the crack propagation is not restricted by the prescribed mesh. The mesh size and orientation impact the computed fracture topology to some extent, but the approach is capable of recovering the theoretical critical energies accurately as long as the length scale parameter is properly set (Tanné et al. 2018;). If one requires explicit properties such as crack openings or exact crack locations for the fracture conductivity or frictional shearing (Fei and Choo 2020;Bryant and Sun 2021), several post-processing techniques have been proposed (Ziaei-Rad et al. 2016;Yoshioka et al. 2020;Yang et al. 2021). Although we simplified the mass conservation in this study, considering the so-called  On the other hand, both the discrete element and lattice element methods use distinct and discontinuous entities (discrete elements and lattice elements) to account for the mechanical response. As a result, cracks can manifest as breakages between the entities (distinct elements or lattice elements). It is straightforward to obtain the explicit crack properties, because cracks are represented explicitly as sharp discontinuities. However, this representation limits the crack propagation to the interfaces of the elements. The fluid flow is also only allowed at the element interfaces and no leak-off to the rock formation is considered in the current implementation, but it is not limited by the approach (Lefort et al. 2020). As for the element size effect, Potyondy and Cundall (2004) point out that it impacts the critical stress for fracture propagation and its roughness needs to be comparable to the actual fracture surface. This issue was addressed by applying the smooth joint model and the model was verified with hydraulic fracture analytical models Damjanac and Cundall (2016).
As far as the crack nucleation is concerned, none of the models presented require a prescribed initial crack to facilitate the nucleation. The variational phase-field model seeks for crack nucleation through the total energy minimization. However, if the problem does not have a clear stress concentration, the minimization algorithm may be trapped in a local minimum and the crack nucleation may not be unique in this situation. This issue around uniqueness was addressed in a recent study by applying the second-order optimality conditions in the constrained minimization Baldelli and Maurini (2021).
For the DEM and the LEM, a crack nucleates "naturally" by breaking the bond that experiences the highest load in the system. However, similarly to the variational phase-field model case, if there is no clear stress concentration in the problem, this highest load may well be the result of element irregularity. In this situation, the nucleation can also be "non-unique". This nucleation uniqueness in the DEM or LEM with varying element size will be a topic of future research.
In the two loading scenarios considered in this study, even though we did not prescribe an initial crack, a stress concentration was present at the open section of the borehole. Combining it with the anisotropic boundary loadings, all the models employed in this study were able to identify a propagation plane orthogonal to the minimum stress orientation. These planes qualitatively match the experimental observations. While the variational phase-field model results show relatively smooth crack surfaces because of the homogeneous mechanical properties, the DEM and the LEM results show more irregular propagation out of the main crack, which were likely induced by the heterogeneous element sizes and placements in the model.
Since geomaterials do possess such irregularities between grain boundaries, simulated failure patterns from DEM and LEM seem more realistic, while the variational phase field still lacks such intergrain boundary interactions in the implementation. And because of these irregularities, one may conclude that the DEM or LEM is more suitable for the simulation of cracking in rock. However, we should note that these irregularities (heterogeneities) in materials are difficult to characterize in reality especially in grain scale. If heterogeneities are characterized and assigned, a continuum approach such as the variational phase-field model would generate an irregular fracture topology. Alternatively, generation of statistically equivalent grain distributions can be used to study crack nucleation and size dependency in rock . A more systematic study on this subject may be for future research.

Conclusions
We have presented three different numerical approaches for crack propagation by fluid pressure. None of the approaches requires prescribed fracture nucleation points or propagation paths. Using fluid percolation experiment performed on rock salt with different triaxial boundary loadings, we have compared three models against experimental results. The three models simulated crack propagation on the plane orthogonal to the minimum stress in accordance with the experimental observations. This predictive capability is important in Fig. 11 The simulation of pressure-driven fluid percolation for a first case with the vertical stress and b second case with the horizontal stress being the lowest principal stress. Crack images are mirrored on the symmetry plane for visualization subsurface integrity analysis where we cannot assume crack nucleation or propagation a priori. In large-scale models, while it is not impossible, computational cost can increase quickly with the discrete approaches. To bridge such gaps between scales, multi-scale modeling framework is a possible future study.