The CryoGrid community model (version 1.0) – a multi-physics toolbox for climate-driven simulations in the terrestrial cryosphere

The CryoGrid community model is a flexible toolbox for simulating the ground thermal regime and the ice–water balance for permafrost and glaciers, extending a well-established suite of permafrost models (CryoGrid  1, 2, and 3). The CryoGrid community model can accommodate a wide variety of application scenarios, which is achieved by fully modular structures through object-oriented programming. Different model components, characterized by their process representations and parameterizations, are realized as classes (i.e., objects) in CryoGrid. Standardized communication protocols between these classes ensure that they can be stacked vertically. For example, the CryoGrid community model features several classes with different complexity for the seasonal snow cover, which can be flexibly combined with a range of classes representing subsurface materials, each with their own set of process representations (e.g., soil with and without water balance, glacier ice).

We present the CryoGrid architecture as well as the model physics and defining equations for the different model classes, focusing on one-dimensional model configurations which can also interact with external heat and water reservoirs. We illustrate the wide variety of simulation capabilities for a site on Svalbard, with point-scale permafrost simulations using, e.g., different soil freezing characteristics, drainage regimes, and snow representations, as well as simulations for glacier mass balance and a shallow water body. The CryoGrid community model is not intended as a static model framework but aims to provide developers with a flexible platform for efficient model development. In this study, we document both basic and advanced model functionalities to provide a baseline for the future development of novel cryosphere models.

1 Introduction

The terrestrial cryosphere is currently undergoing unprecedented changes, including thawing of permafrost, melting of glaciers and ice sheets, and changes in snow cover extent. In the last decade, permafrost temperatures have warmed almost everywhere in the circum-Arctic (Biskaborn et al., 2019), and the melting of excess ground ice can accelerate thawing in a positive feedback loop, leading to the fast transformation of permafrost landscapes through thermokarst (Farquharson et al., 2019; Nitzbon et al., 2020; Turetsky et al., 2019). Glaciers worldwide have been retreating at increasing rates (e.g., Hugonnet et al., 2021; Huss and Hock, 2018) and are important contributors to global sea-level rise. On regional scales, glacier retreat can affect, e.g., freshwater availability, infrastructure, and wildlife (e.g., Kaser et al., 2010).

Cryosphere land surface models are important tools to investigate the sensitivity of the terrestrial cryosphere under complex environmental and climatic conditions. In particular, the models allow us to project climate change impacts and thus answer urgent questions on the future of the cryosphere. As an example, glacier mass balance models are important tools for estimating the response of ice masses to a changing climate. They aid the investigation of the current state of the cryosphere in areas where in situ observations are hard to obtain and can be used to estimate the past and future evolution of glaciers (e.g., Mankoff et al., 2021; Schmidt et al., 2020; van Pelt et al., 2021). Similarly, numerical models are highly important to investigate the current state of permafrost, in particular since permafrost is usually not visible at the Earth's surface and only a limited number of measurement sites exist. Permafrost models provide insights into the evolution of Arctic landscapes and help us understand how these fragile ecosystems respond to natural and human-caused disturbances.

The purpose of land surface models is to describe nature in an adequate way, which means that the models must reproduce observations of the targeted physical parameters. This can be achieved by models of different complexity, from simple semi-empirical models trained by observations to physically based schemes which run independently of observations. Examples of semi-empirical models are degree day melt models for glacier mass balance (Gabbi et al., 2014; Reveillet et al., 2017) and the top of the permafrost table (TTOP) equilibrium model to estimate permafrost temperatures (Smith and Riseborough, 1996), which are both designed for a certain application. In contrast, physically based land surface models can simulate both glacier mass balance and the permafrost thermal regime with the same model framework, relying on universal formulations, such as the surface energy balance and Fourier's law of heat conduction. Over the past decades, land surface models have grown in complexity to incorporate a wide range of processes from various disciplines, such as biophysics, biogeochemistry, hydrology, and ecology (Fisher and Koven, 2020). In theory, continuous improvements over time could eventually lead to a unified “land surface model of everywhere, everything and all times” (Blair et al., 2019), which can reproduce and explain observations of all land surface variables, irrespective of their spatial and temporal scales. In reality, however, complex land surface models feature a large number of model parameters whose variations in space and time are poorly constrained. This severely compromises their advantage over simpler model approaches in many use scenarios, in addition to strongly increased computation demands. Therefore, simple, less process-rich models have significant advantages in many practical applications and are typically employed for high-resolution (e.g., Obu et al., 2019), long-timescale, and/or large-ensemble simulations.

The CryoGrid suite of permafrost models have provided three model categories with an increasing level of complexity to conduct a wide range of permafrost studies. CryoGrid 1 is an equilibrium model to compute mean annual ground temperatures at the top of the permafrost table (TTOP) as the only output (Gisnås et al., 2013), which in particular makes it possible to infer the presence or absence of permafrost. It relies on surface or air temperatures as input, in addition to n-factors to parameterize the effects of snow cover and active layer properties on the seasonal heat exchange. While CryoGrid 1 was used for fine-scale process studies (Gisnås et al., 2014, 2016), its main use was for large-scale mapping of permafrost extent and temperatures, e.g., the generation of a permafrost map for Scandinavia (Gisnås et al., 2017). By using globally available remote sensing and reanalysis data (MODIS land surface temperature, ERA reanalysis) to force CryoGrid 1, permafrost maps on the continental scale could be produced (Westermann et al., 2015). Later, this processing chain was extended to produce 1 km resolution permafrost maps of the Northern Hemisphere (Obu et al., 2019) and Antarctica (Obu et al., 2020). Due to a low number of parameters and an efficient and simple implementation, CryoGrid 1 allowed for large-scale ensemble simulations at 1 km grid cell size so that the effect of small-scale spatial variability of snow depths and ground properties on the thermal regime could be represented statistically. Similar modeling approaches building on analytical formulations for ground temperature and active layer thickness have also been applied at regional scales (e.g., GIPL1 in Alaska, Sazonova and Romanovsky, 2003).

Being an equilibrium model, CryoGrid 1 is generally not well suited for climate change simulations, missing transient processes. For instance, changes of the model forcing impact ground temperatures without time delay in the equilibrium model, when in reality subsurface processes like ground ice melt delay the response of ground temperature. Therefore, the transient model CryoGrid 2 is employed for mapping climate change impacts using similar spatially distributed forcing datasets as CryoGrid 1 (Czekirda et al., 2019; Westermann et al., 2013, 2017). Similar to the GIPL2 model (Jafarov et al., 2012), CryoGrid 2 computes ground temperatures from conductive heat transfer through the ground and the snowpack, validated in permafrost regions of northern Siberia and Norway (Beermann et al., 2017; Langer et al., 2013; Westermann et al., 2011). Furthermore, CryoGrid 2 can be used to compute ground temperatures in deeper layers from field measurements of near-surface ground temperatures (Peter et al., 2023). Due to its relative computational efficiency, CryoGrid 2 has also been adapted to simulate multi-millennial paleo-permafrost evolution, for example during deglaciation on Iceland (Etzelmüller et al., 2020). CryoGrid 2 was also used to model the evolution of permafrost beneath the circum-Arctic continental shelves (Overduin et al., 2019) with model forcing computed for the last 450 kyr from model reconstructions of glaciation, sea level, and air temperature. The implementation of coupled heat and salt diffusion equations in offshore sediments showed mitigation of seabed seasonal freezing and enhancement of the subsea permafrost degradation rate because of the presence of salt (Angelopoulos et al., 2019). Further applications of the model to thermokarst lagoon and coastal settings demonstrated how brine rejection lowers sediment freezing temperature and slows the refreezing of thawed sediments (Angelopoulos et al., 2020, 2021).

CryoGrid 3 is a land surface model that accounts for land–atmosphere coupling by simulating the surface energy balance, similar to the COUP (e.g., Marmy et al., 2013), GEOtop (Rigon et al., 2006; Endrizzi et al., 2014), and SURFEX (Barrere et al., 2017) models. CryoGrid 3 features a representation of excess ground ice so that ground subsidence and thermokarst pond formation upon thaw can be simulated (Westermann et al., 2016). Furthermore, it was used to simulate heat transport in water bodies as well as their impact on the thermal regime and the thaw threshold of the permafrost below (Langer et al., 2016). The implementation of the state-of-the-art snow scheme Crocus (Vionnet et al., 2012) into CryoGrid 3 allowed for transient representation of internal snow properties as well as wind redistribution of snow, which was key to realistically simulate local ground temperature dynamics in snow-rich regions (Zweigel et al., 2021). As about 55 % of permafrost area is covered by boreal forest, CryoGrid 3 was extended by a multi-layer vegetation scheme (Bonan et al., 2018) for the modeling of the thermal and hydrological permafrost conditions under boreal forest cover (Stuenzi et al., 2021a, b). CryoGrid 3 was further extended by a bucket hydrology scheme for unfrozen conditions, as well as lateral transport of water, heat, and snow; this version has been evaluated and applied for different permafrost ecosystems (Martin et al., 2019; Nitzbon et al., 2019). Nitzbon et al. (2020), Nitzbon et al. (2021), and Martin et al. (2021) further demonstrated the applicability of CryoGrid 3 to simulate complex permafrost landscape evolution over a range of spatial (plot to landscape) and temporal (years to centuries) scales. Such spatially distributed realizations of CryoGrid 3 (denoted as “laterally coupled tiles”) aim for a three-dimensional representation of permafrost hydrology, similar to cold-region hydrological models, such as WASIM (Debolskiy et al., 2021), TopoFlow (Schramm et al., 2007), SUTRA-Ice (McKenzie et al., 2007), permaFoam (Orgogozo et al., 2023), PFLOTRAN-Ice (Karra et al., 2014), and Amanzi-ATS (Painter et al., 2016; Jan et al., 2018).

While CryoGrid 1, 2, and 3 are partly based on the same model formulations and process parameterizations, they are essentially different models regarding numerics and code structure. Furthermore, they have been adapted for many different use cases, creating numerous derivatives of slightly different model versions that are not necessarily compatible.

In this study, we present the architecture of a new CryoGrid community model, which comprises most of the functionalities demonstrated in CryoGrid 1–3, while going beyond in many aspects. In particular, the CryoGrid community model is not a single model, but a modular collection of models with different functionalities which can be combined with each other to fit the requirements of a variety of applications. We describe key aspects of the model physics for one-dimensional simulations, especially when going beyond the capabilities of previously documented CryoGrid 1–3. An example is a new glacier mass balance module which extends the capabilities of the CryoGrid community model beyond permafrost. We showcase the possibilities of this new simulation tool with point simulations for Svalbard, as well as benchmark simulations against analytical solutions and reference experiments.

Figure 1 Hierarchy (left to right) of the different types of CryoGrid classes required for multi-physics simulations within the CryoGrid community model, as described in this study. Arrows represent pointers to classes employed in the hierarchical level below. Other hierarchies, potentially comprising additional and different types of classes, are fully possible within the CryoGrid community model and can be implemented in the future, e.g., a set of classes for spatially distributed applications with less customizable, but computationally more efficient, simulation tools (Sect. 4.7). See text. The CryoGrid stratigraphy is depicted in more detail in Fig. 2.

2 CryoGrid community model description

2.1 Architecture and setup

2.1.1 Model concept – modularity through object-oriented programming

The CryoGrid community model is based on an object-oriented programming paradigm implemented in the programming language MATLAB in which “objects” are referred to as “classes”. A class is a defined structure, which consists of a class-specific set of variables, as well as class-specific functions to modify these variables. A variable within a class can once again be a class (more precisely a pointer to another class), typically of a different type, which makes it possible to create a tree-like structure with different hierarchical levels (Fig. 1). Hereby, each class at a given hierarchy level contains pointers to the lower-level classes, with different levels representing different functionalities within the simulation system. In CryoGrid, each functionality level is represented by specific class types for which typically several options (i.e., different classes) are available, allowing customizing and adapting the model setup. Classes of a given type feature mandatory variables and functions with standardized inputs and outputs so that they become interchangeable building blocks of a modular simulation system with intrinsic compatibility. For each class type, the best-fitting class for the particular use case must be selected by the user. Furthermore, each class has specific parameters which must be set by the user and further control the behavior of the model system. A description of all available CryoGrid classes employed in this study is presented in Sect. S2.

Figure 1 depicts the class types and associated hierarchy that are employed to realize the multi-physics simulations described in the remainder of this study. However, it is possible to implement other configurations (e.g., with different class types and hierarchies) within the CryoGrid community model (see Sect. 4.7 for a discussion of such possibilities). The RUN_INFO class is the only mandatory class type, which represents the uppermost level of any class hierarchy. For the multi-physics simulations described in this work, the second level in the class hierarchy uses the class type TILE (Fig. 1). The purpose of a TILE class is to perform a classic one-dimensional model simulation over a pre-defined time period with a pre-defined forcing dataset (as in, e.g., CryoGrid 2 and 3). RUN_INFO classes, on the other hand, organize the model simulations. Depending on which RUN_INFO class is selected, they can, for example, launch only a single TILE class or several TILE classes (e.g., several independent simulations, representing grid cells or ensemble members) either sequentially or in parallel. Additionally, model spin-up can be implemented by sequentially simulating a number of TILE classes. For example, one TILE class can be used for the spin-up phase and another for the target period of the simulations, which is initialized by the model state of the spin-up TILE class (see Sect. 3.1.4 for an accelerated spin-up procedure using a sequence of TILE classes). A description of all RUN_INFO classes employed in this study is provided in Sect. S2.

The multi-physics simulations described in this work all employ the TILE class TILE_1D_standard (see Sect. S2 for more details) that performs a full model simulation, from model initialization to the generation of the model output. For this purpose, a range of specialized class types is employed (hierarchy level 3 in Fig. 1), which control different aspects of the simulation, such as the initialization of model state variables (STRAT_STATVAR classes), the model forcing (FORCING classes), and the model output (OUT classes). We do not describe the entire functionality of these class types here (see Sect. S2 for details), but only provide a few examples showcasing the modularity. FORCING classes, once again adhering to the strict protocol of mandatory variables and internal functions, are designed to provide the required model forcing at a specific time step. Different FORCING classes are available: for example, a class simply interpolating the raw model forcing and a class reprojecting the radiative components of the raw model forcing based on slope and aspect. The choice of the OUT class determines what kind of and how model output is stored. For development and testing, a class storing the entire variable space can be used, while users may want to design a purpose-built OUT class which only stores their model variables of interest. There is also an OUT class storing the full model state after the final simulation time step, which can be used to start a new simulation based on that state (i.e., initialize a new TILE class). STRAT_STATVAR classes (abbreviation for “stratigraphy of state variables”, not shown in Fig. 1) are employed to calculate depth profiles of model variables to define the initial state on the model grid. Depending on the class, these can be provided as layers with constant values or by interpolation between values at defined depths.

The backbone of the modularity within TILE_1D_standard is the possibility to define a vertical stack of classes, each employing different model physics and parameterizations within layers of the model domain (see Fig. 2, Sect. 2.1.2, 2.2). In the following, we refer to this vertical stack of classes as “CryoGrid stratigraphy” and the classes within this stack, which encode the model physics, as “stratigraphy classes”. The vertical domain covered by each of the stratigraphy classes (i.e., the stratigraphy of these classes) is again assigned via a purpose-built class (STRAT_CLASSES, denoted as STRATIGRAPHY in Fig. 1). Finally, TILE_1D_standard features a class type controlling lateral interactions with an external environment (denoted as LATERAL in Fig. 1). These LATERAL classes are described in Sect. 2.3.

2.1.2 Multi-physics representation with stratigraphy classes

A one-dimensional model simulation (in TILE_1D_standard, see above) is realized by vertically stacking different stratigraphy classes (Fig. 2) which are each defined by their specific model physics and state variables. Examples of stratigraphy classes are ground columns with and without water balance, water bodies, glaciers, and snow with different levels of process representation (see Sect. 2.2 for details). Each class occupies a certain vertical domain, with the boundary conditions applied to the uppermost and lowermost class (Fig. 3).

Figure 2 Example of a CryoGrid stratigraphy in the tile class TILE_1D_standard, showing the different classes connected by pointers (symbolized by arrows). The CryoGrid stratigraphy consists of stratigraphy classes 1 to 3 coupled by interaction (IA) classes that specify the exchange of heat and mass between pairs of stratigraphy classes. The stratigraphy is realized as a linked list with pointers between classes and interaction classes, and the top and bottom of the list are represented by a dedicated TOP and BOTTOM class (which have no other functionality).

Figure 3 Example of a CryoGrid stratigraphy in the tile class TILE_1D_standard with two stratigraphy classes. Each stratigraphy class has its own state variables and model grid, and energy and water are exchanged between stratigraphy classes (as coded in the interaction class, Fig. 2). At the upper boundary, energy and water are exchanged according to the model forcing in a class-specific way (see Sect. 2.2.2), while a heat flux is typically applied at the lower boundary. Note that water fluxes only occur between stratigraphy classes which account for soil hydrology, i.e., for the transient evolution of soil water contents (see Sect. 2.2.4).

Within each stratigraphy class, CryoGrid computes the time evolution for state variables, such as ground temperature and water–ice contents. We distinguish between prognostic state variables, for which a time derivative is calculated and which are then integrated in time to advance to the next time step, and diagnostic state variables which are not time-integrated. For the prognostic state variables, CryoGrid uses the simple time integration scheme “first-order forward Euler”; i.e., the new model state is computed as the old model state plus time derivatives times the model time step (see Sect. 2.2.9 for details). When the new model state of the prognostic state variables is obtained, diagnostic state variables are calculated from the prognostic variable by constitutional relationships. This calculation is instantaneous, i.e., does not depend on the employed time step. It is possible that a physical property (e.g., temperature) is a prognostic variable in some stratigraphy classes and diagnostic in others. Stability and accuracy of the time integration are ensured by automatically selecting an appropriate time step. The calculation of a suitable time step is not accomplished by well-known stability criteria for the first-order forward Euler scheme, but explicitly takes the physics represented by the individual stratigraphy classes into account to ensure both a stable and accurate simulation (see Sect. 2.2.9).

Interactions between stratigraphy classes are realized by interaction classes (Fig. 2) that compute fluxes across the boundaries between pairs of stratigraphy classes (Fig. 3). Thus, compatibility of two stratigraphy classes (i.e., if they can border each other in the CryoGrid stratigraphy) can be ensured by providing a dedicated interaction class. Note that the compatibility depends on the order of the two classes in the stratigraphy. In particular, interaction classes compute fluxes across boundaries between stratigraphy classes which are required for computing time derivatives of state variables in the prognostic step of the time integration. For example, if two stratigraphy classes have temperature as state variables and the two modules are connected through heat conduction, the interaction class computes the conductive heat flux between the adjacent grid cells of the classes. If two stratigraphy classes feature different state variables, the interaction class must contain the necessary code to compute the correct fluxes for both involved classes. For example, if only one of the classes is hydrologically active, while the sum of water and ice contents is static for the other class (see Sect. 2.2), the interaction class must provide a zero water flux boundary condition for the hydrologically active class to reflect the fact that water flow through the boundary is not possible. During initialization of the CryoGrid stratigraphy, the correct interaction class is automatically selected for each pair of stratigraphy classes.

2.1.3 Dynamic behavior with stratigraphy class triggers

An important feature of stratigraphy classes is their ability to modify or rearrange the CryoGrid stratigraphy itself if a certain condition (referred to as a “trigger”) is met. In particular, a class can remove itself from the CryoGrid stratigraphy or insert a new stratigraphy class above its own position. As an example, a dynamic representation of ponds (using water body and excess ice classes, see Sect. 2.2) can be achieved by such triggers which modify the CryoGrid stratigraphy. When surface water pooling up over initially dry ground reaches a user-defined threshold, a water body class representing the physics of energy transfer within a water body is created and inserted in the CryoGrid stratigraphy. Likewise, if the water depth of a water body drops below that threshold, the water body class is automatically removed. In this process, all state variables are automatically adjusted to ensure mass and energy conservation.

A special situation is the representation of the seasonal snow cover, which again is handled by stratigraphy class triggers creating a snow class (Sect. 2.2.6) upon initial snowfall and removing it when all snow has melted. For the numerical scheme, handling a very shallow initial snow cover poses significant problems, as this results in a small grid cell size and thus very small time steps. Therefore, snow classes are attached and detached in two stages in the CryoGrid community model. After the first snowfall, the snow class is added as a so-called CHILD to the uppermost stratigraphy class (Fig. 4 left); i.e., it is not part of the CryoGrid stratigraphy, but evolves as part of the uppermost stratigraphy class. In this CHILD state, the snow does not grow vertically, but is assigned a fixed water equivalent in the vertical direction, while at the same time covering only a fraction of the uppermost stratigraphy class. This way, the snow volume is correct, but it is possible to assign the snow cover a sufficient thickness to prevent numerical problems. Algorithm 1 shows how the upper boundary condition is applied to both the uppermost stratigraphy class (“class1”) and the class representing the snow cover (“snowclass”). Most importantly, this procedure allows for a complete separation between these two classes: class1 does not contain any information about snowclass, which applies its class-specific upper boundary condition, while snowclass does not contain any class-specific knowledge on class1 (which also applies its class-specific upper boundary condition). This scheme works for any compatible pair of snowclass–class1 so that, for example, the whole range of stratigraphy classes for ground material (Sect. 2.2.2–2.2.5) can be combined with the different stratigraphy classes representing snow (see Sect. 2.2.6).

Algorithm 1 Pseudocode for applying the upper boundary condition to the stratigraphy class “class1” in Fig. 4, for the case when then ground is entirely snow-free and for the case when the snow cover is in the CHILD phase (Fig. 4, left side). Stratigraphy class “class1” represents a snow-free ground column, while stratigraphy class “snowclass” represents a snow column. “class1” has a variable CHILD, with class1.CHILD = snowclass (pointer to class “snowclass”) if snow is in the CHILD phase, and class1.CHILD = 0 if no snow is present. Both “class1” and “snowclass” have a function get_upper_boundary_condition() which applies the class-specific upper boundary condition (e.g., using a different albedo for “snowclass” than for “class1”). The total area of the model domain is A [m²].

if class1.CHILD = 0 then

– apply get_upper_boundary_condition() to “class1”, using the full area A

else

– assign fractional area f_snow to “snowclass” stored in “ground.CHILD”

and fractional area f_ground to “class1”, so that $f_{\text{snow}}+f_{\text{ground}}=\mathrm{1}$

– apply get_upper_boundary_condition() to “class1”, using area f_groundA

– apply get_upper_boundary_condition() to “snowclass” stored in “class1.CHILD”, using area f_snowA

end if

Figure 4 Schematic representation of the CryoGrid stratigraphy when a snow cover is present (represented by dedicated snow classes, Sect. 2.2.6). (a) Snow is initially added as a so-called CHILD to the uppermost subsurface class (class 1); i.e., it is not part of the regular stratigraphy, but is addressed by specific pointers. In the CHILD phase, the snow is assumed to only cover part of the surface area; e.g., the surface energy balance is calculated as a mix of snow-covered and snow-free ground (see Algorithm 1 for details). (b) When the snow water equivalent exceeds a user-defined threshold, the stratigraphy is rearranged and the snow class becomes the uppermost class of the regular stratigraphy. The process is reversed when the snow melts.

As more snow accumulates, the snow class simply expands its aerial coverage until the amount of snow is sufficient to be handled without numerical problems. Then, the snow class and the associated interaction class are rearranged so that the snow class becomes part of the normal CryoGrid stratigraphy (Fig. 4 right) and thus covers its full area. The procedure is mirrored upon snowmelt, with the snow class first becoming a CHILD and finally being removed completely upon completion of melt. The threshold in snow water equivalent at which the snow class becomes a normal part of the CryoGrid stratigraphy should generally be chosen to be small enough that the CHILD phase does not last longer than a few days, thus only having a negligible effect on model results on timescales longer than a few weeks.

2.1.4 Model operation

All CryoGrid simulations are controlled by a parameter file which defines all aspects of the run, such as the definition of the CryoGrid stratigraphy and lateral interactions, the forcing data, the model grid, and the model output format. At this point, the parameter file can either be set up as a spreadsheet (MS Excel or compatible programs) or as a text file in YAML format.

In the parameter files, all classes required for the simulation are defined, in no particular order. Each class is identified by its name and a unique index, which makes it possible to define the same class several times with different parameters. Furthermore, all mandatory parameters specific to each of the classes must be specified. In the hierarchy of the CryoGrid classes (Fig. 1), the classes on the following level are defined as parameters in the classes of the previous level. Upon initialization, the uppermost hierarchy level (defined as the RUN_INFO class with index 1) is read first, which provides the information for reading the classes in the second level, and so on. In this process, the class connections by pointers are established. Note that the standardized class structure with mandatory variables and functions facilitates a generalized initialization routine which does not make use of specific knowledge of the involved CryoGrid classes and their hierarchy.

2.2 Physics and defining equations of stratigraphy classes

At present, there are around 10 stratigraphy classes, each with different defining equations and model physics which generally contain additional parameters and options to customize its behavior. However, the classes share many common parts and features. In the following, we describe the defining equations and parameterizations of the different model components and categories. A description of each stratigraphy class is provided in Sect. S2.

2.2.1 State variables and model grid

All stratigraphy classes (except the equilibrium TTOP model class, Sect. 2.2.2) feature variables for the subsurface properties volumetric mineral content θ_m, volumetric organic content θ_o, volumetric water content θ_w, and volumetric ice content θ_i. A depth stratigraphy of mineral and organic contents is generally provided by the user, which also defines the porosity as $\mathit{\varphi }=\mathrm{1}-{\mathit{\theta }}_{\mathrm{m}}-{\mathit{\theta }}_{\mathrm{o}}$. In some stratigraphy classes, the sum of water and ice contents ${\mathit{\theta }}_{\mathrm{wi}}={\mathit{\theta }}_{\mathrm{w}}+{\mathit{\theta }}_{\mathrm{i}}$ is constant and provided by the user, while it evolves dynamically in others, driven by precipitation, evapotranspiration (Sect. 2.2.4), and potentially lateral runoff (Sect. 2.3.2). Finally, there is an air phase θ_a defined by ${\mathit{\theta }}_{\mathrm{a}}+{\mathit{\theta }}_{\mathrm{wi}}=\mathit{\varphi }$. Each model grid cell has an enthalpy state, e [J m⁻³], which is composed of a sensible and a latent part:

$$\begin{array}{}\text{(1)}& e(T,{\mathit{\theta }}_{\mathrm{w}})=cT-L_{\text{sl}}^{\text{vol}}({\mathit{\theta }}_{\mathrm{wi}}-{\mathit{\theta }}_{\mathrm{w}})\text{,}\end{array}$$

where c [$\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{3}}$] is the volumetric heat capacity, T [^∘C] the temperature, and $L_{\mathrm{sl}}^{\mathrm{vol}}$ the volumetric latent heat of water freezing [J m⁻³]. The zero point of the enthalpy is thus defined as T=0 ^∘C and θ_i=0; i.e., the grid cell is at 0 ^∘C, but all water is unfrozen. While $L_{\mathrm{sl}}^{\mathrm{vol}}$ is a constant, c is computed from θ_m, θ_o, and θ_wi and the specific volumetric heat capacities of the mineral, organic, water and ice phases (c_m, c_o, c_w, c_i) as follows.

$$\begin{array}{}\text{(2)}& c=\left\{\begin{array}{ll}{\mathit{\theta }}_{\mathrm{m}}{c}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{o}}{c}_{\mathrm{o}}+{\mathit{\theta }}_{\mathrm{wi}}{c}_{\mathrm{i}}& \text{for}T<\mathrm{0}{}^{\circ }\mathrm{C}\\ {\mathit{\theta }}_{\mathrm{m}}{c}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{o}}{c}_{\mathrm{o}}+{\mathit{\theta }}_{\mathrm{wi}}{c}_{\mathrm{w}}& \text{for}T\ge \mathrm{0}{}^{\circ }\mathrm{C}\end{array}\right.\end{array}$$

Note that the heat capacity of a potential air phase is neglected. While these state variables are employed in the stratigraphy classes described in this study, their use is not mandatory within the CryoGrid community model, and fully valid stratigraphy classes with different sets of state variables can be created. The compatibility with existing stratigraphy classes is ensured by appropriate interaction classes (Sect. 2.1.2), which compute fluxes between classes and, if necessary, convert between different sets of state variables.

The model grid is defined by the user, again using a dedicated GRID class. At present, only one grid class is implemented, in which constant grid cell sizes are specified within a sequence of layers. Typically, the smallest grid cell size is defined for the top layer and the largest for the bottom layer. Other grid classes, e.g., with grid cell sizes increasing logarithmically with depth, could be implemented in a straightforward way.

2.2.2 Coupling to model forcing and boundary conditions of model domain

At the uppermost stratigraphy class of the CryoGrid stratigraphy (Fig. 2), the upper boundary condition is applied, which simulates the coupling to the model forcing. Three different schemes are implemented at this point, broadly providing the functionality of the CryoGrid 1, 2, and 3 models.

Equilibrium TTOP approach. Used within CryoGrid 1, the TTOP approach offers an efficient way to estimate mean annual ground temperature (MAGT) directly from the model forcing. For this purpose, freezing and thawing degree days at the surface (FDD_s and TDD_s) are calculated from the temperature forcing (often using air temperature to approximate surface temperature) and semi-empirical n-factors, which phenomenologically simulate the asymmetry of heat transfer in the ground between freezing and thawing periods:

$$\begin{array}{}\text{(3)}& \text{MAGT}=\left\{\begin{array}{l}\frac{\mathrm{1}}{\mathit{\tau }}(n_{\mathrm{f}}{\text{FDD}}_{\mathrm{s}}+r_k{n}_{\mathrm{t}}{\text{TDD}}_{\mathrm{s}})\\ \text{for}{n}_{\mathrm{f}}{\text{FDD}}_{\mathrm{s}}+r_k{n}_{\mathrm{t}}{\text{TDD}}_{\mathrm{s}}\le \mathrm{0}\\ \frac{\mathrm{1}}{\mathit{\tau }}(\frac{\mathrm{1}}{r_k}{n}_{\mathrm{f}}{\text{FDD}}_{\mathrm{s}}+n_{\mathrm{t}}{\text{TDD}}_{\mathrm{s}})\\ \text{for}{n}_{\mathrm{f}}{\text{FDD}}_{\mathrm{s}}+r_k{n}_{\mathrm{t}}{\text{TDD}}_{\mathrm{s}}>\mathrm{0}\text{,}\end{array}\right.\end{array}$$

with τ the number of days in the period for which the TTOP model is applied. Setting n_f and n_t unlike unity causes a temperature offset between the model forcing and the ground surface related to processes during the frozen season (in particular caused by the insulating snow cover) and the thawed season (e.g., caused by incoming radiation modified by slope and aspect). In the same fashion, r_k causes a temperature offset between the ground surface and the top of the permafrost due to differences in active layer thermal conductivities and thus heat transfer between summer and winter. Detailed derivations of the TTOP equation (Eq. 3) are presented in Romanovsky and Osterkamp (1995) and Westermann et al. (2015). The TTOP approach calculates the ground temperature in equilibrium with the applied model forcing, i.e., the ground temperature that would eventually be reached if the model forcing was repeatedly applied for an infinite time period. The approach is not suited to capture ground temperatures during periods of rapid change, especially when an insulating top layer and ground ice delay the penetration of the surface temperature signal into the ground. For this reason, the TTOP approach should preferably be applied to longer time periods, e.g., for one or several decades of model forcing. Moreover, only a single temperature value is delivered without specifying a depth, or even a depth profile of ground temperatures. However, the TTOP model class can be vertically coupled to stratigraphy classes simulating heat conduction into the ground (Sect. 2.2.3), effectively providing a modified temperature boundary condition (see below) which accounts for temperature offsets caused by the snow cover, active layer dynamics, or exposition. Such model setups are particularly useful for simulating temperature dynamics in deeper layers for long timescales, as they do not need to resolve the seasonal freeze–thaw cycle, making them very efficient computationally.

Temperature boundary condition. Transient simulations with heat-conduction-based models (Sect. 2.2.3) require specification of boundary conditions, which can be either time (t) series of temperature, T_ub(t), or time series of the energy flux into the first model grid cell, F_ub(t) [W m⁻²]. In the CryoGrid community model, it is possible to specify a temperature boundary condition (as in the CryoGrid 2), which for each time t is translated to a heat flux into the first model grid cell as

$$\begin{array}{}\text{(4)}& F_{\text{ub}}\left(t\right)=K_{h,\mathrm{1}}\left(t\right){\displaystyle \frac{T_{\text{ub}}\left(t\right)-T_{\mathrm{1}}\left(t\right)}{\mathrm{\Delta }{z}_{\mathrm{1}}}},\end{array}$$

where K_h,1 [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] denotes the thermal conductivity (Sect. 2.2.3), T₁ the temperature, and Δz₁ the thickness of the uppermost model grid cell. Note that Eq. (4) assumes that T_ub(t) is assigned to a virtual grid cell of thickness Δz₁ above the uppermost grid cell. Air temperatures are most commonly used for T_ub(t), but this can be adapted by selecting (and if necessary modifying) an appropriate FORCING class (Sect. 2.1.1) for the simulations.

Surface energy balance. The energy flux into the first model grid cell, F_ub(t), can also be calculated from the surface energy balance (SEB), largely similar to the implementation in CryoGrid 3. The required forcing data include incoming shortwave and longwave radiation (S_in and L_in) [W m⁻²], precipitation (P_s, solid and P_l, liquid) [mm d⁻¹], and air pressure p [Pa], as well as air temperature T_air, wind speed U [m s⁻¹], and specific humidity q_air [kg water vapor kg⁻¹ air] at height above ground h [m] (e.g., Westermann et al., 2016). F_ub is then calculated from the surface energy balance equation:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{F}_{\text{ub}}\left(t\right)& =S_{\text{in}}\left(t\right)-S_{\text{out}}\left(t\right)+L_{\text{in}}\left(t\right)-L_{\text{out}}\left(t\right)\\ & -Q_h\left(t\right)-Q_e\left(t\right)\text{,}\end{array}\end{array}$$

where S_out and L_out denote outgoing shortwave and longwave radiation (both defined positive), and Q_h and Qe denote the sensible and latent heat flux (both defined positive when cooling the surface). Note that the ground heat flux Q_g is not explicitly represented, but becomes manifest in both an enthalpy change of the uppermost grid cell and a conductive heat flux between the two uppermost grid cells. It can be calculated by the user as Q_g=F_ub (assuming the same sign convention as for sensible and latent heat flux). The outgoing shortwave radiation is computed with the surface albedo α_s.

$$\begin{array}{}\text{(6)}& S_{\text{out}}={\mathit{\alpha }}_{\mathrm{s}}{S}_{\text{in}}\end{array}$$

In most classes, a single broadband (i.e., spectrally averaged) albedo is employed, but some classes resolve several spectral bands for which the albedo can vary. Furthermore, a constant albedo is provided by the user in some classes, while the albedo is parameterized as a function of model state variables in others. The outgoing longwave radiation is computed from Stefan–Boltzmann and Kirchhoff's law as

$$\begin{array}{}\text{(7)}& L_{\text{out}}=\mathit{ϵ}{\mathit{\sigma }}_{\text{sb}}(T+T_{\text{mfw}}{)}^{\mathrm{4}}+(\mathrm{1}-\mathit{ϵ})L_{\text{in}}\text{,}\end{array}$$

with σ_sb the Stefan–Boltzmann constant [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$], T_mfw=273.15 K the freezing temperature of free water, and ϵ [–] the surface emissivity.

The sensible heat flux is computed from air temperature at a defined height above ground h and the temperature of the first grid cell as

$$\begin{array}{}\text{(8)}& Q_h=-{\displaystyle \frac{{\mathit{\rho }}_{\text{air}}{c}_p}{r_{\mathrm{a}}}}(T_{\text{air}}-T_{\mathrm{1}})\text{,}\end{array}$$

with ρ_air air density [kg m⁻²] and c_p the air heat capacity [$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] at constant pressure. The aerodynamic resistance r_a is calculated from Monin–Obukhov similarity theory (Monin and Obukhov, 1954), as in CryoGrid 3, with

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{r}_{\mathrm{a}}& ={\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}U}}\left(\text{ln}{\displaystyle \frac{h}{z_{\mathrm{0}}}}-{\mathit{\psi }}_{\text{M}}\left({\displaystyle \frac{h}{L_{*}}},{\displaystyle \frac{z_{\mathrm{0}}}{L_{*}}}\right)\right)\\ & \cdot \left(\text{ln}{\displaystyle \frac{h}{z_{\mathrm{0}}}}-{\mathit{\psi }}_{\text{H,W}}\left({\displaystyle \frac{h}{L_{*}}},{\displaystyle \frac{z_{\mathrm{0}}}{L_{*}}}\right)\right)\text{.}\end{array}\end{array}$$

Here, U is the wind speed at height h above ground, κ=0.4 the von Kármán constant, z₀ [m] the roughness length (assumed to be equal for heat, water, and momentum), and ψ_M and ψ_H,W integrated atmospheric stability functions equal to the ones used in CryoGrid 3 (Westermann et al., 2016). The Obukhov length L_* [m] is calculated from the sensible and latent heat flux values in the same way as CryoGrid 3 using the flux values computed for the previous time step (see Westermann et al., 2016).

The latent heat flux is calculated as

$$\begin{array}{}\text{(10)}& Q_e=-{\mathit{\rho }}_{\text{air}}{L}_{\text{lg,sg}}{\displaystyle \frac{f}{r_e}}(q_{\text{air}}-q_{\mathrm{1}})\text{,}\end{array}$$

where L_lg,sg [J kg⁻¹] represents the latent heat of evaporation (lg) and sublimation (sg), which are employed for T₁≥0 ^∘C and T₁<0 ^∘C, respectively. Depending on the subsurface class used, different formulations for the reduction factor f [–], the resistance r_e [s m⁻²], and the specific humidity above the surface q₁ are employed. Four schemes can be broadly distinguished.

1. For subsurface classes with unlimited surface water or ice supply (e.g., classes representing snow cover or a water body), r_e=r_a, f=1, and q₁=q_sat(T₁), i.e., the specific humidity at saturation for the temperature of the uppermost grid cell, calculated with the Magnus equation. The evaporation thus corresponds to the potential evaporation.

2. For subsurface classes without water balance (i.e., water plus ice content constant in time), f=1, q₁=q_sat(T₁), and $r_e=r_{\mathrm{a}}+r_{\mathrm{s}}$, with r_s [s m⁻¹] a user-defined surface resistance to evaporation (see Westermann et al., 2016).

3. For subsurface classes with a bucket water scheme (Sect. 2.2.4), the potential evapotranspiration (i.e., r_e=r_a and q₁=q_sat(T₁)) is multiplied by the reduction factor f taking soil water availability into account. For unfrozen ground, f is calculated from the water availability coefficients of soil grid cells i, ${\mathit{\alpha }}_i^{\mathit{\theta }}$, as

   $$\begin{array}{}\text{(11)}& \begin{array}{rl}f& =f_{\text{tr}}\sum _i{\mathit{\alpha }}_i^{\mathit{\theta }}\mathrm{\Delta }{z}_i{e}^{-d_i/d_{\text{tr}}}/\sum _i\mathrm{\Delta }{z}_i{e}^{-d_i/d_{\text{tr}}}\\ & +f_{\text{ev}}\sum _i{\mathit{\alpha }}_i^{\mathit{\theta }}\mathrm{\Delta }{z}_i{e}^{-d_i/d_{\text{ev}}}/\sum _i\mathrm{\Delta }{z}_i{e}^{-d_i/d_{\text{ev}}}\text{,}\end{array}\end{array}$$

   which allows the user to specify the partitioning in transpiration (fraction f_tr) and evaporation (fraction $f_{\text{ev}}=\mathrm{1}-f_{\text{tr}}$) with different characteristic depth d_tr and d_ev. Both evaporation and transpiration are assumed to decay exponentially with depth d_i below the surface (taken positive). Furthermore, the weighting of each grid cell depends on the thickness of grid cell i, Δz_i, and the water availability coefficient calculated with the user-defined field capacity θ_fc,i [–] as ${\mathit{\alpha }}_i^{\mathit{\theta }}=\mathrm{0.25}(\mathrm{1}-\cos (\mathit{\pi }{\mathit{\theta }}_{\mathrm{w},i}/{\mathit{\theta }}_{\mathrm{f}c,i}){)}^{\mathrm{2}}$ for ${\mathit{\theta }}_{\mathrm{f}c,i}\ge {\mathit{\theta }}_{\mathrm{w},i}$ and ${\mathit{\alpha }}_i^{\mathit{\theta }}=\mathrm{1}$ for ${\mathit{\theta }}_{\mathrm{f}c,i}<{\mathit{\theta }}_{\mathrm{w},i}$. When the ground is frozen, sublimation is set to zero (i.e., f=0), which in most real-world cases is not a limitation, as snow cover builds up for which sublimation according to scheme 1 (see above) can occur.

4. In subsurface classes in which soil moisture is governed by the Richards equation, water can flow upwards to compensate for evaporative losses and all evaporated water is hence drawn from the uppermost grid cell. Similar to scheme 3, f is set to $f=\mathrm{0.25}(\mathrm{1}-\cos (\mathit{\pi }{\mathit{\theta }}_{\mathrm{w},\mathrm{1}}/{\mathit{\theta }}_{\mathrm{f}c,\mathrm{1}}){)}^{\mathrm{2}}$, while r_e=r_a is assumed and the specific humidity is set to

   $$\begin{array}{}\text{(12)}& q_{\mathrm{1}}=\text{exp}\left({\displaystyle \frac{{\mathit{\psi }}_{\mathrm{1}}g}{R_{\text{wv}}(T_{\mathrm{1}}+T_{\text{mfw}})}}\right)q_{\text{sat}}\left(T_{\mathrm{1}}\right)\end{array}$$

   (Philip, 1957), where R_wv is the gas constant for water vapor [$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$], g is the gravitational acceleration [m s⁻²], and ψ₁ [m] the matric potential of the uppermost grid cell (Sect. 2.2.3). Note that this scheme can only represent evaporation and should be combined with a dedicated vegetation module, such as the one demonstrated in Stuenzi et al. (2021a), to also represent transpiration.

Lower boundary. The lower boundary condition is applied to the lowermost class in the stratigraphy. In all classes described in this study, a user-defined constant heat flux F_lb [W m⁻²] is added to the lowermost grid cell which corresponds to the geothermal heat flux, Q_geo [W m⁻²], for sufficiently deep model domains. Although not yet implemented, it is possible to create classes with temperature boundary conditions similar to the upper boundary (see above).

2.2.3 Subsurface heat transfer and temperature calculation

Heat conduction. Depending on the selected stratigraphy class, CryoGrid considers heat conduction as well as heat advection as the dominant modes of heat transport in the subsurface. Thus, the change of enthalpy e (see Sect. 2.2.1) is given by the continuity equation

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial e}{\partial t}}=-{\displaystyle \frac{\partial j_{\mathrm{hc}}}{\partial z}}-{\displaystyle \frac{\partial j_{\mathrm{hw}}}{\partial z}}\text{,}\end{array}$$

with z [m] the vertical coordinate, j_hc the flux due to heat conduction, and j_hw the flux due to heat advected by water. Heat conduction is calculated from Fourier’s law as

$$\begin{array}{}\text{(14)}& j_{\mathrm{hc}}=-K_h{\displaystyle \frac{\partial T\left(e\right)}{\partial z}}\text{,}\end{array}$$

with K_h [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] the thermal conductivity and temperature T as a function of e (see “Soil freezing characteristics” below). The flux from heat advection with water flow is calculated as

$$\begin{array}{}\text{(15)}& j_{\mathrm{hw}}=c_{\mathrm{w}}T\left(e\right)j_{\mathrm{w}}\text{,}\end{array}$$

with c_w [$\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{3}}$] the volumetric heat capacity of liquid water and j_w [m s⁻¹] the water flux, which consists of a term for vertical advection and a term for evapotranspiration (Sect. 2.2.4). For the thermal conductivity K_h, different parameterizations in terms of the volumetric contents of water, ice, minerals, organics, and air can be selected by the user. For soil material, the parameterization by Cosenza et al. (2003) implemented in CryoGrid 2 and 3 is available, as is the parameterization used in the Community Land Model (CLM) 4.5 (Oleson et al., 2013). While the former treats all soil constituents as equal and thus functions for single-phase materials, e.g., pure ice or rock, the latter is strictly focused on soils with reasonable porosity values. It first computes thermal conductivities for dry and saturated soil, which are then weighted with the Kersten number to yield the final conductivity value (Johansen, 1973). For snow, the thermal conductivity is computed as a function of snow density, with two parameterizations available, namely the exponential relationship described in Yen (1981) and the quadratic relationship from Sturm et al. (1997).

Soil freezing characteristics. The soil freezing characteristics is a constitutive relationship between soil temperature and the unfrozen water content. In the CryoGrid community model, we generalize this concept to derive soil temperature T and water content θ_w from enthalpy e and water plus ice content θ_wi in the diagnostic step (Sect. 2.1.2). Depending on the subsurface class, either the “free water” freezing characteristic or the soil freezing characteristic described in Painter and Karra (2014) is implemented. In the free water case, the phase change of water occurs at 0 ^∘C, i.e.,

$$\begin{array}{}\text{(16)}& T=\left\{\begin{array}{ll}e/c& \text{for}e\ge \mathrm{0}\\ \mathrm{0}& \text{for}-L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}}\le e<\mathrm{0}\\ (e+L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}})/c& \text{for}e<-L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}}\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(17)}& {\mathit{\theta }}_{\mathrm{w}}=\left\{\begin{array}{ll}{\mathit{\theta }}_{\mathrm{wi}}& \text{for}e\ge \mathrm{0}\\ {\mathit{\theta }}_{\mathrm{wi}}(\mathrm{1}+e/(L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}}\left)\right)& \text{for}-L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}}\le e<\mathrm{0}\\ \mathrm{0}& \text{for}e<-L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}}.\end{array}\right.\end{array}$$

Here, $L_{\mathrm{s}l}^{\text{vol}}$ denotes the volumetric latent heat of freezing [J m⁻³] and θ_wi the sum of the volumetric water and ice contents.

For the soil freezing characteristic by Painter and Karra (2014), the free water case functions are employed for e≥0. For e<0, unique functions θ_w(T,θ_wi) relating liquid water content to T and θ_wi exist, which we use to calculate e(T,θ_wi) according to Eq. (1). For e(T,θ_wi), lookup tables are compiled, which allow efficiently evaluating the inverse function T(e,θ_wi) (and θ_w(e,θ_wi), combining T(e,θ_wi) and θ_w(T,θ_wi)) in each time step, thus computing the diagnostic variables temperature and volumetric water content from the prognostic variable enthalpy. θ_w(T,θ_wi) is calculated using the matric potential ψ [m], which also governs water flow in subsurface classes based on the Richards equation (see Sect. 2.2.4). First, the matric potential in unfrozen state, ψ₀, is evaluated with the van Genuchten–Mualem model (van Genuchten, 1980) as

$$\begin{array}{}\text{(18)}& {\mathit{\psi }}_{\mathrm{0}}=-\mathrm{1}/\mathit{\alpha }{\left(({\mathit{\theta }}_{\mathrm{wi}}/\mathit{\varphi }{)}^{-\mathrm{1}/m}-\mathrm{1}\right)}^{\mathrm{1}/n}\text{,}\end{array}$$

with α [m⁻¹] and n [–] soil-type-specific parameters, $m=\mathrm{1}-\mathrm{1}/n$, and assuming no residual water. While ψ=ψ₀ for unfrozen soil (T≥0 ^∘C), the matric potential for freezing soil (T<0 ^∘C) is calculated as

$$\begin{array}{}\text{(19)}& \mathit{\psi }={\mathit{\psi }}_{\mathrm{0}}+\mathit{\beta }{\displaystyle \frac{L_{\text{sl}}^{\text{vol}}}{g{\mathit{\rho }}_{\mathrm{w}}}}{\displaystyle \frac{T-T_{\text{mfw}}}{T_{\text{mfw}}}}\text{.}\end{array}$$

Here, T is in Kelvin, g [m s⁻²] is the gravitational acceleration, ρ_w is water density [kg m⁻³], and β is the ratio of ice–liquid to liquid–air surface tensions for non-colloidal soil, set to 2.2 as suggested in Painter and Karra (2014). The water content is finally calculated as

$$\begin{array}{}\text{(20)}& {\mathit{\theta }}_{\mathrm{w}}=\mathit{\varphi }{\left(\mathrm{1}+(-\mathit{\alpha }\mathit{\psi }\right)}^n{)}^{-m}\text{.}\end{array}$$

The values of α and n are determined by the soil type, and users can define an unlimited number of layers with different soil types (limited by the vertical resolution of the model grid). However, only a limited number of different soil types is possible within a stratigraphy class due to the need for lookup tables, which are specific for combinations of α and n. Currently, four soil types (sand, silt, clay, peat) are implemented to provide users with a convenient interface, but it is possible to change the α and n values associated with each of them so that other soil types can also be realized.

2.2.4 Water balance

In the CryoGrid community model, three schemes to compute the time dynamics of soil water contents are available, namely (1) no flow (i.e., constant water plus ice contents), (2) a “bucket” scheme with only downward vertical water flow driven by gravity, and (3) vertical water flow governed by the Richards equation. For schemes 2 and 3, the hydrological boundary conditions at the top of the soil column, such as rainfall input P_l, snowmelt, and evapotranspiration (related to the latent heat fluxes, Sect. 2.2.2), drive the time dynamics of the soil water content. Therefore, these only work in conjunction with the surface energy balance as the upper boundary condition (Sect. 2.2), while scheme 1 can be applied for both temperature and surface energy balance boundary conditions, as in CryoGrid 2 (Westermann et al., 2013) and the initial version of CryoGrid 3 (as in Westermann et al., 2016). Within the soil domain, the time dynamics of the sum of water and ice contents is governed by water fluxes j_w according to the continuity equation:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{wi}}}{\partial t}}=-{\displaystyle \frac{\partial j_{\mathrm{w}}}{\partial z}}\text{.}\end{array}$$

The three water balance schemes differ in their representation of j_w, which generally consists of vertical water fluxes $j_{\mathrm{w}}^{\mathrm{v}}$ and fluxes due to evapotranspiration $j_{\mathrm{w}}^{\text{ET}}$ (or evaporation $j_{\mathrm{w}}^{\text{E}}$ and transpiration $j_{\mathrm{w}}^{\text{T}}$).

1. For the no flow scheme, j_w=0; i.e., the sum of water and ice contents is fixed for each grid cell (and thus only determined by the initialization) and not affected by rainfall, snowmelt, and evaporation. The no flow scheme therefore needs to be combined with either the temperature boundary condition or the surface energy balance with scheme 2 for evaporation (Sect. 2.2.2).

2. In the bucket scheme, the water in a grid cell is either immobile and bound to the soil matrix, or it flows downwards driven by gravity. The threshold between the two regimes is the user-defined field capacity, θ_fc. In the unsaturated domain, the vertical water flux is hence given by

   $$\begin{array}{}\text{(22)}& j_{\mathrm{w}}^{\mathrm{v}}=\left\{\begin{array}{ll}-K_{\mathrm{w}}& \text{for}{\mathit{\theta }}_{\mathrm{w}}>{\mathit{\theta }}_{\mathrm{fc}}\\ \mathrm{0}& \text{for}{\mathit{\theta }}_{\mathrm{w}}\le {\mathit{\theta }}_{\mathrm{fc}},\end{array}\right.\end{array}$$

   with K_w the hydraulic conductivity [m s⁻¹]. It is not the goal of the bucket scheme to reproduce the exact time dynamics of the infiltration event, and the hydraulic conductivity is broadly set to $K_{\mathrm{w}}=K_{\text{w,sat}}{\mathit{\theta }}_{\mathrm{w}}/\mathit{\varphi }$, with K_w,sat [m s⁻¹] the saturated hydraulic conductivity specified by the user. This in particular prevents or slows infiltration in ice-saturated ground in spring, when the water content is low. No vertical water flux occurs in the saturated domain unless water losses due to evapotranspiration must be compensated for, i.e., $j_{\mathrm{w}}^{\mathrm{v}}=max(-K_{\mathrm{w}},j_{\mathrm{w}}^{\text{ET}})$. In subsurface classes representing soil, the bucket scheme is combined with scheme 3 for evapotranspiration (Sect. 2.2.2). The water flux due to evapotranspiration from grid cell i is calculated from the latent heat flux Q_e as

   $$\begin{array}{}\text{(23)}& j_{\mathrm{w},i}^{\text{ET}}=-{\displaystyle \frac{f_i}{{\sum }_i{f}_i}}{\displaystyle \frac{Q_e}{L_{\text{lg}}{\mathit{\rho }}_{\mathrm{w}}}}\text{,}\end{array}$$

   with f_i calculated from the water availability coefficients ${\mathit{\alpha }}_i^{\mathit{\theta }}$ (see scheme 3 in Sect. 2.2.2 for the other variables) as

   $$\begin{array}{}\text{(24)}& f_i=f_{\text{tr}}{\mathit{\alpha }}_i^{\mathit{\theta }}\mathrm{\Delta }{z}_i{e}^{-d_i/d_{\text{tr}}}+f_{\text{ev}}{\mathit{\alpha }}_i^{\mathit{\theta }}\mathrm{\Delta }{z}_i{e}^{-d_i/d_{\text{ev}}}\text{.}\end{array}$$

   In essence, this ensures that a water flux corresponding to the weight of a grid cell in the calculation of the latent heat flux is extracted, taking water availability and exponential damping with depth into account. Note that CryoGrid 3 features a different bucket scheme, which does not treat soil water as a prognostic variable but redistributes water in the bucket after each time step (Nitzbon et al., 2019). As an example, a rain event leads to an instant increase of the water level in CryoGrid 3, while an infiltration front penetrating downwards with time is simulated by the subsurface classes available in the CryoGrid community model.

3. For water flow governed by the Richards equation (Richards, 1931), movement of water in unsaturated soils through vertical gradients of the matric and gravitational potentials is accounted for, in addition to gravity-driven flow in the saturated domain. In this water balance representation in the CryoGrid community model, evaporation is drawn from the uppermost grid cell, i.e., $j_{\mathrm{w},\mathrm{1}}^{\text{E}}=Q_e/\left(L_{\text{lg}}{\mathit{\rho }}_{\mathrm{w}}\right)$. If transpiration is considered (e.g., by the canopy scheme described in Stuenzi et al., 2021a), a grid cell weighting similar to Eq. (23) is used to compute the transpiration flux from each cell. The vertical water fluxes are calculated according to the Richards equation,

   $$\begin{array}{}\text{(25)}& j_{\mathrm{w}}^{\mathrm{v}}=-K_{\mathrm{w}}\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+\mathrm{1}\right)\text{,}\end{array}$$

   using the matric potential ψ, which also accounts for soil freezing (see Sect. 2.2.3). For the hydraulic conductivity, we use the classic formulation by van Genuchten (1980),

   $$\begin{array}{}\text{(26)}& \begin{array}{rl}{K}_{\mathrm{w}}& =K_{\text{w,sat}}{I}_{\text{ice}}({\mathit{\theta }}_{\mathrm{w}}/\mathit{\varphi }{)}^{\mathrm{0.5}}\\ & \cdot {\left(\mathrm{1}-{\left(\mathrm{1}-({\mathit{\theta }}_{\mathrm{w}}/\mathit{\varphi }{)}^{n/(n+\mathrm{1})}\right)}^{(n-\mathrm{1})/n}\right)}^{\mathrm{2}}\text{,}\end{array}\end{array}$$

   with an additional ice impedance factor $I_{\text{ice}}={\mathrm{10}}^{-\mathrm{\Omega }{\mathit{\theta }}_i/{\mathit{\theta }}_{\mathrm{wi}}}$ (defined as in Hansson et al., 2004) to account for the blocking of water-filled pores by ice (see Sect. 3.1.3). In some subsurface classes, the permeability of the subsurface material, k_w [m²], needs to be specified instead of the saturated hydraulic conductivity, which is calculated according to $K_{\text{w,sat}}=k_{\mathrm{w}}/\left({\mathit{\eta }}_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}g\right)$. Here, η_w is the temperature-dependent dynamic viscosity of water derived as ${\mathit{\eta }}_w=A\text{exp}(B/T+CT+D{T}^{\mathrm{2}})$, with T in Kelvin and coefficients $A,B,C$, and D as defined in Reid et al. (1987).

2.2.5 Excess ground ice

The CryoGrid community model comprises a subsurface class to simulate melting of excess ground ice, which results in the subsidence of the ground surface. It is based on the bucket water scheme with a freezing characteristic and surface energy balance (Sect. 2.2.2 to 2.2.4), and in most respects it is similar to the CryoGrid 3 excess scheme (Westermann et al., 2016) based on Lee et al. (2014). However, freezing and melting of excess water–ice are treated differently than pore water–ice that is contained in the sediment matrix. While the pore water–ice freezes and melts according to the soil freezing characteristic, the excess water–ice portion is always treated as free water; i.e., it undergoes phase change at T=0 ^∘C (see Sect. 2.2.3). Two additional state variables θ_χi and θ_χw denote the volumetric fractions of excess ice and water so that ${\mathit{\theta }}_{\mathrm{m}}+{\mathit{\theta }}_{\mathrm{o}}+{\mathit{\theta }}_{\mathrm{w}}+{\mathit{\theta }}_{\mathrm{i}}+{\mathit{\theta }}_{\mathit{\chi }\mathrm{i}}+{\mathit{\theta }}_{\mathit{\chi }\mathrm{w}}+{\mathit{\theta }}_{\mathrm{a}}=\mathrm{1}$. The initial excess ice content is specified by the user, and the excess ice fraction in a grid cell is unchanged (i.e., neither increases nor decreases) as long as its temperature is below 0 ^∘C.

Once excess ice melts, the excess water is mobilized and transported by the hydrology scheme, with an additional vertical water flux term $j_{\mathit{\chi }\mathrm{w}}^{\mathrm{v}}=K_{\mathrm{w}}$ directed upwards. This excess water is first routed between the excess water variables of adjacent grid cells, with grid cell thickness changing accordingly (i.e., shrinking for net outflow, expanding for net inflow). If excess water exists in an unsaturated grid cell (i.e., it contains a nonzero air content), water is moved from the excess water to the water phase, reducing the air content and leading to the grid cell thickness shrinking. In the uppermost grid cell, the excess water variable can be regarded as water pooling up above the surface, either due to melted excess ice routed upwards or from rainfall and melted snow. This excess water can either evaporate, be routed away laterally (Sect. 2.3), or evolve into a pond or lake represented by a water body class (see Sect. 2.2.7). The latter two depend on the user-defined model setup, which specifies what happens when the excess water in the first grid cell exceeds a threshold depth.

In the user interface, the amount of excess ice in a subsurface grid cell is specified as a fraction (χ) relative to the amount of soil without excess ice; i.e., χ=1 corresponds to a cell consisting of 50 % soil and 50 % excess ice.

2.2.6 Snow cover

Three stratigraphy classes representing snow are currently available in the CryoGrid community model, which all employ heat conduction and the free water freezing characteristic (Sect. 2.2.3) to calculate snow temperatures. New snow is added to the first grid cell, which is split into two cells when the “ice depth” (defined as θ_iΔz, equivalent to snow water equivalent for dry snow) exceeds 1.5 times a user-defined target value (with the lower cell containing the target ice depth and the upper cell the remaining part). When the ice depth in a grid cell decreases below half the target ice depth, it is merged with the grid cell below. Meltwater becomes mobile when the volumetric water content exceeds a user-defined field capacity (provided as a fraction of the porosity of the ice matrix, 1−θ_i) and can flow downwards, but also laterally if a corresponding lateral interaction class (Sect. 2.3) is selected. The dynamic interaction of snow classes with other stratigraphy classes, including their creation upon snowfall and removal when all snow has melted, is described in Sect. 2.1.3. Of the three snow classes, class (a) can be combined with subsurface classes with temperature boundary condition and classes (b) and (c) with subsurface classes with surface energy balance (Sect. 2.2.2).

• (a)

  Constant snow density, temperature boundary condition, and degree-day-based melt model. In this snow class, new snow is added with a user-defined constant density, which could, for example, be derived from field observations. Temperature calculations inside the snowpack rely on a temperature boundary condition, heat conduction, and the free water freeze curve. For snowmelt, a degree-day-based melt model is employed using a melt factor calculated from latitude and day of year (as in Obu et al., 2019). The product of day length and solar culmination angle (i.e., the highest sun angle above the horizon on a given day of year) is used as a measure of snowmelt activity, which scales the degree day melt factor between confining values of 0.002 and 0.012 $\mathrm{m}({}^{\circ }\mathrm{C}\mathrm{d}{)}^{-\mathrm{1}}$ water equivalent. The snowmelt is assigned to the uppermost grid cell from which meltwater is removed once it exceeds the pore space, without infiltrating the snowpack.

• (b)

  Constant snow density, surface energy balance, and snow hydrology (bucket scheme). This snow class largely follows the snow parameterization of the CryoGrid 3 model, as described in detail in Westermann et al. (2016). As for (a), snowfall is added with prescribed density, but the surface energy balance is used as the upper boundary condition with a transient albedo that decreases from a maximum value for fresh snow to a minimum value for old snow, with decrease rates depending on whether the snow is dry or wet (as inferred from the liquid water content of the first snow grid cell, see Sect. 2.2.3). Furthermore, shortwave radiation penetrates the snowpack following de Beer's law with a defined extinction coefficient, and sublimation or resublimation derived from the latent heat flux is extracted from or added to the uppermost grid cell. The snow hydrology follows the bucket scheme (Sect. 2.2.4), with water from both rainfall and snowmelt percolating downwards when the water content exceeds the field capacity. In the simple snow cover module, refreezing of meltwater is the only process that can alter the density of a snow layer.

• (c)

  Snow microphysics, surface energy balance, and snow hydrology (bucket scheme). Introduced within CryoGrid 3 by Zweigel et al. (2021), this snow class is based on the Crocus snow scheme (Vionnet et al., 2012), including transient snow grain property and density evolution. The defining equations and parameterizations are largely identical to the ones described in Vionnet et al. (2012), so we only provide a brief description, concentrating on aspects treated differently. As for the previous class, the energy transfer at the upper boundary is prescribed according to the surface energy balance, but relies on spectrally resolved calculation of albedo and shortwave penetration and absorption in each grid cell. Snowfall is added with properties (density, grain size, dendricity, and sphericity) derived from air temperature and wind speed, and the snow within each grid cell evolves and metamorphoses based on internal temperature gradient and water content. In particular, snow compaction due to the weight of overlying snowpack is derived by computing the snow viscosity, which is computed for each snow layer, parameterized as a function of snow density, temperature, liquid water content, and snow grain size. The class also accounts for the impact of wind drift on snow grain properties and density, which in particular leads to compaction and density increase of the uppermost snow layers. If the uppermost snow grid cell features liquid water, both evaporation and sublimation are calculated and their fractions linearly interpolated between θ_w=0 (all sublimation) and ${\mathit{\theta }}_{\mathrm{w}}={\mathit{\theta }}_{\text{evap}}^{\text{snow}}$ (all evaporation). The threshold ${\mathit{\theta }}_{\mathrm{w}}={\mathit{\theta }}_{\text{evap}}^{\text{snow}}$ is set to twice the field capacity, but this should be revisited in future studies. For evaporation, the corresponding amount of water is extracted, while the same happens for the ice phase for sublimation. The original Crocus setup described in Vionnet et al. (2012) is associated with a variety of model parameters, some of which Royer et al. (2021) suggested revising to better reproduce snowpack characteristics in the Arctic. These are in particular related to the wind speed dependence of the new snow density and the compaction dynamics due to wind drift, as well as parameterization of the snow thermal conductivity. In the snow class, it is possible to choose between the parameter sets for the original (according to Vionnet et al., 2012) and the “Arctic Crocus” (according to Royer et al., 2021) or to independently adjust the parameters in question. With this, the performance of either scheme within the CryoGrid community model can be evaluated against observations (see Sect. 3.2).

While the most process-rich scheme (c) is generally expected to deliver a superior performance, it is also more sensitive to biases in the model forcing, especially the wind speed, which strongly impacts snow density. In some cases, it is therefore preferable to employ the simpler schemes (a) or (b), especially if field measurements constraining the snow density are available.

2.2.7 Water bodies

In water bodies, heat transfer from the surface is strongly different between the ice-free and ice-covered seasons. The key characteristics of this seasonal asymmetry were conceptualized in CryoGrid 3 (see Westermann et al., 2016) for the highly relevant case of shallow water bodies. During the ice-free season, the water column is assumed to be fully mixed due to wind action, while a stable, temperature-driven stratification forms below an ice cover, with water at 0 ^∘C at the ice interface less dense than the warmer water at deeper layers. Within this water column, heat conduction is the main pathway of energy transfer with the relatively small thermal conductivity of water (K_h,w=0.57 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) severely restricting energy losses of the ground below the water column. The CryoGrid community model provides a water body class based on the CryoGrid 3 model physics. In fact, the two seasonal regimes are implemented as two separate stratigraphy classes, which mutually create and destroy each other upon defined conditions (see below). In the “ice-free class”, the entire water column is simply represented by a single grid cell which assumes well-mixed conditions. The surface energy balance is applied at the upper boundary, and shortwave radiation penetrates the water column with a bulk (i.e., not spectrally resolved) absorption coefficient. Both rain and snowfall are added to this grid cell with their respective enthalpy e, leading to a change in both temperature and the grid cell thickness (and thus the water level). When the enthalpy reaches e=0 (which is ensured by the time-stepping scheme, Sect. 2.2.9), an ice cover forms and the ice-free class is exchanged by the ice-covered class. In this process, all state variables are split to the pre-defined model grid so that the surface energy balance is now applied to the uppermost grid cell, which can subsequently freeze according to the free water freezing characteristics. While energy transfer in the water and ice column is by means of heat conduction, grid cells are reordered after each time step, with fully frozen cells (i.e., the ice cover) always on top and the unfrozen cells arranged by their temperature-dependent densities (according to Kell, 1975). When all ice has melted, i.e., e≥0 for all grid cells, the grid cells are merged into a single grid cell and the ice-free class resumes. Note that the FLake water body scheme presented in CryoGrid 3 (Langer et al., 2016) is not yet available in the CryoGrid community model, but will be implemented in the future as an additional stratigraphy class.

2.2.8 Glaciers

The CryoGrid community model contains a glacier class, which consists of layers of pure ice (using the free water freezing characteristic) with a user-defined constant ice thickness. Energy transfer is governed by heat conduction, with the surface energy balance as the upper boundary condition. The scheme is usually coupled to snow scheme (b) or (c) (Sect. 2.2.5), which allows the buildup of a seasonal snow layer for simulations of the ablation area or, if run over longer periods, the buildup of a firn layer for simulations of the accumulation zone. The densification scheme in Crocus (Vionnet et al., 2012, Sect. 2.2.5) has been implemented into several models for simulations of snow and firn densification on glaciers (e.g., Cullather et al., 2016; Langen et al., 2017; Verjans et al., 2019), which have been successfully applied for, e.g., the Greenland ice sheet, Antarctica, and Icelandic glaciers (e.g., Agosta et al., 2019; Fettweis et al., 2017; Schmidt et al., 2017). Similar applications are conceivable in the CryoGrid community model, with the glacier class representing ice and the Crocus-based snow class firn and seasonal snow cover.

Water cannot infiltrate the ice, so any water which does not refreeze during the time step will run off instantaneously if there is no snow on the glacier surface. Likewise, if snow is present, liquid water in the snow class will build up above the glacier ice, where it can eventually refreeze or run off, depending on the selected lateral interaction class.

If additional ice is added to the surface grid cell, either from refreezing of rainwater or deposition, mass is advected downwards to ensure constant water equivalent water thickness. Similarly, if mass is removed from the module by runoff, evaporation, or sublimation, mass is advected up, with the lowest model layer receiving additional mass from an infinite ice reservoir below the model domain. This reservoir is assumed to have the same temperature as the lowest model layer. The movement of mass within the model column is accompanied by a vertical transfer of sensible heat. A similar approach has previously been used in other glacier models which do not account for glacier flow (e.g., Langen et al., 2017) in order to prevent glacier areas with highly negative mass balance from disappearing during spin-up. However, the glacier class can also be employed without this option and coupled to a subsurface class representing subglacial sediments. In this case, the glacier can completely melt away, exposing the ground below, which, for example, offers the possibility to study glacier–permafrost interactions (Sect. 4.6).

2.2.9 Numerical implementation

For prognostic variables like enthalpy and water content, we use the integral form of the respective continuity equations (Eqs. 13 for heat and 21 for water). For a scalar volume-normalized quantity s (e.g., e for enthalpy and θ_w for water contents), the time change within a volume V is explained by fluxes j_s across surface Ω_V of V, with normal vector n:

$$\begin{array}{}\text{(27)}& \underset{V}{\int }{\displaystyle \frac{\partial s}{\partial t}}\mathrm{d}V+\underset{{\mathrm{\Omega }}_V}{\oint }{\mathit{j}}_{\mathrm{s}}\cdot \mathit{n}\mathrm{d}{\mathrm{\Omega }}_V=\mathrm{0}\text{.}\end{array}$$

The numerical implementation is based on finite differences with grid cells (index i increasing downwards, vertical thickness Δz [m], and area A [m²]) within which s is considered constant. Furthermore, for one-dimensional simulations, only vertical fluxes through the upper and lower boundary of the grid cell have to be considered so that the continuity equation for grid cell i simplifies to

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\partial }{\partial t}}{S}_i:={\displaystyle \frac{\partial }{\partial t}}\left(A\mathrm{\Delta }{z}_i{s}_i\right)=-j_{\mathrm{s}}^{i-\mathrm{1}\leftrightarrow i}+j_{\mathrm{s}}^{i\leftrightarrow i+\mathrm{1}}\text{,}\end{array}$$

i.e., the time derivative of the bulk quantity S_i is simply obtained from the fluxes $j_{\mathrm{s}}^{i-\mathrm{1}\leftrightarrow i}$ (defined positive when directed upwards) across the interfaces between grid cells i−1/i and i/i+1. For the numerical implementation in the CryoGrid community model, it is therefore of practical advantage to use the extensive bulk quantities as state variables and not the volume-normalized quantities (which are employed as direct model state variables in, e.g., CryoGrid 2 and 3) for which the defining differential equations in the previous sections are provided. In the CryoGrid community model, each stratigraphy class covers an explicit area A [m²], and the model state variables for mineral, organic, water, and ice contents become volumes: ϕ_m=AΔzθ_m [m³], ϕ_o=AΔzθ_o [m³], ϕ_w=AΔzθ_w [m³], ϕ_i=AΔzθ_i [m³]. Likewise, the bulk value for the enthalpy for each grid cell is used, E=AΔz e [J]. For time integration of Eq. (28), we use a simple first-order forward Euler scheme as in CryoGrid 3 (Westermann et al., 2016), i.e.,

$$\begin{array}{}\text{(29)}& S_i(t+\mathrm{\Delta }t)=S_i\left(t\right)+\mathrm{\Delta }t(-j_{\mathrm{s}}^{i-\mathrm{1}\leftrightarrow i}+j_{\mathrm{s}}^{i\leftrightarrow i+\mathrm{1}})\text{.}\end{array}$$

Stability and accuracy are guaranteed by selecting small enough time steps with conditions specifically designed for each state variable, the particular requirements of the model physics of each stratigraphy class, and the typical orders of magnitude (and timescales of change) of model forcing and parameters. For enthalpy, for example, a maximum change of volume-normalized enthalpy e between grid cells is defined by the user, and time steps Δt are adjusted to not exceed this value for any grid cell. Therefore, small time steps are generally required if large fluxes occur, slowing down computation. In a similar manner, changes in soil water content between time steps can be limited, while it is also possible to prevent “overfilling” of a grid cell (so that the water content exceeds the pore space) by limiting the time step accordingly. Another example is the water body class (Sect. 2.2.7) for which the time step calculation guarantees that the condition e=0 (which triggers the switch between ice-free and ice-covered water body classes, Sect. 2.2.7) is exactly met. In addition, a maximum time step can be defined to satisfy the CFL (Courant–Friedrichs–Lewy) condition (Courant et al., 1928) for the parameters (e.g., thermal conductivities and heat capacities) and grid cell sizes of the simulation setup. The fluxes of both heat and water have the general form

$$\begin{array}{}\text{(30)}& j_{\mathrm{s}}=-{\mathit{\kappa }}_{\mathrm{s}}\left(z\right){\displaystyle \frac{\partial \mathit{\sigma }\left(z\right)}{\partial z}}\text{.}\end{array}$$

Conductivities κ_s are defined for individual grid cells, while j_s in the finite-difference scheme is expressed in terms of fluxes between grid cells i−1 and i, $j_{\mathrm{s}}^{i-\mathrm{1}\leftrightarrow i}$. For this reason, the effective conductivities governing the flux between grid cells are calculated as series of the resistances $\mathrm{1}/{\mathit{\kappa }}_{\mathrm{s},i-\mathrm{1}}$ and $\mathrm{1}/{\mathit{\kappa }}_{\mathrm{s},i}$ using half of the grid cell thicknesses Δz:

$$\begin{array}{}\text{(31)}& j_{\mathrm{s}}^{i-\mathrm{1}\leftrightarrow i}={\displaystyle \frac{\mathrm{2}{\mathit{\kappa }}_{\mathrm{s},i-\mathrm{1}}{\mathit{\kappa }}_{\mathrm{s},i}}{{\mathit{\kappa }}_{\mathrm{s},i-\mathrm{1}}\mathrm{\Delta }{z}_i+{\mathit{\kappa }}_{\mathrm{s},i}\mathrm{\Delta }{z}_{i-\mathrm{1}}}}({\mathit{\sigma }}_i-{\mathit{\sigma }}_{i-\mathrm{1}})\text{.}\end{array}$$

In the uppermost and lowermost grid cells of the CryoGrid stratigraphy (which are in different stratigraphy classes if the stratigraphy consists of more than one class), the fluxes derived at the upper and lower boundaries (e.g., F_ub and F_lb for heat) are added. To connect adjacent stratigraphy classes, corresponding fluxes $j_{\mathrm{s}}^{N_u\leftrightarrow {\mathrm{1}}_l}$ are calculated by the interaction class (Sect. 2.1.2), with grid cell N_u being the lowermost grid cell of the upper stratigraphy class, while grid cell 1_l is the uppermost grid cell of the lower stratigraphy class. $j_{\mathrm{s}}^{N_u\leftrightarrow {\mathrm{1}}_l}$ depends on the state variables and model physics of both classes involved. It can be of the same form as Eq. (31), but also $j_{\mathrm{s}}^{N_u\leftrightarrow {\mathrm{1}}_l}=\mathrm{0}$, for example the water flux between a stratigraphy class with water balance and one without water balance.

2.3 Lateral interactions with an external environment

With CryoGrid 3, several studies were presented that simulate lateral exchange of energy and matter (heat, water or snow), either with external reservoirs for single-tile simulations (Martin et al., 2019; Langer et al., 2016) or between different CryoGrid stratigraphies for three-dimensional multi-tile configurations (Martin et al., 2021; Nitzbon et al., 2019, 2020, 2021; Zweigel et al., 2021). In the CryoGrid community model, we extend these possibilities by providing a standardized interface to implement a variety of lateral interactions, which are compatible with the CryoGrid stratigraphy consisting of a stack of classes (Fig. 5). This functionality is accomplished by two further types of classes, “lateral classes” and “lateral interaction classes”, both of which are selected in the TILE_1D_standard class. The choice of the lateral class determines whether the CryoGrid stratigraphy interacts with external (and static) reservoirs (LATERAL_1D) or whether several CryoGrid stratigraphies interact with each other (LATERAL_3D), which corresponds to the “laterally coupled tiling” demonstrated in CryoGrid 3. In this study, we focus on interactions with external reservoirs, while laterally coupled tiling with LATERAL_3D will be described in a separate study in the future. Other than for the time integration in the vertical CryoGrid stratigraphy, lateral fluxes are added or subtracted after a fixed, user-determined interaction time step (which is a parameter in the lateral class), so they are not part of the time integration scheme of the regular stratigraphy. This is largely due to the computation-related requirements of laterally coupled tiling (which requires parallel computing, see Nitzbon et al., 2019), but constant time steps Δt^lat are also employed for one-dimensional simulations with external reservoirs to ensure consistency. In general, lateral fluxes exchanged with external reservoirs should be small compared to the corresponding vertical fluxes within the CryoGrid stratigraphies, which means that interaction time steps can be significantly longer than the typical time steps required for the vertical integration. The lateral interaction time step must be selected by the user, seeking a balance between runtime and model accuracy and stability. Typical lateral time steps for the applications presented in Sect. 3.2 are between half an hour and 6 h.

The lateral class sets up the environment for the lateral interaction classes, which represent the actual model physics of the lateral interactions with external reservoirs. Lateral interaction classes can be combined with each other, with the lateral class calling them one after the other in the order provided by the user (as a list in TILE_1D_standard).

Figure 5 Schematic representation of CryoGrid stratigraphy (for the example of fully developed snow cover) interacting with a lateral interaction (IA) class (see text). Note that it is specific to each stratigraphy class if and how it is modified by the lateral interaction class. In this example, stratigraphy class 2 is unaffected, while the snow class and class 1 are modified by the lateral interaction class.

2.3.1 Lateral coupling to heat reservoir

The CryoGrid stratigraphy can be laterally coupled to an external heat reservoir, which can be used to mimic the thermal impact of infrastructure, water bodies (similar to Langer et al., 2016), or more general adjacent areas with strongly different ground thermal regimes, e.g., at the edge of a permafrost-underlain peat plateau (Martin et al., 2021). The heat reservoir is characterized by a user-defined, constant temperature T^lat and lateral distance d^lat [m] from the CryoGrid stratigraphy. Furthermore, a lower and an upper elevation for the heat reservoir must be provided, which makes it possible to confine the effect of the heat reservoir to a part of the stratigraphy. Furthermore, several heat reservoirs with different temperatures as well as upper and lower elevations can be combined, which achieves a similar effect as the coupling to a temperature stratigraphy (Langer et al., 2016). If a grid cell i is located between the lower and upper bounds of the heat reservoir, the lateral heat flux is calculated as

$$\begin{array}{}\text{(32)}& j_{\text{hc,i}}^{\text{lat}}=-K_{\text{h,i}}{\displaystyle \frac{T_i-T^{\text{lat}}}{d^{\text{lat}}}}\text{,}\end{array}$$

with K_h,i the thermal conductivity of grid cell i and T_i its temperature. The change in bulk enthalpy in grid cell i over time step Δt^lat is given by

$$\begin{array}{}\text{(33)}& \mathrm{\Delta }{E}_i=\mathrm{\Delta }{t}^{\text{lat}}\mathrm{\Delta }{z}_i{l}_{\mathrm{c}}^{\text{lat}}{j}_{\text{hc,i}}^{\text{lat}}\text{,}\end{array}$$

with $l_{\mathrm{c}}^{\mathrm{lat}}$ [m] the lateral contact length and $l_{\mathrm{c}}^{\text{lat}}\mathrm{\Delta }{z}_i$ [m²] the cross section through which the lateral heat flux occurs.

2.3.2 Lateral water transport

Surface water removal. When snow melts or rain falls on a saturated soil column, water will either pool up on the surface or be lost as surface runoff. The exact way to treat this surface water depends on the particular stratigraphy class, with most classes removing surface water (but storing it in a state variable), while, for example, the excess ice class (Sect. 2.2.5) can explicitly represent a surface water pool. The CryoGrid community model provides a lateral interaction class that constantly removes any surface water and stores the accumulated surface runoff, irrespective of the way surface water is treated in the uppermost stratigraphy class. This in particular makes it possible to generate an unbroken time series of surface runoff, which comprises both snow classes (i.e., snowmelt runoff) and subsurface classes during the snow-free season.

Overland flow. For the excess ice class (Sect. 2.2.5), standing surface water can be represented as excess water in the uppermost grid cell, with water depth given as d_w=θ_χwΔz [m]. Instead of removing surface water like with the previous interaction class, surface water can also be removed by overland flow governed by the Gauckler–Manning equation (Gauckler, 1867; Manning, 1891) as

$$\begin{array}{}\text{(34)}& j_{\text{w,1}}^{\text{lat}}=-G{d}_{\mathrm{w}}^{\mathrm{2}/\mathrm{3}}({\mathit{\delta }}^{\text{lat}}{)}^{\mathrm{1}/\mathrm{2}}\text{,}\end{array}$$

with G [${\mathrm{m}}^{\mathrm{1}/\mathrm{3}}$] the Gauckler–Manning coefficient and δ^lat [–] the local gradient which drives the surface flow. The change in bulk water content in grid cell 1 over time step Δt^lat is given by

$$\begin{array}{}\text{(35)}& \mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{w},\mathrm{1}}=\mathrm{\Delta }{t}^{\text{lat}}{d}_{\mathrm{w}}{l}_{\mathrm{c}}^{\mathrm{lat}}{j}_{\mathrm{w},\mathrm{1}}^{\text{lat}}\text{,}\end{array}$$

with $l_{\mathrm{c}}^{\mathrm{lat}}{d}_{\mathrm{w}}$ [m²] the cross section through which the lateral water flux occurs ($l_{\mathrm{c}}^{\text{lat}}$ is the lateral contact length, see above). Accordingly, the change in bulk enthalpy due to water advection is given by

$$\begin{array}{}\text{(36)}& \mathrm{\Delta }{E}_{\mathrm{1}}=\mathrm{\Delta }{t}^{\text{lat}}{d}_{\mathrm{w}}{l}_{\mathrm{c}}^{\mathrm{lat}}{j}_{\mathrm{w},\mathrm{1}}^{\text{lat}}{c}_{\mathrm{w}}{T}_{\mathrm{1}}\text{.}\end{array}$$

Seepage face. For stratigraphy classes with water balance (Sect. 2.2.4, schemes 2 and 3), a seepage face lateral boundary condition is implemented as a lateral interaction class, leading to drainage in the saturated domain of the soil column. The lateral interaction class first determines the elevation of the water table, z_wt, and then removes a lateral water flux,

$$\begin{array}{}\text{(37)}& j_{\mathrm{w},i}^{\text{lat}}=-K_{\mathrm{w},i}{\displaystyle \frac{z_{\text{wt}}-z_i}{d^{\text{lat}}}},\end{array}$$

for grid cells i with elevations z_i<z_wt. Here, K_w,i is the hydraulic conductivity and d_lat the lateral distance to the seepage face, which determines the strength of the drainage. In unsaturated grid cells, no outflow occurs. Note that the water table elevation is tracked across stratigraphy classes; e.g., when the water table is located in a water body class, it also governs the outflow from the subsurface class below. The seepage face always leads to outflow, and it is possible to define the upper and lower elevations for the domain through which outflow occurs. The change in bulk water content in grid cell i over time step Δt^lat is given by

$$\begin{array}{}\text{(38)}& \mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{w},i}=\mathrm{\Delta }{t}^{\text{lat}}\mathrm{\Delta }{z}_i{l}_{\mathrm{c}}^{\text{lat}}{j}_{\mathrm{w},i}^{\text{lat}},\end{array}$$

with $l_{\mathrm{c}}^{\text{lat}}\mathrm{\Delta }{z}_i$ [m²] the cross section through which the lateral water flux occurs. Accordingly, the change in bulk enthalpy due to water advection is given by

$$\begin{array}{}\text{(39)}& \mathrm{\Delta }{E}_i=\mathrm{\Delta }{t}^{\text{lat}}\mathrm{\Delta }{z}_i{l}_{\mathrm{c}}^{\text{lat}}{j}_{\mathrm{w},i}^{\text{lat}}{c}_{\mathrm{w}}{T}_i\text{.}\end{array}$$

Water reservoir. Similar to seepage flow, hydrological coupling to an external water reservoir is possible, located at elevation z^lat and lateral distance d^lat. The lateral water flux is calculated as

$$\begin{array}{}\text{(40)}& j_{\mathrm{w},i}^{\text{lat}}=-K_{\mathrm{w},i}{\displaystyle \frac{z_i-z^{\text{lat}}}{d^{\text{lat}}}}\text{,}\end{array}$$

with parameters as described above. As for seepage flow, water flow is restricted to the saturated zone, but inflow can occur if z_i<z^lat. This lateral interaction class can therefore be used to keep the soil water table at a certain level, with in- and outflow depending on rainfall and evapotranspiration. While outflowing water has the temperature of the respective grid cell, inflowing water can optionally be assigned a (constant) reservoir temperature T^lat, which is taken into account in terms of heat advection through water. In this case, the change in bulk enthalpy is given by

$$\begin{array}{}\text{(41)}& \mathrm{\Delta }{E}_i=\mathrm{\Delta }{t}^{\text{lat}}\mathrm{\Delta }{z}_i{l}_{\mathrm{c}}^{\text{lat}}{j}_{\mathrm{w},i}^{\text{lat}}{c}_{\mathrm{w}}{T}^{\text{lat}}\text{,}\end{array}$$

while the change in bulk water content Θ_w,i is the same as for the seepage face (see Eq. 38).

3 Results

3.1 Benchmarking of selected model components

3.1.1 Step change in upper boundary temperature

To document the basic numerical performance of the model framework, we model the temperature response of an infinite half-space to a step change in temperature at the upper boundary. An analytical solution to this problem is available (Carslaw and Jaeger, 1959) and is given by

$$\begin{array}{}\text{(42)}& T(x,t)=T_{\text{init}}+(T_{\text{ub}}-T_{\text{init}})\text{erfc}\left(\sqrt{{\displaystyle \frac{c{x}^{\mathrm{2}}}{\mathrm{4}{K}_{h}t}}}\right)\text{,}\end{array}$$

where T(x,t) [^∘C] is the temperature at time t [s] and depth x [m] below the surface, and T_init [^∘C] and T_ub [^∘C] are the initial domain temperature and temperature applied at the upper boundary, respectively. c and K_h are the heat capacity [$\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$] and thermal conductivity [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] of the medium, and erfc is the complementary error function. The formulation is valid for a homogenous material with no variation in thermal properties in space and time, as well as without phase change.

In CryoGrid we simulate the response using a stratigraphy class with zero heat flux at the lower boundary and temperature boundary condition (Sect. 2.2.2) at the upper boundary. A 100 m deep model domain is selected, discretized with a spacing of 0.1 m and initialized with a constant temperature of T_init=1 ^∘C. At t=0, an upper boundary temperature of 10 ^∘C is applied, and the simulation is allowed to run for 5 years.

Figure 6 shows the model and analytical responses for selected time steps. Using a maximum change in internal energy per time step of 50 kJ m⁻³ (which is the default CryoGrid setting used in all following simulations), the maximum absolute difference between the analytical and numerical simulations is less than 0.003 ^∘C. The maximum error is reduced to less than 0.0003 ^∘C by decreasing the maximum energy step to 10 kJ m⁻³.

Figure 6 Comparison of analytical and numerical simulations of a step change in temperature at the upper boundary of a homogeneous half-space. Results are shown for the initial condition (t=0 s) and for 24 h, 100 d, and 1, 2, and 5 years after the step change. Numerical results are plotted for every fourth grid cell.

3.1.2 Stefan problem

In 1889, the Slovene physicist Josef Stefan published an analytical solution to the one-sided freezing problem and applied it to sea-ice formation in the Arctic Ocean (Stefan, 1891). The solution was reformulated and adapted for modeling the thawing of frozen soils by Nixon and McRoberts (1973). To obtain an analytical solution, the following simplifying assumptions were introduced: (1) the heat capacity of the medium is neglected, (2) the temperature in the medium below the freezing front does not affect the rate of thaw, (3) the temperature distribution in the thawed medium is linear, and (4) a constant surface temperature is applied. With these assumptions, the thaw depth as a function of time is given by (Nixon and McRoberts, 1973)

$$\begin{array}{}\text{(43)}& d_{\text{thaw}}=\sqrt{{\displaystyle \frac{\mathrm{2}{K}_h{T}_{\text{ub}}}{L_{\text{sl}}^{\text{vol}}{\mathit{\theta }}_{\mathrm{wi}}}}}\sqrt{t}\text{,}\end{array}$$

where d_thaw [m] is the thaw depth at time t [s], K_h [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] the thermal conductivity of the medium (considered constant with temperature), T_ub [^∘C] the upper boundary temperature, $L_{\text{sl}}^{\text{vol}}$ [J m⁻³] the volumetric latent heat of fusion, and θ_wi [–] the initial volumetric water (for freezing) or ice (for thawing) contents.

This analytical solution to the problem of top-down thawing of frozen soils, which is often referred to as the Stefan solution, has been used here to benchmark the CryoGrid model for a situation with phase change. We used a stratigraphy class (GROUND_freeW_ubT) with a free water soil freezing characteristic (Sect. 2.2.3), temperature boundary condition (Sect. 2.2.2), and zero heat flux at the lower boundary. The domain was initialized with an initial temperature of −0.02 ^∘C, which is a good approximation to satisfy assumption 2. The volumetric ice content was chosen to be 0.3 and the mineral content 0.7. The upper boundary temperature was fixed at 1 ^∘C to minimize the effect of the heat capacity of the medium (to align with assumption 3). The domain was discretized with a grid spacing of 0.01 m from the surface to 3 m depth and 0.1 from 3 m depth to the lower boundary at 100 m depth. The model was run with a constant surface temperature for 5 years, and the numerical and analytical results are compared in Fig. 7.

Figure 7 Comparison of analytical and numerical simulations of one-sided thawing of a frozen soil column. The soil temperature is initialized at −0.02 ^∘C, and an upper boundary temperature of 1 ^∘C is applied for a period of 5 years. The domain shown is discretized with a node spacing of 0.01 m.

The numerical and analytical solutions are nearly identical, with the numerically derived thaw depths slightly shallower than the ones derived from the analytical solution. The fact that the numerical formulation accounts for the heat capacity of the medium and thus storage of sensible heat, while the analytical solution does not, accounts for the very small offset observed towards the end of the period.

3.1.3 Mizoguchi (1990) experiment on cryosuction

The Mizoguchi (1990) experiment is a classic model benchmark for cryosuction, i.e., the redistribution of soil water during soil freezing. A number of unfrozen, 0.2 m long soil columns with constant soil water content are exposed to one-sided freezing, while the other side is kept insulated. After 12, 24, and 50 h, the combined content of water and ice is determined gravimetrically for one of the columns so that the water redistribution over time can be followed. The experimental conditions including suggestions for model setup are presented in detail in Hansson et al. (2004). We modified the upper and lower boundary conditions of the subsurface class with water balance according to Richards equation (Sect. 2.2.4) accordingly, setting the heat flux at the lower boundary (0.2 m) to zero. The heat flux at the upper boundary is set proportional to the temperature difference between the first grid cell and the coolant (at −6 ^∘C). As in Hansson et al. (2004), we tried two settings for the proportionality coefficient (i.e., the convective heat transfer coefficient), one constant at 28 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ (denoted as the “linear heat transfer scenario”) and one decreasing nonlinearly from 40 (above 0 ^∘C) to 10 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ (below −4 ^∘C) as a function of the first grid cell temperature squared (denoted as the “nonlinear heat transfer scenario”). The soil porosity was set to 0.535 with zero organic content, and we use an initial water content of 0.345 to match the initial state (0 h, see Fig. 8) depicted in Hansson et al. (2004). Similar to Painter (2011), the soil permeability is set to $\mathrm{3.25}\times {\mathrm{10}}^{-\mathrm{13}}$ m², which corresponds to a saturated hydrological conductivity of $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ (Hansson et al., 2004) for a water viscosity at room temperature. For the silty soil employed in the experiment, we used the van Genuchten parameters α=1.11 m⁻¹ and n=1.48 (Hansson et al., 2004). For the thermal conductivity, we use the parameterization of Cosenza et al. (2003), with a thermal conductivity of the mineral fraction of 2.0 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. This roughly reproduces the measured frozen and thawed thermal conductivities for volumetric water contents of 0.4 and 0.3, while thermal conductivities for lower water and ice contents are overestimated (note that such low water contents are not reached in the simulation).

Figure 8 shows the best fit to the measurements, which is achieved with the nonlinear heat transfer scenario. Except for small deviations near the freeze front for the 12 and 24 h states, CryoGrid manages to reproduce the measurements very well. Using the linear heat transfer scenario (Sect. S3, Fig. S1) leads to only small changes in the results, with an almost equally good fit. Likewise, changing the thermal conductivity parameterization to the one used in CLM4.5 modified the results only marginally (not shown). Most importantly, the parameterization of the hydrological conductivity had a pronounced influence on the results, especially the choice of the additional ice impedance factor I_ice (Sect. 2.2.4, Eq. 26). The best fit (Fig. 8) is achieved for Ω=5 (i.e., the hydrological conductivity is decreased by a factor of 10⁻⁵ when water contents approach zero). With the default factor of Ω=7 (Dall'Amico et al., 2011), water redistribution is considerably weaker, resulting in a notably worse fit to the measurements (Sect. S3, Fig. S2). While this should be investigated in more detail, we have set Ω=5 in all simulations with the Richards equation in this study.

Figure 8 Simulated sum of volumetric water and ice content (blue lines) vs. measurements (crosses, digitized from Hansson et al., 2004) for the Mizoguchi (1990) experiment for 0, 12, 24, and 50 h freezing time; nonlinear heat transfer scenario, ice impedance factor calculated with Ω=5 (see Eq. 26).

3.1.4 Accelerated spin-up to reach steady-state temperature profile

In most real-world examples of thermal simulations, the initial temperature profile is not known and must be estimated by a model spin-up, i.e., by running the model for a certain time period (denoted the spin-up period) until the simulated temperature profile becomes independent of the initial profile (i.e., it is only determined by the model forcing applied at the upper and lower boundaries). While this state is usually reached within a few years in the uppermost meters, it takes much longer for deeper layers. For climate change simulations at a centennial timescale with a typical model domain depth of 100 m (e.g., Westermann et al., 2016), a spin-up of several hundred years can be necessary, thus requiring significant additional computation time. However, reliable spin-up model forcing this long back in time is often not available, so it is approximated by repeatedly looping a shorter period, often the first part of the regular model forcing. For this practically relevant case, the CryoGrid community model offers the possibility to considerably accelerate model spin-up by estimating the equilibrium temperature profile for the spin-up period, broadly following the procedure outlined in Westermann et al. (2013, 2017) for the CryoGrid 2 model:

1. use the TTOP model (Sect. 2.2.2) to obtain a first estimate for the temperature at the top of the permafrost and bottom of the seasonally frozen layer, T_TOP,1;

2. estimate a plausible maximum thaw depth (or maximum freeze depth in case of seasonally frozen ground), apply T_TOP,1 above, and compute a steady-state temperature profile below with the lower boundary heat flux and thermal conductivities computed for the respective temperatures;

3. run the model twice for the selected spin-up period and store freezing and thawing degree days (FDD and TDD) for each model grid cell for the second iteration;

4. determine the maximum thaw depth (freeze depth in case of seasonally frozen ground), as well as the temperature at the top of the permafrost (top of permanently unfrozen layer for seasonally frozen ground), T_TOP,2, from FDD and TDD;

5. repeat step 2 with T_TOP,2 and the maximum thaw depth (maximum depth of seasonal frost) calculated in the previous step;

6. run the model once for the selected spin-up period prior to starting the regular simulations.

As only steps 3 and 6 require significant computation, the total runtime of the accelerated spin-up approximately corresponds to 3 times the normal model runtime for the spin-up period. Figure 9 displays an example of the accelerated spin-up procedure using a period of 10 years. For comparison, a classic spin-up is performed, starting with a constant initial temperature of −8 ^∘C throughout the entire profile (which is clearly too cold, considering the model forcing). While the classic spin-up takes 200 to 400 years (i.e., 20 to 40 iterations of the 10-year period) to reach a reasonable approximation of the targeted equilibrium state, the accelerated spin-up requires only 30 model years (three iterations of spin-up period) to reach a similar performance. Figure 9 shows that the accelerated spin-up is not exact, with the temperature at 100 m depth about 0.1 ^∘C warmer than the classic spin-up after 400 years, which is still cooling slightly at this point. The reason for this is that the thermal conductivity in the uppermost permafrost layers still fluctuates annually with temperatures, giving rise to a small additional thermal offset (Romanovsky and Osterkamp, 1995). However, this small deviation of the accelerated spin-up is not relevant for real-world applications, as the true temperature profile does not correspond to the equilibrium state for the spin-up period anyway, but is determined by the real climate conditions in the past. Therefore, the accelerated spin-up procedure in the CryoGrid community model provides a compromise between computation time and initialization accuracy that is fully sufficient for many applications.

Figure 9 Comparison of accelerated and classic spin-up for a 10-year period from 31 July 1980 to 31 July 1990; simulated temperature profile on 31 July for accelerated spin-up (black, achieved with 30 simulation years) and classic spin-up (blue) after 0, 30, 100, 200, and 400 simulation years. Setting: Svalbard forcing (see Sect. 3.2) and ground stratigraphy, Crocus-based snow class, 100 m deep model domain with geothermal heat flux of 50 mW m⁻² at the lower boundary.

3.2 Example simulations for Svalbard

To demonstrate different configurations of modular CryoGrid stratigraphies, we perform simulations for the well-studied Bayelva permafrost research site outside Ny-Ålesund, Svalbard (Boike et al., 2018). A soil and climate monitoring station established in 1998 provides records of active layer temperatures, meteorological variables, and snow depth. The Bayelva site is situated on top of a hill, within a short distance (ca. 1 km) to the Brøggerbreen glacier. The unglaciated coastal areas are underlain with continuous permafrost extending to depths of about 100 m and with active layer depths of 1–2 m. The landscape is characterized by patterned permafrost ground and sparse vegetation coverage. The soils on the hill range from silty loam to silty clay, while there is coarser-grained material (silty sands to gravel) in the surrounding area. The site includes a 9 m deep borehole with hourly measurements of permafrost temperatures dating back to 2009, featuring a mean permafrost temperature of −2.8 ^∘C. In the active layer and in the top of permafrost the annual mean temperatures increased by 1.8 ^∘C per decade for the period 1998–2017 (Boike et al., 2018). The soil and air temperature trends show the largest temperature increase in winter (Maturilli et al., 2015). At the Bayelva site, the snow cover buildup starts in September–October and maximum snow depths can reach up to 1.5 m. The timing of spring snowmelt can vary by several weeks and typically starts between May and June. Data from the Bayelva site have been used for validation of various Earth system modeling studies (Chadburn et al., 2017; Ekici et al., 2014, 2015).

We prepare a common forcing time series for all simulations by downscaling surface fields from the ERA-5 reanalysis (Hersbach et al., 2018; closest model grid cell) using the Bayelva measurements as a reference. For each day of year, a linear regression is computed between reanalysis data and available measurements within a window of 20 d before and after each day of the year. This linear regression is subsequently employed to correct the reanalysis data. This procedure is in particular able to account for a potential seasonality of the bias of the reanalysis data, but relies on a sufficient correlation between the two time series (Westermann et al., 2015, 2016). The downscaling is applied to air temperature and incoming longwave radiation for which the reanalysis is strongly cold-biased during the summer period. The specific humidity is corrected accordingly for the change in temperature. Incoming shortwave radiation and wind speed are used without downscaling, as their average values generally agree well with measurements, while the correlation was lower due to their strong dependence on the timing of synoptic events, which is poorly captured by the reanalysis. Furthermore, precipitation is not corrected as there are no reliable measurements, especially for snowfall. Instead, we perform simulations for different multiplication factors for snowfall to evaluate the sensitivity to this important variable. All simulations are performed for the period 1981 to 2018 using the accelerated spin-up procedure (Sect. 3.1.4) for the first 10 years to initialize the ground temperature profile.

In the following, we first establish a model baseline configuration that can reproduce measurements reasonably well (Sect. 3.2.1) and then change one model component at a time, while the others are kept as in the reference (Sect. 3.2.2). We emphasize that only the model parameters of the reference simulation are optimized to match the measurements, while better-fitting configurations likely exist for the other simulations.

3.2.1 Reference simulations

We use the model setup of Zweigel et al. (2021) for the wider area around the Bayelva station as a starting point to optimize model settings specific for the location of the measurement station. A stratigraphy class with surface energy balance (Sect. 2.2.2), the soil freezing characteristic (Sect. 2.2.3) for soil type silt and bucket water balance (Sect. 2.2.4) is used to represent sediments in the uppermost 5 m. Bedrock is assumed below, represented by a less process-rich stratigraphy class with the free water freezing characteristic (Sect. 2.2.3) and constant sum of soil water and ice contents (Sect. 2.2.4). As the Bayelva station is located on a small hill, lateral drainage with the lateral interaction class representing a seepage face (Sect. 2.3.2) is assumed. The seasonal snow is represented by the Crocus-based snow class (described under (c) in Sect. 2.2.6), with parameters for wind compaction, initial snow density, and shortwave albedo adapted to fit measured ground surface temperatures, snow depths, and pre-melt snow densities of 300–350 kg m⁻³ which is a rough average for measurements conducted in 2011 to 2015. Figure 10 shows simulated annual average temperatures at 1.3 m depth for snowfall multiplication factors of 0.9, 1.0, and 1.1. For the majority of the years, the measurements fall within this range, which indicates that at least the general magnitude of snowfall is captured well by the model forcing. Figures 11 and 12 display the full time series of ground temperatures near the ground surface and at 1.3 m depth for the measurements and the simulations with a snowfall multiplication factor of 1.1. The simulations can describe the seasonality of measured temperatures well for most of the years, but winter temperatures are considerably too warm for some years, e.g., 2000–2001, 2002–2003, and 2009–2010. In these years, however, the measured snow depths at the Bayelva station are lower than in the simulations for most of the snow season (Fig. 13), which is the likely reason for the overestimated ground temperatures. In general, simulated snow depths increase smoothly throughout the winter season in the simulations, while the measurements bear evidence of individual snowfall events which lead to strong increases in snow depth.

Figure 10 Simulated (lines) and measured (crosses) annual average temperatures at 1.3 m depth for the reference simulation (see text). Years with measurements available for less than 90 % of the time are not shown. Black line: 100 % snowfall; broken blue line: 90 % snowfall; solid blue line: 110 % snowfall (used as a comparison to other model configurations, Sect. 3.2.2).

Figure 11 Simulated (black) and measured (red) ground temperatures at 4 cm depth for the Bayelva soil and climate station, Ny-Ålesund, Svalbard (data from Boike et al., 2018). Reference simulations with 110 % snowfall; see text. The model setup uses the surface energy balance as upper boundary, a subsurface module with soil freezing characteristic and lateral drainage with seepage face, and the Crocus-based snow model with optimized parameters.

Figure 12 Simulated (black) vs. measured (red) ground temperatures at 133 cm depth for the Bayelva soil and climate station, Ny-Ålesund, Svalbard (data from Boike et al., 2018). Reference simulations with 110 % snowfall; see text. The model setup uses the surface energy balance as upper boundary, a subsurface module with soil freezing characteristic and lateral drainage with seepage face, and the Crocus-based snow model with optimized parameters.

Figure 13 Simulated (black) and measured (red) snow depths at the Bayelva soil and climate station, Ny-Ålesund, Svalbard (data from Boike et al., 2018). Reference simulations with 110 % snowfall. See text.

3.2.2 Simulations with different CryoGrid stratigraphies

We use the warmest reference simulation with 110 % snowfall as a baseline for which simulated ground temperatures are close to the thaw threshold in some years. In all model configurations displayed in this section, snowfall is hence increased to 110 % as well. We then exchange one model aspect at a time, leaving all others unchanged so that differences in simulation results can be attributed to the model aspect in question. This is either accomplished by selecting a different stratigraphy class or by changing model parameters. It is important to note that a poor performance in reproducing measurements cannot necessarily be attributed to the investigated model aspect, as only the reference simulation is optimized for all model parts together. We also emphasize that the results of the model comparison are specific to the particular conditions at the Bayelva site and cannot necessarily be generalized. Therefore, the main purpose of this comparison is to showcase the effect of various setups within the modular CryoGrid community model and not to evaluate their performance at the Bayelva site in a strict sense.

Subsurface representation and water balance. Figure 14 displays simulated ground temperatures for different soil freezing characteristics, from silt for the reference simulation to sand and finally free water freezing for which the phase change of water occurs at 0 ^∘C. The differences in simulated ground temperatures are relatively small, which suggests that the choice of the soil freezing characteristic is not critical. At least in some situations, it is therefore possible to speed up simulations by using free water freezing instead of the computationally more demanding simulations with the soil water freezing characteristic.

Figure 14 Simulated annual average temperatures at 1.3 m depth. Blue, crosses: reference simulation; green (line only): as the reference simulation, but with the soil freezing characteristic for sand instead of silt; red, triangles: soil water freezing as free water.

Secondly, we investigate the effect of different lateral drainage conditions when using the bucket water scheme, with the reference simulation representing well-drained ground conditions due to the applied seepage face (Fig. 15), while excess water pooling up at the surface is removed instantly. The same configuration is also simulated with the excess ice class (Sect. 2.2.5), which allows water to pool up at the surface (and potentially refreeze), being removed laterally by overland flow (Sect. 2.3.2). Note that no excess ground ice is added to the ground itself, so model physics and parameters are identical to the reference simulation, except for the representation of surface water. While the differences to the reference simulation are generally small, the simulated ground temperatures are both warmer and colder in individual years. The main reason for the differences is the representations of surface ice layers, which regularly occur in fall and early winter in W Svalbard. If the ground is frozen but snow-free during a rain event (and saturated with ice and/or water), the rainwater is instantly removed in the reference simulation, while it can freeze and form a surface ice layer in the excess ice class. This not only changes the freeze–thaw dynamics of the uppermost ground layer, but also delays the onset of thawing in spring, as the melting surface ice layer confines the temperature of the uppermost grid cell to 0 ^∘C. The impact on the ground thermal regime of these two effects strongly depends on the conditions of the individual year, which explains the both negative and positive differences in the simulated mean annual ground temperatures. The third model configuration is once again based on the exact setup of the reference simulation, but without lateral inflow or outflow of water, representing a classic, one-dimensional simulation without interactions with external reservoirs. In the fourth configuration (denoted “wet”), inflow of water from a reservoir at the height of the ground surface ensures permanently wet ground conditions, mimicking a wetland located in a depression. While there is little difference between the one-dimensional and the wet simulation case, they are both significantly warmer than the well-drained reference simulation (Fig. 15), and a partly unfrozen zone above the permafrost develops in some of the years after 2006, corresponding to permafrost degradation. This can be explained by the higher soil water and thus latent heat content of the soil, which prevents complete freeze-back before an insulating winter snow cover forms. The small differences between the 1D and the wet simulation case are explained by the generally high summer precipitation at the Bayelva site, which often exceeds evapotranspiration, thus keeping the ground wet. Furthermore, in the bucket scheme employed for this comparison, water is not accessible to evapotranspiration below a certain depth, which is determined by the parameters evaporation and transpiration depth in the model. While this generally does not occur for a shallow active layer, a permanently saturated layer forms for the deep active layer in the Bayelva simulations, which requires a long time to freeze in fall and winter. When lateral drainage is applied, as in the reference simulations, the soil water content is held at field capacity instead so that the freeze-back occurs faster, thus also leading to colder ground temperatures and stable permafrost after 2006.

Figure 15 Simulated annual average temperatures at 1.3 m depth. Blue, solid line, crosses: reference simulation for drained soil conditions; blue, dashed line (line only): as the reference simulation, but representing surface water and ice (with excess ice module, Sect. 2.2.5) and lateral overland flow (Sect. 2.3.2); green (line only): no drainage; red, triangles: wet soil with influx of water.

The simulations not only deliver ground temperatures like the CryoGrid 2 and 3 models, but it is also possible to evaluate the components of the water balance, consisting of precipitation P, evapotranspiration ET, runoff (both from surface and subsurface) R, and storage S in the form of snow and soil water–ice, i.e., $P=\text{ET}+R+S$. In particular, CryoGrid stratigraphies which allow for standing surface water and removal by overland flow are well suited for this task, as all lateral runoff is handled in a physically based way (and not simply by removing excess water instantly). Figure 16 displays the accumulated water balance for 3 example years for this model configuration (dashed blue line in Fig. 15), showing the seasonality of the storage term corresponding to the winter snow cover, which is largely converted to runoff during and after snowmelt. At the same time, evapotranspiration is largest in July, corresponding to snow-free surface conditions and large amounts of incoming radiation. Smaller amounts of water equivalent are lost to sublimation during the snow-covered season, which is clearly visible during the first winter of 2009–2010. Note that there are also smaller changes in storage during the snow-free season, which corresponds to changes in soil water storage.

Figure 16 Simulations of accumulated water balance for simulations with surface water representation (i.e., excess ice class), lateral drainage by a seepage face and overland flow, and 110 % snowfall (dashed blue in Fig. 15) for 3 example years. P (thick black line): precipitation (sum of rainfall and snowfall); ET: evapotranspiration; R: runoff (sum of surface and subsurface runoff); S: storage (soil and snow).

We further investigate the influence of the soil water balance on simulation results by considering different model representations (Fig. 17). The reference simulation with the bucket water scheme is first compared to a simulation with constant soil water content set to the field capacity of the reference simulation, with evapotranspiration calculated independently of the soil water content by setting a constant model parameter for the “surface resistance against evapotranspiration” (as in the scheme used for Westermann et al., 2016). The simulated soil temperatures are very similar to the reference simulation, which reflects the overall similar water contents within the active layer. In addition, the reference simulation is compared to the hydrological scheme driven by the Richards equation, including redistribution of soil water during the freezing process. While this results in slightly higher ground temperatures, they are still reasonably similar to the reference simulation (Fig. 17). Figure 18 displays the simulated volumetric water contents for the three different model representations of the soil water balance, showing significant differences during the unfrozen season. While the soil water content does not change in the “constant water content” case (middle panel), approximately the upper half-meter (controlled by the user-defined characteristic depths of evaporation and transpiration, Sect. 2.2.4) dries out as a result of evapotranspiration for the bucket water scheme (left panel), while strong rainfall events periodically reset the water content to the field capacity. At the same time, water accumulates near the bottom of the active layer, which is removed laterally by the applied seepage face module. Finally, the entire active layer remains unsaturated in the simulations with the Richards equation, as water is redistributed upwards to compensate for evaporation (right panel) during most of the summer season. For this reason, no lateral subsurface runoff occurs (see Sect. 2.3.2), although the soil water content is overall higher than in the reference simulation, especially in fall. The simulated thaw depths are more or less equal for the bucket water and the Richards equation schemes, while they are slightly deeper for the scheme with constant soil water contents. We point out that the simulations with the Richards equation yield a different assessment of permafrost stability in the wider Bayelva area, despite the relatively similar ground temperatures: while permafrost in the reference simulation is only stable when lateral drainage occurs (Fig. 15), it is also stable for the one-dimensional case if the soil water balance is governed by the Richards equation.

Figure 17 Simulated annual average temperatures at 1.3 m depth. Blue, crosses: reference simulation (bucket scheme for water, including lateral subsurface drainage); green (line only): constant soil water plus ice content; red, triangles: Richards equation (no lateral subsurface drainage occurs).

Figure 18 Simulated profiles of volumetric water contents for 2 example years for the three water balance schemes used in Fig. 17. (a) Reference simulation (bucket scheme for water, including lateral drainage); (b) constant soil water plus ice content; (c) Richards equation (no lateral subsurface drainage occurs).

An example of the effect of cryosuction in the Richards equation scheme is displayed in Fig. 19, which shows the sum of the volumetric water and ice contents, θ_wi, for 1 example year. The accumulation of ground ice below and above the freeze front is clearly visible in this year; this in particular leads to an ice-rich layer at the bottom of the active layer, which can modify thaw behavior in the following years (Schuh et al., 2017). Furthermore, cryosuction can replenish ground ice that has melted during summers with deep active layers and is lost to the surface due to evaporation (Fig. 18, right panel). The previous maximum of the active layer thickness is clearly visible in Fig. 19, and ice accumulation over several years with a shallower active layer could replenish the lost ground ice at least partly. In the Bayelva simulations, this does not really happen as the active layer is overall deepening during the simulation period so that ground ice formed at the bottom of the active layer due to cryosuction eventually melts out.

Figure 19 Simulated sum of volumetric water and ice contents using the Richards equation for a 1-year period in 2000–2001 (red in Fig. 17); white line: freeze front. Cryosuction upon soil freezing leads to ice accumulation above and below the freeze front, as marked by the two orange arrows. In this example period, the effect of cryosuction is particularly pronounced, while it is of smaller magnitude in other periods.

Snow cover representation. Figure 20 displays the simulation results for different snow cover classes, in particular the class with snow microphysics representation (following the Crocus-based snow model) and the snow class with a constant snow density (as in the classic CryoGrid 2 and 3 models). Both snow classes feature a bucket scheme for water infiltration, including refreezing of meltwater percolating through the snowpack. While the parameters in the Crocus-based snow class are specifically adapted for the reference simulation, two of the parameter settings presented in the literature produce ground temperatures that are too high and permafrost that is not thermally stable. These are the original Crocus parameter set (Vionnet et al., 2012), which produces snow densities that are too low and thus too much insulation of the ground, and the “Arctic” setting (Royer et al., 2021). This produces reasonable snow densities, but the employed parameterization for the snow thermal conductivity (after Sturm et al., 1997) yields significantly lower values than in the original Crocus configuration (using the parameterization after Yen, 1981). As a result, the simulated temperatures are relatively similar for the two settings, despite significant differences in the simulated snow densities and thus snow depths. Note that the reference simulation uses the parameter set of the “Arctic” configuration, but the Yen parameterization for the snow thermal conductivity. We point out that snow accumulation at the Bayelva site is subject to strong wind redistribution of snow, which likely also changes over time, as the measurement site is located on a hill and thus exposed to snow ablation, but is also surrounded by a protective fence providing wind resistance and thus promoting deposition. The results must therefore be regarded with caution and should not be generalized unless backed up by simulations from other sites.

The simulations with constant snow densities are generally colder than in the reference, although overall low snow densities of 250 and 275 kg m⁻³ were used (compared to end-of-season measurements of 300–350 kg m⁻³), which represent a seasonal and depth average for the entire snowpack. In fall and early winter this likely overestimates the snow density and thus thermal conductivity of the fresh snowpack so that the ground can refreeze faster. In the annual balance, this offsets the overly strong insulation in the later parts of the winter, when the true snow densities are higher in the majority of the snowpack.

Figure 20 Simulated annual average temperatures at 1.3 m depth. Blue, crosses: reference simulation (Crocus-based snow scheme with adapted parameters); green solid (line only): constant density snow scheme with initial density of 250 kg m⁻³; green dashed (line only): constant density snow scheme with initial density of 275 kg m⁻³; red solid, triangles: Crocus-based snow scheme with original parameters from Vionnet et al. (2012); red dotted (line only): Crocus-based snow scheme with the “Arctic” configuration (Royer et al., 2021).

Model type and upper boundary condition. In the reference simulation, energy exchange at the upper boundary is driven by the surface energy balance, similar to the original CryoGrid 3 model. While most stratigraphy classes in the CryoGrid community model employ the surface energy balance, simpler approaches, such as a temperature boundary condition (as in CryoGrid 2) and the semi-empirical TTOP model (as in CryoGrid 1), are also available. Figure 21 shows a comparison of different model types, with a CryoGrid 2-like model setup with air temperature as the upper boundary condition for both the ground and the snow surface, as well as snow ablation derived from a degree-day melt model (Sect. 2.2.6). While the model does not account for the water balance and is generally a much less process-rich approach compared to the reference simulation, it delivers a similar representation of annual average temperatures for the Bayelva site. Interestingly, even the simple TTOP approach can capture some of the warming trend observed in the last 2 decades, although the n-factors are not adapted to represent interannual differences of snow depths. However, as a stand-alone, this equilibrium approach should strictly speaking only be applied to longer time periods (dashed green line in Fig. 21), as it cannot account for the time lag of ground temperature changes due to heat conduction and possibly ice melt. This problem can be partly overcome by using the TTOP model as the upper boundary condition for a stratigraphy class simulating heat conduction (dotted green line in Fig. 21). In this example application, the TTOP model computes mean annual ground temperatures for annual time slices, which is used as temperature boundary condition for the heat conduction model (Sect. 2.2.2). While the active layer dynamics is not resolved, this hybrid approach can efficiently simulate the temperature dynamics in deeper layers (similar to Myhra et al., 2017), which is particularly interesting to speed up centennial- and millennial-scale simulations. Although the TTOP model is a simplistic approach associated with significant uncertainty, it is straightforward to train and validate with more process-rich model setups in the CryoGrid community model.

Figure 21 Simulated annual average temperatures at 1.3 m depth. Blue, crosses: reference simulation (upper boundary condition surface energy balance); red, triangles: temperature boundary condition model (air temperature as upper boundary condition, degree-day-based snowmelt model); green (lines only): TTOP model with time-constant n-factors n_f=0.5, n_t=1, r_k=0.8 (dashed: entire period; solid: 2-year intervals; dotted: 1-year interval as the upper boundary condition for CryoGrid stratigraphy as in the reference simulation).

Glaciers and water bodies. The CryoGrid community model provides stratigraphy classes for land surface types other than ground and snow, in particular glaciers and water bodies. Figure 22 displays simulated ground temperatures for a glacier and a shallow (0.8 m deep) water body when driven with the same forcing data as the reference simulation. The water body simulations clearly show that average ground temperatures are positive throughout the entire time, which is due to the asymmetry of heat transfer between summer and winter (Westermann et al., 2016). This suggests that there is a permafrost-free zone below water bodies of that depth around Ny-Ålesund, with the water body not freezing to the bottom in winter. Figure 23 shows the simulated lake ice dynamics for 2 example years, with significantly thicker ice cover during a year with low overall low snow depths. However, even in this year, the water body does not freeze to the bottom, which explains the high ground temperatures and permafrost-free conditions. Conversely, annual average temperatures are significantly lower for the glacier than for the reference simulation, as the energy input during summer is consumed by ice melt, rather than warming of the subsurface material. Furthermore, the summer meltwater fully drains so that freezing temperatures in fall and early winter instantly penetrate the subsurface, other than for soil freezing whereby the water in the active layer must first refreeze, before the ground below can cool out further. For the employed stratigraphy class, the annual ice melt is automatically compensated for by mass from an infinite ice layer below the model domain. Figure 24 shows the annual ice melt for a theoretical glacier located at the Bayelva site. The ice melt follows the observed mass balance from nearby glaciers Austre Brøggerbreen and Midtre Lovènbreen well, which have an increasingly negative mass balance trend since the 1980s with the highest mass loss year in 2016.

Figure 22 Simulated annual average temperatures at 1.3 m depth. Blue, crosses: reference simulation for soil material; red, triangles: 0.8 m deep water body overlaying soil material (as for the reference simulation); green (line only): glacier.

Figure 23 Simulated ice dynamics (volumetric ice content) for the water body simulation (red in Fig. 22) for a winter with low and high snow depths. As snowmelt water is assumed to pass through the water body ice in the model, the ice rises in elevation during the melt period, which is particularly evident in 2008–2009.

Figure 24 Simulated annual net ice melt (in meters of water equivalent) with the glacier module for the Bayelva forcing dataset.

Figure 25 displays a comparison of the runtimes of the model configurations presented in this section, which are between 150 and 300 s per model year for most of the cases. The largest differences in runtime are caused by the boundary condition and/or model type, with simulations of the surface energy balance requiring a considerably longer runtime compared to simulations in which a temperature boundary condition is applied. The shortest runtimes are achieved by the even simpler TTOP model configurations which do not simulate heat conduction at all and only compute a single equilibrium temperature. Among the model configurations using the surface energy balance as the upper boundary condition, schemes with simpler model physics (e.g., constant water plus ice contents, constant snow density) have shorter runtimes than the more process-rich model configurations, but not by a large margin.

Figure 25 Comparison of runtimes for the different CryoGrid configurations presented in Sect. 3.2, all using a single core of an AMD EPYC 7742 64-core processor (2.25 GHz base frequency). Bold font: setting for the reference run; the color coding and line properties are as in the individual figures. Note that the runtime includes the computational overhead for model initialization and preprocessing of model forcing, which in particular for “TTOP 2 years” and “TTOP full period” dominate the overall runtime.

4 Discussion and outlook

The modular concept of the CryoGrid community model was designed to unlock added value of model development conducted in several research projects with a strongly different research focus, such as Nunataryuk, PermaRisk, and ESA Permafrost_CCI. Instead of building specialized model tools for a narrow set of applications, the community model aims for a library of different model parts that can be flexibly combined by users, while providing developers with an advanced starting point to design new model functionalities. As such, the goal of the CryoGrid community model is to become a dynamic platform for continued model development, not a static model applied by its user community to solve the same classes of problems.

4.1 Differences to CryoGrid 1–3 model implementations

CryoGrid 1. The same TTOP equilibrium model that forms the core of CryoGrid 1 is implemented by a dedicated stratigraphy class in the CryoGrid community model. However, most CryoGrid 1 applications relied on additional parameterizations for the model parameters. In particular, the winter n-factor (Eq. 3) has been calculated from snow depth (e.g., Gisnås et al., 2013), snowfall (Westermann et al., 2015), or snow depth and winter temperature (Obu et al., 2019). While it is straightforward to add these parameterizations to the CryoGrid community model, they have all been developed and adapted for specific input datasets. We therefore leave it to concrete future applications to complement the “raw” TTOP model currently implemented with additional parameterizations for its model parameters.

CryoGrid 2. The main features of CryoGrid 2 are a temperature boundary condition, a simple snow model with constant snow density, and a numerical scheme optimized for computational efficiency. The first two features are retained in the corresponding stratigraphy classes of the community model, but the numerical scheme has been changed to the simple first-order forward Euler scheme with adaptive time step, which is not as computationally efficient as previous CryoGrid 2 implementations. In the future, computationally more efficient implementations of CryoGrid 2 could be added by means of further TILE classes which support different numerical schemes (see also Sect. 4.7).

CryoGrid 3. While most stratigraphy classes in the community model follow the general design of CryoGrid 3, there are some notable differences: in CryoGrid 3, the main prognostic variable integrated in time is the intensive variable temperature, while extensive variables, such as enthalpy, are employed in the community model. This has several advantageous aspects, in particular that extensive quantities are intrinsically conserved when integrated in time. For the same reason, it is straightforward to combine and split grid cells when extensive variables are used as state variables. Using enthalpy as state variable is particularly convenient as both the apparent heat capacity and the thermal conductivity behave smoothly as functions of enthalpy, while the apparent heat capacity has a significant discontinuity at the freezing point of water when expressed as a function of temperature. For numerical integration, this is a challenging aspect, so using enthalpy as state variable is especially favorable for the simple time integration scheme employed in the community model. Another difference between CryoGrid 3 and the stratigraphy classes of the community model is that soil water is treated as a prognostic variable (i.e., time derivatives are calculated to advance the variable in time). In CryoGrid 3, an instantaneous bucket scheme is employed instead, which computes the steady-state soil water profile after each time step, i.e., the water profile that would eventually be reached for the particular amounts of water added or extracted during the time step (which means that soil water is treated as a diagnostic variable). The time dynamics of an infiltration event are therefore only captured in the community model. We emphasize that this is not a hard-coded feature of the community model, but only applies to the stratigraphy classes described in this study. A final difference between the two models is the use of a different parameterization for the soil freezing characteristic in the community model (see Sect. 2.2.3). For most thermal simulations, we expect the impact of this change to be small (see Sect. 3.2.2 on the effect of different freezing characteristics), but the parameterization employed in the community model in particular offers a better performance for combined thermal and hydrological simulations with the Richards equation. The interested reader is referred to the discussion in Painter and Karra (2014) on this topic.

4.2 CryoGrid 1–3 functionalities lacking in the community model

The CryoGrid community model provides the majority of the simulation capabilities that have been demonstrated within the CryoGrid 1 to 3 models, while at the same time providing modularity and a joint operation framework that makes additional functionality, such as the accelerated spin-up (Sect. 3.1.4), available to all model configurations. However, there are a few simulation tools published within the CryoGrid 1–3 models that are at this point not available within the CryoGrid community model. Firstly, this concerns the two-dimensional heat conduction model CryoGrid2D (Myhra et al., 2017) for which a temperature boundary condition must be specified at a number of nodes at the edge of the model domain, which is not possible with the currently available FORCING classes. Moreover, it is not directly compatible with the inherently one-dimensional stratigraphy classes described in this work. Nevertheless, the general framework of the CryoGrid community model in its present form contains structures by which CryoGrid2D could be incorporated. The most direct way is to implement dedicated TILE, FORCING, and OUT classes (see Fig. 1) that are geared to the requirements of two-dimensional simulations. While not directly compatible with the simulation tools presented in this work, functionalities enabled at the RUN_INFO level (Fig. 1), e.g., the accelerated spin-up procedure, could directly be employed.

A second functionality lacking in the community model is the FLake water body model (Mironov et al., 2005), which has been demonstrated within CryoGrid 3 (Langer et al., 2016). Most notably, FLake computes a wind-speed- and temperature-dependent stratification of the unfrozen water column, while the currently implemented stratigraphy class assumes completely mixed conditions (i.e., the entire water body has the same temperature) at all times, similar to the water body representation in Westermann et al. (2016). It is clear that the simple representation is only adequate for shallow water bodies, while FLake must be implemented to simulate a realistic stratification of the water column in deeper water bodies. While the simple water body representation appears to be a sufficient approximation for many permafrost applications with shallow ponds and lakes (Westermann et al., 2016), it is exactly the purpose of the modular community model that such assumptions, though seemingly logical, can be verified by comparing to simulations with a more process-rich scheme. Therefore, implementing a stratigraphy class based on FLake in the community model should become a priority for future model development.

Finally, freezing of saline ground, including diffusion of salt driven by gradients in brine concentrations (Angelopoulos et al., 2019), is a model functionality demonstrated within CryoGrid 2 but not yet included in the community model. Being a one-dimensional setup, it can be implemented in the community model as a stratigraphy class in a straightforward way. Coupling to existing stratigraphy classes can be facilitated by dedicated interaction classes that compute the fluxes of energy, water, and salt between the class interfaces. We point out that freezing of saline ground in CryoGrid 2 so far has only been demonstrated for saturated conditions and without salt advection due to water flow. Within the community model, these limitations can likely be overcome, as a stratigraphy class for saline ground can be built on top of existing stratigraphy classes with water balance representations. In a similar fashion, the existing water body class could be supplemented with a salinity representation so that freezing of saline water bodies can be simulated. Due to the possibility to vertically connect classes, simulations of subaquatic permafrost below saline lagoons (e.g., Angelopoulos et al., 2020) and shallow nearshore zones could be accomplished within a single model framework which simulates ice formation and brine exclusion in the water column in conjunction with salt diffusion and ground freezing in the ground column below.

4.3 Lateral interactions

The CryoGrid community model offers a standardized way to implement lateral interactions of a one-dimensional model column with its environment, which makes it possible to account for the influence of topography and terrain heterogeneity at least in a phenomenological way. This is first and foremost useful to drive the water balance in point-scale simulations to reality. Many wetlands, for example, are located in terrain depressions and receive inflow from surrounding areas, which is the very reason for their existence. Model schemes with one-dimensional water balance cannot reproduce soil moisture conditions in such ecosystems, which can be a significant problem, for example, in carbon cycle simulations. In the CryoGrid community model, lateral inflow can be realized through a dedicated lateral interaction class which connects the model column to a water reservoir at a defined elevation. Likewise, continuous drainage of locations on slopes and hills can be simulated by using lateral interaction classes, which simulate drainage through either a seepage face or a lower-lying water reservoir. The Bayelva simulations for different drainage regimes (Fig. 15) clearly demonstrate the large impact of the drainage regime on ground temperatures, with permafrost thermally stable when water can drain, but degrading when the ground is held wet through inflow of water. In addition to the water balance, lateral fluxes of heat may play a role in permafrost areas, for example by thermally stabilizing permafrost below water bodies (Langer et al., 2016) or destabilizing permafrost in the vicinity of infrastructure elements (Schneider von Deimling et al., 2021). To simulate such effects, a lateral interaction class can be applied in the community model to connect the CryoGrid stratigraphy to an external heat reservoir, following the same standardized way of operation as for lateral water fluxes.

A major asset of the CryoGrid 3 model is the possibility to simulate lateral interactions between different 1D model columns (denoted “tiles”, corresponding to a TILE class in the community model), which facilitated simulations of complicated landscape evolution in ice-rich permafrost areas (Martin et al., 2021; Nitzbon et al., 2019, 2020, 2021). As of now, the CryoGrid community model offers the same possibilities, which will be described in detail in a future study. In short, three-dimensional lateral interaction classes are used, which rely on exactly the same structures as the one-dimensional interaction classes presented in this study. Users can therefore first conduct one-dimensional simulations coupled to a static environment, e.g., a water reservoir, which can reveal many of the effects of lateral fluxes in a semi-quantitative way. Such simulations subsequently offer an excellent starting point to set up more complicated, three-dimensional simulations with laterally coupled tiles, simply by adding information on the type of interaction and the relative positions of the different tiles to the parameter files, while the properties and model parameters related to the one-dimensional CryoGrid stratigraphies, such as ground properties and forcing data, can remain unchanged.

4.4 Model structure and computational aspects

The CryoGrid community model generally follows the concept of “loose coupling”, keeping the degree of knowledge that each class has of the others at a minimum. Lateral interaction classes, for example, are not specifically assembled for a certain stratigraphy class but work for a number of different stratigraphy classes. Therefore, the lateral interaction classes do not contain any code representing specific knowledge of the stratigraphy classes that they can be used with. Instead, each stratigraphy class first processes its internal information to a defined common interface format that is then passed to and processed by the lateral interaction class. This concept is applied throughout the entire model framework; while new code needs to comply with the established interfaces, it can immediately exploit much of the existing functionality. Furthermore, code additions can be made without modifying the source code of the existing classes, which is important to facilitate parallel model development by several independent developers without interference. Thus, the CryoGrid community model provides an integrated platform that is designed to accelerate and simplify the development of new simulation tools for the terrestrial cryosphere. As an example, its modular setup provides a clear pathway for the development of new stratigraphy classes, with several defined steps that can be implemented, verified, and tested independently of one another.

1. Implement the stratigraphy class itself and test it with simulations of a CryoGrid stratigraphy consisting only of this one class.

2. Implement interaction classes between the new stratigraphy class and existing stratigraphy classes that the new class shall be compatible with; test with simulations of CryoGrid stratigraphies consisting of all supported class combinations.

3. Implement lateral interaction classes that the new stratigraphy class shall be compatible with; this can mean the implementation of an entirely new lateral interaction class or implementing the functionality to make existing lateral interaction classes compatible with the new stratigraphy class.

This procedure accelerates the development of complex model configurations, as developers in the first step can fully concentrate on the novel model physics, while the capacities of the new model functionality are amplified in the second and third steps by interfacing to already existing stratigraphy and lateral interaction classes.

In land surface modeling, such modular approaches are increasingly adopted to deal with the increasing process complexity (Fisher and Koven, 2020). As an example, the version 5.0 of the Community Land Model (CLM) offers the possibility to exchange several model components and parameterizations (Lawrence et al., 2019), in particular the soil hydrology scheme and the snow model. As CLM is designed as the land component of Earth system models, e.g., CESM2 (Danabasoglu et al., 2020) and NorESM (Seland et al., 2020), all model setups feature full land–atmosphere coupling, while “simple” schemes without surface energy balance are not provided. The CryoGrid community model, on the other hand, is more used as a stand-alone tool with broader application focus and therefore includes such simple options, in particular the TTOP equilibrium model and heat conduction models with temperature boundary condition, which are largely inherited from CryoGrid 1 and 2.

The time integration in the CryoGrid community model (see Sect. 2.2.9) does not rely on an established partial differential equation solver (as CryoGrid 2, Westermann et al., 2013), but rather employs a simple forward Euler solver with explicit time step control. A main reason for this is to keep the code more readable, with the goal to lower the barriers for further model development, also by non-expert users. Moreover, many established solver packages require the model grids to be static, and retaining full control over time integration makes it easier to include processes requiring constant changes to the model grid, such as excess ice melt and compaction of the seasonal snow cover. An obvious disadvantage of the simple time integration scheme is the relatively slow computation (see Fig. 25 for runtimes), which may be prohibitive for some applications. Another limiting factor for runtime is the use of MATLAB, which is an interpreter language leading to slower runtimes. In the future, it may be possible to use automatic code translators available within MATLAB (MATLAB coder) or provided by independent developers (e.g., Paulsen et al., 2016) to produce a faster C++ version of CryoGrid, in particular since C++ objects feature many similarities with the MATLAB classes employed in the community model.

4.5 Multi-model simulations

While evaluations of multi-model ensembles are increasingly common in cryospheric studies (e.g., van Pelt et al., 2021), most applications assume a fixed model structure. In this typical case, the evaluation of the uncertainty of simulation results is hampered by the fact that the structural uncertainty (i.e., due to imperfect model physics) tends to be lumped together with the parametric uncertainty (i.e., due to poorly constrained model parameters). As demonstrated in Sect. 3.2, a number of different model schemes can be efficiently simulated in the CryoGrid community model, which allows for a more thorough qualitative and quantitative assessment of model uncertainty. In this study, we have restricted ourselves to a comparison of ground temperatures simulated with different model configurations, but we have not attempted to identify which model configurations best explain the field measurements. In the following sections, we discuss techniques to evaluate both parametric and structural model uncertainty from observations, which could be implemented in the CryoGrid community model in the future.

To quantify parametric uncertainty in a model with fixed structure, one typically assumes that the impact of unknown model parameters, including the uncertainties due to model forcing and initial conditions, can be represented by parametric error models. This is commonly achieved probabilistically through the framework of Bayesian inference (Berliner, 2003). The goal then becomes to estimate the “posterior” distribution of the parameters, which quantifies parametric uncertainty in the model after accounting for the likelihood of the data under various plausible parameter settings (Wikle and Berliner, 2007). Since analytical forms of the posterior are not usually available, it is typically calculated approximately using numerical sampling methods (Berliner, 2003; Gelman et al., 1995). In the setting of geoscientific modeling, wherein each forward evaluation of the model is often computationally demanding, approximate inference is commonly performed using data assimilation (DA) methods that efficiently fuse uncertain information from physical models with observations (Carrassi et al., 2018). These methods allow physical models to be fitted to observations (corresponding to a model calibration), while evaluating the uncertainty through a probabilistic inversion that corresponds to the first level of inference in the Bayesian hierarchy (MacKay, 2003). For cryospheric applications, ensemble-based data assimilation schemes based on the particle and ensemble Kalman filter (Carrassi et al., 2018) have been shown to be especially promising (e.g., Aalstad et al., 2018; Fiddes et al., 2019; Alonso-González et al., 2021).

One potential issue with this kind of model fitting is that all uncertainty is attributed to the model parameters, while one effectively assumes the model structure itself to be perfect. Since all models are ”wrong” to some degree (Box, 1976), this assumption is never completely justified, potentially resulting in biases in the fitted model parameters that compensate for the imperfect model. An example of this dilemma is the Bayelva simulations for different water balance schemes and drainage conditions (Figs. 15, 17): if only the bucket scheme with lateral drainage was available, one could, for example, adjust the strength of the lateral drainage to fit the temperature observations at the Bayelva station and conclude that permafrost is only thermally stable for well-drained conditions (Fig. 15). Likewise, one would infer from this model configuration that permafrost may have already started to degrade in areas with poor drainage for which no observations are available (the Bayelva station is located on a hill where good drainage can indeed be expected). However, the simulations with the Richards equation suggest that permafrost could also be stable if no lateral drainage occurs (Fig. 17), which provides a qualitatively different assessment of permafrost stability in the wider area around the Bayelva station. In principle, these ambiguities can be resolved by considering additional observations, for example soil water contents (Boike et al., 2018). In this study, we have not explored and resolved this issue in further detail, but the example from the Bayelva area showcases a problem inherent in many cryospheric simulations.

Running multiple models with different representations of physical processes alongside one another can help to identify potential ambiguities arising from this structural uncertainty. With the CryoGrid community model, adopting such a multi-physics ensemble of models becomes relatively straightforward in practice due to the common interface (e.g., for the model forcing). Furthermore, differences in simulation results can be attributed to a specific stratigraphy class or model parameterization, as it is possible to exchange stratigraphy classes one by one and keep the rest of the stratigraphy unchanged (Sect. 3.2). In the future, the CryoGrid community model could also be used to compare different models objectively and quantitatively, which corresponds to the second level of inference in the Bayesian hierarchy (MacKay, 2003). Here, the key quantity is the “model evidence” (also called marginal likelihood), which is the probability of the observations, given a particular model. The model evidence serves as a normalizing constant when fitting a model at the first level of inference. At the second level of inference, the model evidence is used as a likelihood in Bayes' rule to update the prior probabilities of each model so as to effectively constrain model structural uncertainty using the available data. As such, the model evidence framework can be used for both Bayesian model selection and averaging. This framework is widely used in physics (Feroz et al., 2009) and machine learning (MacKay, 2003), and it is slowly gaining traction in geoscience (Schöniger et al., 2014). Recent work has shown that ensemble-based data assimilation schemes are also well suited for obtaining robust estimates of the model evidence (Carrassi et al., 2017). As such, a promising future avenue in CryoGrid will be the implementation of ensemble-based data assimilation for both cryospheric model fitting (e.g., Aalstad et al., 2018; Alonso-González et al., 2021; Fiddes et al., 2019) and model comparison (Carrassi et al., 2017).

The modular structure of the CryoGrid community model offers the possibility to estimate effective parameters in simplified model schemes by assimilating results from high-fidelity models with richer process representations. This vision is in line with the new Earth system modeling blueprint presented by Schneider et al. (2017), which proposes using ensemble-based data assimilation to let coarse-scale models learn from both observations and targeted high-resolution models. For example, transient heat-conduction-based models could be used to compile climate-dependent parameterizations for n-factors in the simple TTOP model, which can then be used to map large areas at a lower computational cost. In such a scheme, the simple model is essentially used to spatially extrapolate the output of the more sophisticated model in a quasi-physical fashion. Similar schemes could be devised to achieve long-term (e.g., multi-millennial) simulations, e.g., by parameterizing the net effect of the surface energy balance and the seasonal snow cover so as to force the heat conduction model by monthly or even annual averages (as in Overduin et al., 2019; Etzelmüller et al., 2020, 2022). The application of machine-learning-based emulators also opens up the possibility for full-fledged Bayesian inference of the (emulated) posterior via traditional sampling methods, which could provide a more robust and scalable framework for inverse problems (Cleary et al., 2021), thereby unifying machine learning and data assimilation under the Bayesian framework (Geer, 2021).

4.6 Possible applications of the CryoGrid community model

The CryoGrid community model is designed first and foremost for the same applications as the CryoGrid 1–3 models, which were mainly focused on the permafrost thermal regime. In addition, the community model offers the possibility to control the lateral drainage regime of the model domain (Sects. 2.3, 4.3), a capability essential for hydrological applications that is rather limited in conventional 1D models. In particular, this allows conducting simulations with a coupled energy and water balance for situations in which a one-dimensional representation, especially of the water balance, is not appropriate. On the one hand, this concerns permafrost in wetlands that receive net inflow of water from surrounding areas, which also impacts the ground thermal regime (red line Fig. 15). Although one-dimensional model configurations cannot fully capture the time dynamics of the inflow, scenario runs with different strengths of the inflow (e.g., by modifying the distance to the water reservoir, Sect. 2.3.2) could be conducted to evaluate the resilience of wetland permafrost in a changing climate in a semi-quantitative way.

On the other hand, net outflow of water typically occurs for permafrost in sloping terrain or on top of hills, which influences both the thermal regime (blue line Fig. 15) and the water balance (Fig. 16) of the ground. This provides new applications, especially in mountain environments where the complex topography with strong elevation gradients and high spatial variability of the surface energy balance impacts the ground thermo-hydrological regime and thus permafrost distribution (Haeberli et al., 2010; Haeberli, 2013). Mountain permafrost has gained increasing interest in recent years, with emphasis on permafrost degradation and its impact on rock-wall stability (e.g., Gruber and Haeberli, 2007), as well as hydrology (Gao et al., 2021; Yang et al., 2019). Previous studies have used CryoGrid 2 and 3 to simulate thermal processes in steep rock walls (Legay et al., 2021; Magnin et al., 2017; Schmidt et al., 2021), and such applications are also possible with the CryoGrid community model. In addition, thanks to its representations of the surface energy balance, snow cover dynamics, and lateral drainage, the CryoGrid community model is an efficient tool to capture the effect of local topographic controls on the ground thermal regime, making it applicable in most mountain permafrost settings.

These capacities can be combined with simulations of massive ice bodies, such as glaciers and ice patches, with a glacier class (Sect. 2.2.8) implemented in the CryoGrid community model. In mountain environments, this offers the possibility to simulate the complex climate-driven interactions between glaciers and permafrost, for example the development of permafrost upon glacier retreat. While not explicitly simulating the transition, the thermal simulations for Ny-Ålesund (Sect. 3.2.2) showcase the expected changes in the ground thermal regime. Despite the strong net ice melt (Fig. 24), the simulated mean annual temperature of the glacier is much colder than for the permafrost ground (Fig. 22), suggesting that glacier tongues at this location at least partly feature cold-based ice, so the area would directly transition to a permafrost setting upon glacier retreat. Combined with stratigraphy classes representing relevant ground materials (especially with excess ice, Sect. 2.2.7), the climate-driven evolution of features associated with glacier–permafrost transitions, such as debris-covered glaciers, buried glacial ice, ice-cored moraines, and ice-rich frozen debris bodies, could be simulated. The rapidly changing elements of the mountain cryosphere are in particular expected to have profound impacts on stream flow, water resources, and ecosystems, which not only affect the glacial and periglacial realms, but also adjacent valley and lowland areas. With the possibility to jointly evaluate ground thermal regime and water balance, the CryoGrid community model can thus provide simulation capacities for a variety of climate change impact studies in mountain environments.

4.7 Directions and priorities for future model development

In general, the future directions of development of the CryoGrid community model are decided by the users and developers in a bottom-up approach. Here, we identify a few priorities and ideas that are motivated by the results of this work (Sect. 3) and previous modeling efforts with CryoGrid 1–3. Firstly, the sensitivity experiments with different snow schemes (Sect. 3.2) show that the snow representations within CryoGrid need to be improved, as this is one of the most important sources of uncertainty in simulating the ground thermal regime (Fig. 20). With a starting point in the snow class based on the Crocus snow model (Vionnet et al., 2012), both simpler and more process-rich representations are conceivable. In the currently implemented Crocus-based snow class, the snow density, which is the critical control for heat flow through the snowpack, has a strong dependence on wind speed, which can be heavily biased in model forcing datasets, especially in complex terrain (e.g., Delhasse et al., 2020). Therefore, a snow stratigraphy class in which both initial snow density and the strength of wind compaction are controlled by time-invariant empirical parameters (and not by wind speed) could be useful for many applications, as they may provide an improved performance compared to snow schemes with constant snow density (Sect. 2.2.6, Fig. 20) if the quality of the model forcing is poor. The Crocus snow model is undergoing continuous development (e.g., Royer et al., 2021), and many improvements can be ported to CryoGrid. We see two areas in which the CryoGrid community model itself could be a useful tool to catalyze snow model improvements. First, the impact of different near-surface temperature regimes on snow properties could be investigated with relative ease; the drainage regime, for example, changes the dynamics of ground freezing during winter (Fig. 15), which in turn modifies the thermal gradients and thus snow microphysics in the snowpack. Furthermore, the snowpack on Arctic water bodies is remarkably different from surrounding land areas (e.g., Langer et al., 2016), which again can be simulated with relative ease in CryoGrid by selecting a stratigraphy class for water bodies. A comparison with observations of snowpack properties for areas with different subsurface properties can therefore reveal shortcomings of model parameterizations and eventually contribute to improvements. Secondly, explicit redistribution of snow between CryoGrid tiles due to wind drift has been demonstrated by Zweigel et al. (2021), which also has pronounced effects on the stratigraphy of snowpack properties. While lateral interactions between tiles will be described in a later study, ablation of snow from exposed locations due to wind drift could be implemented as a lateral interaction class (Sect. 4.3) through the scheme proposed by Zweigel et al. (2021). Such simulations could again be compared to observations of snowpack properties and to the dynamics of ablation events at wind-exposed locations, such as small hills and ridges, so that both model parameters and (if necessary) the parameterizations for wind drift can be modified.

A second priority for model development is including representations for vegetated surfaces, which is currently a key shortcoming in the CryoGrid community model. A multi-layer canopy model presented in Stuenzi et al. (2021a, b) is already available as a stratigraphy class, but this is designed to represent forests and not short vegetation, such as mosses, sedges, and shrubs. An intermediate-complexity representation of vegetation, e.g., following a big-leaf approach (e.g., Sellers et al., 1992), has significant potential to improve the surface energy and water balance in simulations for vegetated tundra areas. Further development could focus on the interplay between the snow cover and different vegetation types, which has the potential to improve simulations of the ground thermal regime in areas experiencing shrubification (e.g., Sturm et al., 2001). Another focus is the implementation of carbon cycling schemes of different levels of complexity, from dedicated peatland models (e.g., Chaudhary et al., 2020; Frolking et al., 2010) to novel biogeochemistry schemes with explicit representations of microbial population dynamics (e.g., Chadburn et al., 2020). Here, it is desirable to retain the modular concept of the CryoGrid community model so that it is possible to select and test different carbon cycle schemes that can be flexibly combined with stratigraphy classes. To be able to simulate the buildup of organic soils, stratigraphy classes must be amended with an adaptive grid, which would also facilitate including sedimentation, erosion, and possibly weathering in CryoGrid simulations.

At this point, most of the simulation tools in the CryoGrid community model are designed for process-scale simulations, as their computational runtime can be substantial when simulating many grid cells. A conservative benchmark for simulation runtimes (Fig. 25) is about 200 simulated years in 24 h of computation (i.e., wall-clock time), although this can vary to some extent depending on the employed stratigraphy classes and model forcing. To conduct a 1900 to 2100 climate change simulation for the entire permafrost region at 1^∘ resolution would require about 50 000 to 100 000 CPU hours, which is perfectly feasible on modern high performance clusters, but requires significant organizational and possibly financial efforts. It is therefore desirable to implement more efficient simulation tools so that regional-, continental-, and global-scale simulations become more routine tasks.

Clustering techniques provide excellent possibilities to achieve such speedups by determining “typical representations”, given the variability of model parameters. The simulations are then only performed for these typical representations and their output is assigned to all grid cells with similar model parameters, which is justified in the light of model uncertainty. Using k-means clustering, Fiddes and Gruber (2012), Fiddes et al. (2019), and Fiddes et al. (2022) achieved a runtime speedup of 3 to 4 orders of magnitude for simulations in complex topography by clustering for terrain parameters, such as slope, aspect, and elevation. This concept could be extended to land cover maps, which are one of the input layers to model-based permafrost maps (e.g., Obu et al., 2019). Furthermore, clustering could be applied to the model forcing itself so that grid cells with similar forcing are simulated by a single representation. Such clustering would be fully compatible with the modular structure of the CryoGrid stratigraphy and could be integrated on the RUN_INFO level. With a runtime speedup of a factor of 100, for example, global climate change simulations at 1^∘ resolution would only require 500 to 1000 CPU hours, which would make them much more accessible to many users compared to a simulation without clustering.

Secondly, the runtime of the simulation tools could be reduced to facilitate application to many grid cells and/or over long timeframes. This can be achieved by simplifying the model physics, optimizing the model code, and/or using more efficient numerical integration schemes, e.g., Runge–Kutta schemes as employed in some of the CryoGrid 2 realizations (e.g., Etzelmüller et al., 2020). The latter two conflict with the modular setup of the CryoGrid stratigraphy in the TILE_1D_standard class, which is based on a common numerical integration scheme and standardized model structures to facilitate modularity. Therefore, we suggest adding further TILE classes (Fig. 1) custom-made to accommodate one or several efficient model tools. While these would not offer modularity and thus not be compatible with the stratigraphy classes presented in this study, there is ample opportunity for interplay and thus added value. On the one hand, faster models can obviously use parameterizations and model parts described in Sect. 2, even if they need to be optimized for speed. On the other hand, standardized frameworks could be implemented that allow training fast models with semi-empirical elements or unresolved processes on the process-rich but slow models of the modular CryoGrid stratigraphy (see Sect. 4.5). This could be achieved through ensemble-based data assimilation (Aalstad et al., 2018; Fiddes et al., 2019; Alonso-González et al., 2021), machine learning-based emulators (Fer et al., 2018; Dagon et al., 2020), or perhaps most promisingly a combination thereof (Bocquet et al., 2020; Brajard et al., 2020).