Impact of soil heterogeneity and lateral heat fluxes on soil temperature simulations in a permafrost-affected soil

Soil properties vary within centimeters, which is not captured by state-of-the-art land-surface models due to their kilometer-scale grid. This mismatch can lead to systematic errors when simulating the exchange of energy, water, and greenhouse gases between the land and atmosphere – collectively referred to as “aggregation error”.

To quantify the potential aggregation error of soil temperature, we developed the two-dimensional pedon-scale geophysical soil model DynSoM-2D, which has a spatial resolution of 10 cm. We applied DynSoM-2D at a permafrost-affected, non-sorted circle site using three different setups: (i) a homogeneous soil profile representing a typical land surface model, compiled by averaging the heterogeneous soil inputs; (ii) the actual heterogeneous soil profile of a typical non-sorted circle; and (iii) the heterogeneous soil profile including lateral heat fluxes.

Our results show that DynSoM simulates warmer soil temperatures when heterogeneous soil properties are considered, with this warming becoming even more pronounced and consistent across the domain when lateral heat fluxes are included. By aggregating grid cells, we traced the aggregation error back to the spatial distribution of organic matter, which nonlinearly alters soil thermal and hydrological properties, leading to the observed differences between simulations. In our case, the heterogeneity-induced warming led to a deepening of the active layer and an extension of the snow-free period, both of which can strongly alter ecosystem dynamics, while having only a minor effect on soil-atmosphere heat exchange on an annual basis.

1 Introduction

Permafrost-affected ecosystems are unique due to their extreme climatic conditions, which result in a short vegetation season and characteristic soil patterns (Siewert et al., 2021; Schädel et al., 2024). The annual freeze-thaw cycle drives distinct three-dimensional soil movements, leading to pronounced surface and subsurface heterogeneity (Ping et al., 2015; Beer, 2016). These movements, collectively known as cryoturbation, give rise to typical patterned ground features such as non-sorted circles. These circles, which contain a mix of soil materials, span several meters in diameter, while their internal texture variability occurs on the centimeter scale (Ping et al., 2015; Siewert et al., 2021).

Soil texture and organic matter content play a crucial role in determining thermal and hydrological properties, thereby influencing processes such as water flow and the formation of segregated ice (Bockheim, 2007; Nicolsky et al., 2008). In turn, these processes alter particle distributions through frost heave and sediment transport, further modifying soil texture. This bidirectional relationship between soil texture and dynamics is a hallmark of cryoturbation. It is particularly important when considering soil organic matter (OM), which is redistributed by cryoturbation from the topsoil into deeper layers (Ping et al., 2008, 2015). As a result, permafrost soils – especially turbel subsoils – often exhibit peaks in OM content at depth, a pattern that contrasts with the typical vertical decline seen in other soil types (Batjes, 1996; Jobbágy and Jackson, 2000; Harden et al., 2012; Gentsch et al., 2015; Beer et al., 2022). Since OM has distinct thermal and hydrological characteristics, its redistribution feeds back into the same processes that drive cryoturbation (Nicolsky et al., 2008; He et al., 2021; Zhu et al., 2019), forming a tightly coupled system of biotic and abiotic interactions. This system is a defining feature of permafrost soils and is deeply rooted in, and perpetuates, small-scale heterogeneity (Nicolsky et al., 2008; Ping et al., 2015).

Beyond internal soil dynamics, heterogeneity also influences the exchange of energy, water, and greenhouse gases between the land surface and atmosphere (Lawrence et al., 2015; Bonan, 2019; Schädel et al., 2024). Variations in surface and subsurface properties produce spatial patterns in heat fluxes that affect surface temperatures (Aalto et al., 2018). For instance, the texture of the soil and the topography of the land determine the seasonal snow cover and the height of vegetation, both of which vary with the season (Ekici et al., 2014; Aalto et al., 2018). Furthermore, the properties of the soil determine hydrological processes that affect the latent heat flux. Snow insulates the soil, with the effectiveness of this insulation increasing with snow depth. Snow also exhibits high radiative reflectance, which is why the presence and amount of snow have a significant impact on the surface energy balance (Ekici et al., 2014). Vegetation exerts a significant influence on exchange rates through various physiological processes, including water uptake and transpiration, which regulate the availability of water in the environment. Additionally, vegetation's unique capacity to absorb and reflect energy contributes to the regulation of energy balance within the ecosystem. The assimilation and release of carbon (C) is another crucial aspect of vegetation's impact on the environment, playing a pivotal role in the global C cycle (Miralles et al., 2025). As outlined in the work of Bonan (2019), the intricate relationship between vegetation and C dynamics is a crucial factor in understanding and mitigating environmental impacts. The latter encompasses not only the direct C emission resulting from respiration, but also the indirect influence on the C dynamics within the soil. This indirect influence is characterised by the conversion of C into soil OM and the subsequent impact on the microbial community within the soil. The function of this community is to facilitate the decomposition of OM, a process that ultimately leads to the production of greenhouse gases, namely carbon dioxide (CO₂) and methane (CH₄) (Bonan, 2019).

Soil heterogeneity, therefore, has implications not only for site-specific soil processes but also for land–atmosphere interactions on regional and global scales (Beer, 2008; Ping et al., 2015; Nicolsky et al., 2008; Nitzbon et al., 2021). This is particularly relevant for modeling current and future ecosystem dynamics (Hagemann et al., 2016; Aas et al., 2019; Cai et al., 2020; Martin et al., 2021; Schädel et al., 2024). However, Arctic soil heterogeneity is generally not represented in state-of-the-art land surface models (LSMs), as the features occur at sub-grid scales (Beer, 2016; Burke et al., 2020; Schädel et al., 2024). Soil texture is typically aggregated across kilometer-scale grid cells, dramatically reducing heterogeneity (Rastetter et al., 1992). This can lead to distortions in model outputs, since small-scale processes may interact in non-linear ways that cannot be resolved explicitly. This issue is described as aggregation error (Rastetter et al., 1992). In addition to aggregation, there is uncertainty in the parameterization of small-scale permafrost processes, which further limits the pace of improvement in LSMs for Arctic systems (Burke et al., 2020; Schädel et al., 2024).

To tackle this uncertainty, an increasing number of studies explores on the effect of permafrost and permafrost-affected soil representation in models. They range from site-level to global scales, e.g., Nicolsky et al. (2008), Westermann et al. (2016), Aas et al. (2019), Nitzbon et al. (2021), Smith et al. (2022), and De Vrese et al. (2023), and cover various questions and approaches from general plausible representations of permafrost and their implications (De Vrese et al., 2023) over commonly used assumptions and simplifications of permafrost-specific processes (Cai et al., 2020; Gao and Coon, 2022), to an increase of model resolution (Nicolsky et al., 2008; Aas et al., 2019; Aga et al., 2023). Latter allows to account for specific processes in detail, e.g., cryoturbation (Nicolsky et al., 2008) or ice segregation (Aga et al., 2023), but it also incurs high computational costs, making such resolution impractical in LSMs.

To bridge this gap, a common compromise is tiling, where grid cells are subdivided into tiles that interact with one another (Aas et al., 2019; Cai et al., 2020). It enables the execution of a model on a coarser grid, which reduces computational demands, yet implicitly accounts for grid cell-internal heterogeneity and associated processes such as lateral fluxes of energy and water, as well as the movement of mass, e.g., redistribution of snow or soil particle displacement. In the context of permafrost-affected landscapes, reasonable tilings are characterised by the separation of the rim from the centre, or the plateau from the mire, within a polygonal tundra site (Aas et al., 2019; Martin et al., 2021). These tilings differ in terms of soil texture and related soil properties, as well as in vegetation cover, height, and elevation. Alternatively, a classification based on soil ice content can be employed (Cai et al., 2020).

Despite tiling, and especially multi-scale tiling, where tiling works on different sub-grid scales, has been shown to improve the representation of permafrost-specific dynamics such as permafrost degradation (Aas et al., 2019; Nitzbon et al., 2021), it is still an approximation and lacks information about the actual spatial pattern of a certain landscape. For example, tiling schemes cannot distinguish between coarse and fine patterns if the total proportion of properties within a grid cell is the same. This limitation is critical for accurately computing lateral fluxes, which occur at the transitions between contrasting zones. Furthermore, the highest tiling resolution used so far is – to our knowledge – about one meter (Nitzbon et al., 2021), still an order of magnitude coarser than the actual centimeter-scale heterogeneity of soil texture. Consequently, the potential of tiling to address the fundamental questions surrounding the influence of soil property aggregation within a grid cell on simulation outcomes, and the role of actual spatial distribution of soil properties in this process, remains unanswered. Furthermore, the impact of resulting lateral fluxes on aggregation error, which are the foci of this study, remains to be elucidated.

Given the high small-scale heterogeneity and its effects on soil temperature dynamics in permafrost-affected ecosystems, we ask the following questions in this paper:

• What is the aggregation error in soil temperature due to sub-grid scale soil heterogeneity in non-sorted circles?

• What is the impact of heterogeneity-induced lateral heat fluxes on this aggregation error?

To address these questions, we apply a two-dimensional soil temperature model with high spatial resolution to simulate a non-sorted circle. We analyse how small-scale soil heterogeneity and heterogeneity-induced lateral heat fluxes affect surface and subsurface temperatures, surface energy fluxes, snow depth, and active layer thickness. We discuss our findings in the context of land surface modeling in general and the simulation of permafrost-affected ecosystems in particular.

2 Methods

2.1 Soil Temperature Model

To address our research questions, we apply a standard vertical soil temperature model and extend its heat conduction scheme to a second spatial dimension. The third dimension is neglected to reduce computational cost and because high-resolution three-dimensional soil data are not available. To mitigate the effect of this simplification, we adopt a radially symmetric model setup (see Sect. 2.2).

The vertical scheme of the model essentially couples soil thermal and hydrological processes through (i) phase change of soil water (Ekici et al., 2014) and (ii) changes in soil thermal properties due to soil moisture (Bonan, 2019, chaps. 5 & 8). The model also represents snow dynamics after Ekici et al. (2014). Building on this vertical framework, we extended the heat conduction scheme to a second (horizontal) dimension, which is detailed described in Sect. 2.1.1. The resulting two-dimensional, Dynamic pedon-scale Soil Model (DynSoM-2D) simulates geophysical variables such as soil temperature, and soil water and ice content, as well as sensible and latent surface heat fluxes with sub-daily temporal resolution. The spatial resolution is user-defined: column widths are uniform, while layer thickness increases exponentially with depth in the standard configuration. Near-surface layers begin at centimeter scale and gradually increase to several meters at depth to ensure numerical stability. For our simulations, we chose a set-up of 10 columns with a width of 10 cm each, adding to a total horizontal extension of 1 m, and layer thicknesses of 10 cm within the uppermost meter, followed by layers with increasing thicknesses of 1, 3, 5, 10, and 30 m, adding to a total depth of 50 m.

DynSoM-2D requires information on geographic location, vegetation cover and soil texture, which can be recalculated at finer or coarser spatial resolution. The required soil texture information is the organic matter content and the fractions of clay and silt as part of the mineral soil, from which the sand fraction is calculated. Other (potentially) measurable soil information, such as bulk density or soil porosity (Θ_pore) are not needed, but calculated from given the soil texture information. Bulk density is calculated based on equations by Martín et al. (2017) and Leonaviciute (2000), and other site-specific parameters are calculated through so-called pedotransfer functions, which estimate physical parameters, i.e. thermodynamic and hydrological parameters, from measurable soil properties. In the case of soil porosity, which is one of the vanGenuchten parameters (Van Genuchten, 1980), DynSoM-2D applies the pedotransfer functions from Wösten et al. (1999), which give values between 0.36 (silty clay) and 0.45 (silty loam). The estimated physical parameters are then used for further calculations.

In addition to this constant site information, DynSoM-2D requires a site-specific atmospheric forcing as an upper boundary condition. Specifically, DynSoM-2D uses precipitation, which assumed to be snow, if the air temperature is below 0 °C, and serves as upper boundary for the soil hydrological scheme (cf. Appendix A2), air temperature (only for the rain-snow separation), incoming long-wave and short-wave radiation (to calculate the energy balance at the surface, which is used as upper boundary for the soil heat conduction scheme, which is described in the following section), surface pressure, specific humidity, water vapor pressure and wind velocity.

In the following, we present the heat conduction scheme that is used in DynSoM-2D, and the calculation of the surface energy balance, with a special focus on their links to soil texture and its heterogeneity.

2.1.1 Heat Conduction Scheme

We assume that the soil temperature dynamics are governed by heat conduction and water phase change. The classical one-dimensional heat conduction equation is extended to a second (horizontal) dimension (Eq. 1). Instead of solving the two-dimensional heat conduction equation in an implicit scheme, we solve the vertical and horizontal components successively using a classical implicit one-dimensional finite difference method according to Bonan (2019), Sect. 5.2 and 5.3. This approach is justified by the short time step of one hour (Bonan, 2019). In addition, the latent heat of fusion during phase change of water is considered following the approach of Ekici et al. (2014).

$$\begin{array}{}\text{(1)}& {\displaystyle }{c}_{\mathrm{v}}{\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{\Delta }t}}={\mathit{\kappa }}_x{\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{2}}T}{\mathrm{\Delta }{x}^{\mathrm{2}}}}+{\mathit{\kappa }}_z{\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{2}}T}{\mathrm{\Delta }{z}^{\mathrm{2}}}}+\mathrm{\Delta }{T}_{\mathrm{pc}},\text{with}\text{(2)}& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{pc}}=E_{\mathrm{pc}}{\mathit{\rho }}_{\mathrm{i}}{\displaystyle \frac{\mathrm{\Delta }{\mathrm{\Theta }}_i}{\mathrm{\Delta }t}}\end{array}$$

with ΔT representing the change in soil temperature in K in time (Δt, in s) or space (Δz or Δx respectively, in m), as well as due to phase change (ΔT_pc). c_v is the volumetric heat capacity in J m⁻³ K⁻¹, and κ_z and κ_x are the heat conductivities in different directions in W m⁻¹ K⁻¹ (cf. Appendix A1). ΔT_pc is given by the specific latent heat for fusion (E_pc in J kg⁻¹), as well as the change in ice content (kg kg⁻¹) within the given time step ($\frac{\mathrm{\Delta }{\mathrm{\Theta }}_i}{\mathrm{\Delta }t}$) and ice density (ρ_i) in kg m⁻³.

At the bottom boundary, a zero-flux condition is assumed due to the large layer thickness and minimal thermal gradient. At the surface, a prescribed skin temperature (see Sect. 2.1.2) serves as the upper boundary condition. Also for the lateral boundaries we assume a zero flux and define the spatial domain accordingly (see Sect. 2.2).

Soil temperature and heat fluxes are governed by the thermal conductivity (κ) and volumetric heat capacity (c_v), both of which depend on soil texture and composition. Typical values for various soil components are provided in Table 1.

To estimate the effective thermal conductivity for each grid cell, DynSoM-2D uses the Johansen method (Farouki, 1981), which interpolates between dry and saturated conductivity values. These endpoints, κ_dry and κ_pore, are calculated as volumetric-weighted averages based on the thermal conductivities of the constituent materials. For κ_dry, this includes the mineral fractions (clay, silt, sand) and organic matter. For κ_pore, the pore-filling components – air, liquid water, and ice – are considered. The volumetric heat capacity for each grid cell is similarly computed as a weighted average of the component-specific heat capacities, using the volume fractions of each soil constituent.

Table 1 Thermal properties for solid, liquid and ice components of soil originally taken from De Vries and Van Wijk (1963).

2.1.2 Surface Energy Balance & Surface Temperature Calculation

The surface temperature, which is used as the Dirichlet boundary condition for the soil temperature calculation, is determined by the energy exchange between the atmosphere and the soil, i.e. the surface energy balance (Eq. 3). DynSoM-2D calculates the net radiation at the surface (R_n) from incoming solar radiation (S_r), which depends on date, time, and location and is adjusted to the prevailing weather conditions, partly reflected by the surface albedo (α, Eq. 4), as well as incoming and outgoing longwave radiation (L_in, L_out), which is determined by the temperatures, e.g., of the surface and clouds. The net radiation is balanced by the sensible heat flux (H), the latent heat flux (λE), and the soil heat flux (G). Turbulent surface fluxes are parameterized using Monin-Obukhov Similarity Theory (MOST) with a shared atmospheric reference state. We note that the assumptions underlying MOST, in particular horizontal homogeneity of the surface layer, are not strictly fulfilled at the horizontal scale of individual soil patches on the order of decimeters represented in DynSoM-2D. MOST is therefore not interpreted here as a physically valid description of patch-scale turbulence, but is employed as a pragmatic closure to parameterize surface-atmosphere exchange.

All soil patches interact with identical atmospheric conditions at the reference height (2 m, taken from CRUNCEP data; Viovy, 2018), and differences in turbulent fluxes arise exclusively from differences in surface temperature and moisture states. The resulting patch-scale fluxes should thus be interpreted as diagnostic contributions to the area-averaged surface flux rather than as independent local turbulent fluxes.

In the model configuration used in this study, only bare soil is considered and topographic differences between columns are neglected, implying uniform surface height and roughness. Snow is added to the soil scheme but does not affect surface height in this configuration. For other model configurations that introduce roughness differences due to vegetation or topography, DynSoM-2D couples MOST with a roughness sublayer parameterization following Harman and Finnigan (2007, 2008).

$$\begin{array}{}\text{(3)}& R_{\mathrm{n}}=(\mathrm{1}-\mathit{\alpha })S_{\mathrm{r}}+L_{\mathrm{in}}-L_{\mathrm{out}}=H+\mathit{\lambda }E+G\end{array}$$

with α individually calculated for the visible (VIS) and near-infrared (NIR) wavelengths and depending on dry soil albedo (α_dry) and relative soil moisture, if soil is not snow covered (cf. Appendix A2).

$$\begin{array}{}\text{(4a)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{soil}}^{\mathrm{VIS}}={\mathit{\alpha }}_{\mathrm{dry}}^{\mathrm{VIS}}-\mathrm{0.35}\cdot {\displaystyle \frac{{\mathrm{\Theta }}_{\mathrm{w}}}{{\mathrm{\Theta }}_{\mathrm{pore}}}},\mathrm{with}\text{(4b)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dry}}^{\mathrm{VIS}}=\mathrm{0.4},\mathrm{and}\text{(4c)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{soil}}^{\mathrm{NIR}}={\mathit{\alpha }}_{\mathrm{soil}}^{\mathrm{VIS}}+\mathrm{0.05}\end{array}$$

Based on this balance, DynSoM-2D determines the heat in the soil using a skin layer approach (Bonan, 2019, Sect. 7.3), which concentrates the energy balance between the atmosphere and the soil on a theoretical skin layer that is part of the soil profile. The heat is recalculated to a skin layer/surface temperature, which is then used as an upper boundary condition for the soil temperature calculation (Eq. 1).

The surface energy balance is strongly influenced by vegetation and snow, or the presence of a pond, which are possible surface covers in DynSoM-2D that change the surface albedo (α) and thermal properties.

Note, that the term surface in the DynSoM context always refers to the surface that is directly interacting with the atmosphere.

For snow, DynSoM-2D applies a minimum threshold of 3 cm for inclusion in energy and heat transfer calculations, following Ekici et al. (2014). This threshold accounts for the patchy nature of early and late snow cover. Insulation effects depend on total snow depth, which is discretised into 3 cm layers. Each layer is included in the vertical heat conduction calculations once it is fully filled. Snow layers are treated analogously to soil layers for vertical conduction, but are excluded from the lateral conduction scheme. This is based on the assumption that lateral temperature gradients in snow are minimal and that the thermal conductivity of snow is very low. Albedo changes from bare soil albedo (α_soil), which is depending on soil wetness, to a constant snow albedo (α_snow) as soon as the snow depth exceeds the first 3 cm. As albedo values vary with the wavelength of incoming radiation, DynSoM-2D also applies different albedo values for radiation within the visible (VIS) and near-infrared (NIR) spectrum, following Bonan (2019) (Sect. 3.2) (cf. Sect. 2.1 and Appendix A2).

2.2 Model Setup and Protocol

For this theoretical study, we initialised soil texture and organic matter profiles based on a representative example of a non-sorted circle described by Gentsch et al. (2015) (Supplement, Fig. S7, profile CH-E), located near Cherskii, East Siberia, Russia (68.7° N, 161.4° E). We resampled this data to a 10 cm × 10 cm grid over a spatial domain of 1 m in both vertical and horizontal directions. Horizontally, the 1 m transect spans from the center of a non-sorted circle to its rim, allowing us to assume zero lateral fluxes at both boundaries – i.e., at the circle center and the midpoint to the adjacent circle. Vertically, the shown and analysed domain extends to a depth of 1 m, which reaches the parent rock layer at this site, whereby the actual simulation domain is 50 m depth (cf. Sect. 2.1).

In addition, we averaged soil state variables laterally to obtain a reference homogeneous soil. The final set up is represented in Fig. 1 and Table 2.

Figure 1 Schematic representation of soil input matrices for homogeneous soil (left), and heterogeneous soil (right). Values for clay, silt and sand, as well as for organic material can be found in Table 2. Red dashed lines separate topsoil (0–20 cm depth) from subsoil (20–80 cm depth) from deep soil (80–100 cm depth). Please note, that the actual model set-up has also deeper layers up to 50 m depth, which are not shown here. As deeper layers are assumed to be part of parental rock, too, their clay, silt, and sand values are equal to the C-horizon input.

Table 2 Soil texture information for homogeneous soil (left) and heterogeneous soil (right). O-A-B, A-B, B, B-OA, C, A, & OA refer to horizon classifications, whereas top/sub/deep in brackets refer to the position of the specific horizon grid cell in the soil input, i.e. topsoil, subsoil, and deepsoil (cf. Fig. 1). Values refer to shares of clay, silt, sand, and OM in total soil, whereas values in brackets show shares of clay, silt, and sand from mineral soil only.

We excluded pond formation and any vegetation to isolate the effects of soil texture heterogeneity. However, snow cover was included, given its essential role in permafrost energy dynamics. Atmospheric forcing was obtained from the CRUNCEP Version 7 dataset (Viovy, 2018), using data specific to the study site.

We ran three different model simulations:

• M0 – using the homogeneous soil input. This represents the standard one-dimensional soil model without pedon-scale heterogeneity and serves as the reference scenario.

• M1 – using heterogeneous soil input, but assuming no lateral heat fluxes. The model structure remains one-dimensional but with fine-scale horizontal resolution.

• M2 – using heterogeneous soil input with lateral heat transport enabled, thus extending the model to a fully two-dimensional setup.

To initialise the soil state variables, we performed a 50-year spin-up for each setup individually, repeating the 14-year atmospheric forcing. From the final run, we took the penultimate year for our analyses, because we needed three additional days from the years before and after the year being analysed to support the one-week running mean that we used to smooth our simulation results (unless otherwise stated). We chose to smooth our output in this way to avoid overly fuzzy near-surface results, which are strongly influenced by atmospheric forcing when the output is not smoothed, but to preserve the current atmospheric forcing signal, which is strongly superimposed when another averaging method is used, e.g., a 10-year mean.

Our analysis focuses on soil temperature differences among the three model configurations, specifically the effects of sub-grid soil heterogeneity. We distinguish between the inter-model range – i.e., differences in mean results across M0, M1, and M2 –, which may be also presented by the absolute differences (Eq. 5, ΔX) between the heterogeneous soil models (M1, M2) and the homogeneous soil model (M0), and the intra-model range – i.e., spatial variation among columns in M1 and M2.

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }X=X^{\mathrm{het}}-X^{\mathrm{hom}}\end{array}$$

3 Results

3.1 Soil Parameters

Simulated soil porosity differs slightly between the homogeneous and heterogeneous soil setups, primarily due to the spatial variability of organic matter (OM) content in the latter (Table 3). This variability has a pronounced effect on thermal properties. In general, the horizon-wise averaged heat conductivities (κ) differ between the heterogeneous and the homogeneous soil, because the relationship between heat conductivity and organic matter (OM) content (Θ_OM) is not linear. OM-rich patches within the heterogeneous soil lower the mean heat conductivity more than an even distribution of OM within the homogeneous soil, even if the total amount of OM is the same. Conversely, the averaged heat capacity (c_v) is higher in the heterogeneous soil. This is linked to slightly increased porosity, which enhances total water content under saturated conditions. Since water has a higher heat capacity than mineral soil components, this leads to an overall increase in soil heat capacity. In contrast, the calculated surface albedo is less affected. Visible wavelength albedo (α^VIS) varies between 0.19 and 0.34 (M0), 0.21 and 0.28 (M1), and 0.18 and 0.25 (M2) respectively, during the snow-free period depending on the relative soil water content of the uppermost layer (α^NIR varies similarly; cf. Eq. 4 and Fig. B2c in Appendix B).

Table 3 Simulated values for porosity (Θ_pore) and resulting values for heat conductivity (κ), and heat capacity (c_v) for unfrozen, saturated soils. The values are spatially averaged for topsoil (0–20 cm depth), subsoil (20–80 cm depth), and deep soil (80–100 cm depth) separately (cf. Fig. 1), as well as for the entire soil (0–100 cm depth). Additional simulated soil properties can be found in Table B1.

3.2 Temporal Pattern

The surface heat flux to the atmosphere (F_surf), defined here as the sum of the sensible heat flux (H) and the latent heat flux (λE, both see Eq. 3) ranges from −0.32 W m⁻² in winter (all models) to 0.45 W m⁻² (M0) to 0.51 W m⁻² (M2) in summer (Fig. 2a). We find almost no differences between models within the first half of the year and also at the end of the year. The intra-model ranges, i.e. the range among columns within the heterogeneous soil models M1 and M2, are also small. Notable differences emerge only during summer and autumn when the soil is still snow-free. From July onwards, M0 emits significantly less heat on average than M1 and M2, which changes in mid-August until the mean curves converge again in September. This pattern corresponds to differences in (upper) soil moisture content, which influences surface albedo values (Sect. 3.1; cf. Eq. 4 and Fig. B2). However, the inter-model range, i.e. the difference between the model means, does not exceed 0.1 W m⁻², indicating a difference of about 20 % in heat fluxes to the atmosphere due to soil heterogeneity. In contrast, the intra-model ranges of heat fluxes increase to 0.6 W m⁻² (M1) and 0.9 W m⁻² (M2).

Despite these temporarily wider inter- and intra-model ranges, the total simulated annual heat budget differs by less than 1 W m⁻² between all model versions (M0: 13.7 W m⁻², M1: 14.1 W m⁻², M2: 13.5 W m⁻²).

Figure 2 Single-year simulation results (left; M0: black, M1: purple, M2: green) and deviations from reference model M0 (right; M1: purple, M2: green; Eq. 5) for surface heat fluxes (top row), surface temperatures (middle rows), and soil temperatures (bottom row). Lines represent spatial mean, whereas shaded areas show the intra-model ranges, i.e. the ranges among columns for M1 and M2. The full simulated timeseries can be found in Appendix B (Fig. B1). (a) Simulated surface heat flux (F_surf), which is defined as sum of sensible heat flux (H) and latent heat flux (λE) at the surface, which directly interacts with the atmosphere, i.e. soil surface or snow. (b) Differences of surface heat fluxes (ΔF_surf). (c) Simulated surface temperature (T_surf, solid lines), whereby surface refers to the surface that is interacting with the atmosphere directly, as well as soil temperature (T_soil, all depths up to 1 m depth, dashed lines), and atmospheric forcing temperature ($T_{\mathrm{atm}}^{\mathrm{forcing}}$, dotted line) for comparison. (d) Differences in surface temperature (ΔT_surf, solid lines) and soil temperature (ΔT_soil, all depths up to 1m depth, dashed lines). (e) Simulated soil temperatures in topsoil (0–20 cm depth; T_top, solid lines), subsoil (20–80 cm depth; T_sub, dashed lines), and deep soil (80–100 cm depth; T_deep, dotted lines). (f) Differences in topsoil (ΔT_top, solid lines), subsoil (ΔT_sub, dashed lines), and deep soil temperatures (ΔT_deep, dotted lines).

The simulated seasonal variations in surface and soil temperatures (T_surf, T_soil) across all model versions closely follow the atmospheric temperature seasonality ($T_{\mathrm{atm}}^{\mathrm{forcing}}$) (Fig. 2c, Sect. 2.2), but the annual amplitude is reduced from 67.6 °C ($T_{\mathrm{atm}}^{\mathrm{forcing}}$) to less than 40 °C at the surface (M0: 38.5 °C, M1: 39.4 °C, M2: 38.8 °C), and around 25 °C in the soil (M0: 25.5 °C, M1: 25.3 °C, M2: 25.4 °C).

Similar to the surface heat flux results, surface and soil temperature results do not differ much during the first half of the year, when the soil is frozen and snow-covered (Fig. B2). From July to September, and again in November/December, the mean temperature results deviate among models. During these periods, the homogeneous soil model M0 simulates a cooler surface than the heterogeneous soil models M1 and M2 (Fig. 2c and d). However, during the last week of August (DOY 236-244), we observe a cooling effect due to soil heterogeneity, which is strong enough to affect the monthly mean temperatures of the surface and topsoil (Table 4).

The inter-model range of the simulated surface temperatures does not exceed 2.4 °C on daily basis, and is therefore smaller than the individual intra-model ranges of M1 and M2 that increase up to 5.1 °C (M1), and 4.9 °C (M2) respectively. Interestingly the lateral fluxes, which lower the intra-model range in M2 compared to M1, increase the overall difference between M2 and M0 (annual mean: +0.25 °C) compared to M1 and M0 slightly (annual mean: +0.19 °C; see also Fig. 2d). We find a similar pattern for the soil temperature ranges, where either the inter-model range (maximum 2.8 °C) and the intra-model range of M1 (4.5 °C) are similar to the ranges of the surface temperature results, whereas the intra-model range of M2 is much smaller (0.9 °C). Both the inter- and intra-model ranges of simulated temperatures decrease when comparing monthly means, but are still present for annual means, which differ by 0.2 °C (T_surf) and 0.2 °C (T_soil), respectively (Table 4).

Table 4 Monthly and annual mean of surface temperature and soil temperatures in different depths (topsoil: 0–20 cm; subsoil: 20–80 cm; deep soil: 80–100 cm) in °C spatially averaged over all columns (upper row), as well as intra-model range (min/max; lower row), where present.

Compared to the simulated surface temperatures (T_surf), the simulated soil temperatures (T_soil) differ only slightly, due to the fact that they are not only horizontally averaged over all ten columns, but also represent the entire upper soil down to 1 m depth (Fig. 2c, Table 4). We observe the annual cycle of temperatures at all depths, with the amplitudes decreasing with depth, as the signal from the atmospheric forcing is attenuated by each soil layer above (Fig. 2e and f). Simulated topsoil temperatures (T_top, 0–20 cm depth) range from −27.7 to 8.7 °C, subsoil temperatures (T_sub, 20–80 cm depth) range from −23.9 to 0.8 °C, and deep soil temperatures (T_deep, 80–100 cm depth) range from −20.6 to −0.6 °C, with inter- and intra-model ranges varying between seasons.

We find only minor deviations between models (0.2 °C on average until June) and a decreasing trend in the intra-model range during the first half of the year, when the soil is still frozen (Fig. B2e). However, the intra-model ranges for M1 are notable (topsoil: 0.7 °C, subsoil: 0.8 °C, deep soil: 0.9 °C), and are certainly greater than the intra-model ranges for M2 (topsoil: 0.0 °C, subsoil: 0.1 °C, deep soil: 0.0 °C).

During summer, the inter-model ranges increase up to 1.2 °C in the topsoil and 0.2 °C in the subsoil, and the intra-model ranges also increase. Particularly striking is the intra-model range in the topsoil of M1, which exceeds 4 °C on a daily basis and 2 °C on a monthly basis across the simulated distance of 1m, and which is partly reflected in the relatively high inter-model ranges between M1 and M0 (daily: max. 1.2; monthly: max. 0.4 °C) and M2 (daily: max. 1.0; monthly: max. 0.4 °C).

After freezing, i.e. when the soil is fully frozen up to 1 m depth, the simulated monthly mean temperatures differ by up to 1.5 °C on average among all models across all depths, even on a monthly basis. However, the inter-model ranges are still mostly surpassed by the intra-model ranges. The spatial variation is around 2 °C across all depths in M1 (topsoil: 2.2 °C, subsoil: 2.4 °C, deep soil: 1.8 °C), but is significantly lower in M2. The inter-model range of M2 is only high in topsoil in October (1.0 °C), but drops to below 0.5 °C in November and December for both the subsoil and deep soil. Notably, M2 results are consistently warmer than M0 results, whereas the intra-model range of M1 and the respective differences show that some columns are warmer and some are colder than M0, but the mean is warmer (Fig. 2f).

3.3 Spatial Pattern

Simulated monthly mean temperatures in June do not differ significantly when comparing individual grid cells (Fig. 3), which aligns with our temporal analysis since the soil is still frozen. Only two columns simulated by M1 are slightly warmer than M0, but this is balanced out by lateral fluxes in M2. Soil temperature differences begin to evolve in July, with the two heterogeneous-soil models M1 and M2 simulating slightly warmer temperatures than the homogeneous-soil model M0 in the uppermost soil layers, which start thawing in mid-July (DOY 194). Thawing is uniform within all model variants and across all columns, reaching 30 cm depth by the end of the month. We find that the maximum absolute warming of individual grid cells is equal in M1 and M2 compared to M0 (max. +0.8 °C), but the affected area is larger in M2, where 27 grid cells are warmer than in M0 compared to 21 grid cells in M1. This leads to a minor but significant increase in average warming by 0.02 °C across the entire domain compared to M0.

Figure 3 Mean differences in simulated soil temperature (Eq. 5) per grid cell and max. active layer depth (ALD) for the months June–November. Upper row: Differences between M1 and M0. Lower row: Differences between M2 and M0. Colors represent differences in soil temperature. Thick lines show maximum active layer depth for M0 (solid, black) and M1/M2 (dashed, purple/green). Thin dotted lines separate topsoil from subsoil from deep soil. Full year data can be found in Fig. B3 in Appendix B.

The spatial enlargement of heterogeneity-induced effects by lateral fluxes is even stronger in August. M1 shows a patchy pattern of warming and cooling, especially in the topsoil, compared to M0, whereas M2's response is more laterally consistent, reflecting differences in soil moisture content that affect surface albedo (Eq. 4 and Fig. B2). However, spatial mean temperatures in M0 and M1 are almost the same, whereas M2 simulates a cooler topsoil above a warmer subsoil (Table 4, see also Fig. 2f). Notably, the warmer patches in M1 compared to M0 are associated with grid cells that have a relatively high organic matter (OM) content (Fig. 1), while cooler patches are associated with grid cells with a lower OM content. This is most evident at the surface, but also in the deeper soil. Warming influences soil thawing and thus the thickness of the active layer. M0 simulates a maximum active layer thickness of 40 cm, while the active layer thickness varies between 40 and 60 cm in M1 and between 50 and 60 cm in M2. M1 simulates the deepest thawing in columns with deep OM-rich layers, and that this thawing is extended laterally in M2.

The temperature differences between M1 and M0 and M2 and M0 increase further in September, when the air temperature is already cooling at the surface (Fig. 2c and Table 4). M1 simulates a column-wise patchy pattern of warming and cooling (max. ±1.2 °C, +0.2 °C on average) up to 40 cm depth, and a warming of up to 0.4 °C (+0.2 °C on average) below, whereas M2 simulates a strong warming up to 1.0 °C in the upper layers and a moderate warming (+0.3 °C on average) below compared to M0 (Fig. 3). Despite these larger differences in monthly mean temperatures, the maximum active layer depths remain the same in all model versions as in August.

Soil cooling continues in October, which leads to a retreat of the maximum active layer depth in M1 and M2 from a maximum of 60 to 50 cm. However, this is still 10 cm deeper than the deepest thawing in M0. Besides, we find only minor differences between soil temperatures in topsoil, i.e. until 20 cm depth, simulated by M0 and M1 or M2 respectively. Except of a warmer hot spot in the uppermost layer of column 9 and 10, where the OM-rich O-horizon is located (+1.7 °C in M1, +1.5 °C in M2). Below the topsoil, there is a 20 cm-thick layer where no temperature differences are observed among models or columns. This is the so-called zero curtain, where soil temperature has reached 0 °C and is freezing without further cooling until frozen. Since this happens in all simulations, temperature differences vanish. Below this freezing zone, M1 and M2 simulate a significantly warmer soil than M0 (M1: +0.2 °C on average; M2: +0.3 °C on average).

Soil freezing progresses until the soil is completely frozen after the first week of November, leading to the largest temperature differences between models on a monthly basis. During this week the soil freezes from both sides (Fig. 3), i.e. from the surface, which is strongly cooled by the atmosphere and from colder frozen layers below (Fig. 2 and Table 4). The M0 soil is already frozen after DOY 309, whereas M1 has a continuous unfrozen layer only until DOY 308, but remaining unfrozen cells until DOY 312, showing the uneven occurrence of cooling and freezing within a heterogeneous soil. On a monthly basis, M1 simulates soil temperatures ranging from −0.7 to +2.1 °C compared to the corresponding M0 grid cells (+0.6 °C on average). M2 has a continuous unfrozen layer until DOY 311, as laterally transported heat generally delays freezing, but then freezes completely within only two days, as cooling is more uniform due to lateral fluxes. Nevertheless, M2 grid cells are up to 1.9 °C warmer than M0 grid cells (average +1.4 °C; Fig. 3 and Table 4).

3.4 Simulated Soil Temperature Relation to Soil Texture and Soil Organic Matter Content

To explore the influence of soil texture in general and organic matter (OM) content (Θ_OM) specifically on simulated soil temperatures in more detail, we aggregated all grid cells based on the information of the heterogeneous soil input data, i.e. O, A, B, OA and C as shown in Fig. 1, as well as based on their difference in OM content of the heterogeneous and homogeneous soil input data (Table 2), and compare simulated soil temperature differences (ΔT_soil, Eq. 5).

Figure 4 Differences in simulated daily soil temperature per grid cell between M1 and M0 (purple), as well as M2 and M0 (green) calculated by Eq. (5). Grid cells are aggregated by horizon classes of the heterogeneous soil input (left), as well as by differences in OM content (Θ_OM) between the heterogeneous and the homogeneous soil inputs (right; values indicate the number of grid cells that fall in the respective category). The aggregation based on actual differences in OM content can be found in Fig. B4. Boxes represent quartile range of data and whiskers show the 0.1–0.9-quantile range. Circles show data median, whereas stars indicate data mean.

We find that grid cells in the heterogeneous O-horizon simulate higher temperatures than the corresponding grid cells in the homogeneous soil (OAB horizon, see Fig. 1). However, this relationship is weaker or not detectable for the other grid cells (Fig. 4, left).

By aggregating the grid cells by their difference in OM content in the soil input (Fig. 1), which can be higher in the heterogeneous soil compared to the homogeneous soil, equal, or lower, depending on the certain grid cell within its layer, we find that grid cells with a higher OM content tend to be warmer than grid cells with a lower OM content (Fig. 4, right). This relationship is clearer when comparing simulations from M1 and M0, which show the direct effect of OM on soil temperature simulations, than when comparing soil temperature simulations from M2 and M0, where lateral heat fluxes compensate for heterogeneity-induced differences between grid cells. However, differences are not significant, and we cannot derive any clear mathematical relation between soil temperature shifts and OM differences. In particular, for all data sets the median of the differences is much lower than the mean, which is mostly in the upper quartile of the differences. This indicates that there are fewer days with very large differences that substantially increase the mean, and many days with small differences between simulated temperatures.

4 Discussion

We developed the high-resolution geophysical soil model DynSoM-2D and run the model at a typical non-sorted circle profile in three different variants: with aggregated soil properties as reference (M0), considering heterogeneous soil, but without lateral fluxes (M1), and considering both heterogeneous soil and lateral heat fluxes (M2).

We analysed time series of simulated daily mean surface heat fluxes, surface and soil temperature, as well as temperature differences of monthly mean temperatures for each grid cell between these model setups. Additionally, we explored the effect of the homogeneous versus heterogeneous soil representation in general, and the impact of organic matter (OM) content on soil temperature in particular.

Over the whole season, differences among models, i.e. the inter-model differences, are small when comparing mean annual values, and even smaller, when comparing median values. This is reasonable because all three model versions were run with the same forcing, with only the resolution of soil properties differing. In addition to small inter-model differences, we found comparatively large variations between columns of the heterogeneous setup, the so-called intra-model differences of M1 and M2, which in most cases exceed the inter-model differences.

Generally, both, inter- and intra-model, differences are consistent with findings from previous studies, where differences of 0.5 to 2 °C (Beer, 2016), up to 2 °C (Aas et al., 2019), or below 2.5 °C (Smith et al., 2022) are found between tiles, grid cells or pixels, respectively. Notably, those reported values are closer to our wider intra-model ranges, because they actually represent the differences between center and rim specifically or different pixels with varying soil parameters, which correspond to the (outer edge) heterogeneous columns of our models M1 and M2, and thus to the intra-model ranges.

Differences in simulated temperatures and heat fluxes usually occur during summer, when the soil is snow-free and partly unfrozen. In general, soil heterogeneity has a warming effect in our simulations. Although we find temporarily wide deviations among simulated surface heat fluxes on a daily basis among models, but also across columns within models, the annual results differ little, because columns of the heterogeneous soil models that absorb most heat, also emit most heat (Fig. 2a and b). Thus, heterogeneity effects on opposite processes tend to cancel each other out to some extent.

The overall warming of the soil is slightly stronger and certainly more consistent across the columns when considering lateral fluxes. This shows that lateral fluxes can be important for the overall site-level soil temperature dynamics (Figs. 2 and 3). Lateral fluxes reduce the intra-model range of soil temperature and extend the intra-model range of heat fluxes (Fig. 2b and d). This is due to lateral heat transport from surface grid cells that absorb most heat, or towards such grid cells that emit most heat, and explains, why the annual heat budget actually differs more between M0 and M1 than between M0 and M2, whereas the opposite is true for daily heat fluxes.

The difference in simulated temperatures is caused by differences in the thermal properties of the soil that alter the distribution of the heat flux (Bonan, 2019). We traced these differences back to the spatial distribution of organic matter (OM), which affects pore size distribution. Pore space can be either filled with air, which would lead to negative temperature changes in highly porous grid cells due to less heat storage, or water (or ice), which would lead to positive temperature changes in highly porous grid cells due to a much higher heat storage capacity. Our model tends to be very moist and close to saturation, which is why we observe a warming in OM-rich grid cells.

This link between OM and temperature changes is clearer when comparing M1 and M0 than when comparing M2 and M0 (Fig. 4), because lateral fluxes in M2 redistribute heat among columns, thus offsetting heterogeneity-induced temperature differences. However, even for M1, we cannot find a clear relation between the relative difference of OM of a certain grid cell and the resulting difference in simulated temperature. This is due to the non-linear influence of OM on thermal properties, and possibly due to the impact of water content on thermal properties.

Neither for M2, nor for M1 we could detect a clear relationship between the soil input data, except for the O-horizon grid cells, and the change in simulated temperature for individual grid cells. This is, because the actual spatial location of a particular grid cell is important. For example, an A-horizon grid cell of heterogeneous soil will have a higher OM content when compared to an AB-horizon grid cell of homogeneous soil (e.g., 10–20 cm depth, see Fig. 1 and Table 2 for comparison), but a lower OM content when compared to an (homogeneous soil) OAB horizon grid cell (0–10 cm depth). In the first case, the A-horizon grid cell will simulate a warmer temperature than the corresponding AB-horizon grid cell, whereas in the second case it will simulate a cooler temperature than the corresponding OAB-horizon grid cell. This influence of the actual spatial location of certain grid cells, or rather the influence of the actual distribution of OM in the soil, questions the sufficiency of tiling approaches to account for sub-grid scale heterogeneity. Although tiling is a necessary step forward to account for sub-grid scale heterogeneity (Aas et al., 2019; Cai et al., 2020; Martin et al., 2021), the continued neglect of the actual spatial distribution of heterogeneous patches may reduce its utility.

Snow covers the surface and by that overlays any heterogeneity effects on surface albedo and insulates soil from the atmosphere, thus equalizing simulated surface heat fluxes. The only time when snow impacts surface heat fluxes, which are simulated by the three different models, differently is at the end of the snow melt period. Snow melts unevenly on heterogeneous surfaces, which is why some columns of the heterogeneous soil are snow-free earlier than others, and also earlier (or later) than respective columns of the homogeneous soil. Since snow has a very different albedo from (bare) soil, the timing of a snow-free surface majorly affects surface heat fluxes. However, this effect is almost balanced out by the fact that only some columns of the heterogeneous soil show an earlier snow-free surface, while others are snow-covered longer than the homogeneous soil columns. The reason for the uneven snow melt is a change in the thermal properties of the (sub-)surface grid cells due to the heterogeneous soil textures. This change alters heat transport, which can either warm cells and lead to an earlier (and potentially faster) snow melt, or cool cells, leading to a later snow melt (Figs. 2a and b and B2a and b). Since land surface models underestimate the speed of snow melt and hence have soil temperature biases during the snow melt period Ekici et al. (2014), considering this effect of soil heterogeneity on snow melt could help improve future models.

Ice has a higher thermal conductivity than water (Table 1). Therefore, heat is better transported through frozen soil than through unfrozen soil. This happens throughout the year in the deeper, perennially frozen layers, which is why the temperature of the deep soil layers differs least between the M1 and M2 columns, and also between the models. However, it also causes the differences in simulated temperatures of upper layers, which thaw during the summer, to disappear (Fig. 2e and f). This is particularly noticeable in early winter, when large differences due to uneven freezing disappear from November to January (Figs. 3 and B3). As the soil cools, the grid cells reach the freezing temperature of 0 °C at slightly different times due to their previous temperature differences and their specific heat capacity and conductivity, both of which depend on the soil texture. They then begin to freeze, whereby the duration of freezing depends (i) on the required energy for phase change, and (ii) on the allowance of additional lateral heat fluxes. The latter delays the time of initial freezing within a given soil layer, as heat is primarily transported between the columns until all have cooled to 0 °C before the first grid cell freezes, but also shortens the duration of the layer's freezing, as all the grid cells have reached freezing temperature by the time the first cell freezes.

Due to the high heat release of freezing and the resulting time shift of further cooling, temperature differences that evolved during summer are intensified and lead to very high temperature ranges after freezing. Notably, we do not see a similar increase of deviations during thawing, because the temperature differences before thawing are small, and heat is transported comparably fast through frozen soil, so that the soil thaws relatively even among columns, but also among models. Thus, thawing happens rather even, whereas freezing happens very uneven within the heterogeneous soil and builds frozen patches that turn into cooler spots.

In terms of an aggregation error, our findings demonstrate that there is a notable shift towards lower soil temperatures when neglecting sub-grid scale soil heterogeneity and the actual spatial distribution of soil properties, which cannot be resolved by tiling. We admit that the shift towards higher temperatures when considering small-scale heterogeneity is caused by a high soil moisture, and that simulations with lower soil moisture may show the opposite, but the temperature shift, i.e. the aggregation error, itself is substantial. In addition, the temperature change is larger when comparing M2 and M0 than when comparing M1 and M0, showing an amplification of the error when considering heterogeneity-induced lateral fluxes in the small-scale model, which will likely be consistent for any aggregation error. Notably, the detected aggregation error does not affect the annual surface heat budget much, meaning that the heterogeneous soil models simulate a warmer soil without absorbing (or emitting) significantly more energy on annual basis. This finding may reduce the impact of the detected aggregation error on regional or global simulations, as it propagates only marginally into the atmosphere on an annual basis.

However, and although a warming of 0.23 °C on average between M0 and M2 (+0.11 °C on average between M1 and M0) seems minor at first sight, it has great effects on this specific permafrost-affected site by deepening soil thawing and speeding up snow melt, which may impact the dynamics of the entire ecosystem largely (Chadburn et al., 2015b, a; Beer, 2016; Beer et al., 2018; Aas et al., 2019; Burke et al., 2020; Rixen et al., 2022; Schädel et al., 2024). In addition, the bias of individual grid cells up to 5 °C in summer will affect temperature-dependent soil processes such as heterotrophic respiration. Since the temperature dependency of biochemical reactions is highly non-linear (Langridge, 1963; Conant et al., 2011; Bonan, 2019), ultimate effects can be even more amplified and call for more studies including carbon (C) and nutrient cycles.

Simulated soil thaws up to 20 cm deeper in summer when accounting for pedon-scale heterogeneity, which is in line with previous studies (Beer, 2016; Nitzbon et al., 2021), and also the duration of thawing is extended by three days (M1) or four days (M2) in fall (Figs. 2, 3 and B2). Since the faster snow melt extends the snow-free time by up to five days in spring, snow melt and thawing effects may together prolong the simulated vegetation growing season by more than one week, when accounting for sub-grid scale heterogeneity in land-surface models. Additionally, the deeper thawing may allow plants to root deeper and liberate nutrients, which enhances plant growth, as well as it may promote decomposition of deeper soil organic matter. All these effects impact the simulations of ecosystem dynamics and thus may affect its simulated C balance significantly in consequence of the aggregation error due to neglecting sub-grid scale soil heterogeneity (Burke et al., 2012, 2020; Schädel et al., 2024).

Although we included lateral heat fluxes, we did not include lateral water fluxes in the model yet, which affect soil temperature directly through heat advection, but also indirectly due to redistribution of water and ice that largely influence thermal properties (Cai et al., 2020).

The neglect of heat advection (Gao and Coon, 2022) directly affects soil temperature because heat advection may be as strong or even stronger than heat conduction (Nicolsky et al., 2008; Bonan, 2019). In case heat conduction and heat advection work in the same direction, they would enhance their effects and by that mitigate soil temperature differences among columns faster. But in case they flow in opposite directions, they would attenuate the effect of each other and thus potentially diminish temperature differences among columns (Nicolsky et al., 2008; Aas et al., 2019).

Regarding the build of ground ice lateral water fluxes may reinforce the spatial occurrence of ice by flowing towards the freezing front and thus delivering more water that can freeze, or mitigate it, because the delivered water may transport heat towards the freezing front (Nicolsky et al., 2008; Aas et al., 2019; Fu et al., 2022).

Using only one test site reduces our ability to upscale our results to regional or even global scales. We cannot be certain that soil warming due to sub-grid scale soil heterogeneity will remain, if the experiment is replicated at another site or if the spatial resolution is increased. Also, the potential effects on ecosystem dynamics and the associated C budget may vary between sites and change with scale. However, the facts that (1) the aggregation error is consistent in its direction over the entire simulation period of more than ten years and (2) lateral fluxes actually amplify, rather than attenuate, the aggregation error of sub-grid scale soil heterogeneity in terms of soil temperature gives us an indication that the heterogeneity-induced soil temperature change, here a warming, as well as the potential effects on ecosystem dynamics, are likely to be persistent when the resolution is increased or the setup is transferred to another site.

An increase in resolution would cause an even broader averaging over soil properties and since we found that the actual distribution of heterogeneous variables such as OM content plays an important role and the relationship between variables and simulated soil temperature is highly non-linear, the risk of underestimating the aggregation error is high. Also, tiling may only partially account for this aggregation error by neglecting the true distribution of state variables.

When transferring the setup to another site, it would still imply a (linear) averaging over soil properties, which will cause the observed aggregation error due to the non-linearity of relationships. However, the direction of the aggregation error may change with site.

Another factor that limits the interpretability of our results is the use of Monin–Obukhov Similarity Theory (MOST) to derive surface fluxes from atmospheric variables (cf. Sect. 2.1.2). As noted above, MOST was applied as a pragmatic closure to parameterize surface–atmosphere exchange. This choice is consistent with common tile-based land surface modeling approaches (Bonan, 2019, chap. 7) in which subgrid surface heterogeneity is represented while atmospheric feedbacks at the scale of individual tiles are intentionally not resolved. However, this configuration represents a highly idealized modeling framework, and the assumptions underlying MOST are not strictly valid at the spatial scale of our simulations.

Despite this limitation, all soil patches in our simulations are driven by the same atmospheric forcing. Consequently, differences in the simulated turbulent fluxes can be attributed solely to variations in surface temperature and soil moisture conditions. The fluxes calculated at the patch scale should therefore not be interpreted as independent, physically local turbulent fluxes capable of driving atmospheric boundary-layer processes. Rather, they should be understood as diagnostic contributions to the area-averaged surface flux. Interpreted in this way, the results indicate that surface heterogeneity can influence the distribution of turbulent flux contributions at very small spatial scales and should therefore be regarded as an idealized representation of how subgrid soil heterogeneity may influence surface energy partitioning.

Finally, we would like to turn our perspective from our very extreme example with a cm-scale resolution back to LSMs, which usually run on scales of (hundreds of) kilometers. Our results clearly show that with decreasing model resolution, heterogeneity becomes more obvious, and heterogeneity induced lateral fluxes will become more important. Thus, we would like to emphasize that as soon as model developers think about decreasing model resolution in order to improve their simulation results, they have also to think about emerging heterogeneity and induced lateral fluxes and acknowledge their potential influence on simulation results, which may increase simulation uncertainty, especially in combination with the lack of resolution-suitable sub-surface soil data.

5 Conclusion

Simulation experiments with the high-resolution soil model DynSoM-2D show a notable – in our study negative – aggregation error in soil temperature, which affects the simulated depth of the active layer and the duration of the snow-free period in summer. Since the general temperature shift is more pronounced and more consistent when considering heterogeneity-induced lateral heat fluxes, we assume that it will likely remain at other sites and under changed environmental conditions, although the direction of the aggregation error may differ. Temperature changes are most pronounced in the near-surface soil, which make them especially relevant for ecosystem dynamics, but disappear with depth and appear to be less relevant for energy exchange with the atmosphere.

We traced the detected temperature shift back to the actual spatial distribution of organic matter, which nonlinearly alters soil thermal and hydrological properties. Due to this non-linear relationships and the importance of the actual organic matter distribution, current approaches to account for sub-grid scale heterogeneity, such as tiling, may still underestimate the aggregation error of soil heterogeneity and its effect on ecosystem dynamics. Thus, further studies are needed to better quantify these effects and estimate their importance on regional to global scale.

Appendix A: Methods

Here, we provide some additional information about the equations and assumptions used in DynSoM-2D. Further information can be also found in the code (Thurner et al., 2025).

A1 Calculation of κ and c_v

Heat conductivity (κ) and heat capacity (c_v) are calculated within the termal_properties.m routine, which follows Bonan (2019), chap. 5.4.

Heat conductivity is calculated for dry soil (κ_dry) and saturated soil (κ_sat) individually, before they are weighted accordingly to the Kersten number (Ke), which has a value of 0–1 depending on relative soil moisture (for unfrozen soils) or ice content, respectively (cf. Bonan, 2019, Eqs. 5.25–5.26).

$$\mathit{\kappa }={\mathit{\kappa }}_{\mathrm{dry}}+({\mathit{\kappa }}_{\mathrm{sat}}-{\mathit{\kappa }}_{\mathrm{dry}})\cdot \mathit{Ke}$$

κ_dry is calculated for mineral soil (κ_min) first (Bonan, 2019, Eq. 5.27) and then modified by OM-content (following Chadburn et al., 2015b):

$$\begin{array}{c}{\displaystyle }{\mathit{\kappa }}_{\min }={\displaystyle \frac{\mathrm{0.135}{\mathit{\rho }}_{\mathrm{b}}+\mathrm{64.7}}{\mathrm{2700}-\mathrm{0.947}{\mathit{\rho }}_{\mathrm{b}}}}\\ {\displaystyle }{\mathit{\kappa }}_{\mathrm{dry}}={{\mathit{\kappa }}_{\min }}^{\mathrm{1}-{\mathrm{\Theta }}_{\mathrm{OM}}}\cdot {{\mathit{\kappa }}_{\mathrm{OM}}}^{{\mathrm{\Theta }}_{\mathrm{OM}}}\end{array}$$

with ρ_b being bulk density $(=\mathrm{2700}(\mathrm{1}-{\mathrm{\Theta }}_{\mathrm{sat}}\left)\right)$, according to Bonan (2019), and κ_OM being the heat conductivity of dry OM (=0.25) taken from De Vries and Van Wijk (1963).

κ_sat is estimated from the individual heat capacities for solids (κ_soild), which is estimated using the Johansen method as described by Farouki (1981), water (κ_water), and ice (κ_ice) individually, that are then weighted accordingly to their actual volume fractions (Bonan, 2019, Eqs. 5.28–5.31):

$$\begin{array}{c}{\displaystyle }{\mathit{\kappa }}_{\mathrm{sat}}={\mathit{\kappa }}_{\mathrm{solid}}^{(\mathrm{1}-{\mathrm{\Theta }}_{\mathrm{sat}})}\cdot {\mathit{\kappa }}_{\mathrm{water}}^{\left(f_u{\mathrm{\Theta }}_{\mathrm{sat}}\right)}\cdot {\mathit{\kappa }}_{\mathrm{ice}}^{\left(\right(\mathrm{1}-f_u\left){\mathrm{\Theta }}_{\mathrm{sat}}\right)}\\ {\displaystyle }{\mathit{\kappa }}_{\mathrm{soild}}={\mathit{\kappa }}_q^q{\mathit{\kappa }}_{\mathrm{0}}^{\mathrm{1}-q}\end{array}$$

with q being the quartz content as fraction of the total soil solids.

Heat capacity is given as weighted fractions (Θ) of all soil components (solids, including mineral and OM compartments, water and ice, and air), with individual heat capacities taken from literature (cf. De Vries and Van Wijk, 1963; Bonan, 2019; Chadburn et al., 2015b).

$$\begin{array}{c}{\displaystyle }{c}_{\mathrm{v}}={\mathrm{\Theta }}_{\mathrm{dry}}\cdot c_{\mathrm{v}}^{\mathrm{dry}}+{\mathrm{\Theta }}_{\mathrm{water}}\cdot c_{\mathrm{v}}^{\mathrm{water}}+{\mathrm{\Theta }}_{\mathrm{ice}}\cdot c_{\mathrm{v}}^{\mathrm{ice}}+{\mathrm{\Theta }}_{\mathrm{air}}\cdot c_{\mathrm{v}}^{\mathrm{air}}\\ {\displaystyle }{c}_{\mathrm{v}}^{\mathrm{dry}}=(\mathrm{1}-{\mathrm{\Theta }}_{\mathrm{OM}})\cdot c_{\mathrm{v}}^{\min }+{\mathrm{\Theta }}_{\mathrm{OM}}\cdot c_{\mathrm{v}}^{\mathrm{OM}}\\ {\displaystyle }{c}_{\mathrm{v}}^{\min }=\mathrm{2.12}\times {\mathrm{10}}^{\mathrm{6}}\\ {\displaystyle }{c}_{\mathrm{v}}^{\mathrm{OM}}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}\\ {\displaystyle }{c}_{\mathrm{v}}^{\mathrm{water}}=\mathrm{4.19}\times {\mathrm{10}}^{\mathrm{6}}\\ {\displaystyle }{c}_{\mathrm{v}}^{\mathrm{ice}}=\mathrm{1.88}\times {\mathrm{10}}^{\mathrm{6}}\\ {\displaystyle }{c}_{\mathrm{v}}^{\mathrm{air}}=\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{3}}\end{array}$$

A2 Hydrological Scheme (incl. Snow Dynamics)

The hydrological scheme used by DynSoM-2D is based on Bonan (2019), chap. 8. The actual formulations can be found in HSM_water.m (Thurner et al., 2025).

HSM_water.m calculates the soil water balance for each layer (this routine is only vertically applied within this study!) based on precipitation as upper boundary condition and the matrix potential of each layer.

Precipitation is provided as forcing input and only modified given the actual temperature, i.e. added as rain, if air temperature is above/equal 0 °C, and added as snow, if air temperature is below 0 °C. Rain enters the system by infiltrating the soil directly or can either be stored as pond water on top of the surface or leaves the system as run-off water. Since we switched off pond formation for our simulations, all excess water leaves as run-off water. Snow is accumulating at the surface and builds additional layers with a thickness of 3 cm each (cf. Sect. 2.1). Snow layers are included into the energy balance and the heat conduction scheme, as soon and as long as they are “filled”, i.e. are 3 cm thick (specific snow properties can be found in Table A1). When snow temperature reaches melting temperature (0 °C), snow thaws and adds to meltwater that is treated similarly to rain, i.e. it can infiltrate soil, form a pond or leave the system as run-off. DynSoM-2D does not consider snow-internal refreezing, but considers snow compaction, i.e. the gravitational increase in density with time (Ekici et al., 2014; Beer et al., 2018) based on formulations by Hedstrom and Pomeroy (1998) and Verseghy (1991).

The matrix potential of each layer is calculated using the predictor-corrector method based on the Richards equation through Van Genuchten (1980) approach. VanGenuchten parameters are estimated using the pedotransfer functions by Wösten et al. (1999), and the matrix potential at saturation and the Clapp and Hornberger (1978) exponent are calculated using Cosby et al. (1984) functions.

Bonan (2019)Bonan (2019)Mellor (1977)Bonan (2019)Verseghy (1991)Hedstrom and Pomeroy (1998)

Table A1 Snow properties.

Appendix B: Results

This section provides some additional simulation results and additional plots, which are not discussed in detail within this manuscript, but may help the reader to understand the model performance.

B1 Additional Soil Parameters

In the pedotransfer functions, the following parameters are considered: van Genuchten parameters (α, m and n), residual porosity, saturated conductivity and from Cosby parameters (Cosby et al., 1984, saturated matrix potential and the b parameter). In addition, the pedotransfer function of the code includes the calculation of the field capacity matrix potential, field capacity effective saturation and field capacity moisture (all based in van Genuchten parameters), as well as initial porosity water and ice based on the calculated porosity and input fractions. Mean of the most relevant soil parameters used for our study for the homogeneous and the heterogeneous soil can be found in Table B1.

Table B1 Mean values for bulk density (ρ_b), as well as for the van Genuchten parameters (α, m, and n) for the homogeneous and the heterogeneous soil.

B2 Additional Plots

For better presenting deviations between model simulations, we do not only analyse actual simulation results, but also absolute differences (Eq. 5, ΔX) between the heterogeneous soil models (M1, M2) and the homogeneous soil model (M0).

B2.1 Temporal Pattern

The difference in soil porosity due to the difference between the homogeneous and the heterogeneous soils (Table 3) also impacts the soil water and ice contents (Fig. B2). Although the inter-model range is small, the intra-model ranges are wide, especially at the topsoil, where the heterogeneity is most pronounced, and the amount of organic material (OM), which has an extremely high porosity, differs most. In contrast to soil temperature simulations, lateral heat fluxes do not affect the intra-model ranges here, but they impact the inter-model range, i.e. the deviation of M2 from M0 compared to M1 to M0, because the warmer soil of M2 contains more liquid water and less ice during the summer.

Figure B1 Long-term simulation results (left; M0: black, M1: purple, M2: green) and deviations from reference model M0 (right; M1: purple, M2: green; Eq. 5) for surface heat fluxes (top row), surface temperatures (middle rows), and soil temperatures (bottom row) as extension to Fig. 2. For figure clarity, this figure shows only spatial mean values. (a) Simulated surface heat flux (F_surf), which is defined as sum of sensible heat flux (H) and latent heat flux (λE). (b) Differences of surface heat fluxes (ΔF_surf). (c) Simulated surface temperature (T_surf, solid lines), as well as soil temperature (T_soil, all depths, dashed lines), and atmospheric forcing temperature ($T_{\mathrm{atm}}^{\mathrm{forcing}}$, dotted line) for comparison. (d) Differences in surface temperature (ΔT_surf, solid lines) and soil temperature (ΔT_soil, all depths, dashed lines). (e) Simulated soil temperatures in topsoil (0–20 cm depth; T_top, solid lines), subsoil (20–80 cm depth; T_sub, dashed lines), and deep soil (80–100 cm depth; T_deep, dotted lines). (f) Differences in topsoil (ΔT_top, solid lines), subsoil (ΔT_sub, dashed lines), and deep soil temperatures (ΔT_deep, dotted lines).

Figure B2 Simulation results (left; M0: black, M1: purple, M2: green) and deviations from reference model M0 (right; M1: purple, M2: green; Eq. 5) for snow depth (top row), soil water content (kg water per kg dry soil; middle rows), and soil ice content (kg ice per kg dry soil; bottom row). Lines represent spatial mean, whereas shaded areas show the intra-model ranges, i.e. the ranges among columns for M1 and M2. (a) Simulated snow depth. Dashed line indicates 3 cm-threshold for snow effects on surface albedo and soil insulation. (b) Differences in simulated snow depth. (c) Simulated soil water content in topsoil (0–20 cm depth; ${\mathrm{\Theta }}_{\mathrm{w}}^{\mathrm{top}}$, solid lines), subsoil (20–80 cm depth; ${\mathrm{\Theta }}_{\mathrm{w}}^{\mathrm{sub}}$, dashed lines), and deep soil (80–100 cm depth; ${\mathrm{\Theta }}_{\mathrm{w}}^{\mathrm{deep}}$, dotted lines). (d) Differences in topsoil ($\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{w}}^{\mathrm{top}}$, solid lines), subsoil ($\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{w}}^{\mathrm{sub}}$, dashed lines), and deep soil water contents ($\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{w}}^{\mathrm{deep}}$, dotted lines). Black dashed line represents mean porosity (${\mathrm{\Theta }}_{\mathrm{pore}}^{\mathrm{soil}}$, Table 3) as indicator for an upper limit for Θ_w. (e) Simulated soil ice content in topsoil (0–20 cm depth; ${\mathrm{\Theta }}_{\mathrm{i}}^{\mathrm{top}}$, solid lines), subsoil (20–80 cm depth; ${\mathrm{\Theta }}_{\mathrm{i}}^{\mathrm{sub}}$, dashed lines), and deep soil (80–100 cm depth; ${\mathrm{\Theta }}_{\mathrm{i}}^{\mathrm{deep}}$, dotted lines). Black dashed line represents mean porosity (${\mathrm{\Theta }}_{\mathrm{pore}}^{\mathrm{soil}}$, Table 3) as indicator for an upper limit for Θ_i. (f) Differences in topsoil ($\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{i}}^{\mathrm{top}}$, solid lines), subsoil ($\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{i}}^{\mathrm{sub}}$, dashed lines), and deep soil ice contents ($\mathrm{\Delta }{\mathrm{\Theta }}_{\mathrm{i}}^{\mathrm{deep}}$, dotted lines).

B2.2 Spatial Pattern

Figure B3 Mean differences in simulated soil temperature per grid cell (Eq. 5). Upper row: Differences between M1 and M0. Lower row: Differences between M2 and M0. Colors represent differences in soil temperature. Thin dotted lines separate topsoil from subsoil from deep soil.

B2.3 Relation to Organic Material Content

Figure B4 Differences in simulated daily soil temperature per grid cell between M1 and M0 (purple), as well as M2 and M0 (green) calculated by Eq. (5). Grid cells are aggregated by relative changes in OM content between the heterogeneous and the homogeneous soil inputs ($\mathit{\delta }{\mathrm{\Theta }}_{\mathrm{OM}}=\frac{{\mathrm{\Theta }}_{\mathrm{OM}}^{\mathrm{het}}-{\mathrm{\Theta }}_{\mathrm{OM}}^{\mathrm{hom}}}{{\mathrm{\Theta }}_{\mathrm{OM}}^{\mathrm{hom}}}$). Boxes represent quartile range of data and whiskers show the 0.1–0.9-quantile range. Circles show data median, whereas stars indicate data mean.