Deformation of salt structures by ice-sheet loading: insights into the controlling parameters from numerical modelling

Subsurface salt flow is driven by differential loading, which is typically caused by tectonics or sedimentation. During glaciations, the weight of an ice sheet represents another source of differential loading. In salt-bearing basins affected by Pleistocene glaciations, such as the Central European Basin System, ice loading has been postulated as a trigger of young deformation at salt structures. Here, we present finite-element simulations (ABAQUS) with models based on a simplified 50-km long and 10-km-deep two-dimensional geological cross-section of a salt diapir subject to the load of a 300-m-thick ice sheet. The focus of our study is to evaluate the sensitivity of the model to material parameters, including linear and non-linear viscosity of the salt rocks and different elasticities. A spatially and temporarily variable pressure was applied to simulate ice loading. An ice advance towards the diapir causes lateral salt flow into the diapir and diapiric rise. Complete ice coverage leads to downward displacement of the diapir. After unloading, displacements are largely restored. The modelled displacements do not exceed few metres and are always larger in models with linear viscosity than in those with non-linear viscosity. Considering the low stresses caused by ice-sheet loading and the long time-scale, the application of linear viscosity seems appropriate. The elastic parameters also have a strong impact, with lower Young's moduli leading to larger deformation. The impact of both the viscosity and the elasticity highlights the importance of a careful parameter choice in numerical modelling, especially when aiming to replicate any real-world observations.


Introduction
Salt structures, which comprise a variety of salt rocks, are a common feature in sedimentary basins worldwide (e.g., Jackson and Talbot 1986;Baldschuhn et al. 1996;Warren 2006;Hudec and Jackson 2007;Fossen 2016). Salt rocks are mechanically weak and deform by viscous flow when subject to stress at geological time scales. The physical processes that contribute to subsurface salt flow include linear solution-precipitation creep and nonlinear dislocation creep (Van Keken et al. 1993;Urai et al. 2008). Salt flow is triggered by differential loading, which causes salt to flow from areas of higher load to areas of lower load. Differential loading is typically related to regional tectonic stresses, gravitational loading by sedimentation or erosion, or tilting of the salt-bearing basin (Jackson and Talbot 1986;Hudec and Jackson 2007). Classic examples of gravitational differential loading by sedimentation are the filling of supra-salt mini-basins (e.g., Waldron and Rygel 2005;Brandes et al. 2012) or the progradation of sediment wedges above a salt layer (e.g., Wu et al. 1990;Vendeville 2005). Salt flow may be induced by small differential loads and does not require the overlying deposits to be denser than the salt (Jackson and Talbot 1986;Cohen and Hardy 1996;Gemmer et al. 2004Gemmer et al. , 2005Rowan 2019). The salt-flow rate, the efficiency of salt expulsion from the source layer and the resulting geometries depend primarily on the salt viscosity, the thickness of the source layer and the rate of progradation of the overlying sedimentary systems (e.g., Cohen and Hardy 1996;Gemmer et al. 2004Gemmer et al. , 2005Vendeville 2005;Albertz and Ings 2012).
The weight of an ice sheet during a large-scale glaciation can be another source of gravitational differential loading that may trigger subsurface salt flow (Liszkowski 1993;Sirocko et al. 2008). Conceptual models suggest that salt rise and uplift above a salt structure will be amplified in ice-marginal settings, where the load of the ice sheet is applied onto the salt source layer (Fig. 1a). In contrast, the rise of salt structures will be hindered or even reversed to subsidence if an ice sheet transgresses a salt structure (Fig. 1b;Liszkowski 1993;Sirocko et al. 2008). In general, these conceptual models have been confirmed by numerical simulations (Kock et al. 2012;Kock 2013;Lang et al. 2014;Baumann et al. 2022;Zill et al. 2022). Notably, the weight of an ice sheet is sufficient to trigger subsurface salt flow although the load is low and transient compared with loading by tectonic or sedimentary processes (Lang et al. 2014).
The Central European Basin System (CEBS) is a key example of a formerly glaciated salt-bearing basin (Fig. 1c). The infill of the CEBS is characterised by a multitude of salt structures that comprise Permian Zechstein salt and display a variety geometries (Fig. 1c, d;Doornenbal and Stevenson 2010;Warsitzka et al. 2019). Initial salt movement in the CEBS began in the Early Triassic and continued in several phases until the Late Cenozoic (Jaritz 1973;Kockel 2003;Warsitzka et al. 2019). During the Pleistocene, large parts of the CEBS were repeatedly transgressed by major ice sheets. In the German part of the CEBS, for example, every salt structure was at least twice covered by the ice sheets, while salt structures farther to the northeast were covered by ice sheets even more often ( Fig. 1c; Ehlers et al. 2011;Lang et al. 2014Lang et al. , 2018Winsemann et al. 2020). Pleistocene or Holocene salt movements have been inferred for a number of salt structures in this area. However, it remains unresolved if Pleistocene and Holocene salt movements are the continuation of pre-Pleistocene processes controlled by the regional stress field or if ice-sheet loading triggered a reactivation (Sirocko et al. 2008;Al Hseinat et al. 2016;Hardt et al. 2021).
The response of the salt rocks within a salt structure to an external load is controlled by a range of parameters. Fig. 1 a, b Conceptual model of the response of a salt diapir to icesheet loading (modified after Lang et al. 2014). a An ice advance towards a salt diapir causes lateral salt flow from the salt-source layer into the diapir, leading to the rise of the diapir and uplift. b When the ice sheet transgresses the diapir, the diapir is pushed downwards and broadens, leading to subsidence above the diapir. c Map of the salt structures (blue) and maximum ice-sheet extents in the German part of the Central European Basin System. Salt structures are from Doornenbal and Stevenson (2010), ice-sheet extents are compiled from Ehlers et al. (2011), Lang et al. (2018) and Winsemann et al. (2020). d Schematic geologic cross-section of northeastern Germany, showing the characteristic architecture of the basin fills, including several salt structures (modified after DEKORP-BASIN Research Group 1999). The location of the cross-section is indicated by the black line (A-A') in (c) Important controlling factors include the material properties of the salt and the supra-salt rocks, the size and geometry of the salt structure and the magnitude of the load. Previous studies of the response of salt structures to ice-sheet loading indicate that larger displacements are mainly related to more mobile salt, larger amounts of salt rocks involved, larger loads and longer loading phases (Kock et al. 2012;Kock 2013;Lang et al. 2014;Zill et al. 2022).
The major challenge with respect to a better understanding of the response of salt structures to ice-sheet loading and the modelling of such processes is the multitude of controlling factors. The ice load depends on the ice-sheet thickness, the ice-advance and retreat rates and the duration of ice coverage. However, these factors are highly variable in space and time and are typically poorly constrained for specific locations (e.g., Lambeck et al. 2010;Hughes et al. 2016;Lang et al. 2018). In contrast, the mechanical behaviour of the involved geological materials is generally better constrained. However, different approaches exist to describe the rheologies, resulting in a large number of constitutive equations and material parameters (e.g., Van Keken et al. 1993;Ter Heege et al. 2005;Urai et al. 2008). The description of the involved materials is further complicated by the heterogeneity and complex geometries of real rocks.
The aim of this study is to investigate the controlling factors of the ice load-driven deformation of salt structures. The focus is on the impact of different rheologies of the geological materials on the model results. Parameters such as the viscosity of the salt rocks, the elasticity of the suprasalt rocks and the geometry and internal architecture of the salt structure were systematically varied to determine their respective impact. For the salt rocks, both linear and nonlinear viscosities were tested. Furthermore, the elasticity of the salt and non-salt rocks was varied. Salt structures typically display complex internal architectures that comprise evaporitic and non-evaporitic strata. The presence of internal strata with contrasting rheologies may affect the mobility of the entire salt structure. Therefore, salt structures with an internal architecture including layers of very low viscosity (representing potassium salt) or very high viscosity (representing anhydrite), respectively, were modelled to test how the internal architecture and different rheologies affect the results. This study is deliberately intended not to replicate any real-world scenarios although the model setup is inspired by salt diapirs in the CEBS. Instead, the model results will help to evaluate deformation structures and displacements observed at salt structures in nature in the context of ice-sheet loading as a trigger mechanism. The results will also yield important implications for the longterm safety of waste repositories build into salt structures. Such safety assessments need to cover the next 100,000 to 1,000,000 years, making the consideration of the impact of future glaciations necessary. Assessments of the long-term safety of repositories rely heavily on numerical simulations and an appropriate parameter choice, based on natural examples, is of paramount importance. The insights gained from this parameter study can also be applied in more general studies of the evolution of salt structures. Furthermore, the results contribute to a better understanding of recent salttectonic deformation observed in formerly glaciated basins.

Mechanical behaviour of rock salt
Salt rocks are mechanically weak and are commonly regarded as behaving like viscous fluids when deforming at geological time scales and strain rates (e.g., Hudec and Jackson 2007;Fossen 2016). The viscosity of natural rock salt is affected by the salt mineralogy, temperature, water content and grain size and the viscosity varies over several orders of magnitude (Van Keken et al. 1993;Weijermars et al. 1993;Mukherjee et al. 2010). The strain rates in deforming salt bodies typically range from 2 × 10 -16 to 8 × 10 -11 s −1 , corresponding to salt-flow rates of 1 × 10 -2 to 2 mm × a −1 Warren 2006). The differential stresses that drive salt-tectonic deformation typically range from 0.5 to 3 MPa with higher differential stresses of up to 5 MPa occurring at shallow depths (Carter et al. 1982;Spiers and Carter 1998;Urai 2005, 2007;Li et al. 2012;Henneberg et al. 2018Henneberg et al. , 2022Kneuker et al. 2018).
The physical mechanisms responsible for the viscous flow of rock salt are dislocation creep and solution-precipitation creep. During subsurface salt deformation and salt-tectonic processes, both mechanisms contribute to salt flow. The dominant deformation mechanism depends on the stress magnitude, confining pressure, temperature, fluid content and pressure, grain size, impurities and anisotropy (Van Keken et al. 1993;Ter Heege et al. 2005;Warren 2006;Urai et al. 2008;Kneuker et al. 2018).
Dislocation creep is caused by defects in the crystal lattice and the development of subgrains that lead to dislocations. Dislocation creep is a non-linear process, where the strain rate ̇ is related to the differential stress using a powerlaw equation (Fig. 2;Carter and Hansen 1983;Van Keken et al. 1993;Hunsche and Hampel 1999;Urai et al. 2008): The differential stress Δ is defined as Δ = 1 − 3 where 1 and 3 are the maximum and minimum principal stresses, respectively. The value of the power-law exponent n typically ranges from 1 to 7. The factor A summarises the viscous behaviour of the salt and is defined by a materialdependent parameter A 0 , the specific activation energy Q, the universal gas constant R (8.314 J mol −1 K−1 ) and the temperature T. The temperature strongly affects the strain rate for power-law creep (Fig. 2). Values of the material parameter A 0 are mostly derived from creep experiments. Solution-precipitation creep-also called pressuresolution or diffusion creep-is related to sliding along the crystal-grain boundaries and fluid-enhanced dynamic recrystallisation, involving mass transfer at grain boundaries. Dissolution occurs at grain boundaries under high stress, followed by diffusion through the grain-boundary fluid and recrystallisation in areas of lower stress (Carter and Hansen 1983;Spiers et al. 1990;Urai et al. 2008;Li and Urai 2016). For solution-precipitation creep, the strain rate is linearly related to the differential stress ( Fig. 2; Spiers et al. 1990;Li and Urai 2016): The factor B summarises the viscosity of the salt, B 0 is a material-dependent parameter and D is the mean grain size. Solution-precipitation creep is strongly grain-size dependant (Fig. 2). Solution-precipitation creep is a fluid-assisted process mostly active in wet rock salt with a brine content of 10-20 ppm (Urai et al. 2008). Values for the material parameter B 0 and the activation energy Q are derived empirically (e.g., Spiers et al. 1990). However, pressure-solution creep is hardly observed in experiments due to the limited time, low water content and high stresses (Li and Urai 2016;Bérest et al. 2019).
In geological models designed to reconstruct or predict the evolution of salt structures at large temporal (10 3 to 10 8 a) and spatial (10 3 to 10 6 m ) scales salt creep is usually simplified to a linear (Newtonian) viscosity (Gemmer et al. 2004(Gemmer et al. , 2005Albertz and Ings 2012;Nikolinakou et al. 2012Nikolinakou et al. , 2014Rowan et al. 2019;Granado et al. 2021). The constitutive equation for linear viscous materials is The simplification of applying a constant linear viscosity can be justified as linear solution-precipitation creep is considered as the dominant mechanism at the time scales, stress and strain rates relevant for geological processes (Urai et al. 1986(Urai et al. , 2008Weijermars et al. 1993). The linear viscosity η basically summarises all factors that define solution-precipitation creep. Estimates of the viscosity η of rock salt typically range from 10 16 to 10 20 Pa s (Van Keken et al. 1993;Weijermars et al. 1993;Mukherjee et al. 2010;Rowan et al. 2019).

Model set-up
Finite-element simulations of the interaction between salt diapirs and ice sheets were conducted using the commercial software ABAQUS (Version: ABAQUS 2020). The 2D ̇= Δ Fig. 2 Plot of the differential stress versus strain rate for non-linear dislocation creep, linear solution-precipitation creep and linear viscosity. Data for the plots are compiled from Spiers et al. (1990), , Hunsche and Hampel (1999) and Li et al. (2012) models were based on a 50 km long and 10 km deep schematic geological cross-section (Figs. 3,4). The models were meshed using 2D triangular plain-strain elements. The initial edge length of the elements was set to 10 m in the section representing the salt diapir and increased to 50 and 400 m at the sides and the base of the section, respectively. The nodes at the model sides were fixed in horizontal (x) direction and free in vertical (y) direction, while the nodes at the model base were fixed in vertical and horizontal direction.

Tested configurations
A number of model configurations were tested to evaluate the impact of different parameters on the modelling results (Table 1). The tested parameters include the geometry of the geological cross-section, the rheology of the materials, the location of the ice margin and the duration of the different model steps.

Geometry of the modelled section
The internal geometry of the modelled section was modified from Lang et al. (2014). The original section represents a simplified geological cross-section from a salt structure in the Subhercynian Basin at the southern margin of the Central European Basin System ( Fig. 1d; Baldschuhn et al. 1996;Brandes et al. 2012). The selected example salt structure actually forms part of a salt wall with a length of ~ 70 km (Brandes et al. 2012). However, in the modelled cross-section, it appears as a two-dimensional salt diapir. The geological cross-section represents a symmetrical salt diapir, supra-salt rocks, including the infill of symmetrical rim synclines, and basement rocks (Fig. 4). The supra-salt rocks comprise three different stratigraphic units that represent rocks of the Triassic Buntsandstein, Muschelkalk and Keuper Groups and the Cenozoic infill of the rim synclines ( Fig. 4; Tables 2, 3). The base of the diapir is 1500 m deep. The top of the diapir has a depth of 50 m. The diapir has a width of 1000 m. The rim synclines are 6 km wide and up to 300 m deep. The basement is a single unit that begins at a depth of 1500 m and extends to the base of the model at a depth of 10000 m (Fig. 4). The first model series was run using sections with a 250-m-thick salt source layer, which extended laterally from the diapir to the margins of the section (Fig. 4a, b). Internally, the salt structure in the first model series was homogeneous. In the second model series, the salt layer pinches-out 1000 m from the centre of the section (Fig. 4c). The uppermost ~ 50 m of the diapir was defined as caprock. Internally, two steeply dipping, up to 25-m-thick layers were included. These layers could be attributed different properties to represent either beds of potassium salt or anhydrite (Table 4).

Rheology of the model materials
The material parameters were defined to resemble the rheologies of natural sedimentary rocks. The supra-salt and basement rocks were treated as elastic, while the salt rocks were treated as viscoelastic.
The elastic behaviour of the different materials was defined by their Young's Moduli (E) and Poisson's ratios (ν). Two different sets of these parameters were tested (Table 2). Parameter set "LA" was taken from previous models of Lang et al. (2014). Parameter set "RK" was based on parameters that have been applied numerical simulations in the context of the long-term safety of radioactive waste repositories (e.g., Liu et al. 2017;Bertrams et al. 2020).
To simulate the viscous behaviour of salt rocks both linear and non-linear creep was implemented. For non-linear viscosity, the so-called "BGRa" creep law was applied, which describes the temperature-and stress-dependent steady-state creep of salt rocks (Hunsche and Schulze 1994; Hunsche  and Hampel 1999). The BGRa creep law was developed to predict the long-term temporal and spatial development of stress and deformation in underground workings (Bräuer et al. 2011). The BGRa creep law describes power-law creep with an exponent of n = 5, a pre-exponential material parameter A 0 = 2.08*10 -6 s −1 and an activation energy Q = 54000 J mol −1 (Hunsche and Schulze 1994;Hunsche and Hampel 1999). To account for the temperature dependency, temperatures of 300, 325, 350 and 400 K were tested. A temperature of 325 K was chosen as representing the reference temperature. The temperature gradient in the subsurface and temperature changes during a glacial cycle were not considered. It should be noted that the BGRa creep law was defined for specific, homogeneous rock-salt lithologies. For other salt-rock lithologies, the equation needs to be modified by an additional pre-exponential factor, which typically ranges from 1/64 to 16 and defines the so-called creeping class (Hunsche and Hampel 1999;Hunsche et al. 2003;Kock et al. 2012;Liu et al. 2017). For linear viscosity, viscosities ranging from 10 16 to 10 20 Pa*s were tested (Tables 1, 4). A viscosity of 10 18 Pa*s was chosen as representing the reference viscosity, which is in line with previous numerical models of salt tectonics (e.g., Gemmer et al. 2004Gemmer et al. , 2005Albertz and Ings 2012;Nikolinakou et al. 2012;Granado et al. 2021). In models that included internal layers within the salt structure, linear viscosities of 10 13 and 10 23 Pa*s were used for potassium salt and anhydrite, respectively.

Ice-sheet loading
Loading by the weight of an ice sheet was simulated as a pressure load that was applied at the surface of the model. The magnitude of the pressure load was ~ 2.65 MPa, corresponding to the weight of a 300-m-thick ice sheet. Such a thickness can be expected within a distance of 5-50 km from the margin of an ice sheet (Benn and Hutton 2010;Beaud et al. 2016). The surface profile of the ice sheet was not considered for simplicity and to avoid additional horizontal stresses in the upper units of the model (cf., Andersen et al. 2005). To simulate the advance and retreat of an ice sheet across the geological section, the surface of the model was subdivided into four segments that could be independently loaded and unloaded (Fig. 4a). During the ice-sheet advance, the surface segments were sequentially loaded, beginning at the right-hand side of the model. Complete ice coverage was represented by a step where all segments of the surface were loaded. During ice-sheet retreat the segments were sequentially unloaded, beginning at the left-hand side of the model. Loading and unloading of the individual surface segments was instantaneously. In the last step of each run no load was applied, representing post-glacial conditions. Two different configurations of the ice-sheet extent were tested. In models with complete ice coverage, all four surface sections were loaded during the step representing the maximum ice-sheet extent (Fig. 4a). In models with partial ice coverage, only the two surface sections at the right-hand side of the model (sections A and B) were loaded during the step of the maximum ice-sheet extent. In this configuration, the ice margin was located next to the diapir, but the ice sheet did not transgress the diapir (Fig. 4a).
The duration of the modelling steps that represented the ice advance and retreat was 150 a for each step. As the surface sections are 12.5 km long, this corresponds to advance and retreat rates of 75 m*a −1 . The duration of the glacial maximum step was set to 5000 or 15000 a, respectively. The durations and advance/retreat rates are within the range of those estimated for Pleistocene continental ice sheets (e.g., Ehlers 1990;Lüthgens et al. 2011;Hughes et al. 2016;Lang et al. 2018). The load-free post-glacial step had a duration of 10000 or 20000 a, respectively. After these durations, rates of displacement observed in the models became very low.

Results
The models show the deformation of the salt diapir and overlying strata during the advance and retreat of an ice sheet (Figs. 5, 6). Deformation in the sub-salt basement rocks is restricted to some vertical compression by the weight of the ice sheet (Fig. 5). The initial advance of the ice sheet commonly leads to some uplift at the top of the diapir (Fig. 5a). Uplift is caused by lateral salt flow from the source layer into the diapir due to the applied differential load by partial loading of the model surface (steps GA-1 and GA-2; Fig. 5a, c). Some uplift may also be related to the increased stress in the supra-salt rocks, potentially causing horizontal extension that exerts additional stress on the salt diapir. The uplift increases for longer periods of partial loading and lower viscosity salt. The observed maximum uplift was ~ 1.4 m. In models with high-viscosity salt, the initial uplift due to salt flow is very small and exceeded by downwards displacement of the whole model surface due to the elasticity of the material (Fig. 7a).
When the salt diapir is transgressed by the ice sheet, the top of the diapir is displaced downwards, leading to subsidence at the surface above the diapir (steps GA-3 and Gmax; Fig. 6b, c). The downward displacement is accompanied by some broadening of the diapir (Fig. 6c). Depending on the chosen parameters, the downwards displacement may continue for the complete loading phase, reach a maximum or may be followed by another reversal of the movement direction (Figs. 6, 7). The maximum downwards displacement represents an equilibrium state that is controlled by the viscosity of the salt and the elasticity of the surrounding materials. The observed maximum subsidence ranged from -0.6 to -11.4 m. In the lower part of the diapir, some upwards displacement occurs during complete loading, which is related to inflow from the salt-source layer on both sides (Fig. 6a, c). The retreat of the ice sheet leads to another change of the movement directions (Figs. 5, 6, 7, 8; steps GR-1, GR-2, GR-3 and IG) and the previous subsidence at the surface is largely compensated. The compensation of the previous deformation is related to the load applied to the salt source layer by the retreating ice sheet and the elasticity of the supra-salt rocks, which restore the original shape after unloading. Lateral salt flow due to the ice load applied onto the source layer only has a minor impact that can be observed when the ice margin is located at the centre of the model (step GR-2). The extracted curves display a change in the uplift rate during this step (Fig. 7).
The remaining deformation at the end of the runs depends on the choice of the material parameters. Some runs end with a stationary state, while in other runs a stationary state is not yet attained, even if the duration of the load-free step is doubled to 20000 a (Fig. 9a). In both cases the displacements at the end of the runs may be positive (i.e., uplift), negative (i.e., subsidence) or zero, ranging from 0.5 to − 0.4 m.
Models with only partial ice coverage during the maximum ice extent display a different behaviour. The salt diapir is free to rise due to lateral salt flow from the source layer as long as the ice margin is located next to the diapir (Fig. 9b). The uplift at the top of the diapir ranged from 0.5 to 2.8 m and was limited by the elasticity of the materials due to the lack of plasticity.

Linear and non-linear viscosity
The choice of a linear or non-linear viscous model for the salt rheology and the selected values has a large impact on the modelling results. In simulations with a non-linear viscous salt, the attained maximum displacements and the displacement rates are always distinctly lower and the largest part of the deformation can actually be attributed to the elasticity of the materials (Figs.7a, b, 8a). After unloading, the deformation is not fully restored by the elasticity of the overlying rocks, typically causing some remnant subsidence above the diapir. Different values were tested for the pre-exponential parameter defining non-linear creep, corresponding to different temperatures. The differences in the resulting displacements and displacement rates are, however, quite small (Fig. 7b). In configurations with partial ice coverage, the downwards displacement of the complete model surface due to the elasticity clearly exceeds the rise of the diapir due to lateral salt inflow, leading to net-subsidence during the loading phase (Fig. 9b).
In simulations with a linear viscous salt, the magnitude of the viscosity affects the maximum displacements and the shapes of the time-displacement curves. Lower viscosities cause higher maximum displacements (Fig. 7b). At the lowest tested viscosity, the sense of the displacement changes during the maximum loading phase as viscous and elastic deformation attain an equilibrium state. After unloading, displacements are restored by the elasticity of the supra-salt rocks. Remnant displacements at the end of the model runs were very low. In some models, the compensation was not completed at the end of the runs (Fig. 7). Even an extension of the unloaded phase to 20000 a did not allow for a complete compensation (Fig. 9a).
In addition to the displacement, the differential stress within the salt was extracted for three reference points. Reference points are located in the upper part of the diapir (centre of the section, 200 m below surface), the lower part of the diapir (centre of the section, 1100 m below surface) and the salt-source layer (1000 m right from the centre, 1370 m below surface). For both linear and non-linear salt rheologies, the differential stress rapidly changes when the distribution of the load changes. Strong differences between linear and non-linear viscosities occur in the evolution of the differential stress during longer periods of constant loading (Fig. 10). Peaks of the differential stress are attained extremely rapidly during the first two phases of ice advance (GA-1 and GA-2). The peaks are generally lower in models applying linear viscosities than in those applying non-linear viscosities, as could be expected from the relation of the stress and the strain rate (Fig. 2). The attained peaks are followed by exponential-style decreases of the differential stress during further ice advance and complete ice coverage (GA-3 and Gmax). In simulations with a non-linear viscous salt, the differential stress decreases to minimum values that are clearly above zero (Fig. 10b). In contrast, linear viscous salt displays very rapid drops of the differential stress to values close to zero and the values remain very low until the loading conditions change during ice retreat (Fig. 10a). Fig. 7 Temporal evolution of the vertical displacement at the crest of the salt diapir. The diapir has a laterally extensive salt-source layer (cf., Fig. 4b).
Steps of the ice advance (GA-1, GA-2, GA-3) and retreat (GR-1, GR-2, GR-3) have a duration of 150 a per step. The glacial maximum (Gmax) and postglacial (IG) steps have durations of 5000 and 10000 a, respectively. a Salt diapir with source layer: comparison of linear viscous (10 18 Pa*s) and non-linear viscous (BGRa at 325 K) salt and elastic parameter sets "La" and "RK". For comparison, the results of purely elastic models without viscous salt are also shown. b Salt diapir with source layer: comparison of different linear and non-linear viscous salt rheologies. The elastic parameter set was "LA" for all displayed curves The first two phases of ice retreat (GR-1 and GR-2) again trigger a rapid rise of the differential stress, which is again followed by exponential-style decrease during further ice retreat and the load-free phase (GR-3 and IG). The minimum values attained in models with non-linear salt remain clearly above zero even after 10000 a without load (Fig. 10b), while for linear viscous salt the differential stress decreases to near zero (Fig. 10a).

Elastic parameters
The choice of the elastic parameters strongly affects the results of the simulations. Lower Young's moduli (i.e., higher elasticity) generally lead to higher displacements (Figs. 7a, b, 8a). During the maximum loading phases, the subsidence above the diapir is larger and it takes longer to attain the maximum values as the viscosity (load and time dependent) works against the elasticity (load dependent/time independent) in the subsidence phase. For a linear viscosity of 10 18 Pa s, a maximum subsidence of − 2.5 m is attained after 1050 a for parameter set "LA" (BM-A11), while for parameter set "RK" (BM-A8) a maximum subsidence of − 11.4 m is attained after 4800 a (Fig. 7a).
After the removal of the load, models with higher elasticity (i.e., lower Young's moduli) require more time to return to the original shape. The displacement may not be fully compensated even 20000 a after the removal of the load (Fig. 9a). For a linear viscosity of 10 18 Pa*s the remaining displacements after 20000 a are + 0.04 m for parameter set "La" and + 0.47 m for parameter set "RK", respectively. The larger remaining displacements show that materials with higher elasticities take longer to restore the displacement against the viscosity.

Geometry of the model section
The geometry controls the distribution of the different materials in the model sections and has thus a major impact on the results (Figs. 7, 8). The presence of an extensive saltsource layer allows for lateral salt flow even if the distance between the diapir and the ice load is large. In models with a more complex geometry, there is no uplift during step GA-1 Fig. 8 Temporal evolution of the vertical displacement at the crest of a salt diapir without source layer (cf., Fig. 4c). Steps of the ice advance (GA-1, GA-2, GA-3) and retreat (GR-1, GR-2, GR-3) have a duration of 150 a per step. The glacial maximum (Gmax) and post-glacial (IG) steps have durations of 5000 and 10000 a, respectively. a Salt diapir without source layer, no internal layers. Results for linear and non-linear viscous salt and elastic parameter sets "LA" and "RK" are compared. b Salt diapir without internal layers compared to salt structure with low-viscosity layers, representing potassium salt, and high-viscosity layers, representing anhydrite 1 3 due to the lack of a salt-source layer. Some uplift may occur when the ice margin is located next to the diapir during step GA-2 (Figs. 5b, 6a, 8a, b). The subsidence above the diapir during full ice coverage is also lower than in models with a laterally extensive salt-source layer. The presence of internal layers of contrasting rheologies within the diapir has some impact on the results (Fig. 8). In models with higher-viscosity internal layers, the maximum displacements and the rates of viscous movement are lower than in models without internal layers. In contrast, the effect of lower-viscosity internal layers on the displacement is very small (Fig. 8b). In one run, the depth to the top of the diapir was increased to 250 m (Run BM-A22; Fig. 7b). The increase in thickness above the diapir reduced the maximum subsidence and the uplift, compared with the maxima attained in simulations with a 50-m-deep diapir.

Discussion
The model results are discussed with regard to the impact of the tested rheologies, focusing on the impact of linear vs non-linear viscous behaviour of the salt rocks and of the elasticity of the supra-salt rocks. Furthermore, we will explore which additional parameters should be included in future models to attain more realistic results.
Qualitatively, the results of all model runs are similar and concur with those of previous conceptual (Liskowski 1993;Sirocko et al. 2008) and numerical models (Kock et al. 2012;Kock 2013;Lang et al. 2014). Quantitatively, the model results were controlled by the spatial and temporal change of the loading conditions and the viscosity and elasticity of the materials. During partial loading of the model surface (i.e., ice advance and retreat), salt flow from ice-covered to ice-free sections occurred in response to differential loading (Fig. 6a, c). The load applied onto the salt-source layer caused lateral salt flow into the diapir, resulting in uplift above the diapir (Fig. 11a). Lateral salt flow and uplift of the diapir were clearly observed in models applying a linear viscosity. If non-linear viscosity is applied, the uplift was very low. The maximum observed uplift was 2.8 m for linear viscous salt and 0.5 m for non-linear viscous salt, respectively. As no plasticity was included in the models, the maximum uplift was controlled by the elasticity of the supra-salt rocks (Fig. 9b).
During complete ice coverage, subsidence above the diapir occurred due to downwards displacement and broadening of the diapir (Figs. 6b, c, 11b). In the lower part of the diapir, some upwards displacement occurred due to lateral inflow Fig. 9 Plots of the temporal evolution of the vertical displacement at the crest of the diapir. a Salt diapir with source layer. The durations of the glacial maximum (Gmax) and postglacial (IG) steps have been extended to 15000 and 20000 a, respectively. Results for linear and non-linear viscous salt and elastic parameter sets "La" and "RK" are compared. b Salt diapir with source layer in an ice-marginal position. Results for linear and non-linear viscous salt and elastic parameter sets "La" and "RK" are compared. The glacial maximum (Gmax) and post-glacial (IG) steps have durations 10000 a each 1 3 from the source layer from both sides. In the central part of the diapir, the displacement remained very low (Fig. 6b). Such different directions and magnitudes of displacement indicate complex stress and strain patterns within a salt diapir that should be considered when aiming for predictive models.
The downwards displacement of the upper part of diapir and subsidence at the surface during complete ice coverage may seem counterintuitive, as the pressure load applied by the ice sheet is the same everywhere along the model section. The salt movement can thus not be driven by a stress gradient due to differential loading. Instead, deformation is controlled by different ratios between vertical and horizontal stress in elastic (supra-salt rocks) and viscoelastic (salt) materials. In elastic materials, the horizontal stress is related to the vertical stress by the Poisson's ratio. The vertical stress is thus always larger than the horizontal stress. In viscoelastic materials, the horizontal stress approaches the vertical stress after a relatively short period of loading. The length of this period depends on the viscosity (< 10 a for the parameters involved in this study). Therefore, the horizontal stress becomes larger in the viscoelastic salt than in the elastic supra-salt rocks, allowing for a broadening of the salt diapir at the expense of the supra-salt rocks. The maximum subsidence above the diapir ranged from -0.6 m to -11.4 m, depending on the viscosity of the salt and the elasticity of the supra-salt rocks (Fig. 7a). Larger deformation is favoured by lower values of Young's modulus in the supra-salt rocks and lower viscosity of the salt (Fig. 11b). The viscosity (load and time dependent) controls the subsidence rate as it works against instant (load dependent/time independent) elastic deformation. The restoration of the deformation after unloading is controlled by the elastic properties of the suprasalt rocks. The maximum possible deformation in the model is controlled by the ratio of the viscosity of the salt and the elasticity of the supra-salt rocks. Due to the complexity of the model, however, the absolute impact of these factors cannot be defined.
The obtained displacements may seem low in relation to the dimensions of the model section and the minimum element-edge lengths. However, the finite-element method does not require the mesh resolution to be in the range of the expected displacements to obtain valid results (Fagan 1992;Zienkiewicz and Taylor 2005). Instead, element-edge lengths can be much larger than the expected displacements (cf., Hampel et al. 2019). Fig. 10 Plots of the temporal evolution of the differential stress, comparing the effects of linear and non-linear salt rheology. a Evolution of the differential stress for linear viscous salt (viscosity η = 10. 18 Pa*s; elasticity "LA"; run BM-A11). Data are plotted for the upper part of the salt diapir (centre of the section, 200 m below surface), the lower part of the salt diapir (centre of the section, 1,100 m below surface) and the salt-source layer (1000 m right from the centre, 1,370 m below surface). b Evolution of the differential stress for non-linear viscous salt (BGRa at 325 K; elasticity "LA"; run BM-A15). Data are plotted for the upper part of the salt diapir (centre of the section, 200 m below surface), the lower part of the salt diapir (centre of the section, 1100 m below surface) and the salt-source layer (1000 m right from the centre, 1370 m below surface) 1 3

The rheology of salt rocks
The deformation of salt rocks is dominated by non-linear dislocation creep at high stress, while linear solutionprecipitation creep dominates under low stress conditions (Hunsche and Schulze 1994;Urai et al. 1986;Schléder and Urai 2005;Lüdeling et al. 2022). The rheology of rock salt that deforms at low stress and strain rates, however, is not yet properly validated by experimental data and an area of ongoing research (Li et al. 2012;Li and Urai 2016;Bérest et al. 2019;Lüdeling et al. 2022). At geologically relevant stress levels, temperatures and time scales, salt rheology is at the transition between linear and non-linear viscosity (van Keken et al. 1993;Herchen et al. 2018;Bérest et al. 2019). Geological studies on the long-term evolution of salt-bearing basins or entire salt structures generally apply a linear viscosity, commonly using the same viscosity for the entire salt body (Gemmer et al. 2004(Gemmer et al. , 2005Nikolinakou et al. 2012Nikolinakou et al. , 2014. Very few studies have used stratified salt units with different viscosities (e.g., Albertz and Ings 2012) or compared linear and non-linear viscosities (e.g., Li et al. 2012;Herchen et al. 2018;Granado et al. 2021). Estimates of the viscosity of natural salt bodies span several orders of magnitude (Weijermars et al. 1993;Mukherjee et al. 2010;Rowan et al. 2019) and can thus not provide the input parameters required for higher-resolution numerical models.
Laboratory studies, which generally are conducted in short time, at high differential stress and with dry salt-rock samples, indicate that the rheology of salt is represented by a power-law viscosity (Bérest et al. 2019;Lüdeling et al. 2022). An extensive database of the relevant material parameters for different salt-rock lithologies exists (e.g., Hunsche and Hampel 1999;Hunsche et al. 2003;Bräuer et al. 2011), which are widely used as input for high-resolution numerical models (e.g., Kock et al. 2012;Liu et al. 2017;Bertrams et al. 2020). In contrast, few experimental data exist on the behaviour of salt-rock samples at low stress levels and strain rates. Compilations of stress versus strain rate plots using available experimentally derived data suggest that creep at differential stresses below 5 MPa is linear, while creep at differential stresses above 5 MPa follows a power-law with an exponent between 3 and 7 (Herchen et al. 2018;Bérest et al. 2019;Lüdeling et al. 2022). The magnitude of the stress at the transition is temperature dependant, requiring higher stresses at lower temperatures (Herchen et al. 2018).
In the presented models, the displacements are always larger in models with linear viscous salt than in those applying non-linear viscosity (Figs. 7,8,9). The differential stress caused by the 300-m-thick ice sheet within the salt unit is generally less than 2 MPa (Fig. 10). At this magnitude of differential stress, a linear viscous behaviour of the salt seems thus realistic (Herchen et al. 2018;Bérest et al. 2019;Zill et al. 2022). The simplification of salt rheology to linear viscous behaviour is appropriate in geological models focused on the long-term evolution of salt structures and their interaction with the supra-salt rocks and earth-surface processes. Considering salt rocks as linear viscous is an established practice in geological models that has been applied in a large number of studies, successfully reconstructing the geometries observed in natural salt structures and their supra-salt rocks (e.g., Gemmer et al. 2004Gemmer et al. , 2005Albertz and Ings 2012;Nikolinakou et al. 2012Nikolinakou et al. , 2014Rowan et al. 2019;Granado et al. 2021). More sophisticated rheological models are required in studies aiming at the internal behaviour of salt structures, especially with regard to the stability of underground workings, caverns or waste repositories. Complex stress and deformation patterns occur within natural salt structures, locally causing much higher stresses. In such models, the deformation of salt rocks is commonly decomposed into the contributions from linear solution-precipitation and non-linear dislocation creep (cf., Spiers et al. 1990). The application of such combined descriptions of salt rheology commonly applied in numerical modelling (e.g., Spiers and Carter 1998;Urai and Spiers 2007;Baumann et al. 2022;Zill et al. 2022). Zill et al. (2022) also modelled the response of salt structures to ice-sheet loading. Their results showed that a model combining linear and non-linear salt creep results in ~ 5 times more uplift of the salt structure Fig. 11 Summary of the model results. a An ice advance towards the salt diapir causes surface uplift above the diapir due to lateral salt flow and diapiric rise. b Ice coverage of the salt diapir causes subsidence above the salt diapir due to falling and broadening of the diapir than a model applying only the non-linear BGRa creep (Zill et al. 2022).
Natural salt structures display complex internal architectures, comprising a variety of salt rocks and non-salt rocks with contrasting mechanical properties (Hudec and Jackson 2007;Rowan et al. 2019). Even the incorporation of two relatively thin internal layers significantly affected the model results (Fig. 8b). Regional studies aimed at specific salt structures, therefore, need to better resolve the internal architecture to make realistic predictions of ice-load induced deformation. A rich, experimentally derived database of fitting parameters exists for power-law creep of different salt lithologies at high differential stresses (e.g., Hunsche et al. 2003;Bräuer et al. 2011;Liu et al. 2017). In contrast, linear viscous behaviour of different salt lithologies is less well constrained, as experimental data on creep at low differential stresses are rare (Bérest et al. 2019;Lüdeling et al. 2022). Assigning the appropriate linear viscosities when modelling heterogeneous salt comprising various lithologies remains thus challenging.

The impact of elasticity
The presented models demonstrate that the choice of elastic parameters may strongly impact the model results (Fig. 7a). Usually, the elastic behaviour of materials receives far less consideration than viscous or creep behaviour. In models concerned with long-term processes and non-transient loads, where the elasticity has a minor impact, this may be appropriate. The strong impact of the elasticity in the presented models is caused by the rapid and repeated redistribution of the ice load at the model surface. Young's moduli in parameter set "LA" are higher than in parameter set "RK" (Tables 2, 3). The ratio is about 3:1 for the Palaeozoic and Mesozoic rocks and 10:1 for the Cenozoic rocks. Lower values of Young's modulus allow for larger elastic displacements and vice versa. Parameter sets "LA" (Lang et al. 2014) and "RK" (Liu et al. 2017;Bertrams et al. 2020) are both compiled from several studies on the mechanical properties of the various lithostratigraphic units in the Central European Basin System (Tables 2, 3). The wide spread of the values can be explained by the heterogeneous nature of the lithostratigraphic units. The applied values are, however, in the range for generic geological materials (e.g., Turcotte and Schubert 2002). The results demonstrate that a careful choice of the elastic parameters for numerical modelling studies is required, if transient loads are involved.

Towards more realistic models
As in every numerical model, a number of simplifications was required for feasibility and focus on the main questions. The geometry of the model section is relatively complex to depict the architecture around existing salt structures (Figs. 1d, 4). However, with the exception of the Cenozoic infill of the rim synclines, the material parameters of the supra-salt rocks are very similar (Tables 2, 3). Interactions between real salt structures and ice sheets may be further complicated due to 3D effects, especially if salt structures are strongly asymmetric or have an elongated shape. The yield strength of the modelled rocks was deliberately ignored to avoid plastic deformation, which would have hampered the evaluation of the relative importance of the viscosity and elasticity. As shown in earlier models, including plastic behaviour affects the displacements remaining at the end of a model run but does not change the style of the deformation (cf., Lang et al. 2014). Especially poorly consolidated, mechanically weak rocks, as the infill of the rim synclines can be expected to deform plastically. In contrast, the load of an ice sheet may be insufficient to trigger permanent deformation in well-consolidated, mechanically strong rocks. If deformation due to ice-sheet loading occurs in such rocks, it is generally restricted to pre-existing weaknesses (Hampel et al. 2009(Hampel et al. , 2010Brandes et al. 2015). To obtain more realistic magnitudes of the displacements, however, the yield strength of the materials needs to be considered. Furthermore, previous models (Lang et al. 2014) displayed intense plastic deformation in the weak uppermost model unit at the locations of the ice margins. Deformation of near-surface unconsolidated deposits by the weight of an ice sheet is a well-known process (e.g., Andersen et al. 2005;Aber & Ber 2007), but outside the scope of this study. Furthermore, the deformation patterns in the uppermost model units at the ice margin observed by Lang et al. (2014) were partly unrealistic and did not resemble any natural deformation structures (cf., Andersen et al. 2005;Aber and Ber 2007).
Both dislocation creep and solution-precipitation creep of salt rock are temperature-dependant processes (Fig. 2;Van Keken et al. 1993;Urai et al. 2008). Models including the temperature as a boundary condition need to consider variations due to the geothermal gradient and the thermal conductivity of the different lithologies. Furthermore, the surface temperature and, with some time lag, the temperature in the subsurface dramatically change during the course of a full glacial cycle, culminating in the spatially and temporarily variable formation of permafrost. Existing models show that the temperature change during a glacial cycle may impact the behaviour of a salt structure (Kock et al. 2012).
Studies of ice-load induced deformation have primarily focussed on the effect of long-wavelength lithospheric flexure and isostatic rebound caused by glacial loading and postglacial unloading (e.g., Stewart et al. 2000;Hampel et al. 2009Hampel et al. , 2010Brandes et al. 2015). Glacial isostatic adjustment is characterised by subsidence below the ice sheet, which is compensated by post-glacial uplift. A forebulge is formed in a certain distance from the ice margin, 1 3 which is uplifted during glacial loading and subsides during unloading (Walcott 1970;Stewart et al. 2000). Glacial isostatic adjustment is not included in the presented models, as such long-wavelength processes act at very different spatial scales and thus require a very different model set-up. However, glacial isostatic adjustment may affect ice-load driven salt movement by basin-scale tilting of the salt-source layer and by stress changes due to lithospheric flexure. Tilting of the salt-source layer causes an elevation head that can trigger downslope gliding of the salt and overlying rocks (Hudec and Jackson 2007;Brun and Fort 2011). Although gliding may occur at low slopes, gliding at slopes of less than 1° requires basin widths, i.e., lateral extents of the salt-source layer, of 200-400 km (Brun and Fort 2011). Based on the models of the glacial isostatic depression beneath the Fennoscandian ice sheet (Talbot 1999;Stewart et al. 2000;Steffen and Wu 2011) the large-scale tilt would generally be lower than 0.1°. Løtveit et al. (2019) modelled the tilt in the forebulge area, suggesting that maximum values of ~ 0.2° may occur. The potential tilting of the salt-source layer by glacial isostatic adjustment and the lateral extent of the salt-source layer are thus probably too low to have a significant effect on the deformation. In model set-ups lacking a laterally extensive salt-source layer, there should be no effect at all. Stress changes induced by lithospheric flexure due to glacial isostatic adjustment probably have a rather low impact on the deformation of salt structures. Ice load-induced deformation of salt structures mainly occurs when the salt structure is in an ice-marginal position. In contrast, the highest stresses due to glacial isostatic adjustment occur either beneath the centre of the ice shield or in the area of the forebulge, tens to hundreds of kilometres outside the ice margin (Stewart et al. 2000;Hampel et al. 2009Hampel et al. , 2010Brandes et al. 2015). Tilting and stress changes induced by glacial isostatic adjustment may, however, be included in future models to test their impact. For the present study, the focus is on the effects of the material properties and tilting or stress changes due to glacial isostatic adjustment were thus not considered.
This study was intended primarily as a parameter study and was not aimed at reproducing any real-world situations. The deformation observed in the models is generally lower than measurements or estimates from field observations, which commonly indicate several tens of metres of displacement (cf., Sirocko et al. 2002;Stackebrandt 2005Stackebrandt , 2016Müller and Obst 2008;Al Hseinat et al. 2016;Hardt et al. 2021). The discrepancy between models and field observations was already observed in a previous modelling study by Lang et al. (2014). The impact of ice-sheet loading on the salt structures during the Pleistocene glaciations has probably been overestimated. When analysing the Pleistocene or Holocene salt movement, it is therefore important to consider how the studied salt structures behave under the respective geological conditions, especially with regard to the stress field, basin dynamics and sedimentation patterns. Displacements triggered solely by ice-sheet loading are generally low and probably insufficient to explain the larger deformation inferred from some field observations.

Conclusions
Our new modelling results highlight the impact of the rheology on the resulting deformation, both with regard to the viscosity and elasticity. In general, our results support previous conceptual (Liszkowski 1993;Sirocko et al. 2008) and numerical models (Kock et al. 2012;Kock 2013;Lang et al. 2014;Zill et al. 2022). In all conducted simulations, the load of the ice sheet triggers lateral and vertical salt flow, finally resulting in uplift or subsidence above the salt diapir in the centre of the model section. The location of the ice margin relative to the diapir exerts the first-order control on the style of the response of the diapir to ice-sheet loading. Uplift occurs if the load is only applied to the surface above the salt-source layer, causing lateral salt flow into the diapir. If load is applied to the surface above the top of the diapir, the diapir is pushed down and subsidence occurs. The downwards movement of the diapir is accompanied by a slight broadening at the expense of the surrounding supra-salt rocks. The maximum downwards displacement of the diapir is limited by the elasticity of the supra-salt rocks. After the removal of the ice load, the attained displacements are largely compensated by the elasticity of the supra-salt rocks, as no plasticity was included.
The directions and magnitudes of the displacements depend on the model geometry, the viscosity and the elasticity of the modelled materials. The results highlight the importance of a cautious parameter choice for both the viscous and elastic behaviour of the modelled materials. Model sections with extensive salt-source layers generally allow for larger displacements than those with a lateral pinch-out of the salt. The presence of more competent layers within the salt diapir also reduces the displacement.
The viscosity of the salt has a major impact on the model results. Linear viscous salt is more mobile than power-law viscous salt and thus allows for larger displacements. Furthermore, the range of the resulting displacements is far wider, if linear viscosity is applied. The use of a linear viscous salt rheology seems appropriate for the boundary conditions applied in this model, where the differential stress caused by a 300-metres-thick ice sheet is low, the time scale for the deformation is several thousands of years and the focus of the study is on the behaviour of the entire salt diapir and the interactions with the supra-salt rocks and earthsurface processes.
The choice of the elastic parameters of the modelled materials has a large impact on model results that may even exceed the impact of the viscosity. More elastic materials (i.e., lower Young's modulus) allow for larger displacements in all modelling steps. The role of elasticity is commonly neglected in large-scale models of geological processes. The results clearly show that a careful choice of the elastic parameters for numerical modelling studies is required, if transient loads are involved.