Error‐controlled implicit time integration of elasto‐visco‐plastic constitutive models for rock salt

The mechanical behavior of rock salt is rate‐dependent at different time scales. Using caverns in rock salt formations for renewable energy storage implies that the underground structures are subjected to both short‐term and long‐term loads, increasing robustness and flexibility requirements for numerical simulators used to assess the safety of such structures. So far, explicit time integration with model‐specific heuristics for time‐step size determination dominate in application studies. In this paper, the suitability of error‐controlled adaptive time‐integration schemes of the diagonally implicit Runge‐Kutta type is investigated in comparison with the Backward Euler and Crank‐Nicolson schemes when applied to the integration of typical elasto‐visco‐plastic constitutive models of rock salt. The comparison is made both for monotonic and for cyclic loads as well as taking account of thermo‐mechanical coupling. Analyses of the time‐integration errors and the time step‐size evolution show the suitability of the integration scheme for these material models. The automatic adjustment of the time‐step size was found to be robust across all material models and boundary conditions as well as for non‐isothermal situations for a single algorithmic parameter set.


INTRODUCTION
A wide range of time-integration schemes, both implicit and explicit, have been implemented into software packages for the numerical simulation quasi-static and dynamic problems in geotechnics. 1,2 Because a balance has to be sought between using shorter step sizes to increase accuracy and decreasing simulation times by using larger step sizes, 3,4 researchers have paid much attention to the numerical stability and accuracy of different integration schemes depending on the chosen step size. [5][6][7][8] The choice of a suitable time-integration scheme depends to a certain degree on the problem modelled. This article is concerned with the integration of complex inelastic constitutive models for rock salt.
The mechanical behavior of rock salt is rate-dependent, causing time-dependent stress and strain fields associated with primary, secondary, or tertiary creep or stress-relaxation processes developing at different time scales. 9 These processes temperature. [10][11][12][13][14] Time-integration schemes must ensure that the nonlinear and path-dependent stress-strain behavior of the material is captured with a desired accuracy over the course of a numerical simulation in order to provide meaningful results. For that purpose, several means to adjust the time-step size have been applied. In practical approaches, the maximum time-step size for creep problems in numerical simulations has been assessed by limiting the maximum increment of variables representing visco-plastic strain to a certain fraction of the total accumulated strain at all integration points of a system. Further criteria for the time step have generally been imposed, allowing the determination of the time-step change between successive time intervals. 15 Many of these criteria are material model-specific or are justified empirically rather than strictly mathematically and have been implemented in a range of software packages for the control of computational stability and accuracy. 16,17 Explicit algorithms pose stronger requirements in terms of a maximum time-step size to satisfy stability conditions. 18 As inertial effects become practically irrelevant during typical creep simulations, an estimate for the maximum creep time step for numerical accuracy can be determined as a material-specific characteristic time expressed as the ratio of the material viscosity to the shear modulus. 17 Depending on the constitutive formulation, the formulation of the criteria varies, such that for different material domains, additional minimum criteria need to be introduced. A universal integration scheme would, however, be preferable. Finally, also material failure is considered in modern constitutive models for geotechnical materials where hardening or softening behavior can be observerd. 13,[19][20][21] The time scale of various creep stages and failure observed in rock salt differ by orders of magnitude, requiring additional considerations for time adaptivity.
Quasi-static problems with significant nonlinearity and path dependence to be simulated for extended time periods can benefit from implicit integration approaches due to unconditional stability, the possibility for larger time steps, and equilibrium as well as constitutive iterations. Due to the different time scales involved in the process, eg, elastic deformation, primary as well as secondary creep, and finally plasticity as well as damage, an adaptive time integration scheme is desirable also for implicit approaches for many of the reasons outlined above for their explicit counterparts. [22][23][24] Cyclic boundary conditions varying on different time scales encountered in practice may only serve to emphasize this need further.
Time-integration approaches can be characterized by their local truncation error. This error and its propagation for a possibly quite large number of successive time steps is of interest when the analyst wants to establish a certain accuracy level of the simulation. None of the empirical methods listed above offer direct control over the error. Ideally, the engineer could specify a desired level of accuracy in terms of a well-defined error measure and the integration algorithm automatically adapts the time-step size accordingly.
One such class of methods that allows the specification an upper bound for the local truncation error and determines the appropriate time-step size automatically is the class of diagonally implicit Runge-Kutta (DIRK) schemes. 25 Its algorithmic structure in combination with the application of a ultilevel-Newton method for global (equilibrium) and local (constitutive) levels allow for an efficient integration into existing implicit codes. Moreover, the embedded step-size control techniques can be applied efficiently as the local error estimation is computed with negligible computational overhead. 26 The properties of DIRK methods have been thoroughly studied in recent years. Wide-ranging applications of DIRK methods include fluid dynamics, chemical reaction kinetics and mechanisms, applications in medical technology, semiconductor design, plasticity, structral analysis, neutron kinetics, porous media physics, gas transmission networks, and others, cf. 25,27 and references therein. Applications in solid mechanical simulations also extend to rate-dependent materials as studied in this paper. A description of DIRK algorithms for constitutive relations of Norton-type creep models can be found in Mahnken and Stein. 3 Two numerical examples are presented to demonstrate how an adaptive determination of step-sizes improves the accuracy of the solutions. Further application on the classical disc-with-hole problem with visco-elastic material models were studied in Ellsiepen and Hartmann. 26 These algorithms were then used for more complex engineering problems involving porous media, eg, an analysis of a car seat conducted with the software PANDAS. 28 Besides visco-elasticity, recently rate-dependent crystal plasticity simulations further demonstrate the benefits of DIRK methods in terms of accuracy and computational cost. 29 The numerical simulation of the cyclic operation of energy storage caverns in rock salt is characterized by a number of different time scales. 30 Whenever a load cycle occurs, a primary (transient) creep phase is initiated while secondary (steady-state) creep is always prevalent and dominates the long-term response. In the vicinity of the cavern wall, tertiary creep or rather plastic deformation may occur based on a rate-independent rheological submodel. Aside from that, the simulations may run over extended periods of time and include idling phases, such that overall a large number of time steps will be accumulated. The optimal time steps required in each period can differ by orders of magnitude and are 1110 are furthermore dependent on the current state of the material, eg, the acting stress, the presence of damage, or the current difficult to determine a priori without numerous trial simulations. Error-controlled time-integration schemes are thus very attractive for such applications.
The Lubby2 31 and Minkley 32 formulations capture salient features of the constitutive behavior of rock salt as described above and are frequently used in practice. Both have yet to be combined with error-controlled DIRK schemes. In this paper, we close this gap and integrate DIRK schemes with a multilevel-Newton method on the material point level for testing purposes. Testing and analyzing the errors of global and internal variables of both models is based on simple shear and triaxial compression tests under relevant stress conditions. This study is intended for prototyping the method prior to implementation into the open-source finite-element (FE) software OpenGeoSys, 33,34 which is currently based on backward Euler time integration in simulations of structures in rock salt. 35,36 This code is mainly used for multiphysics simulations, which is why we also test the time integration behavior in a thermo-mechanically coupled setting. Thermal strains as well as temperature-dependence of the visco-elastic response of the material are relevant in many subsurface applications in salt rock formations, such as nuclear waste disposal and energy storage.
The paper is organized as follows: Section 2 revisits the algorithms of diagonally implicit Runge Kutta schemes. Section 3 describes the investigated constitutive models for rock salt. Section 4 introduces the numerical implementation of the constitutive integration scheme. In Section 5, the performance of the time integration scheme for the two material models is analyzed using simple shear, triaxial, and cyclic loading tests results for isothermal and non-isothermal conditions. Finally, Section 6 presents some conclusions.

General concept
In a standard continuum mechanical setting, inelastic material models are usually defined by a local differential-algebraic equation (DAE) system. In a FE context, this DAE system has to be solved in each integration point for each global equilibrium iteration. Here, we summarize this local DAE system in a residual form by writing where the dependence on the driving total strain measure provided by the (time-dependent) displacement field is included via the explicit dependence of Equation For a numerical solution, a suitable time-discretization scheme is required. Most FE codes use explicit Runge-Kutta or (implicit) backward Euler schemes which have the main advantage of being easy to implement into software. Here, we employ a DIRK scheme as outlined in Ellsiepen and Hartmann. 26 Consider, for that purpose, a slight reformulation of Equation (1) to a general ODE Given a solution z n at time t n , the solution in the next time step t n+1 = t n + Δt n is formally approximated by introducing n s stages with stage rates . Z and a set of corresponding weight coefficients b as follows: These stage rates follow from the definition of the ODE (2) by substituting stage times T and stage solutions Z as follows: .
The stage times are determined based on a set of coefficients c as T = t n + c Δt n , while the stage solutions follow from By setting a n s i = b i and a i = 0 for i < , one arrives at computationally efficient so-called stiffly accurate SDIRK schemes which have the additional advantage that the final stage solution Z n s coincides with the sought solution z n+1 . 26 With these definitions, the stage solution follows directly from a starting value S and This allows expressing the state rate as The DAE system (1) can now be solved at each stage as follows: By comparing Equation (8) with Equation (1), it becomes clear that the implementation into codes relying hitherto on implicit backward Euler time integration is straightforward in that additional infrastructure is merely needed for managing (computational) stages in addition to (physical) states.
Equation (8) is solved by Newton-Raphson iterations, given an estimate of the strain increment from the current equilibrium iteration. Once these constitutive iterations have converged, the algorithmically consistent tangent  required for the global equilibrium Newton-Raphson iterations is obtained by standard procedure 26,35,[37][38][39] and the next equilibrium iteration provides a new estimate of the strain increment to the constitutive level. Since the present tests operate on material point level only, the equilibrium iterations are simply solving r eq = ( ) − ext = 0, where boundary conditions can be specified in terms of coordinates of either the stress tensor ext or the strain tensor . The consistent tangent  = d derived from this formulation and alluded to above is the one that would likewise be required for the assembly in an actual FE simulation.

Embedded methods for error estimation
A backward Euler scheme is, for example, a time-integration scheme of first order while the Crank-Nicolson scheme is of second order. The order p of a specific SDIRK method is determined by its n s stages. As demonstrated previously, 26 this allows for an efficient way of controlling the integration error. For that purpose, a second approximationẑ n+1 of the sought solution is calculated based on the already determined stage rates . Z and a second set of coefficientsb chosen to constitute a method of order p − 1 as follows: so that an error estimate can be defined based on an upper-bound theorem (see previous studies 26,40 for details) by comparing the results from the two schemes of order p and p − 1 as follows: In an FE setting, several error measures, usually for the primary nodal variables and for internal variables at the integration points, may be of interest. Here, as the implementation of the scheme is at the integration point level, we chose strain at the integration points to represent the primary variables. Internal variables are represented by the state vector z with its enties denoted by z k where the first entry is always stress, z 1 = , see Equation (1). The errors of primary and internal variables were calculated from distinct relations following 26,41 : have been used. This dimension dependence will be addressed by variable-specific tolerances in the sequel. Note in passing that relative tolerances could easily be included. 25,26,41 Based on the maximum of these error measures, e m = max(e , e k z ), a suitable time-step adaptation criterion can now be defined 26 as follows: with typical values of 0.8 ≤ safety ≤ 0.9, 0.2 ≤ min ≤ 0.5, and 2 ≤ max ≤ 3. 26,41 Another constraint on the time-step size exists because the time-integration algorithm was implemented such that it treats time points which are located at the transition between load steps as must points, ie, the time steps are decreased such that these points are hit precisely.

Ellsiepen's rule
In this study, we chose the following coefficients yielding a special case of SDIRK schemes commonly called Ellsiepen's rule. 26 This is an embedded scheme with two stages (n s = 2), where the main scheme is of order p = 2 and A-, S-, and L-stable, while the embedded scheme of orderp = 1 is A-and S-stable. A so-called Butcher tableau can be used to collect the coefficients c , b , and a i introduced above and reads for Ellsiepen's scheme where additionally, the coefficientsb of the embedded scheme have been included. Being of order 2 and not higher saves computational resources. 26 For details on the determination of the value of in the above schemes, see Kennedy and Carpenter. 25 In essence, this means that Equation (1), or rather Equation (8), is solved twice per time step, ie, in two stages, in a fully implicit manner. The updated state vector z is then determined from the stage rates based on the coefficients b (second-order scheme) andb (embedded first-order scheme) using Equation (3). The difference between both state vector updates allows the estimation of the error by means of Equation (10).
For comparison and based on the notational definitions made in Ellsiepen and Hartmann, 26 we give here the Butcher tableaus of the backward Euler and Crank-Nicolson schemes used later on as well. Note that nob coefficients need to be given as these schemes will be used without any embedded scheme as follows:

CONSTITUTIVE MODELS OF ROCK SALT
A variety of constitutive models of rock salt has been proposed in the past for different applications based on experimental and theoretical studies. 13,14,[42][43][44][45][46][47][48] Lubby2 and Minkley are constitutive models composed of Burgers models extended by nonlinear stress-dependent viscosities, and, in the case of Minkley, dilatancy boundaries. 31,32,49,50 Therefore, both models can represent transient and steady-state creep and Minkley's model also dilatant plastic deformation.

The modified Lubby2 model
The Lubby2 model 31,50 is very common in the geotechnical analysis of rock salt behavior. Here, it is slightly modified for general three-dimensional stress states and its numerical implementation into OpenGeoSys based on the Kelvin mapping. 51 The Lubby2 model is derived from the Burgers model, and as such, comprises a Maxwell element to represent stationary creep in series with a Kelvin element for transient creep, see Figure 1. This model can be described by the where D is the deviatoric stress, D is the deviatoric strain, e is the volume strain, K is the bulk modulus, G is the shear modulus, and is the viscosity. The subscript of K or M indicates the parameters attribution to either the Kelvin or the Maxwell element, respectively, as shown in Figure 1. The viscosities and the Kelvin shear modulus of the Lubby2 formulation are functions of the current stress state by means of an exponential dependence on the equivalent deviatoric stress eq where m a are material parameters characterizing the stress dependence and 0 is a reference stress (typically, 0 = 1 MPa).
For the non-isothermal simulations, additional temperature dependence is taken into account. First, the stress equation is modified to account for thermal strains based on the material's linear thermal expansion coefficient T , as follows: Second, the elastic properties were allowed to vary linearly with a temperature change with respect to a reference temperature T ref following Böttcher et al 36 and based on data by Sriapai et al 52 : where K M0 is the initial bulk modulus of Maxwell element, G M0 is the initial shear modulus of Maxwell element, and m KT and m GT are temperature dependent coefficients for K M and G M , respectively. The Maxwell viscosity describing the stationary creep rate varied based on an Arrhenius-type relationship in addition to its stress dependency 35 : where M0 is the initial viscosity modulus of Maxwell element, Q is the activation energy for creep, and R is the universal gas constant.

Simplified Minkley model
The Minkley model, commonly used for rock salt in situations where the dilatancy limit is exceeded in areas of practically relevant dimensions, is based on a visco-elasto-plastic rheological model as shown in Figure 2. The structure is similar to the Lubby2 model but extended by a plastic element allowing rate-independent non-associated plastic flow. The original formulation 32 uses a modified Mohr-Coulomb yield surface, while the current implementation in OpenGeoSys 35 rests on a standard Mohr-Coulomb model and related flow rule. As long as the stress-state remains within the yield/dilatancy boundary, ie, F < 0, the material behavior is visco-elastic (primary and secondary creep) and the solution method is analogous to that outlined for the Lubby2 model above. This procedure is also used to compute the trial stress with which the yield criterion is evaluated. Upon yielding, we seek the solution to .
where I is the identity tensor, e P is dilatant plastic volume strain, Peff is the plastic arc length, f S is the spherical part of tensor of G F , and f D is its deviatoric part. The fourth-order spherical and deviatoric projection tensors are introduced as where  s and  D are spherical and deviatoric projection tensors and I is the second-order identity tensor.
In contrast to the LUBBY2 model, the only parameter in the Burgers part of the model depending on the equivalent stress is the Maxwell viscosity with the following relationship: where m and n are material parameters. The cohesion can evolve according to linear and non-linear hardening/softening laws. 13,32,53,54 Here, we considered only simple linear hardening in the early stages of plastic deformation 35 as follows: where c is cohesion, c 0 is initial cohesion, and H 1 is first-order harding coefficient. Therefore, this model is termed "simplified Minkley." For the non-isothermal simulation, the stress equation was modified to account for thermal expansion The temperature dependence of the elastic properties followed Equation (18), and the Maxwell viscosity was extended by a temperature dependence analogous to Equation (19):

NUMERICAL IMPLEMENTATION OF THE TESTING SCHEME
Python is a high-level programming language with a range of libraries for scientific computing, suitable for the rapid prototyping and testing of code. Therefore, we implemented the DIRK scheme within a multilevel-Newton method in Python to facilitate testing prior to the more elaborate implementation into the full FE solver OpenGeoSys. 33 Taking advantage of the concept of object-oriented programming, classes were introduced for applying boundary conditions, material models, and time-integration methods. The present paper will test the Ellsiepen SDIRK time-integration scheme for the material models Lubby2 31 and Minkley 32 and compare the results with the traditional integration schemes backward Euler and Crank-Nicolson in terms of the obtained accuracy. The main architecture and simulation process is illustrated in Figure 3. Here, the error was obtained by comparison against an analytical solution. (b) Triaxial compression testing starting with the application of an isotropic confining stress followed by a strain-controlled compression in the z direction. (c) Cyclic triaxial compression testing, controlling the stress magnitudes in a range roughly representative of a point in the proximity of a cyclically operated gas storage cavern. 36 In cases (b) and (c), the reference results are obtained by numerical simulation with very small time steps.
The detailed description and parameter values of each test will be provided in the subsequent section. The comparison is expressed in terms of absolute and relative errors, eg, in the case of a shear strain, wheréx is the numerical result and x is the reference result (either analytical or obtained by a high-resolution numerical simulation). The focus of the subsequent tests does not lie on mathematical properties of the scheme itself but it is a practical testing for the targeted material models and loading conditions.

Simple shear
Simulations of simple shear tests with the Lubby2 and Minkley models were conducted to evaluate the accuracy of the different time integration schemes. Applied boundaries are shown in Figure 4A. During the entire simulation, normal strains are constrained to zero to enforce isochoric deformations. A 5-or 2-MPa shear stress jump is applied * and subsequently held for 15 or 1500 days in case of the Lubby2 and Minkley models, respectively. The difference in stress magnitude and creep duration stems from the different parameter values chosen here for both models which correspond to different types of rock salt. The material parameters used for Lubby2 are listed in Table 1, while those used in Minkley's model are listed in Table 2. The temperature-dependent parameters for both models are listed in Table 3. The references for the experimental data or parameter identification are given in the table captions. The unrealistically high cohesion value in the Minkley model ensures that no plastic flow sets in and Equation (33) applies. The absolute tolerance for the residuals of the local (constitutive) Newton iteration l (for the respective Euclidean vector norms) was set to 1.0 · 10 −14 , while the absolute tolerance for the global (equilibrium) Newton iteration g was set to 1.0 · 10 −12 . Since these simulations are performed on a material point level, the strains were chosen to represent the global variables taking the role of the displacements in an actual FE computation. The numerical results were compared to the analytical solution 35 given in Equation (33) to quantify the error.
*Numerically, the jump is realized via a linear ramp within 1.0 · 10 −8 days and 20 time steps.    where, x is the shear strain, x is the applied shear stress, and G M , M , G K , and M are stress-dependent material parameters as defined in Sections 3.1 and 3.2.

Constant time-step size
To perform the comparison based on a roughly comparable computational effort, the Crank-Nicolson and Ellsiepen schemes are run with a constant step-size twice as large as that of the Backward Euler scheme, eg, 0.3 or 0.15 days, respectively, during the creep phase of the Lubby2 model. The shear strain and associated errors are obtained for each time-integration scheme. The material model shows an initial immediate elastic response to the stress jump, followed by the expected transition from primary to secondary creep. The latter transition is accompanied by a reduction of the creep rate. The numerical errors are highest where deformation rates are highest and subsequently drop as the material enters the constant-rate stationary creep phase ( Figure 5).
Backward Euler, a first-order scheme, shows a slight underestimation of the strain in the transient phase due to the Taylor series expansion around the current time step ( Figure 5A). Likewise, a Forward Euler scheme originating from a Taylor series expansion around the previous time step overestimates the strain in this example (data not shown). Transitioning from a first-order to the second-order accurate Crank-Nicolson scheme, which can be considered to use an average of the tangents used in Forward and Backward Euler and is hence sometimes called semi-implicit, decreased the error magnitudes ( Figure 5B). This is despite the fact that the time step for Crank-Nicolson was two times larger than that used for the Backward Euler scheme.
Finally, Ellsiepen's scheme, also second-order accurate, is fully implicit in that both stages are solved akin to a backward Euler scheme. This scheme produced the minimal error magnitudes in the case studied among the three time-integration  Figure 6A, where the maximum errors of the Crank-Nicolson and Ellsiepen schemes can be seen to decrease at the second order compared with the first-order decrease observed with the backward Euler scheme.
A completely analogous behavior is seen for the Minkley model (constant time step size of 15 days for backward Euler sheme and 30 days for Crank-Nicolson and Ellsiepen's scheme), cf. Figure 7, albeit at different magnitudes and time scales. Thus, the implementation of all basics schemes was confirmed to behave as expected. We note in the Figure 7C a jump appears in relative and absolute errors. This is because we plot the absolute values of the errors and as change of sign occurs during the simulation.

Error-controlled time-step size
The main advantage of using Ellsiepen's scheme is that the step-size can be controlled by error tolerances associated with both global primary variables (eg, displacements, in our case strain) and local secondary as well as internal variables (eg, stress and inelastic strain).
Error bounds were provided for all variables, ie, strain, stress, internal tensorial variables, and internal scalar variables. The latter are applicable only in case of the Minkley formulation (see Section 3.2). The tolerances for errors associated with different quantities were set to individual levels according to their physical range and chosen units, compare Table 4 Their magnitude was varied in order to demonstrate their effect on the chosen time step and the error evolution.
The results shown in Figure 8 for the Lubby2 model reveal how the step-size was automatically adjusted to satisfy increasingly strict error tolerances. This is evidenced by an increasing number of time steps as the error tolerances decrease.
The second observation that can be made is that the time-step size changes as the simulation progresses. At the beginning of the computation, the adaptive time integration algorithm limited the step-size to accommodate the rapid nonlinear primary creep deformation. Subsequently, the time-step increased as the rock salt model reached stationary creep stage with lower creep rates.

Level
Stress (MPa) Strain Internal Tensor Internal Scalar High 1 · 10 −1 1 · 10 −6 1 · 10 −6 1 · 10 −2 Middle 1 · 10 −2 1 · 10 −7 1 · 10 −7 1 · 10 −3 Low 1 · 10 −3 1 · 10 −8 1 · 10 −8 1 · 10 −4 Note. The different sets "high," "middle," and "low" correspond to different accuracy requirements. All internal variables were of strain type. Examining the error measures more closely as the tolerance settings were changed from "high" to "low," the errors in the shear stress decrease accordingly ( Figure 6B). Note in passing that while the total absolute error remains below the specified tolerance in our examples, this is not the error for which the bound theorems provide the error estimate. The upper bound derived for the time-integration scheme refers to the local truncation error and does not account for error propagation in itself. The error measures computed here are the total errors which include this error propagation. In the present case, these total errors remain below the bounds specified for the local truncation error.
The convergence of the global and local Newton loops during the simple shear simulation of the Lubby2 model using Ellsiepen's scheme is shown in Figure 9. The fast initial loading are completed within the first five steps in a practically elastic fashion. Therefore, both local and global convergence were rapid. In the creep phase, especially the early and more nonlinear creep phase, convergence slows down both globally and locally and the algorithm repeats time steps with adjusted step sizes. Whenever the solution is approached, convergence accelerates. A time step repeated with a smaller step size converges faster.
The results of the Minkley model shown in Figure 10 exhibit a similar behavior. Two differences can still be observed. First, the Minkley model required fewer time steps in the present case due to the different material parameters and the lower applied stress. In combination, these two aspects decreased the degree of nonlinearity and hence allowed the simulator to ensure a specified level of accuracy with larger time steps. Second, error controls were applied to a wider range of tensorial and scalar variables as compared to Lubby2, compare Equations (15) with Equations (20) to (26).

Non-isothermal simulations with error control
To test the response of the algorithm in a thermo-mechanical setting, we increased the temperature very sharply twice during the stationary creep phase, first to a ΔT of 30 K then to a ΔT of 60 K.
The results of Ellsiepen's scheme applied to Lubby2 and Minkley models with abrupt temperature changes are visualized in Figures 11 and 12. The temperature increase has several consequences. The hydrostatic pressure in the material increases sharply due to the isochoric nature of the test. Due to the decrease in shear modulus, a temperature increase leads to a sudden increase in shear strain. A new primary creep phase is not initiated, as the acting stress on the Kelvin element remains constant and the element remains relaxed. However, the stationary creep accelerates due to the decreased viscosity. This creep acceleration in the absence of a new transient creep phase is not significant enough to warrant a decrease in time-step sizes. The observed decrease is instead related to certain must points the algorithm must hit during  the transition between load cases. Here, these are the start and end points of the heating phases. Such a must point also explains the final small time step. The small time step at the beginning of each creep phase is due to another constrained imposed by us which may also be relaxed: every load, a new load case begins with a conservative time step size assuming that the response of the material during the new load phase is not known. One can observe that the error magnitudes are maintained throughout the entire simulation according to the bounds specified for the local truncation error.

Triaxial loading and unloading
Triaxial loading and unloading are common laboratory testing regimes for rock salt. In this section, boundary conditions motivated by such tests will be applied to the numerical model to evaluate the adaptive time-integration algorithm for more complex boundary conditions and material response ( Figure 4B). The boundary conditions were applied as follows: 3. Maintain the confining stress in the x and directions ( 1 = 2 ) and hold the strain in the z direction ( ax ) at the maximum value (stress relaxation phase). 4. Maintain the confining stress in the x and directions ( 1 = 2 ) and apply a strain-controlled axial unloading ramp in the z direction ( ax ) back to the value that was obtained at the end of stage 1.
Since the parameter values for Lubby2 and Minkley differ due to different experimental references, distinct total time profiles and magnitudes of the applied load were chosen as shown in Table 5.
Instead of computing errors with respect to an analytical solution, a numerical reference solution was created using the Crank-Nicolson scheme, with a very small time-step size of 2 · 10 −5 d for Lubby2 and 2 · 10 −4 d for Minkley, respectively. These values were at least three orders of magnitude below the time step sizes selected by the adaptive time-stepping scheme.
Plasticity in the Minkley model was now allowed with a cohesion value of 1.6 MPa. 35,48 Other parameters were kept unchanged, see Tables 1 and 2. Figures 13 and 14 show that for specified error tolerances, the obtained numerical errors for both models are in the same order of magnitude. Despite the physically different model formulations, the different parameter ranges and the different boundary conditions in both simulations, the controlling error tolerances of the algorithm did not have to be adjusted to obtain this similar behavior. This is despite the major difference that in the Minkley model, a rate-independent plastic flow is activated which is not present in the Lubby2 formulation. Hence, the qualitative behavior of the results is very different in both situations: the Minkley model (in our simplified implementation) predicts plastic flow at a differential stress of about 6 MPa, while Lubby2 shows a smoothly varying stress-strain behavior during the loading phase according to its nonlinear formulation. In both models, the stress relaxation phase is captured. The unloading proceeds qualitatively different akin to the behavior observed in the loading phase.
In both cases, the algorithm adjusts the time-step sizes and hits all must points. When prescribing lower (ie, more stringent) error tolerances, the number of time steps increases (ie, the time-step size decreases) and the errors with respect to the numerical reference solution decrease, both for Lubby2 ( Figure 13A,B) and for Minkley ( Figure 14A,B), by one order of magnitude as prescribed by the error tolerances.

Cyclic triaxial loading
The contemporary or intended use of rock salt caverns for renewable energy storage in the form of compressed air energy storage or power-to-gas entails operational protocols with faster cycles during which storage pressures fluctuate more frequently and over wider ranges than storage caverns have been subjected to traditionally.
To obtain stress paths for the current tests which are roughly representative of such a cavern usage scenario, FE simulations of a cyclically operated cavern were performed (similar to Böttcher et al 36 ). Principal stress paths were extracted from a material point close to the cavern wall in the FEM simulation and applied to the present model shown as seen in Figure 4C. The simulated cavern was ellipsoidal with half-axes (diameters) of 70 M × 30 M and located at a depth of roughly 1000 m.
Cyclic pressure loading ( Figure 15) was applied on the wall of the cavern to simulate gas injection and extraction. The resulting differential stress 1 − 3 is also shown in Figure 15. These loads were studied here on a material point level with the time-adaptive simulation scheme.
The cycle duration of the Lubby2 model was set to 4 days and of the Minkley model to 200 days to initiate appreciable creep responses over only a few cycles for algorithmic testing and to cover operations ranging from short-term to long-term loading. For these simulations, the cohesion of the Minkley model was again set to a high value to isolate the creep/stress relaxation response. The reference solution is again obtained by a high-resolution Crank-Nicolson simulation.
The results of the Lubby2 and Minkley models under cyclic loading shown in Figures 16 and 17, respectively, show a qualitatively different behavior due to the different activated process time scales. In Lubby2, the transient creep is very dominant and clear hysteresis loops are visible. The simulations using Minkley's model show a situation in which the stationary creep response dominates. Again, in both cases, the errors (this time of strain) are of the same order of magnitude. One can also observe an accumulation of numerical errors in this example. Tighter error tolerances once again increase the accuracy of the solution by automatically shrinking the time-step size. The errors decrease by one order of magnitude in both cases, as specified by the tolerance settings.

CONCLUSIONS
In this paper, DIRK schemes with error-controlled adaptive time stepping were implemented and tested at a material point level in an application-oriented manner for two practically relevant constitutive models for rock salt. The version of the Minkley model used was a simplified one in the sense that the elaborate hardening/softening laws, usually a feature of this model, 32,53 were not implemented here. This means that the simulated plastic response does not represent certain experimentally observed features. The model can, however, be considered sufficient for the algorithmic testing performed in this article. In simple shear tests, the fully implicit scheme following Ellsiepen's rule outperformed backward-Euler and Crank-Nicolson schemes in terms of accuracy for a fixed time-step size. Time adaptivity was then tested using simple shear as well as triaxial configurations in isothermal and non-isothermal cases. The automatic adjustment of the time-step size was found to be robust in all cases for a single algorithmic parameter set which was retained at standard values found in the literature. In other words, no adjustments had to be made for material models or boundary conditions contrary to other time adaptivity schemes currently in use. The approach could handle primary and secondary creep phases as well as transitions to rate-independent plastic flow without any adjustments. Similarly, temperature dependence did not pose a problem. In a fully coupled FE simulation, the nodal temperatures could, in addition to the nodal displacements, be added to the global vector of unknowns used for controlling the time-stepping scheme. Refining the error tolerances controlling the time-step size determination also had the desired effect in that decreasing the tolerances resulted in decreased errors by the same order of magnitude. The presented algorithm has demonstrated that it can control the accuracy of time-integration algorithms based on parameter sets that are robust against the choice of constitutive model, boundary conditions, and simple thermo-mechanical couplings. For a given constitutive model, it may nevertheless be possible to obtain further improvements in accuracy or time-step size by considering its specific rheological features. For example, Bazant et al 56 combined a strain-softening creep model with an exponential algorithm based on the analytical solution for creep in its rheological elements. In addition to high accuracy, certain properties of the results can be ensured, such as vanishing stress at high tensile strain. While we aimed at presenting a general algorithm suitable for implementation into software independent of the particular model used here, a combination with approaches exploiting the structure of specific constitutive models can be beneficial. Given its performance in the current tests and its suitability for implementation into FE codes currently employing backward Euler schemes, Ellsiepen's DIRK scheme was identified as a promising candidate for expanding time integration procedures of OGS and similar codes.