The linear elastic behavior of crushed salt is meant to play a minor role in the compaction process. However, its formulation is essential for the computational framework. The increase in stiffness with ongoing compaction is described by a generalized Hook’s law (Equation 7). Young’s modulus and Poisson’s ratio are applied as elastic constants. The isotropic linear elastic model is formulated in combination with a porosity dependent evolution of Young’s modulus (Olivella et al., 2023). By applying Equation 8, the Young’s modulus is increasing by decreasing porosity, simulating the material’s stiffness increase with ongoing compaction. The change of Young’s modulus with the change of porosity ( d E / d φ ) is considered to be a constant parameter. In the case of decreasing Young’s modulus with increasing porosity, the minimum value for Young’s modulus limits its decrease. ε ˙ E L = C e σ ˙ (7) E = E 0 + φ − φ 0 d E d φ ≥ E min (8) where C e is the elastic compliance matrix, E is the Young’s modulus, E 0 is the reference value for Young’s modulus evolution, φ is the current porosity, φ 0 is the reference porosity, d E / d φ is the variation of Young’s Modulus with porosity and E min is the minimum value for Young’s modulus.

The fluid assisted diffusional transfer mechanism describes the humidity creep of crushed salt, depending on the influence of moisture and dominating in areas of low stresses and low temperatures. Dissolution of salt will take place in areas of high stress concentration; the salt will then migrates through the liquid phase and precipitates in areas of lower stress. The FADT model is based on extensive studies by Schutjens (1991), Spiers et al. (1986), Spiers et al. (1989) and Spiers and Brzesowsky (1993).

The model formulation is derived on basis of the idealized geometry (Figure 3) and the assumptions of dissolution of salt in grain contacts, diffusive flux of salt through the liquid phase driven by a chemical potential gradient and precipitation of salt in pore space. The detailed derivation of the strain rate formulation in Equation 9 is to be found in Olivella and Gens (2002). The strain rate is decomposed in a volumetric and a deviatoric part, labeled by v and d, respectively. Equation 10 and Equation 11 are called the volumetric and deviatoric viscosities, respecitvely. They compile the dependencies on temperature T , solid volume d 0 3 , the material parameter B (Equation 12) and liquid saturation S l . The auxilliary functions g F A D T v (Equation 13) and g F A D T d , (Equation 14) as well as the functions g and f are geometry dependent. The functions g and f are derived directly from the idealized geometry in Figure 3, where g is the relative stress concentration and f express the relative pore size (Olivella, 1994). ε ˙ F A D T = 1 2 η F A D T d σ ′ − p ′ I + 1 3 η F A D T v p ′ I (9) 1 η F A D T v = 16 B T S l d 0 3 g F A D T v e (10) 1 2 η F A D T d = 16 B T S l d 0 3 g F A D T d e (11) B T = A B R T exp − Q B R T (12) g F A D T v e = g 2 1 + e (13) g F A D T d e = 3 g 2 e 3 / 2 1 + e (14) g = 1 1 − f 2 = d 2 x 2 (15) f = 2 e 3 1 − e 3 / 2 = s 2 d (16) where I is the identity matrix, S l is liquid saturation, d 0 is characteristic grain size, T is temperature, A B is a pre-exponential parameter, Q B is the activation energy, R is the gas constant and e the void ratio.

The dislocation creep mechanism captures deformations like dislocation glide and climb. It refers to the intracrystalline mechanisms (Olivella, 1994). These mechanisms can be described by power law terms and therefore are grouped here. The model is based on the rock salt power law and combined with a geometrical derivation of a volumetric strain rate and a deviatoric strain rate formulation. It is generalized by the use of a viscoplastic approach including a flow rule and a viscosity parameter. G and F in Equation 18 are formulated as functions of stress invariants. The viscosities in Equations 21, 22 compile geometrical and material properties. The fuctions g and f are the same as defined in Equations 15, 16 (Olivella and Gens, 2002). ε ˙ D C = 1 η D C d Φ D C F D C ∂ G D C ∂ σ ′ (17) F D C = G D C = q D C 2 + − p α p 2 (18) Φ D C F D C = F D C n (19) α p = η D C v η D C d 1 1 + n (20) 1 η D C v = A T g D C v e (21) 1 η D C d = A T g D C d e (22) g D C v e = 3 g − 1 n f (23) g D C d e = 1 + g + g 2 3 n − 1 2 g + 1 3 f + 1 g (24) A T = A A ⁡ exp   − Q A R T (25) where F D C is the stress function, G D C is the flow rule, Φ D C is a scalar function, n is the scaling exponent of the rock salt power law, g D C v and g D C d are non-linear functions of void ratio, A A is a pre-exponential parameter and Q A is the activation energy. The deviatoric stress q D C is defined as: q D C = 3 J 2 0.5 (26) with J 2 the second invariant of the deviatoric part of the stress tensor.

The term for grain rearrangement was added for simulating the irreversible processes of grain displacement during compaction taking place mainly at high porosities and under fast loads. This part is added to describe the behaviour of crushed salt as loose aggregate. In its loose state crushed salt is assumed to behave sand-like. Such material behaviour is described by a critical state approach which is based on a yield surface. Due to expansion of the yield surface, bond creation and densification of the crushed salt material is modelled. An important feature is that the hardening of the crushed salt is mainly driven by creep deformations (Equations 27–33) (Olivella and Gens, 2002). The basic equations are the following (Equations 27–33) (Olivella et al., 2023): ε ˙ G R = Γ Φ G R F G R ∂ G G R ∂ σ (27) Φ F G R = F G R m (28) G = F = q G R 2 − δ 2 p 0 p ′ − p ′ 2 (29) Γ = Γ 0 ⁡ exp   − Q G R R T (30) δ = 6 ⁡ sin   ϕ ′ 3 − sin   ϕ ′ (31) d p 0 = p 0 1 + e Χ d ε v o l (32) where Φ G R is the stress function, F G R is the viscoplastic yield function,   are Macauley brackets, G G R is the flow potential, Γ is the fluidity, m is the stress exponent, δ is the slope of the critical state line, p 0 is a hardening parameter, Γ 0 is the initial value of the fluidity, Q G R is the activation energy, ε v is the volumetric strain, Χ is a parameter in the hardening law, and the invariant q G R is described in terms of octahedral stress: q G R = 3 2 τ o c t = 1 2 σ x − σ y 2 + σ y + σ z 2 + σ z + σ x 2 + 6 τ x y 2 + τ y z 2 + τ z x 2 (33)

For developing a crushed salt constitutive model, it is a common approach to derive an idealized geometry and deformation mechanisms from microstructural observations (Olivella, 1994). This approach lacks some generality since grains deform in a non-regular way during the compaction and a constantly regular arrangement cannot be expected. Deformation mechanisms like dislocation creep (grain breakage, plastic deformations) and humidity creep (dissolution and precipitation) dominate under different conditions and act not uniform on a grain (Hansen et al., 2014).

The creep models in CODE_BRIGHT for FADT and DC are based on a geometry dependent formulation, expressed in terms of void ratio (Olivella, 1994). The functions g and f (Equations 15, 16) are related to the characteristic sizes derived from the idealized geometry (Figure 3) and yield the auxiliary functions g F A D T v , g F A D T d , g D C v and g D C d . These functions can be plotted for their dependence on void ratio as shown in Figure 4.

Within the KOMPASS projects, benchmark calculations of long-term triaxial compaction tests on crushed salt were performed using the model implemented in CODE_BRIGHT and presented above. The reproduction of compaction was not satisfactory for both volumetric and deviatoric deformation, as well as their deformation rates (Friedenberg et al., 2023; Friedenberg and Olivella, 2024). Especially, for the creep contribution, which is the most important process in the long-term for porosity and permeability reduction, the accordance between numerical results and experimental data is weak. The example is shown in Section 4.1. Based on this work an idea for the modification of the crushed salt model in CODE_BRIGHT is developed aiming to improve the numerical simulation of crushed salt compaction.

The modification approach addresses the prescription of an idealized geometry for the creep models. To give some flexibility in the handling of the geometrical basis, it is proposed to formulate g D C v and g D C d in terms of phenomenological functions. The combination of a microstructural basis with phenomenological functions allows the inclusion of both microstructural observations and experimental experiences from triaxial compaction tests.

Since the long-term compaction tests considered in the recent work of the KOMPASS projects were executed on dry crushed salt, the FADT mechanism is assumed to have no or at least a minor effect on the compaction behavior. Therefore, this modification approach focuses on the dislocation creep mechanism.

From the shape of g D C v and g D C d as implemented in CODE_BRIGHT, an exponential formulation was found to build the basis for the phenomenological functions: g D C v = exp a * e b − 1 (34) g D C d = exp c * e d (35) where e is the void ratio and a , b , c and d are parameters.

The parameters a , b , c and d do not have a specific physical meaning and cannot be measured during experiments. However, experimental data from standard triaxial compaction tests can be used to calibrate the proposed phenomenological functions and determine the parameters a , b , c and d . Therefore, measured volumetric and deviatoric strain rates must be considered, respectively. An example for the calibration process is shown in the following section.

Initially, the available model was calibrated against the experimental data. The simulation was executed in a THM-coupled approach following the stress history given in Figure 5. Thermal processes are based on Fourier’s law, for the hydraulics one phase flow with a constant gas pressure of P g = 0.1 MPa is assumed and the mechanics are based on the crushed salt model described before. Based on an existing parameter set, the calibration process was performed in combination with a parameter sensitivity to identify the most suitable parameter combination for the simulation of the TUC-V2 test. Detailed information on the sensitivity and calibration procedure can be found in Czaikowski et al. (2020) and Friedenberg et al. (2024). Table 1 gives an overview of the parameters elaborated during calibration. The majority of the parameters are physically based and can be derived from the material itself. Crushed salt is studied for several decades. The significant parameters in the calibration process are part of the dislocation creep (parameter A A ) and the grain rearrangement (parameters Γ 0 , p 0 , Χ ) models. Due to the lately definition of the KOMPASS reference material and the consequent limited experimental database, the derivation of material specific model parameters is still in progress. Therefore, the calibration process was applied with the aim to keep the parameters in a realistic and meaningful range. This means for parameters with physical meaning to vary them within observed ranges for crushed salt material (e.g., the pre-exponential parameter A A describes the creep ability, therefore values for different salt formations are available) and for the other parameters to vary them with respect to the conceptual idea and their dependences. The parameter dependences of various combinations are investigated by performing the sensitivity analysis.

The comparison in Figure 6 shows clear differences between the experimental data and the numerical results. The volumetric compaction is strongly underestimated by the simulation, especially the compaction under the constant mean stress of 20 MPa. The final values for volumetric compaction differ by 3.5%. In contrast, the deviatoric strain is overestimated. The model response to the increase and decrease of deviatoric stress is higher compared to the experimental data. However, the final values for deviatoric strain just differ by 1%.

An improvement due to a broader calibration with reasonable parameter ranges could not be achieved, thus the modification approach was developed.

The modified approach was developed in the way that Equations 34, 35 can be calibrated against experimental data by using the experimental strain rate data transformed to be comparable with the function values of the g D C v and g D C d functions. As a first attempt the volumetric equation g D C v was calibrated. Considering the experimental data in Figure 7A the trend can be approximated by the choice of parameters a and b . For the deviatoric function g D C d calibration was not performed in this first step and default values are taken (Figure 7B). The corresponding values for the parameters a , b , c and d are shown in Table 2.

A new simulation with the modified functions and parameters as shown in Table 2 was performed. Due to the calibration of the modified function g D C v against the experimental data an overall improvement of the numerical results is achieved (Figure 8). The evolution of volumetric strain is reproduced satisfactorily. Now the compaction due to the high mean stress of 20 MPa is captured and the whole volumetric compaction evolution is simulated well. The reproduction of the deviatoric strain evolution shows also improvements. The final value for deviatoric strain is met well. However, the response in deviatoric strain to deviatoric load changes is still overestimated by the model.

A further quantification can be realized by considering the strain rates. Figure 9A show the volumetric strain rates for both simulations, the initial and the modified, compared to the experimental data. The modified simulation shows an improved accordance for the volumetric strain rates during the rapid load changes in the beginning of the test. Up to 300 days the numerical strain rates hardly differ from each other. Figure 9B show the deviatoric strain rates. Only slight differences in the numerical results can be observed for the first 300 days, along with a higher accordance of the modified rates with the data. In both simulations, the strain rates are strongly overestimated at the start of the test. From 300 days on no significant difference between the trend of the numerical results is identifiable.

In general, the results show an improvement for the simulation of crushed salt compaction by using the modified functions g D C v and g D C d . Already with the calibration of the volumetric part, the volumetric compaction is simulated in an adequate manner and the deviatoric strain is reasonably reproduced in its magnitude.

The calibration for the deviatoric part provides some open questions. The experimental data for deviatoric strain shows a large scattering over the function values, leading to difficulties in the calibration process (Figure 7B). The data trend is not straight forward as it is for the volumetric strain. Reasons might be traced back to the execution of the test with short durations of deviatoric stress (10 days) and therefore a wide range of strain rates due to the rapid load changes.

Figure 10 presents simple one-factor-at-a-time variations of c and d . By choosing c = 200 (Figure 10A) a good calibration of the modified g D C d function on experimental data is achieved. However, the corresponding model response in Figure 10C shows a rapid increase in deviatoric strain and a complete overestimation. Setting d = 1.5 the function g D C d is calibrated against the lower bound of the experimental data (Figure 10B). Here, the model response also shows an overestimation of the deviatoric strain. All in all, the default values ( c = 50, d = 2) yield the best reproduction of the deviatoric strain in comparison to the measurements (Figure 10D).

In general, the model response to changes in deviatoric stress is strongly overestimated (Figures 10C, D). In former research work, deviatoric stress and strain represented a minor part. They predominantly focused on mean stress and isotropic compaction conditions (Kröhn et al., 2017). However, in recent crushed salt projects, the need for the investigation of deviatoric stress influence on the compaction was highlighted (Friedenberg et al., 2024).

The presented simulation results indicate an overall need to focus on the simulation of deviatoric stresses with the constitutive model available in CODE_BRIGHT. The consistent strain overestimation in response to the deviatoric load changes might be a general issue but is not part of this modification approach.

The Morris method (Morris, 1991) was used to determine which of the input parameters may have the largest effect within the simulation of the triaxial compaction test TUC-V2. The idea was to use a factor’s screening method to determine the weighted effect of changes in the parameter values. Therefore, the Morris method is based on a “one-factor-at-a-time” approach yielding the global sensitivity by performing a series of local changes at different points in a predefined parametric state-space (Saltelli et al., 2004).

The analysis was carried out by using the Julia programming language (Bezanson et al., 2017). Julia’s in-built capabilities for distributed computing were used to conduct the Morris sensitivity analysis with a simple Monte-Carlo sampling strategy on various workers in parallel. The sampling was performed within the deterministic predefined parametric state-space ( a in [100, 1,500], b in [3, 5], c in [50, 200] and d in [2, 5]) and the Morris means are calculated for each parameter. To derive the weighted effect of each parameter, the respective Morris means are divided by the sum of the means.

As indicators for the performance evaluation, the volumetric strain ( ε v o l ) and deviatoric strain ( ε d e v ) are chosen (Figure 11). The results of these analyses show a high influence of the exponent b on the volumetric strain followed by the pre-exponent parameter a . For the deviatoric strain, the evolution depends strongly on the exponent d and on the pre-exponential parameter c . There is a small influence of the two deviatoric parameter c and d on the volumetric strain and the two volumetric parameters on the deviatoric strain resulting from the formulation of the dislocation creep law (Equations 17–26).

In this analysis, the influence of the main parameters on the respective output quantity is investigated. As shown in the previous section with the Morris method, the volumetric strain output is primarily influenced by the parameters a and b (Equation 34) and the deviatoric strain output is primarily influenced by the parameters c and d (Equation 35).

The analysis was conducted by using an adaptive sparse-grid collocation method (see, e.g., Gates and Bittens (2015)) implemented by the open-source project DistributedSparseGrids.jl (Bittens and Gates, 2023) in the Julia programming language. Hereby, a hierarchical Lagrangian basis enables the adaptive refinement of the parametric model in areas with high effect on the output quantity.

The influence on the volumetric strain for a in [100, 1,500] and b in [2, 5] is shown in Figure 12A. Increasing the exponent b leads to a decrease in the volumetric strain (less compaction), thus the material response is stiffer. The pre-exponential parameter a has a small effect on the volumetric strain output, however, with increasing value of an increase in volumetric strain is observed.

The response of deviatoric strain for c in [50, 200] and d in [2, 3] is shown in Figure 12B. The plot shows nearly no change in deviatoric strain over a wide range of parameters. Only a small number of parameter combinations in the left corner of the plot leads to a sensitive response of the deviatoric strain. The boundaries for this sensitive area can be defined by c > 100 and d < 2.5.

Table 3 presents absolute values for volumetric and deviatoric strains for the interval boundary parameter combinations.