The signi cance of ground freezing is becoming ever more germane as the design of new urban tunneling systems requires more complex geometries and higher bear- ing capacities, which are limited with conventional construction methods. Ground freezing is an advanced construction technique to make the water-saturated subsoil impermeable and temporarily increase its strength and sti ness. This study reports experimental investigations consisting of single-stage and multi-stage creep tests un- der uniaxial loading. The comparison of the di erent loading types reveals the in- uence of the stress-strain history on the rate- and temperature-dependent behavior of frozen granular soils. We extend the constitutive model for frozen soils proposed by Cudmani et al. (2022) to consider stepwise loading and creep by coupling creep time with stress-strain history. Moreover, we simulate element tests and compare the simulations with our own experimental data as well as data from the literature to achieve the rst step in validating the extended model. The good agreement of the numerical and experimental results con rms the constitutive model's ability to capture the main features of the complex mechanical behavior of frozen granular soils for single-stage as well as multi-stage loading under constant temperatures. Keywords: ground freezing, frozen soil, creep, constitutive model, element tests Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 4 1 Introduction Ground freezing is widely used in geotechnical applications to temporarily increase the strength and sti ness of the subsoil and provide water tightness below the groundwater table. Ground freezing has many advantages over other ground stabilization techniques. On the one hand, it has a relatively low environmental impact. On the other hand, owing to its versatility and reliable quality control through temperature measurements, ground freezing is often preferable to other available soil improvement techniques in urban situ- ations with complex boundary conditions (Harris 1995; Andersland and Ladanyi 2003). Many ground freezing applications can be found in the literature, e.g., for the construc- tion of underground stations, excavation of cross passages to connect tunnel tubes and building underpinning (e.g., (Cudmani and Nagelsdiek 2006; Russo et al. 2015; Han et al. 2016; Orth 2018; Pimentel and Anagnostou 2019)). Recently, the intensive use of ground freezing to excavate the cross passages under the Suez Canal (Egypt) impressively demon- strated the importance of ground freezing in modern tunnel constructions (Phillips et al. 2021). In general, frozen soil bodies supporting tunnel excavations exhibit varying stress states and di erent loading types through the construction stages (Andersland and Ladanyi 2003; Orth 2018). The excavation leads to the shearing of the frozen soil and a relatively rapid increase in stress. After completing the excavation, the frozen soil body deforms under a predominantly constant stress state and must support the soil above it until the completion of the tunnel lining. Common excavation techniques, such as partial face ad- vance, result in a stepwise increased loading of the frozen soil body. Russo et al. (2015) reported three long-duration excavation steps during the challenging construction of the Toledo underground station in Naples (Italy). Several excavation stages with di erent lengths of advance induced a stepwise increased loading on the supporting frozen soil over a period of more than six months. Classen et al. (2019) and Zhou et al. (2021) considered the application of ground freezing to support very large tunnel excavations with up to fourteen excavation stages. From a practical point of view, these and other examples clearly highlight the need to improve our fundamental understanding of the in uence of the loading history on the mechanical behavior of frozen granular soils. Eckardt (1979b,a, 1982) extensively investigated the creep behavior of frozen medium sand samples under uniaxial loading conditions. Single-stage and multi-stage creep tests were carried out. In the single-stage creep test, the load was increased monotonically to the desired value and then maintained constant. In the multi-stage creep tests, Eckardt (1979b,a, 1982) applied the same nal load as in the single-stage creep tests, but this was achieved stepwise. In each loading step, creep was allowed during a period before the next load increment was applied. In single-stage as well as multi-stage creep tests, Eckardt observed primary (decreasing creep rates), secondary (nearly constant creep rates), and tertiary creep (increasing creep rates) according to the de nitions proposed by Orth (1986) and Andersland and Ladanyi (2003). Moreover, the di erent strain evolutions resulting from the single-stage and multi-stage creep tests evidenced the in uence of the stress and strain history on the creep behavior. Indeed, Eckardt (1979b,a, 1982) also concluded that during secondary creep, the strain evolution of both test types is similar and thus inde- pendent of the stress-strain history. Nevertheless, the evaluation of creep strain rates for di erent stress histories supporting this conclusion is missing in this study. Therefore, the Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 5 equivalence of single-stage and multi-stage loading regarding the creep behavior of frozen soils has not been conclusively clari ed. Indeed, numerous studies like (Andersland and Akili 1967; Vyalov et al. 1989; Zhou et al. 2020) have contributed to the understanding of the complex mechanical behavior of frozen soils depending on the stress and strain history. However, the in uence of the stress and strain history on the rate-, stress-, and temperature-dependent mechanical behavior of frozen soils is not yet fully understood. According to Zhao et al. (2023) and Cudmani et al. (2022), there are many constitu- tive models to either predict the shear or the creep behavior of frozen soils. Nevertheless, there is a lack of sophisticated models able to capture the rate-, stress-, and temperature- dependent behavior of frozen soils under compressive and tensile loading. In addition, only a few models in the literature take into account the in uence of the con ning pres- sure on the shear resistance and creep behavior of frozen soils. The model proposed by Cudmani et al. (2022) intends to ll this gap and is already validated for predomi- nantly monotonic loading. However, the consideration of stress-strain histories deviating from monotonic loading (e.g., stepwise, loading-unloading) on the mechanical behavior of frozen soils is still missing. In fact, there are many novel constitutive models for frozen soils (e.g., (Ghoreishian Amiri et al. 2016; Xu et al. 2017; Yao et al. 2018)), which neither take into account nor are validated for stepwise loading. This study aims to present an important and impactful contribution to the following ends: First, we introduce a comprehensive experimental program including uniaxial single-stage and multi-stage creep tests to observe, understand, and quantify the main characteristics of the stress-strain history on the mechanical behavior of frozen granular soils. Here, the objective is to extend the existing constitutive model by Cudmani et al. (2022) to con- sider stepwise loading and creep based on our experimental ndings. Subsequently, the focus is placed on the model validation for multi-stage creep using element tests. Here, we compare our numerical results with our own experiments as well as data from the literature. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 6 2 Uniaxial single- and multi-stage creep tests 2.1 Testing material, sample preparation, and testing program Table 1 includes the physical properties, and Figure 1 the grain size distribution curve of the tested frozen sand. Table 1: Physical properties of the tested frozen sand emin/emax[-]d10=d30=d60[mm]s[g/cm3] 0.52/0.82 0.44/0.53/0.64 2.65 0.02 0.06 0.2 0.6 2 6.3 Grain size [mm] 020 40 60 80 100Passing by weight [%] silt sand mediumfine coarsegravel Figure 1: Grain size distribution of the tested frozen sand The tested samples were circular and had a diameter of about 50 mm and a height of about 100 mm. They were prepared with the water-sedimentation method (see (Lade 2016) for further details). Afterward, the samples were additionally compacted on a vibrating ta- ble. Subsequently, the specimens were frozen isotropically in a freezer to a temperature of =4:3C. The preparation method led to a dense state with an average dry density of d= 1:66g=cm3and a water content of w = 0 :20. We cut a few samples into several parts and determined their grain size distribution, water content, and dry density. The results for each part were similar, thus con rming the soil sample homogeneity and reproducibil- ity using the selected preparation method. Nevertheless, according to Lade (2016), it is dicult to achieve a full saturation of the soil samples by using the water-sedimentation method. In this study, the specimens were not fully saturated and had an average sat- uration degree of S r90 %. In fact, studies like (Enokido and Kameta 1987) indicate Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 7 that saturation degrees below 80-85% strongly a ect the essential mechanical behavior of frozen coarse-grained soils. In contrast, the in uence of partial saturation with S r90 % on the shear and creep behavior of frozen sand is relatively small, and its mechanical behavior is comparable to fully saturated conditions. Before testing and in accordance with Orth (1986), the sample end plates were trimmed and smoothed out. Moreover, silicone grease on a Te on layer was used as lubrication between the specimen and the load frame to reduce friction. In addition, the samples were wrapped in a rubber cover to prevent sublimation during the tests. In order to investigate the in uence of the loading history on the creep behavior, a series of uniaxial single-stage and multi-stage creep tests on frozen sand were performed at Zen- trum Geotechnik of the Technical University of Munich. The tests were conducted in a climatic chamber (air temperature cooling of the specimens) using a dead load oedometer test apparatus. Hence, the loading speed during the tests was relatively fast as the weight discs for the load application were applied manually and instantaneously after each other. Moreover, the load application occurred after the samples were frozen, resulting in post- freezing con nement, as described by Nishimura and Wang (2019). Figure 2 illustrates the testing procedure for the two di erent load types. ab a, b: tests with a single creep stage (a) tests with three (multiple) creep stagese1=2.5 %e1=3.5 % (b) Figure 2: Test procedure for a) single-stage creep tests; b) multi-stage creep tests In the single-stage creep test, vertical pressure was monotonically applied, and then the stress was kept constant as the frozen soil samples crept. In contrast, the multi-stage creep test consisted of three loading stages. After the initial loading stage, the second and third load increases occurred after an axial strain 1= 2:5 % and1= 3:5 %, respec- tively. In each step, the load was maintained constant, and the creep deformations of the frozen soil samples were monitored. The axial strains chosen for the load increase in the multi-stage creep tests result from two di erent constraints. On the one hand, the load increase should occur during primary creep (see Andersland and Ladanyi (2003)) and thus at decreasing axial strain rates. Consequently, the selected axial strains at the load increases must be smaller than the expected average axial strain mat the turning point (secondary creep) before tertiary creep begins. Here, literature data from Orth (1986) and our own single-stage creep tests indicated that mwas in the range of 4 6 %. On Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 8 the other hand, a sucient creep time is required during each loading stage to accurately determine and evaluate the axial strain rates. As a compromise between these two con- straints, we chose to increase the load in the multi-stage creep tests at 1= 2:5 % and 1= 3:5 %. The testing program, including six single-stage and three multi-stage creep tests, is sum- marized in Table 2. Table 2: Uniaxial creep tests on frozen sand at =4:3C Test number 1[kN/m2] Single-stage tests ST1 2.5 ST2 3.0 ST3 3.5 ST4 4.0 ST5 4.5 ST6 5.0 Multi-stage tests MT1 2.5/3.0/3.5 MT2 3.0/3.5/4.0 MT3 4.0/4.5/5.0 2.2 Test results The results of the single-stage and multi-stage creep tests are compared in Figure 3 and include the plotted total strain evolution (left side) and the total strain rate evolution (right side) over the testing time. However, for the creep stage, the plotted total strains in Figure 3 are predominantly viscoplastic. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 9 0 500 1000 1500 200000.020.040.060.080.1 1 10 100 1000 10000 10 -4 10 -3 10 -2 10 -1 10 0 (a) 0 500 1000 1500 200000.020.040.060.080.10.12 1 10 100 1000 10000 10-410-310-210-1100 (b) 0 100 200 300 400 50000.020.040.060.080.10.12 1 10 100 1000 10000 10-410-310-210-1100 (c) Figure 3: Results of the single- (symbols) and multi-stage (lines) creep tests at =4:3C: left: evolution of axial strain over time; right: evolution of axial strain rate over time The single-stage creep tests (symbols) in Figure 3 show the well-known creep behavior of frozen soils (e.g., (Andersland and Ladanyi 2003)): On the one hand, the axial strain increases with time. On the other hand, the strain rate rst decreases (primary creep) Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 10 and then increases (tertiary creep) with time. The testing time, at which the minimum axial strain rate _ m(secondary creep) is reached, and the tertiary creep begins, is called lifetime t m, according to Orth (1986) and Cudmani et al. (2022). The following is a brief summary of the fundamental physical mechanisms that lead to the well-known creep be- havior of frozen soils observed in the single-stage creep tests in Figure 3. Andersland and Akili (1967) and Parameswaran (1980) rst reported physical interpretations of the ex- perimentally observed deformation behavior of frozen soils. Based on crystal mechanics, Orth (1985, 1988) comprehensively described and explained the underlying fundamental mechanisms that lead to primary (decreasing creep rates), secondary (minimum creep rate), and tertiary (increasing creep rates) creep. There are two physical processes in the ice that control the well-known rate-, stress- and temperature-dependent frozen soil be- havior: dislocation glide limited by discrete obstacles (thermal activation) and di usion. In addition, the simultaneous initiation and propagation of microcracks in the ice matrix plays an important role. Dislocation glide limited by discrete obstacles and/or di usion leads to ice hardening and strengthens the frozen soil, whereas crack initiation and propa- gation weakens it. Both strengthening and weakening processes occur simultaneously, but are temporally distinct. During primary creep, ice hardening dominates the behavior of the frozen soil, resulting in decreasing creep rates. During secondary creep, i.e. when the minimum creep rate _ mis reached, the frozen soil reaches its maximum creep resistance. Subsequently, cracking of the ice matrix becomes the dominant process rather than ice hardening, resulting in continuously increasing creep rates and unstable frozen soil creep behavior. In principle, the multi-stage creep tests (solid and dashed lines) in Figure 3a)-c) reveal a similar evolution of axial strain over time as the single-stage creep tests. However, the strain evolution resulting from the multi-stage creep tests lies between that of the single-stage tests. As seen on the right side of Figure 3, the load increase in the multi- stage creep tests led to a jump in the strain rate. As expected, during the rst loading stage, the evolution of the axial strain rate is similar in both the multi-stage and the corresponding single-stage creep tests. At the beginning of the second and third creep stage, the strain rate in the multi-stage creep tests lies above that of the corresponding single-stage creep tests for the same testing time. With increasing time, the strain rate approaches that of the single-stage creep tests. Deviations from this behavior pattern occur when the third and last creep stage begins near the lifetime t m(the turning point of the strain evolution), as can be seen in test MT3 in Figure 3c). In this case, the axial strain rate does not decrease anymore with time. Instead, it remains almost constant over a relatively short period of time and then increases, thus indicating the beginning of tertiary creep. Here, the behavior of the frozen soil becomes unstable. We nd that both the deviatoric stress level and the previous stress history in uence the lifetime t m. In contrast, based exclusively on the results in Figure 3, there is no clear evidence that the minimum axial strain rate _ mat the turning point is independent of the loading history. However, in accordance with the principles of crystal mechanics applied to frozen granular soils by Orth (1985, 1988), the physical processes in frozen soils leading to ice hardening, e.g. dislocation glide limited by discrete obstacles (thermal activation) and di usion, re- main largely independent of the loading history. Orth explained physically that these two processes are only in uenced by the deviatoric stress level and temperature. Moreover, the resulting ice hardening process dominates the frozen soil behavior until the minimum Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 11 axial strain rate _ mand corresponding lifetime t mare reached. Thus, the independence of _mfrom the loading history at constant temperature is physically reasonable. In contrast, the simultaneous crack initiation and propagation in the ice may be in uenced to some extent by the loading history. Based on the experimental results in Figure 3 and the consideration of crystal mechanics, _ mis assumed to be dominated only by the current deviatoric stress level at constant temperature. Based on the results of uniaxial single-stage creep tests, Ting et al. (1983) and Orth (1986) concluded that the total axial strain at the turning point is roughly independent of the temperature and the axial stress 1. Figure 4 shows the average axial strain m at the turning point (t = t m) for our tests compared with test results published by Orth (1986), which were obtained using frozen sand similar to that described in Section 2.1. 0 2 4 6 8 10 00.01 0.02 0.03 0.04 0.05 0.06 0.07 Mean value - 20 % Mean value + 20 % Mean value - 4.3 °C Multi-stage loading - 4.3 °C - 5.0 °C (Orth, 1986) - 10.0 °C (Orth, 1986) Single-stage loading Figure 4: Average axial strain mat the turning point (t = t m) for uniaxial single-stage and multi-stage creep tests for di erent axial stresses 1and temperatures According to Orth (1986), the observed scatter of min Figure 4 can be traced back to di erences when setting up the sample end plates before the start of the tests. In addition, slight di erences in the initial density and degree of saturation of the samples could also contribute to this scattering. As can be seen, the results of our single-stage and multi- stage creep tests con rm the ndings of Ting et al. (1983) and Orth (1986) regarding the negligible in uence of the temperature and the axial stress 1onm. In addition, since the total average axial strains mdetermined in the single-stage and multi-stage tests are in the same range (see Figure 4), our preliminary conclusion is that mis almost independent of the load history. Nevertheless, additional investigations with more general loading and creep stages are required to further validate the loading history independence of both _mandm. In addition, it is worth mentioning that in our constitutive model, mis not an input quantity, but results from the integration of the constitutive equations (see following sections). Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 12 3 Conceptual framework to consider the in uence of the loading history on the creep behavior of frozen soils A crucial relationship of the constitutive model for frozen granular soils proposed by Cudmani et al. (2022) is the relationship between the normalized axial strain rate _ 1=_m and the normalized time t =tm. According to Orth (1986) and Cudmani et al. (2022), for single-stage creep and monotonic shearing with constant strain rate, this relationship is independent of the temperature and the stress level and can be described by Equation 1. _1 _m= exp ( ) exp t tmt tm (1) Here, _mis the minimum axial strain rate, t mis the lifetime, and is a material parameter. As shown by the experimental results of the multi-stage creep tests in Section 2.2, in contrast to _ m, the lifetime t mdepends on the loading history. For this reason, Equation 1, which establishes a unique relationship between t mand _m, is unable to describe the creep behavior observed in the multi-stage creep tests. In the following, a time transformation procedure is proposed to capture the behavior observed in the multi-stage creep tests using Equation 1. The proposed procedure is based on the following assumptions derived from the experimental results discussed in Section 2.2. Assumptions a)-f): a) All tested frozen soil specimens had the same initial density, degree of saturation, and amount of frozen pore water. b) The total stress is equal to the creep stress. c) The average axial strain mat the turning point is mostly independent of the temperature, axial stress, and loading history, as shown by the experimental results in Figure 4. d) The minimum axial strain rate _ mdepends only on the temperature and the axial stress. A possible dependence of _ mon the loading history can be disregarded, as indicated by the experimental results in Figure 3. e) The lifetime t mdepends on the temperature, the axial stress, and the loading history. f) There is a unique relationship between the strain rate _ 1, the strain 1, and the stress1. This relationship does not depend on the loading history. Figure 5 illustrates the developed procedure to transform the creep time t and the "life- time" t mdetermined in multi-stage creep tests into an equivalent time and lifetime of the corresponding single-stage creep tests. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 13 S2 S1 M2 b c tcS2 tbM2 = taM1 tmS2 tmM2 tmS1 M1 Required for the determination of tcS2 a Dtbc eM2 (tbM2 ) = eS2 (tcS2 ) eM2 (t)= eS2 (t*,M2) with t*,M2=t- Dtbc >0emS2 » emM2 » emS1 (a) b c tcS2 tbM2 = taM1 tmS2 tmM2 tmS1 M2 S2 S1 M1 = Dtbc Dtbc = tbM2 - tcS2 a == emS1 emM1 emS2 emM2 ecS2 ebM2 eS2 (t*,M2) eM2 (t)= with t*,M2=t- Dtbc >0 (b) Figure 5: Determination of the normalized axial strain rate _ 1=_mover normalized time t=tmfor a multi-stage creep test M consisting of the loading stages M1 and M2 compared to the single-stage creep tests S1 and S2 with S1=M1andS2=M2under a constant temperature The sub-index S denotes single-stage and M multi-stage creep tests for the same tem- perature. For simplicity, the multi-stage creep test M assumably consists of two loading stages, M1 and M2. The stepwise-applied stresses M1andM2are equal to those in the single-stage creep tests S1 and S2, with S1=M1andS2=M2. In addition, t describes the actual testing time. Naturally, the time t, the strain , and the strain rate _ mof test M1 are identical to test S1 during the rst loading stage. Hence, the actual strain rate _S1= _M1, the corresponding minimum strain rate _ S1 m= _M1 m, and lifetime tM1 m= tS1 mare also the same. As can be seen in Figure 5b, it is assumed that the frozen soil sample has not achieved the lifetime during stage M1. At the end of the rst loading stage, M1 at t = tM1 a= tM2 b, the stress increases from M1toM2(see Figure 5a). In accordance with the aforementioned assumptions c)-f), the relationship M2(t) =S2(t) must be ful lled for t tM2 band ttS2 c, respectively. Thus as shown in Figure 5a, the strain versus time curves for M2 and S2 are identical starting from points b and c, but they are shifted by a period of time t bc= tM2 btS2 c. Using the condition M2 tM2 b =S2 tS2 c and assuming that S2 tS2 c is known from the single-stage test, the value of tS2 ccan be determined to calculate the time shift t bc. According to Figure 5a, the transformed time t;M2= ttbcin combination with M2(t) =S2 t;M2 for t;M2>0 describes the creep behavior in the loading stage M2 based on the corresponding single-stage creep test. The transformed lifetime for M2 is tM2 m= tS2 m= tM2 mtbcderived from the lifetime tS2 mof the single-stage test S2. Summarizing, the transformed time t;M2is the time in a single-stage creep test required to achieve the deformation M2 t>tM2 b withS2=M2. In particular, t;M2at the point b of M2 equals t;M2 b= tM2 b tM2 btS2 c = tS2 c. According to the assumption c), both samples M2 and S2 achieve the same strain m=M2 tM2 m =S2 tS2 m and the same minimum strain rate _ M2 tM2 m = _S2 tS2 m at the turning point. In accordance with the Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 14 above-described time transformation procedure, we propose the following equations to calculate the time-dependent evolution of the strain rate _ M2(t) for the loading stage M2 based on the results of the corresponding single-stage creep test S2: _S2 t;M2 _S2 m= exp ( ) exp t;M2 t;M2 mt;M2 t;M2 m (2) t;M2= ttbcwith t>0 and t>tbc (3) tM2 m= tM2 mtbc= tS2 m (4) _M2(t) = _S2 t;M2 (5) In addition, the equations 2-5 can be generalized to describe the strain rate in the loading stage Mj of a multi-stage creep test consisting of n stages based on their corresponding single-stage test Sj. The loading stage number is de ned with the superscript index j and the corresponding time increment with the subscript index i: _Sj t;Mj _Sj m= exp ( ) exp t;Mj t;Mj mt;Mj t;Mj m (6) t;Mj= ttiwith t>0 and t>ti (7) tMj m= tSj m (8) _Mj(t) = _Sj t;Mj (9) Here, t i= tMj 0tSj 0is the time di erence between the beginning of Mjat tMj 0and the single-stage creep time tSj 0, at which Mj tMj 0 equalsSj tSj 0 . In addition, the multi- stage and single-stage creep behavior for the loading stage j = 1 and time increment i = 1 correspond to t = tS1 0= tM1 0, t 1= 0 and t;M1= tS1= t. To sum up, Figure 6 compares both the normalized axial strain rate _ 1=_mover normal- ized time t=tmfor the single-stage creep tests presented in Figure 3, as well as _ 1=_mover the normalized transformed time t;Mj=t;Mj mfor the additional multi-stage creep tests. In accordance with the proposed time transformation procedure, the multi-stage test results are converted to equivalent single-stage tests, and therefore, as expected, both single- and multi-stage test results converge to Equation 1. These ndings provide another strong argument for the load history independence of _ m, and agree well with the physical con- siderations in Section 2.2 from which we originally made this preliminary assumption. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 15 0.001 0.01 0.1 1 10 100 Figure 6: Uniaxial creep tests at =4:3C according to Table 2: Evolution of the normalized axial strain rate _ 1=_mover normalized time t =tmfor the single-stage as well as _1=_mover the normalized transformed time t;Mj=t;Mj mfor the multi-stage creep tests In summary, the evolution of the axial strain rate for multi-stage creep depends on the temperature, the axial stress, and the loading history. The latter in uence only the lifetime t mbut not the minimum axial strain rate _ m. In the following, the existing constitutive model for frozen granular soils proposed by Cudmani et al. (2022) will be extended to take into account the in uence of the previous loading history on the stress- strain behavior based on the conceptual model presented in this section. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 16 4 Extension of the constitutive model to consider multi-stage creep 4.1 Constitutive model for predominantly monotonic loading, according to Cudmani et al. (2022) Cudmani et al. (2022) developed a constitutive model for predominantly monotonic load- ing that takes into account the in uence of the con ning pressure and di erentiates be- tween compressive and tensile strength and creep. In this section, only the basic equations of the model are presented; for further details, see Cudmani et al. (2022). The elastic- viscoplastic model for frozen granular soils is based on a rheological model of the Maxwell type ( _=_e+_v)): _=L: (__v) (10) where _is the stress rate tensor, _is the total strain rate tensor, _vis the viscous strain rate tensor, and Lis the fourth-order isotropic elastic sti ness tensor. The isotropic elastic sti ness tensor can be expressed in terms of Young's modulus E and Poisson's ratio , see details in (Cudmani et al. 2022). The viscous strain rate tensor for a general stress state can be obtained from the following equation: _v=jj_mjjexp ( ) exp t tmt tm s jjsjj(11) where sis the deviatoric stress tensor, and is a dimensionsless material constant. The minimum creep rate jj_mjjcan be obtained by Equation 12. jj_mjj=r 3 2_ expK1 + 273:4+ ln _ cr(p;q; )  ()1 (12) Here,represents the temperature in Celsius, K 1is a material constant, _  is a reference strain rate (_  = 1%=min), while the temperature-dependent uniaxial reference stress  at the reference strain rate, the lifetime t m, and the equivalent uniaxial creep strength cr are de ned in Equations 13-15, respectively.  () = 1() 2(13) tm=r 3 2c=jj_mjj (14) cr(p;q; ) =1 2 B cos  3 + C q + Dp +q B cos  3 + C q + Dp2+ 4Aq2 (15) The parameters 1, 2, and c in Equations 13 and 14 are material constants of the 1D model. In Equation 15 the parameters A, B, C, and D represent material constants of the 3D model, while p and q denote the Roscoe stress invariants, and denotes Lode angle. In this paper, the mechanical sign convention is used, i.e., compressive stress is negative and tensile stress positive. Therefore, the Lode angle for compressive stress is = 3and for tensile stress, it is = 0. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 17 4.2 Generalization of the transformed creep time The constitutive model developed by Cudmani et al. (2022) simulates the time-dependent stress-strain behavior of frozen soils for predominantly monotonic loading. As described in Section 3, the behavior observed in a multi-stage creep test can also be predicted with Equation 1, originally valid for single-stage creep tests, with the help of a time transfor- mation. In this section, this concept is generalized to consider the in uence of multi-stage stress paths on the mechanical behavior of frozen soils. Figure 7 illustrates the proposed extension of the constitutive model to consider multi- stage creep. Here, the loading stage number is de ned with the superscript index j and the corresponding time increment with the subscript index i. sSj+1(t) sSj+1(t) > sSj (t) q = const.Increase of vertical stressWith creep time transformation (extended model)No creep time transformation (Cudmani et al., 2022) sSj (t) tmSj+1tmSj Dt itS i j+1tS i j = t i tS i j - tS i j+1t * i+1t i+1 ||emSj+1|| ||emSj || Figure 7: Generalization of the transformed creep time t i+1for multi-stage loading We assume that the constitutive equations 10 to 15 have been integrated from t 0to t i. The stress state Sj(t), the viscous strain i, and all state variables depending on stress state, temperature , and time t are known at the speci c time tSj i. Since we consider loading under a constant temperature , the time-dependent evolution of the viscous strain rate _v;i(t) (solid line in Figure 7), corresponding to single-stage creep, is de ned by the stress state Sj(t) according to Equations 11 to 15. In particular, the minimum strain ratejj_Sj mjjand lifetime tSj mfor single-stage creep with Sj(t) can be determined with Equations 12 to 15. For increment i+1, the stress increases from Sj(t) to Sj+1(t) and the creep time steps forward from t ito t i+1= ti+ t i. For single-stage creep with Sj+1(t), the viscous strain rate _v;i+1(t) (dashed line in Figure 7), the minimum strain ratejj_Sj+1 mjj, and the lifetime tSj+1 mcan be determined with Equations 11 to 15. Here, the constitutive model by Cudmani et al. (2022), which does not consider the in uence of the loading history on the viscous behavior, predicts a tertiary creep strain rate with ti+1>tSj+1 mfor this loading history in Figure 7. However, this deviates from the strain Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 18 rate evolution experimentally observed in Section 2.2, where no tertiary creep behavior after the load increase was observed, and thus t i+1;exp<tm;exp. In fact, the actual viscous strain rate _v;i+1 t i+1 results from the single-stage creep curve for Sj+1(t) by using a transformed creep time tinstead of the global time t. In accordance with the description of the viscous strain rate under 1D conditions in Section 3, we consider the in uence of the loading history on the creep behavior by including the transformed creep time t i+1in Equation 11: _v;i+1=jj_mjjexp ( ) exp t i+1 t m;i+1t i+1 t m;i+1 s jjsjj(16) t i+1= ti+1 tSj itSj+1 i+1 = tSj i+ t i tSj itSj+1 i = tSj+1 i+ t i (17) t m;i+1= tSj+1 m (18) Here, tSj+1 mis the lifetime during single-stage creep with Sj+1(t). In addition, tSj+1 iis the creep time required to achieve the norm of the current viscous strain jjv(ti)jjin a single-stage creep test with Sj+1(ti). The time tSj+1 ican be determined from the solution of the following equation: jjv(ti)jj=jj_Sj+1 mjjexp ( )ZtSj+1 i t0exp t t m;i+1t t m;i+1 dt (19) In Equation 19, jjv(ti)jjis the norm of the viscous strain tensor at the end of the incre- ment i.jj_Sj+1 mjjand t m;i+1can be determined from Equations 12 to 15 for Sj+1(t). To sum up, the extended constitutive model able to consider multi-stage creep is de ned by Equations 10 and 12 to 19. As suggested by Gudehus et al. (2008), we implemented these constitutive equations in the form of a 'user-de ned material' (UMAT) FORTRAN subroutine. In this manner, the model can be used in commercial Finite Element Analysis (FEA) codes and thus is applicable in geotechnical boundary value problems. Figure 8 presents the numerical integration scheme of the extended constitutive equations for FEA. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 19 2nd - 5th iteration: updated si+1abs( || Dev(t i+1,||em,i || ,tm,i )|| - || Dei+1|| ) / || Dei+1|| < 0.001Provided by the FEA code, among others, at the start of increment i+1: ei, Dei, si, ti+1, Dti, state variables: t* i , || em,i ||, tm,i t* i +1 = ti+1, (tS i j =tS i j+1)No creep time transformationNumerical integration in order to determinate tS i j+1 according to Equation 19 t* i +1 = tS i j+1 + DtiYes si+1 » si and Dei+1»Dev,i+1 , (Deel ,i+1»0) No Calculate ev,i+1 and si+1 according to Equation 16 and 10 and update si+1 after 5th iteration Return of final si+1 and updated state variables to the FEA codeCalculate || em,i+1 || , tm,i+1 with si+1, (1st iteration: si+1 = si) User-defined material (UMAT) User-defined material (UMAT) si+1 ¹ si and Dei+1»Dev,i+1+Deel ,i+1, (Deel ,i+1¹0) 2nd - 5th iteration: updated si+1 Figure 8: Numerical integration scheme of the constitutive equations for FEA to consider multi-stage creep For the calculation of the increment i+1, the FEA code provides, among others, the stress tensor i, the total strain tensor i, the time t i, and the increments of the total strain tensor  iand the time increment t i, respectively. In addition, the creep time t i, the minimum strain rate jj_m;ijj, and lifetime t m;iof the increment i are also available as state variables. Using the incremental values provided by the FEA code, the total strain tensor i+1=i+ iand the time t i+1= ti+ t ican be determined, while the stress tensor i+1has to be updated from the integration of the constitutive model. As can be seen in Figure 8, the UMAT calculation of increment i+1 consists of ve iterations. For the rst iteration, the minimum strain rate jj_m;i+1jjand the lifetime t m;i+1of the increment i+1 are calculated with i+1=i, because the stress tensor i+1has not been updated yet. Subsequently, the norm of the total strain tensor increments jji+1jjprovided by the FEA code is compared with the calculated norm of the viscous strain tensor increments jjv;i+1(ti+1;jj_m;ijj;tm;i)jjto identify a stress change for increment i+1. The stress i+1 is equal to i, if  i+1only consists of a viscous strain rate component ( el;i+1= 0) and has approximately the same value as  v;i+1(ti+1;jj_m;ijj;tm;i). Hence, there would be no stress change between increment i and i+1, and no creep time transformation would be necessary. This strain increments (strain rate) comparison intends to avoid unnecessary computing time caused by the numerical integration to nd tSj+1 ifor increments without stress changes. As can be seen in Figure 8, on the one hand, no stress changes between increments i and i+1 result in the creep time t i+1= ti+1. On the other hand, numerical integration is required to determine tSj+1 iaccording to Equation 19 when i+1is not equal to i. The numerical integration in Equation 19 is solved by Newton's method to approximate tSj+1 i. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 20 Subsequently, the transformed creep time t i+1is the sum of tSj+1 iand the time increment ti. After completing the rst iteration process inside the UMAT subroutine, the viscous strain rate tensor _v;i+1and the resulting stress rate tensor _i+1are calculated accord- ing to Equations 16 and 10. Subsequently, the stress tensor i+1=i+1+_i+1tiis updated and used for the second iteration inside the UMAT to determine jj_m;i+1jjand tm;i+1again. At the end of the iteration process (here: ve iterations), the nal stress tensor i+1and the updated state variables are returned to the FEA code for convergence consideration. Summarizing, the extended constitutive model presents a simple approach without in- troducing more material constants to consider the in uence of the loading history on the creep behavior of frozen soils during steady-state temperature conditions. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 21 5 Model validation This section compares the prediction of the extended constitutive model with experimental data. First, we focus on the shear and creep behavior under single-stage loading. Here, it is important to verify that the novel transformed creep time tin Equation 16 does not in uence the already validated predictive capacity of the model by Cudmani et al. (2022). Subsequently, we analyze the model response for multi-stage loading by comparing experimental and numerical results using our own test data as well as data from the literature. 5.1 Monotonic loading Gudehus and Tamborek (1996) tested a frozen Karlsruhe sand (medium sand). The dry densityd= 1:72 g=cm3, void ratio e = 0 :54, and the water content w = 0 :18 of the frozen soil samples, were similar to the freezing test characteristics in Orth (1986). Therefore, our study adopts the material parameters for frozen Karlsruhe sand determined in (Cudmani et al. 2022), as listed in Table 3. It is worth mentioning that Cudmani et al. (2022) presented their study using the geotechnical sign convention, i.e., compressive stress is positive and tensile stress negative. In this paper, the mechanical sign convention is used (see also Section 4.1). Therefore, the 3D material parameters A, B, C, and D in Table 3 were converted accordingly to ful ll the mechanical sign convention. Table 3: Material constants for frozen Karlsruhe sand after (Cudmani et al. 2022) one-dimensional model three-dimensional model * E c 1 2 K1 A B C D [MPa] [-] [%] [MPa/C] [-] [-] [K] [-] [-] [-] [-] 500 0.3 2.40 3.05 0.59 0.69 3817 2.11 -3.18 3.18 3.33 * converted to the mechanical sign convention Gudehus and Tamborek (1996) conducted con ned compression and creep tests at a tem- perature of =20C under di erent constant traces of stress I =1+2+3. Fig- ure 9 compares the experimental and numerical results. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 22 0 0.05 0.1 0.15 0510 15 20 25 0 0.05 0.1 0.15 0510 15 20 25 (a) Con ned compression test with di erent constant traces of stress I =1+2+3 10 010 110 210 310 410 510 -5 10 -4 10 -3 10 -2 10 -1 10 010 110 210 310 410 510 -5 10 -4 10 -3 10 -2 10 -1 (b) Con ned creep test with di erent deviator stresses q Figure 9: Experimental (symbols) and numerical (solid lines with lled symbols) results of con ned freezing tests at =20C [Data from (Gudehus and Tamborek 1996)] As can be seen, the proposed model continues to realistically predict the increasing shear resistance for increasing I in con ned compression tests under a constant strain rate (Figure 9a). Nevertheless, for very high stress states (I = 24 MPa and I = 30 MPa), the measured compressive peak strength does not increase signi cantly. Gudehus and Tamborek (1996) explained this nding as being due to the pressure-dependent reduction of the freezing point. As a result, this reduction led to ice softening and, simultaneously, a reduction of the shear resistance between the ice matrix and the grain skeleton. These weakening e ects are complex and dicult to calculate numerically. Therefore, the sim- ulation of the con ned compression tests in Figure 9a leads to a slight overestimation of the peak strength for very high stress states. In addition, slight di erences can be observed between the simulations with the model by Cudmani et al. (2022) (dashed lines) and the extended model presented in this study (solid lines) for the con ned compression tests in Figure 9a (right side). The use of the ex- tended model improves the description of the actual deformation behavior as the ultimate shear strength q u=13reaches the peak shear deformation u;1>0:06 more slowly. In contrast, the predicted behavior with the model by Cudmani et al. (2022) is sti er (u;1>0:06) than observed in the experiments. Hence, the model accuracy improves in Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 23 terms of the shear behavior of frozen soil through the proposed equivalent creep time. In addition, the comparison of con ned creep tests with I = 12 MPa and di erent devi- ator stresses q at =20C is depicted in Figure 9b. The evolution of the axial strain rate over time is slightly underestimated by the extended model. In fact, the simulated minimum axial strain rates _ mare lower in comparison to the experiments. As a result, the predicted lifetime t mbecomes longer. As expected, there are no essential di erences between the model response by Cudmani et al. (2022) (dashed lines) and the extended model (solid lines). The introduction of the transformed creep time tin Equation 16 does not change the model response in terms of the predicted creep behavior. In summary, the implementation of the transformed creep time timproves the model response for monotonic shearing without having any noteworthy in uence on the model accuracy for single-stage loaded creep. 5.2 Multi-stage loaded creep 5.2.1 Using our own test data After con rming the model response for single-stage loading, we investigate the validation for multi-stage loading by simulating the multi-stage creep tests described in Section 2.2. According to Cudmani et al. (2022), seven material constants have to be determined for the 1D model. The material parameters determined in (Cudmani et al. 2022) based on Orth's (1986) experimental ndings cannot be used for our analysis. Our own freezing test characteristics are di erent from those in (Orth 1986) mainly due to the di erent dry density of our tested frozen sand samples (for in uences on testing results, see, e.g., (Ting et al. 1983)). Therefore, it is necessary to re-calibrate the model parameters using our own testing data. According to Cudmani et al. (2022), at least two uniaxial creep tests and three uniaxial compression tests are needed to determine the seven 1D-material constants. Thus, we performed additional uniaxial creep and compression tests to avoid using the same experimental data for the model calibration and validation. The most important results of these additional tests are summarized in Table 4. Table 4: Additional testing data to calibrate the 1D-model Uniaxial creep tests with constant loads Test number   1 _m tm [C] [MPa] [%/min] [min] M-7 -10.0 6.5 6.50 E-3 334.2 M-8 -10.0 8.8 9.01 E-2 25.8 Uniaxial compression tests Test number   c _1 [C] [MPa] [%/min] C-1 -10.0 10.4 1.0 C-2 -15.0 13.7 1.0 C-3 -20.0 16.5 1.0 The method to determine the 1D-material constants using experimental data is extensively described in (Cudmani et al. 2022) and (Orth 1986) and will not be repeated here. The Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 24 material parameters used for the 1D-model are presented in Table 5. Table 5: Material constants for the 1D model E c 1 2 K1 [MPa] [-] [%] [MPa/C] [-] [-] [K] 500 0.3 2.40 2.44 0.63 0.69 3817 The results of the uniaxial experimental and numerical multi-stage creep tests are com- pared in Figure 10. Note that the experimental test results of MT1 to MT3 have already been introduced in Section 2, viz. Table 2 and Figure 3. 0 200 400 600 800 1000 00.02 0.04 0.06 0.08 0.1 0.12 Experiment 0 200 400 600 800 1000 00.02 0.04 0.06 0.08 0.1 0.12 Simulation (a) Axial strain evolution over time 1 10 100 1000 10 -3 10 -2 10 -1 10 0 Experiment 1 10 100 1000 10 -3 10 -2 10 -1 10 0 Simulation (b) Axial strain rate evolution over time Figure 10: Experimental and numerical results of uniaxial multi-stage creep tests at =4:3C The model prediction agrees well with the experimental data. The minimum axial strain rate _m;1and the lifetime t mfor the last load step are in accordance, thus con rming that the proposed constitutive model extension is appropriate. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 25 To extend the model validation, we simulated uniaxial multi-stage creep tests at dif- ferent temperatures using data from the literature. This discussion follows in the next section. 5.2.2 Using data from the literature As mentioned in the introduction, Eckardt (1979a,b, 1982) extensively investigated the mechanical behavior of frozen sand with the help of single-stage and multi-stage creep tests under compressive stress states. The granulometric properties (grain size distribution) and state variables (void ratio, degree of saturation, water content, and dry density) of Eckardts tested frozen medium sand were very similar to the freezing tests by Orth (1986), who investigated frozen Karlsruhe sand in the same frost laboratory at the University of Karlsruhe (Germany). Thus, we also adopted the determined material parameters for frozen Karlsruhe sand in Table 3 for the tested sand by Eckardt. The experimental and numerical results of uniaxial multi-stage creep tests for di erent increasing vertical stresses at=10C and=20C are shown in Figure 11. 0 200 400 600 00.020.040.060.080.10.12 1 [MPa] s1= 2.0/3.5/7.0s1= 1.7/3.0/4.0/5.0/6.0s1= 3.0/4.0/5.0 s1= 1.0/2.0/3.0/4.0/5.0/6.0 s1= 0.6/0.9/1.2/1.5/2.1 (a)=10C [Data from (Eckardt 1979b)] 0 100 200 300 400 500 00.020.040.060.080.10.12 s1= 4.0/8.0/12.0s1= 4.0/6.0/8.0/10.0 s1= 1.0/3.0/4.0/5.0/6.01 [MPa] (b)=20C [Data from (Eckardt 1982)]] Figure 11: Evolution of axial strain over time: Experimental (symbols) and numerical (solid lines) results of uniaxial multi-stage creep tests The multi-stage creep tests include three to six load steps N1;j(j = 36) with incre- mental stress increases  1between 0:3 to 4:0 MPa. The vertical stress 1for each load step ranged from 0 :6 to 12:0 MPa. The model accuracy concurs more strongly at =10C (Figure 11a) when compared to =20C (Figure 11b). The deviations between the calculated and measured deforma- tions at=20C occur during the load applications due to sti ness di erences. Here, the simulations overestimate the actual elastic deformations during the load application in comparison with the deformations observed in the experiments. However, the creep deformations (evolution of the predominantly viscous axial strain) for every loading step are in good accordance. For both the experiments and the simulations, the time-dependent axial strain increases Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 26 as a function of the stress state. As shown in Figure 4, the average axial strain mat the turning point (t = t m) for this frozen Karlsruhe sand is around 0 :04 to 0:06. The model realistically captures the turning point of the creep curves in accordance with the experiments in Figure 11, as the axial strains increase strongly for 1>0:06 and creep failure occurs. In fact, there are two simulations in Figure 11, one at =10C with1= 2:0=3:5=7:0 MPa and one at =20C with1= 4:0=8:0=12:0 MPa, which both predict creep failure immediately after the last load increase. The model response concurs with the corre- sponding experiments that also show creep failure shortly after the last load step. The proposed models ability to predict creep failure after a load increase is a powerful feature for safely improving the geotechnical design of frozen soil bodies. The consideration of multi-stage creep provides clear economic bene ts as the lifetime of the frozen soil is ulti- mately longer in comparison to a single-stage lifetime evaluation under the same loading stages. In addition, our proposed model extension guarantees the early identi cation of sudden creep failure. Besides, for very long creep times (here: t >500 h) under low stress states (1= 0:6=0:9=1:2=1:5=2:1 MPa and 1= 1:0=3:0=4:0=5:0=6:0 MPa) the model does not predict creep failure which is also in accordance with the corresponding experiments. Therefore, we assume that the over-engineering of frozen soil bodies by using the model is unlikely. Summarizing, the essential features of the rate-, stress-, and temperature-dependent char- acteristics of frozen soils are well reproduced by the model for multi-stage loaded creep. Moreover, the improved model can precisely predict the rate-dependent evolution of axial strain, both for di erent numbers of load steps N1;jand for varying stress with increasing 1. Here, the model validation includes a very large stress state range 1. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 27 6 Limitations and recommendations for the use of the constitutive model Cudmani et al. (2022) have already discussed several shortcomings of the constitutive model, e.g., upper and lower strain rate boundaries and predominantly monotonic loading. Even though the improved, far-reaching constitutive model, designated by the acronym EVPFROZEN ( elastic- viscoplastic frozen ), tackles some of these issues, we should point out that some limitations remain. Model calibration Cudmani (2006) and Cudmani et al. (2022) have already presented the unambiguous 1D and 3D model calibration procedure for the in total eleven parameters. In particular, eight freezing tests consisting of three uniaxial compression and two creep tests, one uniaxial tensile test and two triaxial compression or creep tests are required. Thus, EVPFROZEN only requires a single, unambiguously determinable set of material parameters to model a wide range of temperatures, stress states, and strain rates. Presented model validation for stepwise loaded uniaxial creep Despite the comprehensive model validation in Section 5.2, it should be noted that based on uniaxial creep tests and their corresponding axial strain measurements, EVPFROZEN cannot be fully validated in terms of all its possible predictive capabilities, e.g. volumet- ric deformations during creep. Therefore, the present model validation for multi-stage creep is naturally limited to a certain extent. However, the scope of our study is the experimental investigation of frozen granular soils under uniaxial multi-stage loading and its constitutive modeling consideration. In this context, studies such as Ting (1983) and Orth (1988) revealed that the uniaxial creep deformations of frozen granular soils are almost purely deviatoric. The comparison between uniaxial and triaxial creep tests by Orth (1986, 2018) and Cudmani et al. (2022) con rms that the creep deformations under con nement are still predominantly deviatoric. In addition, our own experimental results in Figure 3 and Figure 6 also show that stepwise loading does not a ect the essential (creep) strain rate evolution. Consequently, the measured frozen soil creep deformations are still mainly deviatoric. Therefore, it is sucient to experimentally measure and nu- merically compare the axial strain as well as axial strain rate and use these data as a rst step to validate the model for uniaxial multi-stage loading. Mechanical behavior upon unloading Despite the e orts made so far (e.g., (Andersland and Akili 1967; Zhou et al. 2016; Li et al. 2017)), there is a lack of sophisticated experimental data regarding the creep behavior of frozen soils under stepwise unloading or a combination of loading and unloading. We plan to close this knowledge gap and con rm the presented approach accordingly. At this stage, it is recommended to rst perform stepwise unloaded multi-stage creep tests and compare the model response with the experimental results if unloading plays a dominant role in the boundary-value problem under consideration. Strain rate changes during shearing Orth (1986) performed uniaxial compression tests of frozen sand with varying strain rates _1at a constant temperature. The uniaxial compressive strength cfor stepwise increas- ing or decreasing strain rates was not similar to that under a constant strain rate. Orth (1986) argued that changing strain rates strongly in uence the initiation and propagation of cracks in the ice matrix. Hence, uniaxial compression tests with varying strain rates Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 28 were not directly comparable to tests with a constant strain rate due to di erent viscous behaviors. In contrast, based on con ned compression tests with frozen clay under con- stant and varying strain rates at a constant temperature, Wang et al. (2017) suggested a potentially unique stress-strain relationship for a given strain rate. This would be con- sistent with the isotach theory proposed for unfrozen clays ( Suklje 1957). Despite these e orts, the in uence of varying strain rates on the mechanical behavior of frozen soils has not been fully understood. Our currently proposed model does not consider a possible in uence of varying strain rates on the shear behavior of frozen soils. The predicted stress state is uniquely de ned by the current strain and strain rate under a constant tempera- ture. In uence of varying temperatures The mechanical behavior of frozen soils changes not only in cases of varying stress states under constant temperatures but also for varying temperatures at constant stress states. However, there are few experimental studies like (Wang et al. 2017; Yao et al. 2018; Wang et al. 2022) dealing with the in uence of transient temperature conditions on the shear and creep behavior of frozen soils. According to Andersland and Ladanyi (2003), decreas- ing temperatures can result in (partial) refreezing of the developed cracks in the frozen soil. As a result, the creep rates sharply decrease. In contrast, increasing temperatures lead to ice softening, a possible increase in the unfrozen water content, and an overall reduction of shear and creep strength (Wang et al. 2022). The proposed concept may be further adjusted to consider these highly non-linear characteristics. At this stage, we suggest limiting the use of the constitutive model for steady-state temperature conditions. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 29 7 Conclusion This study deals with the in uence of varying stress states on the mechanical behavior of frozen soils at constant temperatures. The comparison between uniaxial creep tests with both single-stage loading and multi-stage loading shows that the minimum axial strain rate _mis mostly independent of the loading history. In contrast, the lifetime t mof the frozen soil depends on the loading history. These ndings are consistent with previous experimental studies found in the literature. Moreover, the introduction of the transformed creep time tenables the conversion of multi-stage creep tests into equivalent single-stage creep tests. Here, the evolution of the normalized axial strain rate _ 1=_mover normalized time t =tmfor single-stage as well as the normalized axial strain rate _ 1=_mover normalized transformed time t;Mj=t;Mj mfor multi- stage creep tests converge and are consequently comparable with each other. Hence, this crucial relationship of the constitutive model for frozen soils proposed by Cudmani et al. (2022) to describe the rate-, stress-, and temperature-dependent mechanical behavior of frozen soils is not only valid for predominantly monotonic (single-stage) loading but also for multi-stage loading. Therefore, we extend the constitutive model by taking into ac- count multi-stage loading based on our proposed determination of the normalized axial strain rate _ 1=_mover normalized transformed time t;Mj=t;Mj m. The extended model, designated by the acronym EVPFROZEN, captures the equivalent creep time after a changing stress state by coupling the transformed creep time tto the previous loading history. Comparing the simulations with our own experiments and com- prehensive data from the literature successfully completes the rst step in validating the extended model. The model validation includes a large number of incremental axial stress changes N 1;step, and incremental axial stress increases  1by simultaneously covering a wide range of total axial stress states 1. Considering the in uence of varying stress states on the shear and creep behavior of frozen soils enables a more realistic and economically optimized design of ground freezing applications due to the resulting prolonged lifetime of the frozen soil. In fact, it is possible to take into account di erent construction stages and techniques and, therefore, di erent stress levels of the frozen soil body. In summary, this study has signi cantly improved and further validated an existing nu- merical tool for assessing the stability and deformations of frozen soils in engineering applications for single-stage and multi-stage loading. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 30 Competing interests The authors declare there are no competing interests. Data availability Some or all data, models, or code that support the ndings of this study are available from the corresponding author upon reasonable request. Funding statement The authors declare no speci c funding for this work. Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 31 Notation Total strain tensor elElastic strain tensor vViscous strain tensor _Total strain rate tensor _elElastic strain rate tensor _vViscous strain rate tensor 1Axial strain _1Axial strain rate mAverage axial strain at the lifetime of the frozen soil _ Reference strain rate _mMinimum creep rate / Minimum axial strain rate _mMinimum creep rate tensor tmLifetime of the frozen soil after which tertiary creep starts  Temperature in Celsius t Actual testing time / Global total time t0Starting creep time tTransformed creep time p Roscoe's invariant representing mean pressure q Roscoe's invariant representing deviator pressure 1Axial stress cUniaxial compression strength crEquivalent uniaxial creep strength after Cudmani et al. (2022) Stress tensor _Stress rate tensor sDeviatoric part of stress tensor ITrace of stress N1Number of load steps in an uniaxial multi-stage creep test LFourth-order isotropic elastic sti ness tensor E Young's modulus Poisson's ratio Lode angle dDry density of the frozen soil sample w Gravimetric water content SrDegree of saturation e Void ratio Schindler et. al., 'Multi-stage creep behavior of frozen granular soils' 32