# Coupled crystal plasticity-cohesive zone modeling of rock salt viscoplasticity

Nour Habib, Saber El Arem, Amine Ammar

► To cite this version:

Nour Habib, Saber El Arem, Amine Ammar. Coupled crystal plasticity-cohesive zone modeling of rock salt viscoplasticity. Finite Elements in Analysis and Design, 2025, 252 (1), 10.1016/j.finel.2025.104438. hal-05247232

HAL Id: hal-05247232

https://hal.science/hal-05247232v1

Submitted on 9 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

# Coupled crystal plasticity-cohesive zone modeling of rock salt viscoplasticity

Nour Habib<sup> ID </sup>*, Saber El Arem, Amine Ammar

LAMPA Arts et Metiers Institute of Technology, 2 boulevard du Ronceray, 49035 Angers cedex 01, BP 93525, France

## ARTICLE INFO

**Keywords:**

Cohesive zone model

Rate dependent traction separation law

Viscoplasticity

Rock salt

## ABSTRACT

Rock salt, owing to its viscoplastic behavior and structural integrity under high pressure, is a promising candidate for safe and large-scale underground energy storage. This study presents a comprehensive numerical framework for modeling the viscoplastic deformation of rock salt, accounting for both intragranular and grain boundary (GB) deformation mechanisms. Intragranular deformation is modeled using a crystal plasticity approach governed by a power-law relation, capturing the activity of crystallographic slip systems. Concurrently, a cohesive zone model (CZM) is introduced to simulate grain boundary sliding (GBS) and opening via a rate-dependent traction-separation law. This modeling strategy enables a detailed analysis of the coupled interplay between crystal plasticity and intergranular decohesion phenomena.

## 1. Introduction

The accelerating transition toward renewable energy necessitates the development of efficient and sustainable storage solutions. Among alternative energy carriers, hydrogen is particularly promising, owing to its high energy density and substantial potential for decarbonizing major industrial sectors. However, the widespread adoption of hydrogen technologies fundamentally depends on the availability of safe, reliable, and scalable storage infrastructures. In this context, underground storage within rock salt formations is increasingly recognized as an optimal solution. This prominence derives from the distinctive physical attributes of rock salt, including its low permeability, inherent self-healing capability, and exceptional mechanical stability under high confining pressures. These features have inspired numerous investigations into the behavior and underlying deformation mechanisms of rock salt [1–4].

When subjected to stress, crystalline materials such as rock salt undergo irreversible deformation primarily through dislocation motion, which distorts the crystal lattice. Rock salt crystallizes in a face-centered cubic (FCC) structure, with sodium and chloride ions occupying the lattice sites. Due to the electrostatic interactions between these ions, deformation is preferentially accommodated along specific crystallographic planes and directions [5]. The dominant slip system in halite is the {110} plane in the 〈110〉 direction [1], while slip along the cubic {100} planes [2] and, at elevated temperatures along octahedral systems [6], has also been observed. In addition to the formation of slip lines on grain surfaces, GBS plays a significant role in the overall deformation of polycrystalline rock salt [7,8]. The efficiency of GBs in either impeding or facilitating dislocation motion is strongly influenced by their crystallographic orientation and impurity content. Grain size strongly influences whether crystal slip plasticity (CSP) or GBS deformation predominates. In fine-grained materials, GBS is favored, as predicted by the Hall-Petch effect, wherein smaller grains increase the relative GB area and thus enhance sliding [9,10].

Experimental investigations into the mechanical response of rock salt have provided the foundation for various constitutive models, most of which adopt a phenomenological, macroscopic perspective. These models generally fall into two categories:

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: nour.habib@ensam.eu (N. Habib).</footnote>
1. Long-term creep models, which account for time-dependent deformation and describe the accumulation of creep strain [11],

2. Viscoplastic models, which incorporate rate-dependent plastic strain and internal variables [12,13].

In this work, we conduct a detailed numerical investigation of the viscoplastic deformation mechanisms in rock salt, emphasizing the interplay between CSP and GBS. By explicitly resolving both intragranular and intergranular phenomena and their dependence on microstructural features, we aim to provide a predictive modeling framework for the reliable assessment of rock salt behavior under the demanding conditions relevant to large-scale hydrogen storage.

We implement a viscoplastic power-law model within the grains to account for CSP, while GBS is represented by a CZM situated at the interfaces between grains. Early numerical studies of ductile polycrystals predominantly focused on intragranular plasticity, often neglecting the effects of GBs. In these early approaches, the computational mesh was typically constructed from regularly arranged single elements, each representing an individual grain [14,15]. Such simplifications yielded largely uniform grain-scale deformation and overlooked the local heterogeneity intrinsic to real polycrystals. Subsequent advancements, such as those by Evers et al. [16], introduced more refined meshing techniques and finite element formulations, allowing for the modeling of dislocation pile-up at GBs. Anand [17] further advanced this line of research by adopting Voronoi tessellations, enabling a more realistic representation of grain morphology and heterogeneity in the simulation of polycrystalline deformation. [18] marked a key development by coupling crystal plasticity to cohesive interface elements to account for GBS and separation phenomena. This approach introduces specialized interface elements. These elements, governed by traction-separation laws, enable a phenomenological description of microstructural processes leading to crack initiation. When the local strength and energy release rate criteria are met, decohesion occurs, allowing cracks to propagate along GBs. This methodology provides a robust framework for analyzing the complex interplay between intragranular deformation and GB effects. Building upon these foundations, later studies have highlighted further complexities. Kamaya et al. [19] demonstrated the critical influence of GB orientation and the associated stress concentrations arising from deformation constraints imposed by neighboring grains. Such factors, when combined with the inclination of boundaries, significantly modulate the local stress field within polycrystalline aggregates. Musienko and Cailletaud [20] further refined the approach by simulating stress corrosion cracking (SCC) in polycrystalline materials using Voronoi tessellations. In their models, cohesive elements of finite thickness represent GBs, and both sliding and opening are identified as dominant modes of failure. The evolution of stress intensity along boundaries is characterized using appropriate fracture mechanics parameters, illustrating the central role of GBS in crack propagation and material degradation. Recent developments have coupled crystal plasticity with GB effects in increasingly complex settings. For instance, [21] employed a crystal plasticity finite element method (CPFEM) to simulate behavior in additively manufactured metals while incorporating experimentally informed GB interactions. Similarly, [22] investigated the interaction between plasticity and intergranular fracture using a bilinear mixed-mode CZM, effectively coupling crystal slip mechanisms with GB decohesion to capture the onset and evolution of intergranular damage.

Given the rate-dependent viscoplastic behavior of rock salt, it is natural to extend this feature to its GBs. Accordingly, the CZM considered in this study is based on a rate-dependent traction-separation law. The computational implementation is organized into two principal components. First, the bulk grain response is described by a crystal plasticity model governed by a viscoplastic power law and incorporating the characteristic slip systems of rock salt. Although the rate-dependent crystal plasticity and CZMs have traditionally been used in metals and composite materials, their underlying physical mechanisms are well suited to capturing the dominant deformation modes in rock salt. Specifically, the intragranular viscoplastic deformation of halite is primarily controlled by crystallographic slip, analogous to that observed in metals, which justifies the use of crystal plasticity theory in this context. Second, the GB interfaces are represented as zero-thickness cohesive zones, modeled as zero-thickness interfaces by duplicating nodes at grain junctions and incorporating cohesive elements that adhere to the rate-dependent traction-separation law. GBS and decohesion observed in rock salt under loading share mechanistic similarities with delamination in composites. Therefore, a cohesive zone approach governed by a rate-dependent traction-separation law is appropriate for modeling these intergranular phenomena. The novelty of the present study lies in the coupling of these two mechanisms within a unified computational framework specifically designed for simulating the viscoplastic behavior of rock salt. This dual-mechanism approach allows for the detailed simulation of the interplay between intragranular deformation and intergranular decohesion, an interaction that is particularly relevant for rock salt under loading but has not been adequately addressed in previous modeling studies. Furthermore, the framework introduces a new quantitative method for characterizing GBS and opening, offering enhanced insight into localized deformation phenomena and their influence on macroscopic mechanical response.

The following sections present, in turn, the detailed constitutive formulation for the bulk (intragranular) response, followed by the interface (GB) model. A numerical experiment on a bicrystal is then introduced to illustrate interface behavior. Finally, we describe the fully coupled implementation for a polycrystalline halite aggregate and investigate the interactions between bulk plasticity and interface mechanisms by systematically varying the ratio of the plastic slip threshold to the cohesive strength.

## 2. Grain bulk constitutive model: Crystal plasticity

The mechanical behavior of both halite single crystals and pure NaCl crystals has been extensively studied by numerous researchers. The viscoplastic behavior of single-crystal rock salt is highly anisotropic due to its crystallographic structure, necessitating a modeling approach that accurately captures its deformation mechanisms. Classical macroscopic viscoplastic models, which typically assume isotropy and homogenize the material's response, are insufficient to capture the essential features of slip-system-dependent plasticity exhibited by NaCl. In contrast, crystal plasticity models explicitly account for the activation of individual slip systems, thereby providing a more accurate and physically grounded description of strain localization, texture evolution, and
<table><caption>**Table 1**
Slip systems of rock salt.</caption><thead><tr><th colspan="2">Dodecahedral Slip Systems</th></tr><tr><th>Slip plane normal</th><th>Slip direction</th></tr></thead><tbody><tr><td>(011)</td><td>[0 1 1]</td></tr><tr><td>(101)</td><td>[1 0 1]</td></tr><tr><td>(110)</td><td>[1 1 0]</td></tr><tr><td>(0 1 1)</td><td>[0 1 1]</td></tr><tr><td>(1 0 1)</td><td>[1 0 1]</td></tr><tr><td>(1 1 0)</td><td>[1 1 0]</td></tr></tbody><thead><tr><th colspan="2">Cubic Slip Systems</th></tr><tr><th>Slip plane normal</th><th>Slip direction</th></tr></thead><tbody><tr><td>(100)</td><td>[0 1 1]</td></tr><tr><td>(100)</td><td>[0 1 1]</td></tr><tr><td>(010)</td><td>[1 0 1]</td></tr><tr><td>(010)</td><td>[1 0 1]</td></tr><tr><td>(001)</td><td>[1 1 0]</td></tr><tr><td>(001)</td><td>[1 1 0]</td></tr></tbody></table>

rate-dependent deformation. This framework is particularly well-suited for simulating the time-dependent, anisotropic mechanical response of halite, offering a robust foundation for high-fidelity numerical simulations and predictive modeling.

In the present study, the bulk grain response is modeled using a phenomenological crystal plasticity framework, wherein plastic deformation is governed by the critical resolved shear stresses (CRSS) associated with crystallographic slip systems. This approach does not explicitly track microstructural features such as dislocation density, but it captures the salient features of plastic flow by resolving slip on a set of representative systems. Specifically, the viscoplastic response of rock salt is described by incorporating slip on six dodecahedral and six cubic systems, which together account for the dominant deformation mechanisms observed experimentally.

## 2.1. Resolved shear stress

Plastic deformation in crystalline materials results from the rearrangement of atoms along specific crystallographic planes and directions in a new stable configuration. This phenomenon, governed by the law of Schmid and Boas, is activated when the resolved shear stress on a given slip system reaches a critical value, known as the CRSS, denoted $\tau_c$. The resolved shear stress depends on the applied Cauchy stress and the Schmid tensor, $\mu_\alpha$.

$$ \mu_{\alpha} = \frac{1}{2} (\mathbf{m}^{\alpha} \otimes \mathbf{n}^{\alpha} + \mathbf{n}^{\alpha} \otimes \mathbf{m}^{\alpha}) \qquad (1) $$

Here, $\mathbf{n}^\alpha$ and $\mathbf{m}^\alpha$ denote the unit normal to the slip plane and the unit vector in the slip direction for slip system $\alpha$, respectively. For rock salt, the slip systems considered in this work are shown in Table 1.

## 2.2. Plastic deformation and finite strain kinematics

The finite deformation of crystalline solids is described by a multiplicative decomposition of the total deformation gradient, $\mathbf{F}$, into an elastic part, $\mathbf{F}^e$, and a plastic part, $\mathbf{F}^p$:

$$ \mathbf{F} = \mathbf{F}^e \mathbf{F}^p \qquad (2) $$

Here, $\mathbf{F}^e$ characterizes reversible lattice stretching and rotation, while $\mathbf{F}^p$ represents the irreversible plastic distortion due to slip.

The evolution of the deformation gradient is governed by the velocity field of the material. If $\mathbf{v} = \frac{d\mathbf{u}}{dt}$ denotes the material velocity, then the spatial velocity gradient is $\mathbf{L} = \nabla \mathbf{v}$. In terms of the deformation gradient, this becomes:

$$ \mathbf{L} = \dot{\mathbf{F}} \mathbf{F}^{-1} \qquad (3) $$

This formalism enables a rigorous description of the kinematics of finite deformations and provides the basis for the subsequent constitutive modeling of crystal plasticity.

## 2.3. Shear rate and plastic flow

The evolution of plastic deformation in crystals is described by the plastic velocity gradient, $\mathbf{L}^p$, which is related to the plastic deformation gradient, $\mathbf{F}^p$, by:

$$ \mathbf{L}^p = \dot{\mathbf{F}}^p (\mathbf{F}^p)^{-1} \qquad (4) $$
This gradient can be expressed as a sum over all active slip systems:

$$L^p = \sum_a \dot{\gamma}^a (\mathbf{m}^a \otimes \mathbf{n}^a) \qquad (5)$$

where $\dot{\gamma}^a$ is the shear rate associated with slip system $\alpha$, and $\mathbf{m}^a$, $\mathbf{n}^a$ are the corresponding slip direction and slip plane normal unit vectors.

The shear rate for each slip system is governed by a viscoplastic power-law relation that incorporates both isotropic and kinematic hardening, as well as a threshold for slip activation:

$$\dot{\gamma}^a = \dot{a} \langle |\tau^a - x^a| - g^a \rangle_+^n \text{sign}(\tau^a - x^a) \qquad (6)$$

Here, $\dot{a}$ is the reference shear rate, $n$ is the stress sensitivity exponent, $\tau^a$ is the resolved shear stress, $x^a$ is the backstress (kinematic hardening), and $g^a$ is the critical resolved shear stress (isotropic hardening). The Macaulay brackets $\langle X \rangle_+ = \max(X, 0)$ ensure that plastic slip only occurs once the yield criterion is met.

## 2.4. Hardening laws

Kinematic hardening, which accounts for the translation of the yield surface in stress space, evolves as function of the shear rate as:

$$\dot{x}^a = c \dot{\gamma}^a - d |\dot{\gamma}^a| x^a \qquad (7)$$

where $c$ and $d$ are material constants.

Isotropic hardening is modeled by allowing the critical resolved shear stress $g^a$ to evolve with the contribution of slip on different slip systems due to interactions with other systems:

$$g^a = g_0^a + H_0 \sum_{\beta} h^{\alpha\beta} [1 - \exp(-g_{\text{sat}}v^{\beta})] \qquad (8)$$

$$\dot{v}^{\beta} = |\dot{\gamma}^{\beta}| \qquad (9)$$

The term $h^{\alpha\beta}$ represents the hardening rate on slip system $\alpha$ due to shearing on slip system $\beta$ [23]. It provides insight into the interactions between dislocations across different slip systems. $h^{\alpha\beta}$ is the latent hardening matrix that accounts for the interaction between different slip systems.

$$h^{\alpha\beta} = \begin{pmatrix} 1 & q & q & \cdots & q \\ q & 1 & q & \cdots & q \\ q & q & 1 & \cdots & q \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ q & q & q & \cdots & 1 \end{pmatrix} \qquad (10)$$

where $q$ is the latent hardening factor, representing the ratio of latent hardening to self-hardening. For coplanar slip systems, $q$ is typically equal to 1. $H_0$ is the initial hardening rate, and $g_{\text{sat}}$ is the saturation value of the slip resistance.

## 2.5. Governing equations

In summary, the constitutive model is completed by enforcing mechanical equilibrium in the absence of body forces and inertia:

$$\text{div } \mathbf{P} = 0 \quad \text{in } \Omega \qquad (11)$$

where $\mathbf{P}$ is the first Piola–Kirchhoff stress tensor. Dirichlet boundary conditions are imposed on $\Gamma_D$:

$$\mathbf{u} = \mathbf{u}_D \quad \text{on } \Gamma_D \qquad (12)$$

While all forces are null on the boundary:

$$\mathbf{T} = 0 \quad \text{on } \Gamma_F \qquad (13)$$

The finite-strain kinematics are described by the multiplicative decomposition:

$$\mathbf{F} = \mathbf{F}^e \mathbf{F}^p \qquad (14)$$

and the evolution of plastic slip is governed by the power-law relation with hardening and threshold effects:

$$\dot{\gamma}^a = \dot{a} \langle |\tau^a - x^a| - g^a \rangle_+^n \text{sign}(\tau^a - x^a) \qquad (15)$$

This framework enables the accurate simulation of time-dependent, anisotropic plasticity in rock salt crystals under complex loading conditions.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Stress distribution near the crack tip, after [28].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Typical traction-separation relationship for a cohesive interface.</figcaption></figure>

## 3. Grain interface constitutive model: Cohesive zones

### 3.1. Phenomenological foundations of the cohesive zone model

Structural imperfections, such as voids and cracks, induce pronounced stress concentrations in their immediate vicinity, which can lead to substantial magnification of local stresses. Inglis [24] was the first to analyze this phenomenon, showing that, in the vicinity of a sharp flaw, the theoretical stress may approach infinity as the material's resistance to rupture is locally diminished. Building on these ideas, Griffith [25] postulated that crack propagation ensues when the energy available for crack advance exceeds the intrinsic resistance of the material. At the atomic level, fracture initiation occurs when the applied mechanical work is sufficient to break the atomic bonds, held together by electrostatic forces between charged particles. The CZM is used to represent the crack propagation process within a structure. This approach introduces cohesive interfaces that simulate the progressive degradation and eventual separation of material, confining the damage to the interface which is also the assumed crack path. Inside the cohesive zone, internal stresses serve as restorative forces, comparable to molecular or atomic attractions, working to maintain the cohesion of the crack surfaces. The magnitude of this force depends on the separation distance between the crack surfaces. In the Dugdale model [26], the cohesive stress is constant and equivalent to the bulk yield stress, as illustrated in the right crack tip in Fig. 1, whereas in the Barenblatt model [27], it varies with the displacement jump at the crack tip as illustrated in the left crack tip in Fig. 1.

The mechanical behavior of the cohesive interface is governed by a constitutive relation linking the traction across the interface to the displacement jump, commonly referred to as the traction-separation law (TSL). As external loading increases, the separation between crack faces, denoted $\delta$, grows until a critical threshold $\delta_c$ is reached. Beyond this point, the cohesive bond is irreversibly broken and the traction decreases, as depicted in Fig. 2. Total decohesion occurs when the traction falls to zero, resulting in the creation of two new surfaces, denoted $\Gamma^+$ and $\Gamma^-$.

### 3.2. The traction-separation law

The displacement jump across the interface is the vector $\delta = \begin{Bmatrix} \delta_1 \\ \delta_2 \\ \delta_3 \end{Bmatrix}$ with :

$$ \delta_i = \delta_i^+ - \delta_i^- \qquad (16) $$

where i denotes the normal (i = 1) and tangential (i = 2, i = 3) components. While Barenblatt originally formulated the traction law as a function of the distance along the crack tip, contemporary formulations typically express the traction as a function of the
separation, $\delta$, across the cohesive interface. Given the small physical extent of the process zone, the exact constitutive law is often established phenomenologically, based on experimental data [29].

The initial application of the CZM to concrete fracture by Hillerborg [30] led to the concept of brittle TSLs, characterized by the immediate onset of damage upon any separation. In contrast, ductile behavior is represented by TSLs exhibiting an initial elastic regime. Various TSL shapes have since been proposed: quadratic [31], exponential [32,33], and bilinear [34] models each capture distinct aspects of fracture behavior. Furthermore, real materials frequently display a combination of opening (Mode I) and sliding (Mode II) fracture. To address such mixed-mode effects, the notion of effective separation is introduced, as developed in the works of Alfano and Crisfield [35], Mi and Crisfield [36], and Camanho and Davila [37]. The TSL is fundamentally linked to the energy release rates G for each mode of fracture, with the area under the traction-separation curve corresponding to the critical energy release rate, $G_c$, required for crack propagation. In each interface element, the relative displacement can thus be decomposed into normal and tangential components to distinguish between opening and sliding mechanisms.

### 3.3. Rate-dependent traction-separation laws

In this study, particular attention is given to the rate dependence of the cohesive response. According to Bažant and Oh [38], the time-dependent nature of fracture arises from two primary sources: the viscoelastic behavior of the bulk material and the rate processes governing the breakage of bonds within the process zone, leading to a softening law for crack opening that is inherently rate-dependent. Corigliano et al. [39,40] further elaborated viscoplastic models that accommodate fracture by introducing rate-dependent traction-separation and damage evolution laws. The work of Allen et al. [29] explored the time-dependent response in viscoelastic fibril structures, while Seidel et al. [41] incorporated damage evolution within a similar framework. In the context of rock mechanics, Aliguer et al. [42] implemented time-dependent CZMs to investigate the homogenized behavior of rock masses, employing a viscoplastic continuum formulation at the interface for efficient numerical integration.

For the present work, we adopt a trapezoidal form of the traction-separation law, as this best reflects the experimentally observed behavior of halite polycrystals [43]. The traction depends on the displacement rate and a scalar damage parameter that governs the progressive degradation of the interface stiffness and ultimate decohesion.

#### 3.3.1. Thermodynamic potential

To derive the traction-separation law at the interface from first principles, we consider the Helmholtz free energy of the cohesive interface. The free energy $\Psi$ is taken as a function of the elastic displacement jump, $\delta_{el}$, and a scalar damage variable, $d$:

$$ \Psi = \Psi(\delta_{el}, d) \qquad (17) $$

The total separation is additively decomposed into elastic and viscoplastic components,

$$ \delta = \delta_{el} + \delta_{vpl} \qquad (18) $$

This choice is justified by the use of a zero-thickness interface, not volumetric cohesive elements, with a separation based on relative displacement jumps. While a deformation gradient is computed at the interface, it is done to extract the local rotation tensor and to transform between the local coordinates at the interface and the global coordinates of the material. The additive decomposition offers an intelligible distinction between the reversible elastic separation and irreversible inelastic mechanisms leading to the decohesion at the GBs. This leads to the energy potential:

$$ \Psi = \frac{1}{2} g(d) K (\delta - \delta_{vpl})^2, \qquad (19) $$

where K is the cohesive stiffness and $g(d) = (1 - d)$ is the degradation function reflecting the damage state.

The time derivative of $\Psi$ yields

$$ \dot{\Psi} = \frac{\partial \Psi}{\partial \delta_{el}} \dot{\delta}_{el} + \frac{\partial \Psi}{\partial d} \dot{d} \qquad (20) $$

which, after substitution, gives

$$ \Psi = (1 - d) K (\delta - \delta_{vpl}) \dot{\delta} - \frac{1}{2} K (\delta - \delta_{vpl})^2 \dot{d} \qquad (21) $$

Imposing the Clausius-Duhem inequality ensures thermodynamic admissibility:

$$ (\tau - (1 - d) K (\delta - \delta_{vpl})) \dot{\delta} + \frac{1}{2} K (\delta - \delta_{vpl})^2 \dot{d} \geq 0 \qquad (22) $$

which leads to the constitutive relation for the cohesive traction,

$$ \tau = \frac{\partial \Psi}{\partial \delta} = (1 - d) K (\delta - \delta_{vpl}) \qquad (23) $$

with the constraints $d \geq 0$ and $d \leq 1$.

The evolution of damage is typically taken as linear between a softening threshold $\delta^s$ and the critical separation at failure $\delta^f$:

$$ d = \frac{\delta^f (\delta_{max} - \delta^s)}{\delta_{max} (\delta^f - \delta^s)} \qquad (24) $$
where $\delta_{max}$ is the maximum separation reached during the loading history, enforcing irreversibility.

The viscoplastic component of the separation is governed by a Perzyna-type overstress model, which introduces rate dependence into the cohesive response. This is particularly suitable for modeling rock salt, whose microstructure exhibits time-dependent deformation even at low stress levels due to crystal plasticity and diffusion-related mechanisms. The Perzyna formulation allows for a smooth transition from interface bonding to sliding, reflecting the gradual mobilization of inelastic slip observed at GBs. This rate-dependent model introduces a plateau in the traction-separation response, capturing the stable phase of inelastic sliding during which the interface accommodates deformation without immediate loss of load-bearing capacity, a feature that shows more accurate results when simulating crack evolution in rock salt. The viscoplastic rate is defined by:

$$ \dot{\delta}_{vpl} = \dot{\gamma} \frac{\partial F}{\partial \tau}, \quad \dot{\gamma} = \frac{\langle F(\tau) \rangle_{+}}{\eta} \qquad (25) $$

with $\langle X \rangle_+ = \max(0, X)$, $\eta$ the viscosity, and $F(\tau) = |\tau| - \tau_c$ the yield function. The standard Karush-Kuhn-Tucker (KKT) conditions apply:

$$ \gamma \ge 0, \quad f(\tau) \le 0, \quad \gamma f(\tau) = 0 \qquad (26) $$

The numerical resolution of this model will be computed by way of a regularization of the perfect plasticity model, based on the assumption that the viscoplastic solution relaxed to the rate-independent solution as $t \to \infty$. This method is detailed in [44]. The period of relaxation is introduced as:

$$ r = \frac{\eta}{K} \qquad (27) $$

So the viscoplastic evolution equation can be rewritten as:

$$ \dot{\delta}_{vpl} = \frac{\tau - \tau_c}{rK} \frac{\tau}{|\tau|} \qquad (28) $$

The initial value problem to obtain the time evolution of the stress is of the form:

$$ \dot{t} + \frac{\tau}{r} = K \dot{\delta} + \frac{\tau_c}{r} \qquad (29) $$

Considering that in a relaxation test the deformation is kept constant, we can draw the solution as:

$$ \tau(t) = (\tau - \tau_c) \exp(-\frac{t}{\tau}) + \tau_c \qquad (30) $$

$\tau_c$ is the rate independent solution, the asymptotic limit toward which the solution converges in a relaxation test.

$$ \tau_{\infty} = \tau_c \qquad (31) $$

This leads to the expression of the rate dependent traction as :

$$ \tau = \frac{K(\delta - \delta_{vpl})}{1 + \frac{\Delta t}{r}} + \frac{\frac{\Delta t}{r}}{1 + \frac{\Delta t}{r}} \tau_c \qquad (32) $$

### 3.3.2. Constitutive equations for the traction-separation law

The traction vector across the interface is defined as $\tau = \begin{cases} \tau_1 \\ \tau_2 \\ \tau_3 \end{cases}$ where i denotes the normal ($i = 1$) and tangential ($i = 2, i = 3$) components. The normal displacement jump $\delta_1$ is decomposed into two parts:

$$ \delta = \delta_{\text{active}} + \delta_{\text{inactive}} \qquad (33) $$

* Active displacement, which corresponds to the positive (opening) part of the displacement:

$$ \delta_{\text{active}} = \langle \delta_1 \rangle_+, \quad \langle X \rangle_+ = \max(X, 0) \qquad (34) $$

* Inactive displacement, which accounts for compressive (closing) behavior with negative component value:

$$ \delta_{\text{active}} = \langle \delta_1 \rangle_+, \quad \langle X \rangle_+ = \max(X, 0) \qquad (35) $$

The complete expression of the rate-dependent TSL constitutive equations with a trapezoidal form (see Fig. 3) is:

$$ \tau_1 = \begin{cases} K_1 (\delta_{\text{active}} - \delta_{\text{active,vpl}}) + K_{\text{inactive}} \delta_{\text{inactive}}, & \delta_{\text{active}}^{\max} \le \delta_1^0 \\ \tau_c + K_{\text{inactive}} \delta_{\text{inactive}}, & \delta_1^0 < \delta_{\text{active}}^{\max} < \delta_{\text{softening}} \\ (1-d)K_1 (\delta_{\text{active}} - \delta_{\text{active,vpl}}) + K_{\text{inactive}} \delta_{\text{inactive}}, & \delta_{\text{softening}} < \delta_{\text{active}}^{\max} < \delta_f \\ K_{\text{inactive}} \delta_{\text{inactive}}, & \delta_{\text{active}}^{\max} \ge \delta_1^f \end{cases} \qquad (36) $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. Schematic representation of the traction-separation law.</figcaption></figure>

for the normal component ($i = 1$). And thus, the negative displacement is not considered in the damage evolution and does not contribute to the degradation of the cohesive elements. For the shearing components with $i = 2$ or $i = 3$, we write:

$$ \tau_i = \begin{cases} K_i (\delta_i - \delta_{i,vpl}), & \delta_i^{\max} \le \delta_i^0 \\ \tau_c, & \delta_i^0 < \delta_i^{\max} < \delta_{softening} \\ (1-d) K_i (\delta_i - \delta_{i,vpl}), & \delta_{softening} < \delta_i^{\max} < \delta_f \\ 0, & \delta_i^{\max} \ge \delta_i^f \end{cases} \qquad (37) $$

The stiffness $K_{inactive}$ corresponds to the response under negative displacement jumps. The complete elimination of interpenetration is not feasible in numerical models, this stiffness is introduced to penalize compressive (inactive) displacements. To effectively limit interpenetration, $K_{inactive}$ is typically chosen to be significantly higher than the normal and tangential stiffnesses associated with active (opening or sliding) jumps.

However, this approach may lead to undesirable numerical effects under cyclic loading, including the accumulation of residual displacement discontinuities. To address this, we introduce a penalty modification to $K_{inactive}$, following the methodology proposed in [45].

$$ K_{inactive} = P K_n \qquad (38) $$

Here, $P$ is a penalty coefficient designed to increase the stiffness response in compression, thereby limiting artificial interpenetration. The penalty function $P$ is defined as a linear function of the compressive (inactive) displacement jump $\delta_{inactive}$:

$$ P = a \delta_{inactive} + c \qquad (39) $$

where $c = 1$ ensures consistency with the active stiffness in the absence of compressive overlap. The slope $a$ is calibrated to achieve a prescribed penalty value $P_{inactive}^c$ at a critical compressive overlap $\delta_{inactive}^c$. The parameter $a$ is thus given by:

$$ a = \frac{P_{inactive}^c - 1}{\delta_{inactive}^c} \qquad (40) $$

This construction yields a penalty stiffness that increases linearly with compressive displacement, thereby enforcing strong resistance to interpenetration. However, to avoid excessive stiffness accumulation under large cyclic compressive loads, which may induce spurious residual discontinuities, a saturation mechanism is introduced. Specifically, the penalty coefficient $P$ is capped at a maximum value of $P_c$:

$$ P = \begin{cases} a \delta_{inactive} + 1, & \text{if } P \le P_c \\ P_c, & \text{otherwise} \end{cases} \qquad (41) $$

This capping strategy ensures both numerical stability and consistency with the underlying physical behavior under cyclic or high-amplitude compressive loading. To maintain mechanical coherence across scales, the maximum penalty value must be selected in accordance with the global material parameters, thereby avoiding significant stiffness mismatches between the grain interiors and GBs. In the present work, the selection of numerical parameters is guided by experimental data obtained from mechanical testing on rock salt specimens, ensuring that the model reflects the intrinsic properties and deformation characteristics of the material.
<table><caption>**Table 2**

Material parameters for the bicrystal example.</caption><thead><tr><th>C<sub>11</sub></th><th>C<sub>12</sub></th><th>C<sub>44</sub></th><th>K<sub>N</sub></th><th>K<sub>S</sub></th><th>τ<sub>c</sub><sup>n</sup></th><th>τ<sub>c</sub><sup>f</sup></th><th>δ<sub>s</sub></th><th>δ<sub>f</sub></th><th>η</th></tr></thead><tbody><tr><td>48.7 × 10<sup>3</sup> MPa</td><td>11.9 × 10<sup>3</sup> MPa</td><td>13.1 × 10<sup>3</sup> MPa</td><td>10<sup>4</sup> MPa/mm</td><td>10<sup>4</sup> MPa/mm</td><td>8 MPa</td><td>8 MPa</td><td>0.02 mm</td><td>0.2 mm</td><td>0.01 MPa s /mm</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Geometry and boundary conditions for the bicrystal shear test.</figcaption></figure>

In real materials and practical engineering applications, fracture processes often involve coupled Mode I (opening) and Mode II (sliding) behavior, wherein the response of one mode influences the other. In compressive stress regimes, such as those relevant to subsurface cavern applications, a frictional coupling mechanism becomes significant. Under such conditions, compressive loading enhances the so-called *inactive* normal traction, which in turn restricts tangential (shear) displacements along the fracture interface. To account for this coupling, the shear yield function is formulated as:

$$F(\tau_s) = \tau_s + \mu\tau_{n,\text{inactive}} - \tau_c^s \quad (42)$$

where $\tau_s$ denotes the tangential traction, $\tau_{n, inactive}$ represents the normal traction arising from compressive (negative) displacement jumps, $\mu$ is the coupling coefficient reflecting frictional resistance, and $\tau_c^s$ is the critical (peak) tangential traction. This formulation implies that under compression, shear failure is inhibited by normal tractions, analogous to frictional contact models. When $\mu = 0$, the shear and normal responses are entirely decoupled, corresponding to independent fracture modes.

### 3.4. Governing equations at the interface

The continuity of tractions across the interface elements is enforced as

$$\tau_i^+ d\Gamma^+ = \tau_i^- d\Gamma^- \quad (43)$$

with the local boundary conditions

$$\sigma \cdot \mathbf{n}^+ = \tau^+ \quad \text{on } \Gamma^+, \quad \sigma \cdot \mathbf{n}^- = \tau^- \quad \text{on } \Gamma^- \quad (44)$$

ensuring equilibrium between the adjacent elements and the interface.

### 3.5. Numerical example: Bicrystal shear test

To illustrate the influence of viscosity in the rate-dependent CZM, we consider a representative numerical experiment involving a bicrystal specimen. The model comprises two rectangular blocks with a total height of 2 mm and width of 5 mm, as shown in Fig. 4. Shear loading is imposed along the x-axis by prescribing a displacement rate $\Delta u_x = 10^{-5}$ mm/s at the top edge ($y = 2$ mm), while the bottom edge ($y = 0$ mm) is fixed with $u_x = 0$. The vertical displacement $u_y$ is constrained to zero along all boundaries to prevent rigid body motion. The material parameters considered in the simulation are summarized in Table 2.

The results, presented in Fig. 5, reveal the sensitivity of the interface response to the viscosity parameter $\eta$. As $\eta$ increases, the model exhibits more pronounced rate-dependent plastic deformation, resulting in higher interface tractions under loading. This behavior arises from the viscoplastic formulation, in which the viscosity parameter governs the rate at which damage and separation evolve once the cohesive strength is exceeded. A higher viscosity delays the onset of the plastic plateau, allowing the interface to temporarily sustain overstress conditions and leading to an increased peak traction. Moreover, the viscous response smooths the traction-separation evolution, mitigating abrupt failure and enhancing numerical stability. Extending this interpretation to polycrystalline aggregates, a higher viscosity at the interfaces would tend to suppress GBS by resisting rapid separation and shear across neighboring grains. Conversely, a lower viscosity accelerates the onset of interface damage and sliding, promoting the concentration of deformation at grain joints and potentially enhancing intergranular damage accumulation.

The interface response is also strongly influenced by the applied loading rate, as demonstrated in Fig. 6. At higher loading rates, the model exhibits higher peak tractions and delayed damage initiation, indicative of rate-dependent strengthening behavior.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Reaction force versus displacement for varying viscosity values.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Reaction force versus displacement for different loading rates.</figcaption></figure>

This effect is also a direct consequence of the viscoplastic formulation, wherein the separation rate drives the evolution of damage and inelastic deformation. Under faster loading, the cohesive interface has less time to accumulate damage, resulting in a higher overstress before significant separation occurs. Consequently, the peak traction increase with loading rate. The interplay between the loading rate and viscosity thus controls not only the magnitude of interfacial forces but also the timing and progression of failure. At the polycrystalline scale, higher loading rates would tend to suppress GBS by limiting time-dependent separation processes at the interfaces. In contrast, lower loading rates provide more time for damage accumulation and sliding to develop at the GBs, particularly in models characterized by low interfacial viscosity.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Duplication of nodes at grain interfaces.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Numerical model schematic: (a) internal mesh, (b) interfaces, (c) grain structure.</figcaption></figure>

# 4. Numerical implementation: Application to halite polycrystal

## 4.1. Microstructural topology

The material under investigation is a polycrystalline aggregate of halite. To realistically capture the natural variability of grain orientations, each grain is assigned a unique orientation, generated through random sets of Euler angles using a quaternion-based approach, as outlined by Diebel [46]. The use of quaternions provides a singularity-free representation of rotations and facilitates uniform sampling, ensuring that the orientation distribution function (ODF) is statistically uniform and the resulting polycrystal exhibits a random texture. The finite element mesh is constructed so that nodes are duplicated at GBs, thereby creating cohesive interfaces between adjacent grains (Fig. 7).

These interfaces, modeled as zero-thickness cohesive zones, enable explicit treatment of intergranular decohesion and crack propagation. As summarized in Fig. 8, the simulation domain consists of bulk elements governed by crystal plasticity constitutive laws, and interface elements described by the traction-separation law.

## 4.2. Model calibration

Model parameters are initially estimated from the literature, then refined by comparison with experimental data. Only the dodecahedral and cubic slip systems are considered, as the octahedral systems are activated exclusively under high temperature and pressure. All numerical implementations and simulations are performed using the open source parallel finite element framework MOOSE [47]. An implicit time integration scheme is employed to solve the coupled viscoplastic crystal plasticity and cohesive zone equations. While adaptive time stepping is not explicitly implemented, the solver automatically reduces the time step in response to convergence issues, particularly during stages involving significant interfacial decohesion at late softening and failure stages.

The elastic properties, crystal plasticity parameters, and cohesive zone parameters used in the simulations are summarized in Tables 3, 4, and 5.
<table><caption>**Table 3**
Elasticity parameters of Rock Salt, taken from [48].</caption><thead><tr><th>C<sub>11</sub></th><th>C<sub>12</sub></th><th>C<sub>44</sub></th></tr></thead><tbody><tr><td>48.7 × 10<sup>3</sup> MPa</td><td>11.9 × 10<sup>3</sup> MPa</td><td>13.1 × 10<sup>3</sup> MPa</td></tr></tbody></table>

<table><caption>**Table 4**
Crystal plasticity parameters.</caption><thead><tr><th>γ<sub>0</sub></th><th>h<sub>0</sub></th><th>g<sub>0</sub><sup>do</sup></th><th>g<sub>0</sub><sup>cub</sup></th><th>g<sub>sat</sub></th><th>q</th><th>m</th><th>c</th><th>d</th></tr></thead><tbody><tr><td>2 × 10<sup>-3</sup> s<sup>-1</sup></td><td>15MPa</td><td>3.5MPa</td><td>17.5MPa</td><td>7MPa</td><td>0.5</td><td>0.5</td><td>6</td><td>4</td></tr></tbody></table>

<table><caption>**Table 5**
Cohesive zone parameters.</caption><thead><tr><th>K<sub>N</sub></th><th>K<sub>N inactive</sub></th><th>K<sub>S</sub></th><th>τ<sub>c</sub><sup>n</sup></th><th>τ<sub>c</sub><sup>s</sup></th><th>η<sub>N</sub></th><th>η<sub>S</sub></th><th>μ</th><th>δ<sub>x</sub></th><th>δ<sub>y</sub></th></tr></thead><tbody><tr><td>10<sup>4</sup> MPa/mm</td><td>10<sup>5</sup> MPa/mm</td><td>5 × 10<sup>3</sup> MPa/mm</td><td>5 MPa</td><td>5 MPa</td><td>0.1 MPa s/mm</td><td>0.1 MPa s/mm</td><td>0.7</td><td>0.08 mm</td><td>0.1 mm</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Model A: tessellation, finite element mesh, and boundary conditions. Average grain diameter = 3.6 mm.</figcaption></figure>

The specimen geometry and boundary conditions are depicted in Fig. 9, it will be called Model A. The modeled domain is a cylinder of height 36 mm and diameter 18 mm, with a cut at $x = 15$ mm, similar geometry to experiments carried out by [43] on rock salt specimen in laboratory setting. The average grain size is set to 3.6 mm, as is the average grain size of the grains in the experimental sample, resulting in a tessellation consisting of 333 grains. The mesh, generated with NEPER [49], consists of 246,923 linear tetrahedral finite elements. Boundary conditions are applied as follows: the base ($z = 0$) is fixed in the vertical direction, while the top ($z = 36$ mm) is subject to a prescribed displacement increment $Δu_z = -10^{-5}$ mm/s to mimic compression, consistent with the laboratory experimental loading conditions too. To eliminate rigid body motions, one node on the bottom boundary is fixed in $x$ and $y$ directions, and a second in the $z$ direction only.

### 4.3. Macroscopic response and model calibration

The macroscopic response, presented in the form of a stress–strain curve, is shown in Fig. 10. The parameters of the CZM were calibrated to reproduce the overall mechanical behavior observed in experimental studies, specifically those reported by [43], and the corresponding dataset is plotted as Lab experimental data.

The penalty parameter associated with the inactive (compressive) displacement jump was set to $K_{N inactive} = 10^5$ MPa/mm, corresponding to a maximum penalty coefficient of $P_c = 10$. This value ensures a sufficiently stiff response under compressive interfacial loading while maintaining coherence with the bulk material stiffness.

Notably, GBS and opening emerge as dominant deformation mechanisms, which are clearly observable at approximately 3% of the total axial deformation in Fig. 11.
<figure><img src="image_13.png" /><figcaption>Fig. 10. Comparison of force-displacement curves: laboratory data versus numerical simulation.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. The deformation at failure, $[F_{zz} - I] = -0.027$.</figcaption></figure>

The equivalent plastic strain is defined as

$$ \epsilon_p = \int_t \sqrt{\frac{2}{3} L_p : L_p} dt, \qquad (45) $$

where $L_p$ is the plastic velocity gradient. Figs. 12–14 show the evolution of the plastic strain, while Figs. 15–17 display the corresponding von Mises stress distributions. Regions with limited GBS are characterized by higher plastic strains. Stress concentrations at GBs depend not only on the orientation of the grains but also on the inclination of the boundary and the grain geometry. Triple junctions consistently exhibit the greatest deformation, as evidenced by sliding and opening along these interfaces (Figs. 18–20).
<figure><img src="image_15.png" /><figcaption>Fig. 12. The plastic strain at 0.8% total strain.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. The plastic strain at 1.4% total strain.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. The plastic strain at 2.7% total strain.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. The von mises stress at 0.8% total strain.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. The von mises stress at 1.4% total strain.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. The von mises stress at 2.7% total strain.</figcaption></figure>

## 4.4. Contributions of grain boundary sliding and opening to overall deformation

A key quantity in the analysis of cohesive interfaces is the surface-averaged tensor, defined by

$$A = \frac{1}{V} \int_{\partial S} \delta_g \otimes n \, ds, \qquad (46)$$

where V is the total volume, ∂S the total GB surface, δg the displacement jump in global coordinates, and n the local normal vector also expressed in global coordinates. This integral provides insights into the interaction between displacement jumps and the boundary normal vectors, enabling a better understanding of interface mechanics. To further analyze the tensor, we decompose it into spherical (volumetric) and deviatoric (shear) components. The spherical part, given by

$$A^{\text{sph}} = \frac{1}{3} \text{tr}(A) \mathbf{I}, \qquad (47)$$
<figure><img src="image_21.png" /><figcaption>Fig. 18. The jump magnitude at the interface at 0.8% total strain.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. The jump magnitude at the interface at 1.4% total strain.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20. The jump magnitude at the interface at 2.7% total strain.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 21. Geometry and boundary conditions for the bicrystal compression test.</figcaption></figure>

quantifies the isotropic expansion or contraction across the interface and is indicative of Mode I (opening) fracture. The deviatoric part,

$$A^{\text{dev}} = A - A^{\text{sph}}, \qquad (48)$$

captures the distortional effects associated with shear-dominated (Modes II and III) fracture mechanisms.

### 4.4.1. Numerical example on bi-crystal

This analysis is first illustrated using a bicrystal subjected to uniaxial compression. The sample geometry and boundary conditions are presented in Fig. 21: a block of height 1 mm and width 2 mm is composed of two grains, with $\Delta u_x = -10^{-5}$ mm/s imposed on the right edge and $u_x = 0$ on the left; $u_y = 0$ at the lower left corner. Material parameters are as in Table 2.

The results, shown in Fig. 22 (blue curve), indicate that the contribution of GBs to axial deformation under compressive loading is null. This is attributed to the exclusion of negative normal displacement jumps, arising from compression, from the displacement-based integration used to compute GB deformation.

A similar simulation conducted under tensile loading, with an imposed displacement rate of $\Delta u_x = +10^{-5}$ mm/s, yields a markedly different response (Fig. 22, orange line). In this case, the quantified deformation arises exclusively from interfacial opening, thereby confirming that the integral effectively captures pure Mode I behavior.

For the shear loading case (see Fig. 23), the contribution of GBS is quantified using the von Mises equivalent of the deviatoric component of the interfacial deformation tensor $A$, as depicted in Fig. 24. Specifically, the equivalent shear measure is defined as:

$$VM_{\text{eq}} = \sqrt{\frac{2}{3} A^{\text{dev}} : A^{\text{dev}}} \qquad (49)$$

This scalar measure captures the accumulation of shear-dominated deformation at the interface and serves as a robust indicator of Mode II sliding activity.
<figure><img src="image_25.png" /><figcaption>Fig. 22. GBS contribution to the total deformation in a compression and tension test.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 23. Geometry and boundary conditions for the bicrystal shear test.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 24. Von Mises equivalent of the GBS contribution under shear loading.</figcaption></figure>

### 4.4.2. Results on the Model A polycrystal

For the polycrystalline aggregate (Model A), the contribution of GBS to the overall axial deformation is illustrated in Fig. 25. At the onset of loading, GB mechanisms play a significant role, accounting for up to 23.5% of the total axial strain. This initial regime is predominantly influenced by the parameters of the CZM, particularly the interfacial stiffnesses $K_N$ and $K_S$. Lower stiffness values result in higher contribution from GB deformation.

As the aggregate begins to yield plastically at the crystal scale, the relative contribution of GB deformation decreases to approximately 15%. However, with continued loading and the initiation of interface softening, the GB contribution rises again, ultimately reaching 30% just prior to macroscopic failure. The experimental results [43] have shown a similar deformation contribution, where the GBS accounts for 20% of the total deformation.
<figure><img src="image_28.png" /><figcaption>Fig. 25. Contribution of GBS to axial deformation in Model A.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 26. Evolution of the volumetric and deviatoric components of the interface integral during loading.</figcaption></figure>

The evolution of the volumetric (spherical) and shear (deviatoric) components of the interface tensor is presented in Fig. 26. Under compressive loading, the deviatoric term dominates, reflecting the prevalence of sliding over opening. The overall response is nonlinear and reflects the cumulative effect of grain geometry and loading path.

### 4.5. Interplay of crystal plasticity and grain boundary sliding

To characterize the competition between intragranular slip and intergranular deformation, we introduce the dimensionless ratio $\chi$, defined as the ratio of the critical resolved shear stress for dodecahedral slip systems, $\tau_{CRSS}$, to the peak shear strength of the cohesive zone, $\tau_c^s$:

$$ \chi = \frac{\tau_{CRSS}}{\tau_c^s} \qquad (50) $$

A higher value of $\chi$ indicates that GBS is favored, as the resistance to interface separation is lower compared to the stress required to initiate slip within the grains. Conversely, a lower $\chi$ implies that intragranular slip (CSP) is more easily activated than interface decohesion. To explore this interplay, we perform a parametric study varying $\tau_{CRSS}$ for both dodecahedral and cubic slip systems,
<caption>**Table 6**</caption>

Critical resolved shear stress, cohesive peak strength, and resulting $\chi$ values.

<table><thead><tr><th>Model</th><th>$g_0^{do}$ (MPa)</th><th>$g_0^{cub}$ (MPa)</th><th>$\tau_c^s$ (MPa)</th><th>$\chi$</th></tr></thead><tbody><tr><td>Model A</td><td>3.5</td><td>17.5</td><td>5</td><td>0.7</td></tr><tr><td>cp1</td><td>7</td><td>35</td><td>5</td><td>1.4</td></tr><tr><td>cp2</td><td>1.75</td><td>8.75</td><td>5</td><td>0.35</td></tr><tr><td>cp3</td><td>17.5</td><td>87.5</td><td>5</td><td>3.5</td></tr><tr><td>cp4</td><td>3.5</td><td>3.5</td><td>5</td><td>0.7</td></tr></tbody></table>

<figure><img src="image_30.png" /><figcaption>**Fig. 27.** Macroscopic stress-strain curves for various $\chi$ values.</figcaption></figure>

<caption>**Table 7**</caption>

Cohesive peak strength and resulting $\chi$ values.

<table><thead><tr><th>Model</th><th>$\tau_c^s$ (MPa)</th><th>$g_0^{do}$ (MPa)</th><th>$\chi$</th></tr></thead><tbody><tr><td>Model A</td><td>5</td><td>3.5</td><td>0.7</td></tr><tr><td>S1</td><td>50 x 10<sup>3</sup></td><td>3.5</td><td>7 x 10<sup>-4</sup></td></tr><tr><td>S2</td><td>2.5</td><td>3.5</td><td>1.4</td></tr><tr><td>S3</td><td>40</td><td>3.5</td><td>0.0875</td></tr><tr><td>S4</td><td>50</td><td>3.5</td><td>0.07</td></tr><tr><td>S5</td><td>3.5</td><td>3.5</td><td>1</td></tr><tr><td>S6</td><td>1.75</td><td>3.5</td><td>2</td></tr></tbody></table>

as reported in Table 6 (we maintain the same ratio between $g_0^{do}$ and $g_0^{cub}$, unless stated otherwise). The macroscopic stress-strain responses and the evolution of GBS contributions are presented in Figs. 27 and 28.

The results indicate that increasing $\tau_{CRSS}$, and consequently the dimensionless ratio $\chi$ (as exemplified by case cp1), delays the onset of intragranular plasticity and increases the macroscopic yield stress. This behavior is expected, as higher external stresses are required to activate slip within the grains. As a result, when crystal plasticity is suppressed, GBS emerges as the dominant deformation mechanism, leading to a sharp increase in its contribution to the overall axial strain. This shift in deformation mode accelerates interfacial damage accumulation and promotes earlier macroscopic failure.

In contrast, for lower $\chi$ values (e.g., cp2), plastic slip is readily activated, distributing deformation within the grains. In cp4, where $\tau_{CRSS}$ is identical for both dodecahedral and cubic systems, all slip systems are activated at very low strain (< 0.005%), resulting in early yielding, considerably reduced GBS contribution, extended strain hardening, and postponed failure.

These trends highlight the critical role of slip system resistance in determining the balance between GB and intragranular deformation mechanisms. Lower slip resistance favors distributed plasticity and delayed failure, while higher slip resistance shifts deformation toward GBS, promoting localization and premature failure.

The effect of cohesive strength is investigated by varying $\tau_c^s$, as summarized in Table 7. The macroscopic response and GBS contributions are depicted in Fig. 29 and Fig. 30, respectively.
<figure><img src="image_31.png" /><figcaption>Fig. 28. GBS sliding contribution to axial deformation for different $\chi$ values.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Fig. 29. Macroscopic stress–strain curves for various cohesive strengths.</figcaption></figure>

The initial slope of the stress–strain curves remains nearly identical across all cases, indicating that early elastic behavior is unaffected by the cohesive strength. Divergence in the GBS contribution emerges only once plasticity is triggered within the grains. Notably, for S2 and S6, the suppression of GBS at the onset of plasticity is minimized, resulting in earlier and more substantial intergranular deformation, which accelerates softening and failure. As the cohesive strength increases beyond a threshold, both the macroscopic response and GBS contribution plateau, indicating a saturation effect (cases S1, S3, S4).

Fig. 31 presents the evolution of the volumetric and deviatoric components of GBS. The shear (deviatoric) component consistently dominates, reflecting the prevalence of sliding mechanisms under compressive loading, while the volumetric (opening) contribution remains comparatively modest.
<figure><img src="image_33.png" /><figcaption>Fig. 30. GBS sliding contribution to axial deformation for different cohesive strengths.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 31. Deviatoric and volumetric components of GBS contribution.</figcaption></figure>

## 4.6. Comparison with different grain sizes

The effect of grain size on the mechanical response of polycrystalline halite is examined by comparing three numerical models with varying average grain diameters: Model A ($\Phi = 3.6$ mm, 333 grains, reference), Model B ($\Phi = 5.6$ mm, 89 grains, Fig. 32), and Model C ($\Phi = 2.2$ mm, 1124 grains, Fig. 33). All material parameters remain the same as prescribed in the Tables 3, 4, and 5 for all three simulations. For reference, all three models exhibit nearly identical macroscopic responses in the absence of CZMs at GBs. Therefore, in this section we mainly focus on the influence of GBS as a function of grain size.

Grain size is a fundamental parameter influencing the mechanical behavior of polycrystalline materials. Traditionally, the Hall-Petch effect describes the yield strength increase with decreasing grain size, attributed to the accumulation of dislocations at GBs which impede plastic deformation. However, at very small grain sizes, this trend often reverses, leading to the inverse Hall-Petch effect, where finer grains result in a reduction of yield strength. This phenomenon is commonly linked to the increased
<figure><img src="image_35.png" /><figcaption>Fig. 32. Model B: average grain diameter 5 mm.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. 33. Model C: average grain diameter 2.2 mm.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Fig. 34. Macroscopic stress–strain response for different grain sizes.</figcaption></figure>

dominance of GB deformation mechanisms. In the present study, the crystal plasticity model used does not incorporate size-
dependent intragranular hardening mechanisms (such as geometrically necessary dislocations), focusing instead on the role of GBs
represented by the pre-described rate-dependent CZM. As a result, the deformation and failure responses are strongly governed by
GBS leading to decohesion. This framework effectively captures the transition from Hall-Petch strengthening to inverse behavior
with decreasing grain size, highlighting the dominant role of GB weakening in ultrafine-grained microstructures.

The macroscopic stress–strain response of the polycrystalline models with varying grain sizes (Fig. 34) reveals a complex interplay
between GB behavior and intragranular deformation. Model B, characterized by grains with the largest average diameter, exhibits
the highest initial stiffness at early loading stage. Higher stress levels are also sustained throughout the loading process, which may
be attributed to the reduced compliance associated with the low GB density. Although Model B demonstrates an initial strength
advantage, it exhibits the earliest onset of softening and failure at approximately 1.5% of total axial strain (see Figs. 35 and 36).
This occurs despite the global contribution of GBS to the overall deformation being relatively modest, approximately 10%, compared
to the other two models (Fig. 37). These observations suggest that the comparatively low GB density in Model B constrains the
localization of damage-inducing interfacial deformation mechanisms, thereby precipitating earlier failure.
<figure><img src="image_38.png" /><figcaption>Fig. 35. Von Mises stress at failure, Model B, $[F_{zz} - I] = -0.014$.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. 36. Von Mises stress at failure, Model C, $[F_{zz} - I] = -0.036$.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Fig. 37. Contribution of GBS to axial deformation for different grain sizes.</figcaption></figure>

Conversely, Model C, with the smallest grains and thus the highest GB area, shows a reduced initial stiffness and stresses, reflecting the increased compliance through GBS and opening (Fig. 37). The greater interfacial surface area allows for more distributed damage-inducing deformation, which delays failure. Importantly, the absence of size-dependent intragranular hardening mechanisms, such as those related to geometrically necessary dislocations, means that the grain size effect arises primarily from GB damage and deformation processes which account to approximately 50% of the total axial deformation, for this model.

These results are consistent with the experimental findings reported by [10], where rock salt specimens with larger grain sizes exhibited deformation dominated by intragranular crystal plasticity (accounting for approximately 90% of the total strain), while specimens with finer grains showed a more balanced response, with GBS contributing up to 50% of the total deformation on their synthetic rock salt samples.

## 5. Conclusion

This work has established a modeling framework to elucidate the viscoplastic deformation mechanisms in polycrystalline rock salt subjected to compressive loading. The approach couples a power-law crystal plasticity model, which resolves the activity of crystallographic slip systems within grains, with a rate-dependent CZM that captures GBS and decohesion. Calibration against
experimental data on halite demonstrates the model's ability to reproduce essential features of the macroscopic response, including the onset of GB-driven deformation and the transition to failure.

The results underscore the decisive influence of gb behavior on the overall mechanical response of rock salt, particularly regarding the interplay between intragranular slip and intergranular separation. Grain size, in particular, modulates both strength and ductility, with finer grains increasing the prevalence of deformation at GBs.

Although the current simulations are focused on relatively short-term loading scenarios, the formulation can be extended to investigate long-term creep by calibrating the rate-dependent parameters using experimental creep data. This micro-mechanical modeling approach provides a mechanistic alternative to traditional macroscopic creep laws by linking grain-scale deformation processes such as crystallographic slip and GBS to the observed time-dependent macroscopic response. As such, the framework offers significant potential for future studies on the long-term performance of rock salt in geotechnical and energy storage applications. While explicit resolution of microstructural interfaces yields valuable insight into local failure processes, its computational expense necessitates the development of homogenization strategies for applications at engineering scale. The present framework thus provides a robust foundation for bridging microscale mechanisms and macroscale performance, offering predictive capability for the design and assessment of salt-based geological storage systems under a range of loading conditions.

### CRediT authorship contribution statement

**Nour Habib:** Writing – review & editing, Writing – original draft, Conceptualization. **Saber El Arem:** Writing – review & editing, Supervision. **Amine Ammar:** Writing – review & editing, Supervision.

### Funding

Nour Habib reports financial support was provided by Angers Loire Métropole, France.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

We would like to thank our colleagues involved in the ANR RockStorHy project for providing the experimental data that supported this study. Financial support from Angers Loire Métropole, France is also gratefully acknowledged.

### Data availability

Data will be made available on request.

### References

[1] P.E. Senseny, F.D. Hansen, J.E. Russell, N.L. Carter, J.W. Handin, Mechanical behaviour of rock salt: phenomenology and micromechanisms, in: International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 29, (4) Elsevier, 1992, pp. 363–378.

[2] W. Skrotzki, O. Frommeyer, P. Haasen, Plasticity of polycrystalline ionic solids, Phys. Status Solidi (A) 66 (1) (1981) 219–228.

[3] N.L. Carter, S.T. Horseman, J.E. Russell, J. Handin, Rheology of rocksalt, J. Struct. Geol. 15 (9–10) (1993) 1257–1271.

[4] D. Picard, A. Dimanov, J.L. Raphanel, Plastic behavior of halite single-crystals at different temperatures and strain rates: New insights from in-situ experiments and full field measures, Mater. Sci. Eng.: A 732 (2018) 284–297.

[5] Peter Anthonie Fokker, The behaviour of salt and salt caverns, 1995.

[6] Neville L. Carter, Hugh C. Heard, Temperature and rate dependent deformation of halite, Am. J. Sci. 269 (3) (1970) 193–249.

[7] Mathieu Bourcier, Alexandre Dimanov, Eva Héripré, Jean Raphanel, Michel Bornert, Guillaume Desbois, Full field investigation of salt deformation at room temperature: cooperation of crystal plasticity and grain sliding, in: Mechanical Behavior of Salt VII, 2012.

[8] Ababacar Gaye, Analyse multiéchelle des mécanismes de déformation du sel gemme par mesures de champs surfaciques et volumiques (Ph.D. thesis), Paris Est, 2015.

[9] Jihui Ding, Grain boundary processes, anelasticity, and test of the effective stress law for semibrittle deformation of synthetic salt-rocks (Ph.D. thesis), Texas A&M University, 2019.

[10] Mathieu Bourcier, Michel Bornert, Alexandre Dimanov, Eva Héripré, JL Raphanel, Multiscale experimental investigation of crystal plasticity and grain boundary sliding in synthetic halite using digital image correlation, J. Geophys. Res.: Solid Earth 118 (2) (2013) 511–526.

[11] Pierre Berest, M. Karimi-Jafari, B. Brouard, Comportement mécanique à très long terme des mines et cavernes dans le sel gemme: loi de norton-hoff ou loi de lemaître? Rev. Française Géotechnique (124) (2008) 45–59.

[12] C.S. Desai, D. Zhang, Viscoplastic model for geologic materials with generalized flow rule, Int. J. Numer. Anal. Methods Geomech. 11 (6) (1987) 603–620.

[13] M. Aubertin, D.E. Gill, B. Ladanyi, An internal variable model for the creep of rocksalt, Rock Mech. Rock Eng. 24 (2) (1991) 81–97.

[14] Curt A. Bronkhorst, Surya R. Kalidindi, Lallit Anand, Polycrystalline plasticity and the evolution of crystallographic texture in FCC metals, Philos. Trans. R. Soc. Lond. Ser. A: Phys. Eng. Sci. 341 (1662) (1992) 443–477.

[15] Andrey Jivkov, Nicholas Stevens, Thomas Marrow, A three-dimensional computational model for intergranular cracking, Comput. Mater. Sci. 38 (2006) 442–453.
[16] L.P. Evers, D.M. Parks, W.A.M. Brekelmans, MGD1100 Geers, Crystal plasticity model with enhanced hardening by geometrically necessary dislocation accumulation, J. Mech. Phys. Solids 50 (11) (2002) 2403-2424.

[17] Lallit Anand, Single-crystal elasto-viscoplasticity: application to texture evolution in polycrystalline metals at large strains, Comput. Methods Appl. Mech. Engrg. 193 (48-51) (2004) 5359-5383.

[18] Y.J. Wei, Lallit Anand, Grain-boundary sliding and separation in polycrystalline metals: application to nanocrystalline fcc metals, J. Mech. Phys. Solids 52 (11) (2004) 2587-2616.

[19] Masayuki Kamaya, Yoshihiro Kawamura, Takayuki Kitamura, Three-dimensional local stress analysis on grain boundaries in polycrystalline material, Int. J. Solids Struct. 44 (10) (2007) 3267-3277.

[20] Andrey Musienko, Georges Cailletaud, Simulation of inter-and transgranular crack propagation in polycrystalline aggregates due to stress corrosion cracking, Acta Mater. 57 (13) (2009) 3840-3855.

[21] Igor Simonovski, Leon Cizelj, Cohesive zone modeling of intergranular cracking in polycrystalline aggregates, Nucl. Eng. Des. 283 (2015) 139-147.

[22] Nicolò Grilli, Edmund Tarleton, Alan C.F. Cocks, Coupling a discrete twin model with cohesive elements to understand twin-induced fracture, Int. J. Fract. 227 (2021) 173-192.

[23] Surya Kalidindi, Curt Bronkhorst, L. Anand, Crystallographic texture evolution in bulk deformation processing of metal, 1990, pp. 281-286.

[24] C. Inglis, Stress in a plate due to the presence of sharp corners and cracks, Trans. Roy. Inst. Nav. Archit. 60 (1913) 219-241.

[25] Alan Arnold Griffith, T. Geoffrey, V.I. Ingram, The phenomena of rupture and flow in solids, Philos. Trans. R Soc. Lond. Ser. A 221 (582-593) (1921) 163-198.

[26] Donald S. Dugdale, Yielding of steel sheets containing slits, J. Mech. Phys. Solids 8 (2) (1960) 100-104.

[27] Grigory Isaakovich Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, Adv. Appl. Mech. 7 (1962) 55-129.

[28] Karl-Heinz Schwalbe, Ingo Scheider, Alfred Cornec, Guidelines for applying cohesive models to the damage behaviour of engineering materials and structures, Springer Science & Business Media, 2012.

[29] David H. Allen, Chad R. Searcy, A micromechanical model for a viscoelastic cohesive zone, Int. J. Fract. 107 (2001) 159-176.

[30] Arne Hillerborg, Mats Modéer, P.-E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, Cem. Concr. Res. 6 (6) (1976) 773-781.

[31] Viggo Tvergaard, Effect of fibre debonding in a whisker-reinforced metal, Mater. Sci. Eng.: A 125 (2) (1990) 203-213.

[32] A. Needleman, An analysis of tensile decohesion along an interface, J. Mech. Phys. Solids 38 (3) (1990) 289-324.

[33] X.-P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (9) (1994) 1397-1434, URL https://www.sciencedirect.com/science/article/pii/0022509694900035.

[34] Philippe H. Geubelle, Jeffrey S. Baylor, Impact-induced delamination of composites: a 2D simulation, Compos. Part B: Eng. 29 (5) (1998) 589-602.

[35] G. Alfano, M. Crisfield, Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues, Internat. J. Numer. Methods Engrg. 50 (7) (2001) 1701-1736.

[36] Y. Mi, M.A. Crisfield, H.B. Hellweg, G.A.O. Davies, J. Owen, Finite element method and progressive failure modelling of composite structures, Comput. Plast.: Fundam. Appl. Part 1 (1997) 239-254.

[37] Pedro P. Camanho, Carlos G. Dávila, Mixed-mode decohesion finite elements for the simulation of delamination in composite materials, Technical Report, 2002.

[38] Zdeněk P. Bažant, Yuan-Neng Li, Cohesive crack with rate-dependent opening and viscoelasticity: I. mathematical model and scaling, Int. J. Fract. 86 (1997) 247-265.

[39] Alberto Corigliano, Michele Ricci, Rate-dependent interface models: formulation and numerical applications, Int. J. Solids Struct. 38 (4) (2001) 547-576.

[40] Alberto Corigliano, Stefano Mariani, Anna Pandolfi, Numerical modeling of rate-dependent debonding processes in composites, Compos. Struct. 61 (1-2) (2003) 39-50.

[41] G.D. Seidel, D.H. Allen, K.L.E. Helms, S.E. Groves, A model for predicting the evolution of damage in viscoelastic particle-reinforced composites, Mech. Mater. 37 (1) (2005) 163-178.

[42] I. Aliguer, I. Carol, Efficient numerical integration of perzyna viscoplasticity, with application to rock slope stability using zero-thickness interface elements, in: ISRM EUROCK, ISRM, 2014, pp. ISRM-EUROCK.

[43] Xinjie Li, Alexandre Dimanov, Michel Bornert, Hakim Gharbi, Simon Hallais, Activation of deformation mechanisms in natural salt rock: A microscale investigation using digital image correlation, in: Mechanical Behavior of Salt XI: Proceedings of the 11th International Conference on the Mechanical Behavior of Salt, CRC Press / Taylor & Francis, London, 2025.

[44] Juan C. Simo, Thomas J.R. Hughes, Computational inelasticity, vol. 7, Springer Science & Business Media, 2006.

[45] A. Rovinelli, M.C. Messner, T.-L. Sham, Initial microstructural model for creep-fatigue damage in grade 91 steel, Technical report, Argonne National Lab. (ANL), Argonne, IL (United States), 2020.

[46] James Diebel, et al., Representing attitude: Euler angles, unit quaternions, and rotation vectors, Matrix 58 (15-16) (2006) 1-35.

[47] Guillaume Giudicelli, Alexander Lindsay, Logan Harbour, Casey Icenhour, Mengnan Li, Joshua E. Hansel, Peter German, Patrick Behne, Oana Marin, Roy H. Stogner, Jason M. Miller, Daniel Schwen, Yaqi Wang, Lynn Munday, Sebastian Schunert, Benjamin W. Spencer, Dewen Yushu, Antonio Recuero, Zachary M. Prince, Max Nezdyur, Tianchen Hu, Yinbin Miao, Yeon Sang Jung, Christopher Matthews, April Novak, Brandon Langley, Timothy Truster, Nuno Nobre, Brian Alger, David Andrš, Fande Kong, Robert Carlsen, Andrew E. Slaughter, John W. Peterson, Derek Gaston, Cody Permann, 3.0 - MOOSE: Enabling massively parallel multiphysics simulations, SoftwareX 26 (2024) 101690, URL https://www.sciencedirect.com/science/article/pii/S235271102400061X.

[48] Philippe Boch, Jean-Claude Ni, et al., Ceramic materials: Processes, properties, and applications, vol. 98, John Wiley & Sons, 2010.

[49] Romain Quey, P.R. Dawson, Fabrice Barbe, Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing, Comput. Methods Appl. Mech. Engrg. 200 (17-20) (2011) 1729-1745.