Simulating ice segregation and thaw consolidation in permafrost environments with the CryoGrid community model

. The ground ice content in cold environments influences the permafrost thermal regime and the thaw trajectories in a warming climate, especially for very ice-rich soils. Despite their importance, the amount and distribution of ground ice are often unknown due to lacking field observations. Hence, modelling the thawing of ice-rich permafrost soils and associated thermokarst is challenging as ground ice content has to be prescribed in the model set-up. In this study, we present a model scheme, which is capable of forming segregated ice during a model spin-up together with associated ground heave. It provides 5 the option to add a constant sedimentation rate throughout the simulation. Besides ice segregation, it can represent thaw consolidation processes and ground subsidence under a warming climate. The computation is based on soil mechanical processes, soil hydrology by Richards equation and soil freezing characteristics. The code is implemented in the CryoGrid community model (version 1.0), a modular land surface model for simulations of the ground thermal regime. The simulation of ice segregation and thaw consolidation with the new model scheme allows us to analyze the evolution 10 of ground ice content in both space and time. To do so, we use climate data from two contrasting permafrost sites to run the simulations. Several influencing factors are identified, which control the formation and thaw of segregated ice. (i) Model results show that high temperature gradients in the soil as well as moist conditions support the formation of segregated ice. (ii) We find that ice segregation increases in fine-grained soils and that especially organic-rich sediments enhance the process. (iii) Applying external loads suppresses ice segregation and speeds up thaw consolidation. (iv) Sedimentation leads to a rise of the 15 ground surface and the formation of an ice-enriched layer whose thickness increases with sedimentation time. We conclude that the new model scheme is a step forward to improve the description of ground ice distributions in permafrost models and can contribute towards the understanding of ice segregation and thaw consolidation in permafrost environments under changing climatic conditions.


Introduction
Permafrost underlies an area of approximately 14 million km 2 , which corresponds to around 15 % of the exposed Northern Hemisphere or 11 % of the global exposed surface (Obu et al., 2019;Obu, 2021).Frozen soils with a high organic content are 1 https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.widespread (Hugelius et al., 2014) and these permafrost environments largely evolved in the last glacial cycle, accumulating massive amounts of soil organic carbon.Especially in Sibiria and Alaska, they are often associated with a high content in excess ground ice (Saito et al., 2020;Schirrmeister et al., 2011).The climatic, sedimentological, geomorphological and ecological conditions in the past govern today's distribution of ground ice (Gilbert et al., 2016).This is especially true for syngenetic permafrost, which builds during sedimentation or formation of organic material in cold environments and is typically ice-rich.
In contrast, epigenetic permafrost forms with a time lag to the sedimentation processes in the ground, and is often characterized by ice-poor soils (Gilbert et al., 2016).Ground ice can be present as pore ice or excess ice, which can occur as relict ice, wedge ice or segregated ice (O'Neill et al., 2019).
Segregated ice forms as discrete ice lenses or layers through cryosuction (French and Shur, 2010) and is typically associated with ground heave as the ice content in the ground increases (Miller, 1972;Taber, 1929).It can build up in epigenetic permafrost when unfrozen soil water from the active layer is attracted towards the freezing front (Guodong, 1983;Mackay, 1983;Taber, 1929), accumulating ice near the permafrost table.Besides, ice segregation can occur at the base of the permafrost, as water is drawn from the soil below as permafrost is forming.These ice-enriched layers are widespread in permafrost environments (French and Shur, 2010) as for example in Canada (O'Neill et al., 2019).However, the evolution of syngenetic permafrost can result in enhanced ice formation.Accumulation of organic material as well as sedimentation in alluvial, eolian or hillslope settings can lead to a rise in the permafrost table and hence a growth of segregated ice (Guodong, 1983;French and Shur, 2010).
In this context, segregated ice can form also together with syngenetic ice wedge growth, forming ice lenses within polygonal permafrost as observed in Yedoma deposits (Schirrmeister et al., 2013).
The ground ice content and its distribution strongly determines the sensitivity of permafrost to thaw (Jorgenson et al., 2010;Nitzbon et al., 2019).If ice-rich layers are present in the soil column, increasing permafrost temperatures and deepening of the active layer can lead to melting of the excess ice, resulting in thermokarst and subsequent ground subsidence (Farquharson et al., 2019;Kokelj and Jorgenson, 2013;Nitzbon et al., 2019).In consequence, substantial geomorphological changes reshape the landscape, manifested in the formation of lakes, thaw slumps, gullies and the transformation of low-centered to high-centered polygons (Kokelj and Jorgenson, 2013;Liljedahl et al., 2016;Nitzbon et al., 2019).These processes could contribute to accelerating the mobilization of permafrost carbon, which may further increase atmospheric carbon concentrations (Miner et al., 2022;Schuur et al., 2008).Furthermore, ground ice controls the hydrologic and mechanical properties of the soil by reducing permeability and increasing the mechanical strength (Painter and Karra, 2014).These parameters control the structural stability of the ground and potentially endanger human infrastructure and settlements upon thawing (Hjort et al., 2022;Schneider von Deimling et al., 2021).
As about 20 % of the permafrost in the Northern Hemisphere is prone to thermokarst processes (Olefeldt et al., 2016), it is highly important to develop models representing the formation and melt of (excess) ground ice, such as segregated ice, to better understand Arctic landscape evolution.This will improve our capabilities to assess soil stability under transient climate conditions, which is essential for the sustainable design of infrastructure in permafrost environments (Dumais and Konrad, 2018;Painter and Karra, 2014).Modelling of ground consolidation requires implementing soil mechanical processes in addition to heat and water transfer (Dumais and Konrad, 2018).Thaw consolidation has been the focus of model development https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.for many decades (Morgenstern and Nixon, 1971;Nixon and Morgenstern, 1973;Sykes et al., 1974;Foriero and Ladanyi, 1995;Dumais and Konrad, 2018), and different approaches have been presented for modelling ice segregation (Fu et al., 2022;Fisher et al., 2020;Lacelle et al., 2022).However, these processes are typically not implemented in land surface models and simulating the long-term evolution of segregated ice has not been performed yet (Fu et al., 2022).
As a consequence, previous models targeting ground ice thaw and thermokarst require excess ice distributions prescribed from field observations, which makes applications at sites without ground ice data challenging.In this study, we demonstrate how distribution of segregated ice in the ground can be obtained through a climate-dependent spin-up procedure, which aims at reproducing the evolution of ground ice stocks.We present a new model approach, which is capable of coherently simulating both ice segregation and thaw consolidation, based on heat and water transfer as well as soil mechanical processes.The model code is implemented within the framework of the modular CryoGrid community model (version 1.0, referred to as "CryoGrid" from here on), a land surface model for permafrost applications (Westermann et al., 2022).The main objectives of this study are the following: -We demonstrate a proof of concept for simulating ice segregation and thaw consolidation with associated ground heave and subsidence.To do so, we run various model scenarios with climate data from two contrasting permafrost sites.
-We evaluate controlling factors for ice segregation by model runs with different climatic conditions (by applying different forcing data sets), the soil type (by using different grain sizes and compositions) and external loads.
-We investigate the capability of the model to simulate long-term ice segregation under a constant sedimentation regime.

Model description
In this work, we extend the capabilities of the CryoGrid community model (Westermann et al., 2022) with a representation of soil mechanical processes.The CryoGrid community model is a modular framework for simulating the permafrost thermal state and water/ice balance, which allows combining different process representations for ground and snow to model a vertical column in the best possible way.Here, we describe a new "stratigraphy class" (Westermann et al., 2022) in CryoGrid called GROUND_freezeC_RichardsEq_seb_pressure, a fully fledged process model, which is applied to a one-dimensional domain in the ground.It is based on the already existing stratigraphy class GROUND_freezeC_RichardsEq_seb (Westermann et al., 2022) and inherits all its functionalities.The new class can be vertically stacked with existing stratigraphy classes, in particular classes representing the seasonal snow cover.In the following, we summarize the main aspects of CryoGrid relevant for this work, as well as the model physics of the base stratigraphy class.A detailed description of all aspects of the model is provided in Westermann et al. (2022) and a summary of the employed stratigraphy classes can be found in Suppl. 1.A soil freezing characteristics is implemented in GROUND_freezeC_RichardsEq_seb and describes the relationship between the ground temperature and the unfrozen water content.We use an improved parametrization for the phase partitioning of water/ice in frozen soils (Painter and Karra, 2014;Stuenzi et al., 2021) to distinguish between liquid water and ice.To do so, we calculate the matric potential for unfrozen conditions ψ 0 [m] and based on that, the matric potential in a frozen state ψ [m], followed by the water content in the soil.A detailed description of the approach can be found in Westermann et al. (2022).

Subsurface properties in
The water balance in GROUND_freezeC_RichardsEq_seb is based on vertical water flow j v w [m s −1 ] controlled by the Richards equation (Richards, 1931): with ψ [m] being the matric potential, z [m] the vertical coordinate and K w [m s −1 ] the hydraulic conductivity.The latter is computed according to Van Genuchten (1980) as well as Hansson et al. (2004) to take into account ground ice blocking waterfilled pores.To do so, a depth stratigraphy of the permeability k w [m 2 ] is defined by the user (see Westermann et al., 2022).
Vertical water fluxes between unsaturated grid cells are initiated through gradients in matric potentials as well as gravitational potentials.In saturated unfrozen regimes, water flow by gravity occurs according to Darcy's law.
In addition to vertical water fluxes, the CryoGrid community model provides the option for lateral water transport, e.g. through overland flow (Westermann et al., 2022), which is induced according to the Gauckler-Manning equation (Gauckler, 1867;Manning et al., 1890).Depending on the Gauckler-Manning coefficient and the local ground surface gradient, standing surface water is removed from the system.Other lateral water fluxes such as subsurface drainage are not included in the model setup.
The seasonal snow cover is implemented by applying the snow class SNOW_crocus2_bucketW_seb, which includes snow microphysics, the surface energy balance and snow hydrology.It is based on the Crocus snow scheme (Vionnet et al., 2012) and described in Westermann et al. (2022) and Zweigel et al. (2021).Excess water from snow melt is transferred to the uppermost grid cell of the new GROUND class which can represent standing surface water by a dedicated pool from which it can either run off laterally through Gauckler-Manning flow or infiltrate into the ground according to soil properties and moisture conditions.

Soil mechanical processes
We build on established concepts for ground consolidation, following Dumais and Konrad (2018).This includes the calculation of (i) the total stress, (ii) the effective stress, (iii) the compaction of the soil and the water fluxes (see section 2.1.2),which are updated in every time step.To do so, a set of additional state variables is necessary, which can be found in table 1.A schematic illustration of the stresses acting on the soil column is depicted in Fig. 1.
External stress is caused by additional loads on the ground surface such as buildings or roads and acts with the same value on all grid cells in the soil column.In contrast, the vertical geostatic stress is caused by the mass of the overburden soil and consequently increases with depth.In our model, we calculate the total stress at the midpoint of each grid cell.Therefore, we consider the weight of all grid cells above in addition to half of the weight of the grid cell itself.The geostatic stress can be determined by dividing the weight of the soil, W [kg m s −2 ], by the area A [m 2 ]: with Θ i [m 3 ] being the volume, and ρ i [kg m −3 ] the density of a grid cell i.As grid cells consist of mineral, organic and water/ice constituents, the density is calculated as the arithmetic average of the constituent densities weighted by their volumetric fractions.(4) To compute the pore water pressure, the buoyancy effect has to be accounted for.When the soil is at or near saturation, the porewater results in a reduction the stress on the soil matrix, as the density of each component is reduced by the density of water.For saturated conditions, we assume that the buoyancy effect is fully effective.If the saturation is less than 50 %, we make the assumption that the total stress is carried exclusively by the soil matrix and the pore water pressure is set to zero.To facilitate a continuous transition of pore water pressure between saturated (100 %) and dry (< 50 %) conditions we apply a linear interpolation.While this is a coarse approximation, neglecting the increase in effective stress where the soil is saturated due to capillarity, it allows us to calculate the effective stress in the soil column avoiding step changes between the two regimes.
(iii) When the soil compresses or relaxes, the porosity changes, which must be equilibrated by either fluxes of air or water.
In this work, we assume that air fluxes occur instantly, i.e. for unsaturated conditions, the soil is compacted to its equilibrium position after each time step.The compressibility of a soil can be expressed with a linear relationship between void ratio e [-] (defined as the ratio of pore volume to grain volume) and the decadic logarithm of the effective stress σ ′ z [Pa] (Murthy, 2002) as illustrated in fig. 2. With increasing effective stress, the void ratio decreases linearly, with steeper curves indicating a higher compressibility of the soil.The slope of the curve is given by the compression index C c [-], which can be expressed as with the initial void ratio e 0 [-] and the residual stress σ ′ 0 [Pa] staying constant during the simulation.The latter is defined as the effective stress, where no consolidation occurs upon thawing in undrained conditions, i.e. the effective stress that can be maintained by the soil skeleton (Nixon, 1973;Dumais and Konrad, 2018).The effective stress is calculated in Eq. ( 4), so that the void ratio can be computed from Eq. ( 5): Knowing the void ratio e [-] of the soil, the porosity ϕ [-] of a grid cell can be calculated by the following equation: With changing porosity, the volume of the affected grid cell has to be adapted.As the area A [m 2 ] is fixed throughout the entire simulation, the volume can be expressed as the thickness of a grid cell d [m]:  (2022) and described in section 2.1.(i) Gradients in matric potential induce a vertical water flow between grid cells.This process is especially active in the upper soil layers, where rainfall, snow melt water and evaporation influence the soil water content.For example, when upper grid cells loose soil water due to evaporation, the matric potentials of the affected grid cell decrease, leading to suction of soil water from the grid cells below and inducing upward water fluxes.(ii) The gravitational potential always leads to water fluxes from the top to the bottom of the ground column, e.g. the infiltration of rainfall or melt water.The highest value occurs in the uppermost grid cell and with progressing depth, the value decreases according to the thickness of the grid cells.For saturated conditions, the gravitational potential in the aquifer matches the hydrostatic potential, i.e. the value at the top of the aquifer and as a result, no water fluxes occur within the aquifer due to the gravitational potential.
It is important to highlight that water fluxes in unsaturated soil do not directly lead to a change in grid cell size in the presented model scheme, as changing water content are replaced by the same volume of air.However, the resulting changes in water content affect the stress level in the soil and therefore the soil mechanical responses.
If the soil is unfrozen and saturated, downward water fluxes occur in case unsaturated grid cells are present in deeper layers, which attract soil water through both gravitational and matric potential similar to GROUND_freezeC_RichardsEq_seb.In case of an aquifer, where the entire soil column below is saturated, no water is exchanged between the saturated grid cells.
When external loads are applied or sedimentation takes place, the additional stress is taken up by the soil water in a first step, leading to excess pore water pressures u e [Pa]: with σ z [Pa] being the total stress and σ ′ z [Pa] being the effective stress, carried by the soil skeleton.Over time, the excess pore water pressure is reduced by consolidation, which results in water fluxes j v w,ue [m s −1 ] away from the affected grid cell The water flux is added to the fluxes calculated based on the Richards equation (eq.1).The consolidation continues until the excess pore water pressure reaches a value of zero, while at the same time the effective stress is increased.As the finite hydraulic conductivity of the soil controls the water fluxes, this process takes a certain amount of time.Therefore, saturated soils cannot respond immediately to changing stresses: Before the soil skeleton can take the additional stress and compress, the water has to flow out of the pores.The opposite effect takes place, when a formerly dry soil becomes saturated and the increased pore water pressures reduce the effective stress on the soil skeleton.
After calculating the water fluxes in the soil column, the change in water content for each grid cell can be derived.For unsaturated conditions, changing water content is replaced by air.For saturated conditions, no air inflow is possible and changes A reduction/increase in water content results in a change in effective stress and thus in compression/swelling of the soil (cf. Eq. 5) and consequently a reduced/increased thickness of the grid cell.Knowing the grid cell thickness, porosity and void ratio can be derived with Eq. ( 7) and Eq. ( 8).

Segregated ice
When the soil freezes, water fluxes are significantly smaller than in unfrozen conditions due to reduced liquid water contents and hydraulic conductivity.However, the remaining soil water is still partly mobile.This is mainly driven by matric potentials (see section 2.1), which reach considerably negative values for ground temperatures below zero degrees, resulting in an attraction of soil water towards the freezing front.In contrast to the stratigraphy class GROUND_freezeC_RichardsEq_seb, the new model scheme enables water fluxes also into saturated grid cells.In this case, the grid cell gradually expands and the void ratio increases.When its void ratio reaches the initial void ratio, the entire weight of the overlying soil is compensated by the matric water redistribution and leads to a stable situation.Upon thawing, excess pore water pressure leads to water fluxes out of the affected grid cell and thaw consolidation takes place.
We highlight that the term "segregated ice" in our model is based on a different definition than the term "excess ice" in previous version of the CryoGrid model.Westermann et al. (2016) and Westermann et al. (2022) distinguish between an incompressible soil phase including pore water/ice, and an excess ice phase, which becomes mobile upon melting and thus leads to a shrinking of the grid cell.In contrast, the presented stratigraphy class in this study defines segregated ice as the ice volume exceeding the pore space, which the soil would have in a compressed state without excess pore water pressures under the applied overburden pressure (illustrated in fig.2).This implies that for a given soil a small additional ice content at great depth is be defined as segregated ice, while the same volumetric ice content could be present as regular pore ice at lower depths.
We note that the new model scheme can only represent segregated ice and that the formation of other forms of excess ice such as wedge ice cannot be accounted for.
Segregated ice that is present at void ratios smaller than the initial void ratio is considered to be bound within the pore structure of the soil and the partitioning of ice and water thus follows the soil freezing characteristics of the given soil type (Westermann et al., 2022).In contrast, any further segregated ice exceeding the initial void ratio is assumed to be "free" water, which freezes at T = 0 °C.

Sedimentation
The model provides the option to include sedimentation by organic and mineral matter in the simulations (stratigraphy class GROUND_freezeC_RichardsEqW_seb_pressure_sedimentation).The user can define the sedimentation rate and the depth at which the material is added.In the current version, the added material is assumed to have the same composition as the already present soil at the time of sedimentation.However, this is not a principle limitation, and sedimentation with user-defined properties could be implemented in a straight-forward way.Furthermore, sedimentation only occurs for positive ground temperatures in the current model version.Water is not added during sedimentation, but the newly created pore space can be filled through water flow from neighboring grid cells or rain/snow melt (section 2.1.2).When the grid cell exceeds a certain target thickness, the grid cell is split into two to maintain a set minimum resolution in the soil column.Intensive variables (e.g.temperature and porosity) are inherited for both grid cells, while extensive parameters (e.g.bulk mineral and organic content, Westermann et al., 2022) are divided between them while maintaining the proportions during the split.

Field sites and forcing data
To demonstrate the capabilities of the new model, we perform sensitivity studies for two different permafrost sites with strongly different climate conditions: Samoylov Island (located in northern Siberia) represents a setting with cold continuous permafrost, while the ancillary site Bayelva (located on Svalbard) is characterized by warm and maritime, but still continuous permafrost.
Samoylov Island is located in northeastern Siberia in the Lena River delta (72°22' N, 126°28' E).Monitoring of soil temperatures and meteorological conditions has been performed at the research site since 1998 (Boike et al., 2013(Boike et al., , 2019)) Boike et al., 2013;Gouttevin et al., 2018;Langer et al., 2011b).Samoylov is situated in continuous permafrost, which reaches depths between 400 and 600 m (Grigoriev, 1960).East of the research station, the landscape is shaped by polygonal tundra with ice wedge polygons dominating the upper meters of the fine-grained and organic-rich soils.
There are both low-and high-centered polygons, and the local microtopography is largely structured by sequences of polygon troughs, rims, and centers (Boike et al., 2013).The ground temperature was measured with -9.0 to -7.9 °C in a depth of 20.75 m between 2007 and 2016 ( GTN-P, 2018).With these records, the field site shows one of the strongest warming of permafrost within 123 globally distributed boreholes (Biskaborn et al., 2019).The ground temperature for validation was taken beneath a polygon center close to the research station (Boike et al., 2013).We use the same forcing data set as Nitzbon et al. ( 2020

Model parameters and initialization
The soil column is described by two different ground classes, which are treated differently during simulation.Below 9 m depth, soil mechanical processes are neglected due to frozen conditions and high total stress during the entire simulation.Therefore, we apply the stratigraphy class GROUND_freeW_seb, which operates without soil mechanical processes and is described in Westermann et al. (2022).Above 9 m depth, we use the new stratigraphy class GROUND_freezeC_RichardsEq_seb_pressure as presented in this study.We chose the snow class SNOW_crocus2_bucketW_seb, based in the Crocus snow model (Vionnet et al., 2012), to represent seasonal snow cover with parameters adapted as in Westermann et al. (2022).
For model validation, we set a soil stratigraphy as described for the center of a low-centred polygon in Holocene deposits in Samoylov in Nitzbon et al. (2020) and Westermann et al. (2016).The key difference is that no excess ice is assumed in the beginning of the simulation as it inherently generates by the new model scheme.To represent the entire system of polygonal tundra, laterally coupled simulations (e.g.Aas et al., 2019;Nitzbon et al., 2019) would be required, but a one-dimensional soil column is sufficient for the purpose of this proof of concept study.
To analyze the performance and sensitivity of our model, we simplify the stratigraphy to reduce uncertainties due to the model setup.Hereby, we let the upper 0.15 m unchanged to account for the insulating effect of the organic-rich layer on top.
The undecomposed organic material features coarse pores with unknown values for the van Genuchten parameters n and α.
Therefore, it is phenomenologically set to the properties of coarse-grained (sandy) soil (table 2), which has a similiar retention characteristic.Between 0.15 and 9 m depth, we set homogeneous soil properties, consisting of soil type, volumetric fractions of mineral and organic content, saturation, initial void ratio, residual stress and compression index.The chosen values for the different settings can be found in table 2. The volumetric fractions of the mineral and organic content are taken from Westermann et al. (2016).The initial void ratio is calculated from the corresponding porosity.The soil is assumed to be fully saturated at the beginning of the simulations, only the surface layer has a saturation of 50 %.As the residual stress is unknown, we choose a value of 1000 Pa for all soils, except for organic-rich layers, which we set to 100 Pa, so that the corresponding thawed void ratio stays meaningful and realistic.A comparable range in values of residual stress has been applied by Dumais and Konrad (2018) and Dumais (2019).We estimate the compression index based on literature value (Gudehus, 1981), as well as the Van Genuchten parameters α and n (Dall' Amico et al., 2011;Gnatowski et al., 2010;Van Genuchten, 1980).It should be noted that the mechanical and hydrological parameters are soil type dependent, but are adjusted independently in the model setup.
Other model parameters are set in accordance with previous model studies in the same area.During the snow-free period, the albedo of the surface is set to 0.2 and an aerodynamic roughness length of 0.001 m is assumed for the ground (Westermann et al., 2016).We use the snow class based on the Crocus snow model (Vionnet et al., 2012;Zweigel et al., 2021), but increase the maximum wind slab density to 500 kg m -3 , in line with Barrere et al. (2017) and Royer et al. (2021).Furthermore, we define a maximum snow depth to account for snow ablation by wind drift as stated in section 2.4 (Westermann et al., 2016).
We  temperature profile corresponding to a geothermal heat flux of 50 mW m −2 .Furthermore, the soil column in the uppermost 9 m is compressed according to the mass of the overlying layers before the start of the simulations.For this pre-compaction, we assume an unfrozen conditions with the saturation given in table 2. Sedimentation is assigned to the grid cell below the organic-rich surface layer to assume a constant thickness of the layer influenced by vegetation.Consequently, the organic-rich surface layer is lifted with sedimentation, but stays unchanged regarding its soil properties.We set the sedimentation rate to 0.002 m per year, active only for positive ground temperatures conditions, which corresponds to an effective sedimentation rate of 0.3 m per 1000 years.While this is a simplified case, not trying to mimic a specific sedimentation process, it allows us to analyze the effect of material deposition on the accumulation and thaw of segregated ice.

Model scenarios
We set up different model scenarios by varying both forcing data and model parameters ( (i) We use two different forcing data sets to simulate conditions of permafrost in a cold environment (Samoylov, Sibiria, model run S-clay) and permafrost closer to the thaw threshold (Bayelva, Svalbard, model run B-clay).We emphasize that the ground stratigraphy designed for Samoylov is not representative for the Bayelva field site, so that modelled ground temperatures cannot be compared to observations.All model scenarios that are forced by the Samoylov forcing include a spin-up by repeating ten times the years 1960 to 1969, while the period 1980-1989 was employed for Bayelva.Furthermore, we include one model run with the Samoylov forcing with rainfall reduced to 50 % to simulate drier conditions (S-clay-rain50).The spin-up of 100 years allows for a stable ground temperature profile, but ice segregation can still continue after the spin-up period.(ii) A comparison of the model runs S-sand, S-silt, S-clay and S-peat show the effect of different soil types.Hereby, we change the soil type in the simplified stratigraphy by adapting the mineral and organic content, the initial porosity, the residual stress, the compression index as well as α and n.In contrast to the other soil types, peat is characterized by a large porosity, low density and large water and organic matter contents.(iii) Further model scenarios include an external load of 5 kPa during the entire simulation (S-clay-load5) and from 1980, when the annual downward thawing in summer deepens (S-clay-load5-1980).It corresponds to an about 40 cm thick gravel pad and allows us to discuss the influence of external loads both on ice segregation and on thaw consolidation.The load is applied instantaneously and does not change the surface properties (as a real gravel pad would certainly do), so that energy and water transfer occur as without the external load.(iv) Model scenarios S-clay-sed100, S-clay-sed500 and S-clay-sed1000 include sedimentation with a rate of 0.002 m per year over spin-up periods of 100, 500 and 1000 years, respectively.
We set a threshold for the snow depth to account for observed snow ablation due to wind drift following the approach of Westermann et al. (2016).As the snow in the polygon centers is protected by the surrounding rims, we set a value of 0.45 m, which is in line with measured snow depths in polygons centers on Samoylov Island (Boike et al., 2013).For validation, we use temperature measurements in a polygon center.However, these measurements are not conducted in the middle of Samoylov Island but close to the edge of the island where the wind drift is stronger and the height of the rims above the polygon center is reduced.Therefore, we reduce the maximum snow depth for the validation run to 0.20 m.Surface water discharges by overland flow are based on the Gauckler-Manning equation (Westermann et al., 2022).As the terrain at the field site on Samoylov Island is flat, we set a small gradient of 0.1 m km -1 for the model runs that are forced with data from this location..In contrast, the Bayelva field site is situated on gently sloping terrain, so that we increased the gradient to 1 m km -1 .

Results
In the following, we present results of the model scenarios described in section 2.4 and summarized in table 3.

Model validation
We compare the validation run S-val to published in situ temperature data from 0.05 m and 0.40 m depth of a polygon center (Boike et al., 2013;Langer et al., 2011a, b;Westermann et al., 2016).Figure 4 shows the mean daily ground temperatures for both measured and modelled data.The ground temperatures are well reproduced by the model and the seasonality of the measured signal is captured well.However, modelled ground temperatures are slightly too warm during summer in both 0.05 m and 0.40 m depth.In fall, modelled ground temperatures decrease earlier than the measurements and the zero curtain is underestimated.The same effect occurs in the simulations made by Westermann et al. (2016).A possible explanation is a delayed start of the snow season in the forcing data (Westermann et al., 2016)  in the field (Boike et al., 2013).This is explained by uncertainties in distinguishing between snow-and rainfall in the reanalysis product (Westermann et al., 2016) which likely causes the cold bias in our simulations during fall.
The results of the validation run for Samoylov and the validation performed by Westermann et al. (2022) for Bayelva suggest that the model system including CryoGrid and model forcing can capture the main characteristics of the ground thermal regime at the two sites.In the following, we use this model setup to perform a sensitivity analysis for ice segregation and thaw consolidation with the model scenarios described in section 2.4.

Formation of segregated ice and thaw consolidation
The model run S-clay using a simplified stratigraphy of soil type clay shows the formation of segregated ice in the relatively cold period until the 1980s, as well as thaw consolidation during the past decades and the future forcing.During the spin-up, mean annual ground temperatures M AGT are modelled with values around -9°C at a depth of 5 m depth.The relatively cold period is then followed by a significant rise of M AGT .At the end of the 21st century, M AGT at 5 m depth reach temperatures close to the 0 °C (Fig. 5).The increasing ground temperatures are also reflected in the depth of the permafrost table.While the active layer is less than 0.70 m during the spin-up, a deepening can be observed since the 1980s (Fig. 5).In the early 2080s, the seasonally frozen ground does not reach the permafrost table anymore.A continuously unfrozen zone with positive ground temperatures throughout the year (talik) develops, indicating the degradation of the permafrost.The influence of the changing climatic conditions is also apparent in the formation and thawing of the ground ice (Fig. 6).
Ground ice accumulates during fall refreezing at the freeze front and the upper part of the permafrost due to differences in matric potential, which draw liquid water from the active layer into the frozen ground.While the model run is initialized with saturated conditions below 15 cm depth, the soil water content adapts to 70-80 % saturation in the active layer in the first year and remains at these values during spin-up.The process of ice segregation continues as long as thermal conditions allow a mobilisation of soil water, dependent on the freezing characteristics of the soil.The increased ice content is preserved at depths to which the active layer of the following years does not reach.As a consequence, S-clay builds up 0.034 m of segregated ice until the 1980s.Simultaneously, the ground surface heaves over the years by 0.044 m.Hence, most of the ground heave can be explained by ice segregation, while only a small amount is related to soil mechanical processes in the active layer, i.e. swelling and shrinking of the soil due to changes in soil water content caused by precipitation and evapotranspiration.In the following decades, warmer climate conditions lead to a deepening of the active layer.The segregated ice formed during spin-up starts to melt from the top to the bottom and thaw consolidation takes place.The resulting melt water flows upwards in the active layer and the soil subsides again by 0.082 m until 2100.As only 0.034 m can be explained by melting of the segregated ice, the remaining subsidence is due to a drying of the soil column for the deeper soil layers at the end of the century (years 2080 to 2100, see supplement B), increasing the effective stress.
When the active layer stabilizes for several years at a specific depth, new segregated ice is formed at the upper part of the permafrost (Fig. 6), which can lead to temporary surface heave even though the ground surface elevation decreases in the long term.However, the changes in ground surface elevation are too large to be only explained by new ice segregation.Highly

Sensitivity towards climate conditions
Model simulations with different forcing data sets suggest that ice segregation is highly dependent on the climatic conditions.
Figure 7 shows the formation of segregated ice for the model runs S-clay (Samoylov forcing data), S-clay-rain50 (Samoylov forcing data with 50 % rainfall) and B-clay (Bayelva forcing data).The reduction of rainfall in S-clay-rain50 leads to less segregated ice during the spin-up on average with 0.025 m compared to 0.034 m in S-clay.The reduced rainfall results in drier conditions in the active layer (see supplement B), so that less water is available for the formation of segregated ice.
Furthermore, the increased negative matric potential in the drier active layer counteracts the matric potential responsible for ice segregation.In model scenario S-clay-rain50, the permafrost does not degrade to the end of the 21st century and no talik develops as less energy is consumed for thawing and melting of soil water under drier conditions.
When the same model set-up is run with forcing data from Bayelva (B-clay), which represents warmer permafrost under moist conditions, we obtain slightly less segregated ice with 0.028 m in the 1908s compared to 0.034 m with forcing data from Samoylov (S-clay).As the results discussed above have shown that moist conditions lead to more segregated ice, one would expect higher values for B-clay.However, the higher ground temperatures and consequently smaller temperature gradients in the ground during fall refreezing (see supplement B) seem to compensate the effect of the soil moisture.For future forcing, when ground temperatures in Samoylov increase as well, B-clay builds up more segregated ice in deeper soil layers as the active layer stabilizes for longer time periods.Varying the climatic forcing shows, that the model results depend on both soil moisture and temperature gradients in the ground.

Sensitivity towards soil type
The soil type has a strong influence on ice segregation and thaw consolidation (Fig. 8).The formation of segregated ice in the model run with the sand layer (S-sand) is negligible and no significant ground heave occurs.Yet, the other soil types form segregated ice and reveal a distinct change in ground surface position.The model run with the silt stratigraphy (S-silt) builds up a maximum of 0.015 m segregated ice until the 1980s and heaves by 0.018 m.The reference run with the clay stratigraphy (S-clay) has a slightly stronger reaction with a build-up of segregated ice of 0.034 m and a ground heave of 0.044 m (see section 3.2).In contrast to these reactions on a centimeter scale, the model run with the peat stratigraphy (S-peat) yields much more ice-rich ground, with a maximum of accumulated segregated ice of 0.784 m and ground heaving by by 0.815 m (Fig. 9).Here, the ground heave is predominantly caused by ice segregation, which is especially effective for peat due to large differences in matric potential and small vertical geostatic stress due to the low density of the peat, while soil swelling accounts for only a small part.
Upon active layer deepening, the segregated ice melts.The results show that a larger portion of the segregated ice melts between 1990 and 2000 in the silt and clay stratigraphy, while this happens between 2000 and 2010 for peat.A possible explanation for this effect is the lower thermal conductivity of the organic material (0.25 W m -1 K -1 ) compared to the mineral phase (3.0 W m -1 K -1 ) as well as the higher amount of segregated ice that consumes energy upon melting.In the future projection, ice segregation occurs in deeper soil layers (at the top of the permafrost) even though it is typically less intense due to higher total stress and smaller temperature gradients.Therefore, segregated ice does not disappear completely for model runs S-silt, S-clay and S-peat, but potentially slows permafrost degradation as energy is needed for melting of the newly formed segregated ice in deeper soil layers.With decreasing particle diameter, the soil builds up more segregated ice under equal conditions.Peat soils form substantially more ground ice than mineral soils.3. The results show a significant formation of segregated ice below the permafrost table under a constant sedimentation rate.The ground surface rises due to both accumulated material and the build-up of segregated ice.With the deepening of the active layer, the segregated ice is melted and thaw consolidation takes place with associated ground subsidence.

Ice segregation and thaw consolidation in CryoGrid
Previous work with CryoGrid (Westermann et al., 2016(Westermann et al., , 2022) ) does not consider soil mechanical processes when modelling the thermal regime of permafrost environments.This corresponds to the assumption of a time-constant soil porosity, so that e.g. the effect of ground loading cannot be represented.However, CryoGrid already provides a stratigraphy class to include excess ice (Westermann et al., 2022).Earlier studies model ground subsidence upon thawing with CryoGrid 3 in the Samoylov area in a 1D setting (Westermann et al., 2016) and simulate the thaw-induced transition from low-to high-center polygons with laterally coupled tiles (Aas et al., 2019;Nitzbon et al., 2019).Here, a fixed amount and stratigraphy of excess ice is defined at the beginning of the simulation, while it is not possible to account for climate-dependent build-up of excess ice (as demonstrated for segregation ice in this study).Upon thawing, excess water is mobilized and routed away without considering changing hydraulic properties of the ground due to the fixed porosity.
In contrast, the new model scheme can accumulate and melt segregated ice depending on defined soil mechanical properties.
The porosity is adjusted in each time step according to stress conditions in the ground and the soil characteristics (initial void ratio, residual stress, compression index).Therefore it evolves in time in response to the model forcing, which can lead to for example saturated ground conditions and thus soil swelling.Furthermore, the ice content decreases with increasing confining pressure due to stronger compaction.
In particular, the new model scheme can form segregated ice according to the climatic forcing and soil properties.Observation sites in permafrost environments are scarce and ground ice contents are poorly constrained (Cai et al., 2020).Therefore, model schemes like those of Nitzbon et al. (2019Nitzbon et al. ( , 2020)); Westermann et al. (2016) are associated with significant uncertainties when applied at sites without data on ground ice content.With a sufficient spin-up, the segregated ice content is not an input parameter anymore, but is calculated by the model itself.To do so, one would need to conduct model simulations with paleoclimatic forcing back to the periods when today's segregated ice layers were actually formed.This way, it could for example be possible to simulate segregated ice stratigraphies, which have build up during colder climate periods and compare them to observations from permafrost cores today.Yet, the validation of the model remains challenging.While ground ice data from the present-day state is available for many field sites, the evolution of ground ice content in the past is poorly known.It could be possible to validate the model by deriving ice contents for specific landforms created by segregated ice, in particular peat plateaus and palsas, by retracing their formation history in the new model.However, for this, laterally coupled simulations (see section 4.4) are likely required (Martin et al., 2021).
Another possibility for validation could be laboratory freezing experiments.An example is the study by Xue et al. (2021), who conducted a one-sided freezing experiment in saturated soil to investigate the relationship of matric potential, unfrozen water content and segregated ice.The experimental conditions could be applied to a model setup and results in matric potential and in water/ice content could be compared.However, several soil parameters, which are needed as model input, are often not given in the experiment, such as initial void ratio, residual stress and compression index, so that a meaningful comparison

Influencing factors for ice segregation and thaw consolidation
The presented model scheme is capable of simulating ice segregation and thaw consolidation, which are both governed by a strong seasonality.Segregated ice forms during fall refreezing at both the downwards penetrating freeze front near the surface, and the upwards going freeze front at the bottom of the active layer.This leads to ice-enriched zones both near the surface and at the bottom of the active layer.When the ground thaws in spring and summer, the segregated ice in the upper parts of the active layer melts rapidly.In contrast, the ice at the bottom of the active layer only melts at the end of the thaw season or might even persist and eventually become part of the permafrost.Over time, this leads to an ice-enriched zone below the permafrost table (Fig. 3).
With the deepening of the active layer, the segregated ice begins to melt and excess pore water pressures initiate the consolidation process.As this process takes a certain amount of time, soil stability can potentially be reduced and excess pore water pressures might cause active layer detachment slides in sloping terrain (Lewkowicz and Harris, 2005).During thaw consolidation, the void ratio is reduced according to the compression curve of the soil and the ground surface subsides.
We highlight that the new model scheme allows us to analyze the ground ice content both in space and time and ground properties such as porosity and ice content are adjusted in each time step according to the stress conditions in the soil column.
We identify several factors, which influence the soil mechanical processes: (i) climate conditions, (ii) soil type and (iii) external ground loading, which will be discussed in the following.
(i) The climatic forcing controls ice segregation with its influence on ground temperatures and therefore the active layer thickness as well as the water supply by rain-and snowfall and evaporation rates.Moist conditions lead to more segregated ice compared to dry conditions, as more water is available for ice segregation (see supplement B) and the increased hydraulic conductivity supports the mobilization of soil water.Furthermore, high temperature gradients in the soil column (see supplement B) enhance ice segregation due to the temperature dependency of the soil matric potential in frozen state (Westermann et al., 2022).This likely increases the formation of segregated ice in continental settings with large differences in ground temperatures between summer and winter season.
(ii) The soil properties have a strong influence on ice segregation.One reason is different hydraulic properties, which govern the matric potential and its dependency on water and ice contents (Westermann et al., 2022).In addition, the soil types feature different freezing characteristics (α and n) controlling the soil unfrozen water content for negative ground temperatures.For sand, the matric potential only becomes strongly negative, when almost all water is frozen.In this state, however, the hydraulic conductivity is already very low, so that little or no water is drawn to the freezing grid cell.As a consequence, the formation of segregated ice is negligible.For the soil types silt, clay and peat, the matric potential already becomes strongly negative, when there is still a significant amount of unfrozen water available.The hydraulic conductivity is hence high enough so that water is drawn to the freeze front, resulting in ice segregation.When segregated ice is melted, the behaviour of the soil is controlled by (iii) External loads suppress the formation of segregated ice dependent on the weight of the external load.This can be explained by an increased total stress in the soil column and consequently an increased effective stress on the soil skeleton, which needs to be overcome by the matric potential differences in order for segregated ice to form.These stress conditions also speed up the consolidation process as the excess water from the melting of segregated and pore ice ice is pressed out of the soil matrix faster due to the higher total stress.This is especially true for the first year after the external load is applied: The load is acting on the top of the permafrost, where the segregated ice accumulated in the past years.Due to the soil characteristics of the clay, mobile soil water occurs next to the ground ice and is pressed out due to higher total stress.In addition, the reduced porosity due to soil consolidation cause an increase in thermal conductivity, allowing the active layer to thaw to greater depths.
Therefore, additional ground ice near the top of the permafrost may melt as well.
The model results highlight that ice segregation and thaw consolidation are influenced by climatic factors, site-specific factors (soil properties) and anthropogenic factors such as deployment of structures causing external loading.

The effects of sedimentation on ice segregation
Including sedimentation in the model runs allows us to analyze the aggradation of segregated ice in syngenetic permafrost, which often contains ice-rich layers and is susceptible to ground subsidence upon thawing.When material is deposited on the ground, either by sedimentation (e.g. in river deltas or at the bottom of slopes) or as organic matter through vegetation growth, the permafrost table moves upwards under consistent climatic conditions.The shift of the freezing front leads to a thickening of the ice-rich soil layers as demonstrated in this study (section 3.6).Therefore, the time period available for the formation governs today's segregated ice content.This is important as the distribution of ground ice preconditions thaw behaviour under a warming climate.A thin layer of segregated ice below the active layer melts rapidly and causes ground subsidence and possibly positive feedbacks to thawing, even for moderate warming.However, ground ice in deeper soil layers needs more intense warming or longer time periods until thawing and its consequences occur.Consequently, it is not only important how much excess ice can be found in the ground, but as well how it is distributed in the soil column.A setup in our new model scheme driven by paleo-climate data and sedimentation rates may be able to form segregated ice in different depths, improving estimates of ground ice content especially in field sites with poor field observations of ground ice stratigraphies.
In contrast to the previous excess ice module of CryoGrid, the new model scheme allows us to evaluate the full "ground ice mass balance" including both formation and melt: Before, ice contents were prescribed in the model setup, so that they were not necessarily in equilibrium with climatic conditions.The new model scheme allows us to simulate the accumulation of ground ice with a model spin-up.One has to pay attention to the fact that ice segregation under sedimentation begins instantaneously after the start of the simulation, when ground temperatures might not be yet in equilibrium.Therefore, a good knowledge of the ground thermal conditions is necessary, which is given in this study by initializing the model with ground temperatures given by Nitzbon et al. (2019), or a separate spin-up for ground temperatures should be performed.Despite this and other challenges such as knowledge of necessary soil parameters and climate forcing for the periods in the past when segregated ice formed, the presented model can in principle produce ice distributions in equilibrium with climatic conditions so that more realistic thaw trajectories are likely computed.

Limitations and uncertainties
In this study, we test the sensitivity of a new model scheme for ice segregation and thaw consolidation towards a range of climatic and environmental factors.However, there are limitations and uncertainties which are important to address in future studies.
In its present form, the model is one-dimensional and hence does not account for complex lateral processes, e.g.lateral water flow due to cryosuction.This approach is feasible for our proof of concept or when simulating field sites with little lateral variability in the soil structure.In order to consider complex lateral processes, laterally coupled tiles as applied in Martin et al. (2021); Nitzbon et al. (2019Nitzbon et al. ( , 2020) ) are necessary.
Our model does not take into account the different densities of water and ice.In dry conditions, this effect does not influence the void ratio of the soil significantly as the ice can expand into free pore space.However, in saturated conditions, the volume change of the phase change of water can change the void ratio of the soil.The associated error accounts for approximately 8 %, so that a result of 0.1 m segregated ice formation would correspond to about 0.11 m including the effect of phase change.Consequently, the volumetric ice content and the corresponding segregated ice and associated surface heave is slightly underestimated in our model.At least phenomenologically, this could be implemented as an independent process, which does not interact with the calculation of soil mechanics and water fluxes.However, as CryoGrid assumes the density of ice to be equal to the density of water in a range of submodules, e.g. for the snow cover (Westermann et al., 2022), this has not been implemented yet for consistency reasons.
The soil mechanical processes implemented in the model consider only primary consolidation, and do not account for long-term creep processes.This is a first order approximation, which is acceptable for simulating frozen soils, for which creep rates are typically low.For long-term simulations with thick sedimentary deposits near or above thawing temperatures, creep processes should be implemented to get a better representation of the deformation of the soil column.Furthermore, we simplify the relationship between void ratio and effective stress by having a constant compression index for each soil type.As a consequence, the deformation of the soil is completely reversible.This implies that a soil, which is compressed e.g. by external loading, expands to its initial state upon unloading.This simplifying assumption also affects the quantification of segregated ice in the model, which is defined as the ground ice volume that exceeds the initial pore volume (inital void ratio) of the soil.
An improved representation of the soil system and formation of segregation ice could therefore be obtained by introducing the concept of a preconsolidation stress, and a recompression index for unloading and reloading branches of the void ratio vs.
effective stress curve.For the model cases presented in this paper, these simplifications are expected to have a minor impact on the modelled ground surface deformation, as processes such as sedimentation and build-up of ground ice result in an increase the last decades of the model scenarios.In potential modelling scenarios including erosional processes, the proper handling of unloading and reloading would be essential.

Possible applications
The presented study is a proof of concept of the new model scheme for the CryoGrid community model.Applying it to different study setups could contribute to the demonstration of natural processes in permafrost soils.
Soils with a high organic content are widespread in the Northern Hemisphere (Hugelius et al., 2014).As these soils are sensitive to changing stress conditions due to a high compression index, it could be beneficial to apply the new model scheme for field sites affected by carbon-rich soils.The accumulation of peat could be simulated with the associated ice segregation in long-term runs.One model target could be the evolution of peat plateaus.When peat mires are exposed to cooling climate conditions e.g. during the Little Ice Age (Kjellman et al., 2018), epigenetic permafrost develops.Due to the soil properties of peat, we find in this study that considerable amounts of segregated ice can be formed (see section 3.4), resulting in strong ground subsidence under a warming climate, in agreement with observations from peat plateau sites (e.g.Martin et al., 2021).
In addition to peat plateaus, the new model scheme could be a step forward to simulate polygonal tundra.To do so, a model for wedge ice formation is required.It could be coupled laterally to the polygon center (see Martin et al., 2021;Nitzbon et al., 2019Nitzbon et al., , 2020) ) represented by the stratigraphy class in this study, which could accomplish the segregated ice formation occurring potentially in the polygon centers.
Furthermore, applying the model to inclined terrain could enable the user to simulate conditions on slopes affected by permafrost.If segregated ice is melted rapidly and excess water is mobilized faster than it can be drained, high pore pressures can develop.This can promote weak layers that can shear off in form of active layer detachment slides (Lewkowicz and Harris, 2005).Analyzing the distribution of stress conditions would be possible with the new model scheme, enabling to detect potential weak layers in the slope.Furthermore, water fluxes along the slope could be initiated by comparing the hydraulic head of the soil columns.
The new model scheme could as well be applied in geotechnical studies, especially as the effect of external loads on the soil parameters such as porosity can be analyzed.Climate change can warm permafrost temperatures and increase active layer thickness or even lead to the development of a talik.The resulting subsidence rates could be simulated with the model approach and contribute to the analysis of the stability of settlements and infrastructure.

Conclusions
In this study, we present a new model scheme, which demonstrates ice segregation as well as thaw consolidation and can CryoGrid are defined by state variables, which are set manually by the user.This includes the volumetric content of the mineral, organic, water and ice components as well as the corresponding porosity.Each stratigraphy class demands a certain set of state variables depending on the implemented model physics and functions.https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.The upper boundary of a stratigraphy class is defined as the interface between the ground surface and the atmosphere.For the base stratigraphy class GROUND_freezeC_RichardsEq_seb, the upper boundary condition is given by the surface energy balance, controlled by the exchange of short-wave radiation and long-wave radiation, latent heat flux and sensible heat flux as well as ground heat flux.To calculate the surface energy balance, the forcing data requires a time series of air temperature, solid and liquid precipitation, wind speed, short-wave and long-wave radiation, humidity and air pressure.The lower boundary is defined by a constant geothermal heat flux.Subsurface heat transfer is based on both heat conduction and heat advection.The calculation of heat conduction follows Fourier's law with the thermal conductivity of the material being the controlling factor.The heat transfer through advection is determined by vertical water fluxes.
https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.While the general CryoGrid framework and all other stratigraphy classes (such as the snow class, see above) remain unchanged, the new stratigraphy class GROUND_freezeC_RichardsEq_seb_pressure extends GROUND_freezeC_RichardsEq_seb with a representation of soil mechanical processes.The new features are presented below in sections 2.1.1 and 2.1.2.

Figure 1 .
Figure 1.Illustration of the stresses acting on the soil column with geostatic stress σgeo, external stress σext and pore water pressure u.Their relationship as well as the computation of the effective stress σ ′ z is described in Eq. (3) and Eq.(4).
(i) The total stress σ z [Pa] is the sum of both external stress σ ext [Pa] and the vertical geostatic stress σ geo [Pa]: https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.(ii) The effective stress σ ′ z [Pa] describes the stress acting on the soil matrix and is defined as the total stress σ z [Pa] minus the pore water pressure u [Pa]: σ ′ z = σ z − u.

Figure 2 .
Figure 2. Relationship between void ratio e and effective stress σ ′ z for thawing soils.If ice segregation results in void ratios larger than the initial void ratio e0, melt water produced by thawing of such grid cells drains without compression of the soil skeleton, as the effective stress stays constant at the residual stress σ ′ 0 .For increasing effective stress, the soil compresses according to the compression index Cc.The combination of void ratio and effective stress, which the soil has in a compressed state without excess pore water pressures under the applied overburden pressure, are indicated with e1 and σ ′ 1 .All ground ice exceeding the available pore space at this stress regime is defined as segregated ice.
https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License. in water content affect the grid cell size.The thickness of a grid cell d [m] can be calculated for saturated conditions from the volume Θ [m 3 ] of the grid cell (being the volume of water, ice, mineral, organic) divided by the area A [m 2 ]:

Figure 3 .
Figure 3. Illustration of the potentials resulting in water fluxes as calculated in the model (not to scale): During winter season, the entire soil column is frozen and no substantial water fluxes occur.When the upper soil layers are unfrozen during summer season, water fluxes in the unsaturated zone are controlled by the gravitational potential Pg and the matric potential ψ.Rainwater or meltwater infiltrates from the top to deeper layers.Most important for ice segregation is the matric potential during fall refreezing (marked in red).In case the downward thawing from the surface during summer reaches (partly) the segregated ice, excess pore water pressures ue develops, leading to thaw consolidation.
potential.Water fluxes into the grid cell continue, as long as (i) the difference in matric potential is large enough to compensate for the weight of the overlying soil and (ii) the hydraulic conductivity of the freezing soil allows water fluxes into the affected grid cell.With increasing ice content in the grid cell, the hydraulic conductivity is strongly reduced, which effectively stops https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.
. Due to the geographic location, Arctic continental climate leads to mean annual air temperatures of below -12 °C.Summer air tempera-https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.tures can reach up to 28 °C, while minimum air temperatures in winter season can fall below -45 °C.In summer, precipitation and evapotranspiration balance each other.Winter season is typically characterized by a thin snow cover with a mean end-ofwinter thickness of 0.3 m ( ) for the long-term thaw susceptibility run in that study.The climatic forcing between 1960 and 2012 is derived from the reanalysis product CRU-NCEP, downscaled for Samoylov with in situ data from the automatic weather station.After 2012, the forcing data is based on the Representative Concentration Pathway (RCP) 8.5 of CCSM4.To better match site measurements, statistical downscaling is performed for incoming long-wave radiation, air temperature and absolute humidity, using linear regression between measurements and CRU-NCEP and CCSM4.Comparing the reference period 1961-1990 to 2080-2100, the forcing data set shows the following trends: An increase in air temperature from -16.7 °C to -8.1 °C, an increase in longwave radiation from 214 W m -2 to 249 W m -2 , a slight decrease in shortwave radiation from 105 W m -2 to 101 W m -2 , a pronounced increase in rainfall from 628 mm to 870 mm and an increase in snowfall from 530 mm to 668 mm.A detailed description of the forcing data set can be found inNitzbon et al. (2020) andWestermann et al. (2016).In addition, we use an ancillary data set from Bayelva, Svalbard (78°55' N,11°50' E) to assess the sensitivity of the model towards climatic forcing.The field site conditions are controlled by the West Spitsbergen Current leading to a relatively mild climate for that latitude.Measured annual mean air temperature in Ny-Ålesund was -3.1 °C between 2010 and 2021, with 2010 having the lowest (-3.6 °C) and 2017 the highest (-2.8 °C) value (Norwegian Meteorological Institute, 2022a).Mean annual precipitation was recorded to be 481 mm between 2000 and 2021 (Norwegian Meteorological Institute, 2022b).It falls predominantly as snow between September and May and the maximum snow depth typically varies between 0.65 and 1.4 m(Boike et al., 2018).The field site is characterized by hilly tundra with soils featuring silty clay to sandy silt with some larger stones and organic carbon concentration being highest in about 1 m depth (> 6 % weight)(Boike et al., 2008(Boike et al., , 2018;;Westermann et al., 2009).Measurements show ground temperatures of -3.0 to -2.6 °C in a depth of 9 m between 2009 and 2016 (GTN-P, 2018) and thus relatively warm permafrost.We use the same model forcing asSchmidt et al. (2021) for the presented long-term runs in that study.Between 1980 and 2019, the model forcing is derived from ERA-Interim reanalysis, while the future forcing is based on CMIP5 projections (RCP 8.5) of CCSM4 using an anomaly approach.For a detailed description of the forcing data set seeSchmidt et al. (2021).Validation performed byWestermann et al. (2022) shows that CryoGrid can capture the ground temperatures at the Bayelva site well.https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.
initialize the ground temperatures as inNitzbon et al. (2019): 0 m depth: 0.0°C; 2 m: -2.0°C; 5 m: -7.0°C; 10 m: -9.0°C; 25 m: -9.0°C; 100 m: -8.0°C; 1100 m: +10.2°C.These values are based on borehole data from 2006 and a steady-state Table2.Stratigraphies used for the model scenarios in this study.While the multi-layer stratigraphy represents a soil column as described inNitzbon et al. (2020) andWestermann et al. (2016) for Samoylov Island, the simplified stratigraphies consist of only one soil type throughout the entire soil column below the organic-rich layer at the top.Below a depth of 9 m, we define a base layer with a volumetric mineral content of 0.7 and a volumetric ice content of 0.3, which is assumed to be static and does not evolve due to soil compression.θm: volumetric fraction of mineral content; θo: volumetric fraction of organic content; S: saturation; e0: initial void ratio before compaction; σ0: residual stress; Cc: compression index; kw: permeability; α: alpha coefficient; n: n coefficient.
as the onset of snow cover in the applied forcing data (based on reanalysis data) can differ by more than a month to the start of the snow season detected by automatic cameras https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.

Figure 4 .
Figure 4. Measured and modelled ground temperatures at a depth of 0.05 m and 0.40 m.Measured temperature data is taken in a polygon center at the Samoylov field site.Results of the model scenario S-val can capture the measured ground temperatures in the two different depths.

Figure 5 .
Figure 5. Mean annual ground temperatures M AGT in 5 m depth, depth of the top of the permafrost and the seasonal frost in the model scenario S-clay.The spin-up covers the years 1860 to 1960.M AGT warm from around -9°C until the 1980s to values close to the thaw threshold at the end of the 21st century.The depth of the permafrost table deepens from around 0.7 m to 5 m, so that the seasonal frost doesn't reach the same depth towards the end of the century and a talik develops.
saturated conditions in the upper soil layers due to temporarily increased precipitation in the future forcing (see supplement B, https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.e.g. years 2070 to 2080) influencing the soil behaviour: The buoyancy effect of the water decreases the effective stress on the soil structure, which results in an temporarily increased void ratio and thus soil swelling.

Figure 6 .
Figure 6.Sum of volumetric water and ice content in the permafrost on the reference date August 31 for the model run S-clay.Values in the thawed layer are not displayed.The spin-up covers the years 1860 to 1960.For scenario setup see table 3.

Figure 7 .
Figure 7. Column-accumulated segregated ice on the reference date August 31 of each year in the simulation, for the model scenarios S-clay, S-clay-rain50 and B-clay.Drier conditions lead to less formation of segregated ice.The moist conditions in B-clay are compensated by smaller temperature gradients in the ground so that similar segregation ice contents are formed during model spin-up.

Figure 8 .
Figure 8. Column-accumulated segregated ice on the reference date August 31 for the model scenarios S-sand, S-silt, S-clay and S-peat.

Figure 9 .
Figure 9. Volumetric water and ice content in the permafrost on the reference date August 31 in model scenario S-peat.Values in the thawed layer are not displayed.The spin-up covers the years 1860 to 1960.The sum of volumetric water and ice content is shown as soil water can still occur below freezing temperatures dependent on the soil type.For scenario setup see table3.The results show a significant formation of segregated ice below the permafrost table and associated ground heave.

Figure 11 .
Figure 11.Volumetric water and ice content on the reference date August 31 for the model run S-clay-sed1000.Values in the thawed layer are not displayed.The sum of volumetric water and ice content is shown as soil water can still occur below freezing temperatures dependent on the soil type.For scenario setup see table3.The results show a significant formation of segregated ice below the permafrost table under a https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.might be challenging.Nevertheless, model validation with laboratory experiments could be conducted in future studies and will likely contribute to model improvements.
https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.thecompression curve, which is defined by the parameters initial void ratio, residual stress and compression index.The model shows that the reaction of a soil to changing climatic conditions is strongly affected by its properties.
in effective stress and thus a movement in the normally consolidated domain.Unloading effects are minor in the presented model cases and result primarily from changes to the water balance and drainage of excess water due to ground ice melt in https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.
lead to improved simulations of ice-rich ground responding to a warming climate.It combines soil mechanical processes with soil hydrology according to Richards equation and soil freezing characteristics.We run the model with forcing data of the permafrost site Samoylov in Sibiria and analyze different influencing factors such as climatic forcing, soil type, external loads and sedimentation.Our main conclusions are: https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.https://doi.org/10.5194/egusphere-2023-41Preprint.Discussion started: 16 January 2023 c Author(s) 2023.CC BY 4.0 License.

Table 1 .
State variables for soil mechanical properties

Table 3 .
Model scenarios